1 Introduction

For Beyond the Standard Model (BSM) scenarios with additional new particles, the decays of these particles determine the experimental signals we would observe at collider experiments. If the new particles have a well separated mass spectrum, long decay chains will occur when a heavy new particle is produced. For decays involving coloured particles, hard quantum chromodynamic (QCD) radiation at each step in the decay chain will alter the structure of the event and therefore the number of events passing experimental selection criteria. The effects of radiation are also important in models with degenerate new particle mass spectra, where decay chains are typically limited to one step. Searches for these compressed spectra scenarios look for events in which hard radiation in the initial-state shower recoils against missing transverse energy in the final state to give an observable signal.Footnote 1 The emission of hard QCD radiation in the decay of new particles could either enhance or reduce this effect and so must be taken into account. Therefore, accurate simulation of hard radiation in the decays of BSM particles is necessary in order to optimise searches for new physics.

Monte Carlo event generators use fixed-order matrix elements combined with parton showers and hadronization models to simulate particle collisions. In the Herwig++ event generator [2, 3], the decays of unstable fundamental particles are treated separately from the hard process which produced them, prior to the parton shower phase, using the narrow width approximation. Decays are generated using the algorithm described in [4], which ensures spin correlations are correctly treated. The parton shower utilises an approximation that resumes the leading collinear and leading-colour soft logarithms [5] and so does not accurately describe QCD radiation in the regions of phase space where the transverse momenta of the emitted partons are high. The Positive Weight Hardest Emission Generator (POWHEG) formalism [6] is one method that allows the simulation of high transverse momentum (hard) radiation to be improved upon by using the real-emission matrix element to produce the hardest emission in the shower. This approach affects both the overall cross sections for inclusive processes and results in local changes to the shapes of distributions sensitive to the hardest emission. In particular, local changes to observables such as jet transverse momenta are important since they can impact on the proportion of events passing selection criteria in new physics searches. Since BSM signals often consist of only a few events, this can in turn result in significant changes to the exclusion bounds that can be set.

The POWHEG formalism has been successfully applied to a wide range of hard production processes, for example [7–22, 25–27], and particle decays [28, 29] in the Standard Model (SM) as well as selected BSM processes [30–34]. Next-to-leading order (NLO) corrections to BSM particle decays have also previously been studied, for example in [35] where the Supersymmetric-QCD correction to the decay \(\tilde{q} \rightarrow q \tilde{\chi }\) was calculated. In this work, we present results from the implementation of the POWHEG method in Herwig++ for a range of decays relevant for new physics searches. A similar approach based on generic spin structures is used to apply a matrix-element correction to hard radiation in particle decays in PYTHIA [36].

The POWHEG formalism will be reviewed in Sect. 2 and in Sect. 3 our implementation of the POWHEG correction will be described in full for the example of top quark decay. In Sect. 4, details of the decay modes implemented will be given. The impact of the correction on the decay of the lightest graviton in the Randall–Sundrum (RS) model [37] will be studied in Sect. 5.1. Finally, results from a selection of decays in the Constrained Minimal Supersymmetric Standard Model (CMSSM) will be presented in Sect. 5.2.

2 POWHEG method

In this section, a brief outline of the POWHEG method is given. Further details can be found in [38].

In the conventional parton shower approach, the inclusive differential cross section for the highest transverse momentum emission from an \(N\)-body process is given by

$$\begin{aligned} d \sigma ^\mathrm{{PS}}=B(\Phi _N) d\Phi _N \left[ \Delta (p_{T \mathrm {min}}) \, + \Delta (p_T) d\Phi _R {\mathcal {P}} \right] . \end{aligned}$$
(1)

Here we are considering a parton shower ordered in terms of the transverse momentum of the emitted parton, \(p_T\). \(\Phi _N\) are the phase space variables of the \(N\)-body leading-order (LO) process and \(B\) is the Born-level matrix element squared, including the relevant flux factor, such that the total LO cross section is \(\sigma ^\mathrm{{LO}} = \int B(\Phi _N) d\Phi _N\). \({\mathcal {P}}\) is the unregularized Altarelli–Parisi splitting kernel and \(\Phi _R\) is a set of variables parameterizing the phase space of the additional radiated parton. The radiative phase space is limited to the region \(p_T (\Phi _R)>p_{T \mathrm {min}}\), where \(p_{T \mathrm {min}}\) is a transverse momentum cut-off introduced to regularize the infra-red (IR) divergences in the splitting kernel. The Sudakov form factor for the parton shower is

$$\begin{aligned} \Delta (p_T) = \exp \left( - \int d\Phi _R {\mathcal {P}} \, \Theta \left( p_T \left( \Phi _R\right) -p_T \right) \right) . \end{aligned}$$
(2)

The square bracket in Eq. 1 integrates to unity which ensures that the total cross section is given by the LO result.

In the POWHEG approach, the inclusive differential cross section for the hardest emission is given by the QCD NLO differential cross section, that is

$$\begin{aligned} d \sigma ^\mathrm{{PO}}=\bar{B}(\Phi _N) d\Phi _N \!\left[ \Delta ^\mathrm{PO}(p_{T \mathrm {min}}) +\Delta ^\mathrm{PO}(p_T) d\Phi _R \frac{R(\Phi _N, \Phi _R)}{B(\Phi _N)}\right] \!,\nonumber \\ \end{aligned}$$
(3)

where \(\bar{B}(\Phi _N) \) is defined by

$$\begin{aligned} \bar{B}(\Phi _N)&= B(\Phi _N) + \left[ V(\Phi _N) + \int C(\Phi _N, \Phi _R)\,d\Phi _R\right] \nonumber \\&+ \int \left[ R(\Phi _N, \Phi _R) d\Phi _R - C(\Phi _N, \Phi _R) d\Phi _R\right] .\nonumber \\ \end{aligned}$$
(4)

The real-emission contribution, \(R(\Phi _N, \Phi _R)\), corresponds to the radiation of an additional parton from the LO interaction and the virtual contribution, \(V(\Phi _N)\), comes from the 1-loop correction to the LO process. \(C(\Phi _N, \Phi _R)\) is a counter term with the same singular behaviour as the real and virtual contributions and is introduced to ensure the two square brackets in Eq. 4 are separately finite. The Sudakov form factor appearing in Eq. 3 is

$$\begin{aligned} \Delta ^\mathrm{PO}(p_T) = \exp \left( - \int d\Phi _R \frac{R(\Phi _N, \Phi _R)}{B(\Phi _N)} \, \Theta \left( p_T \left( \Phi _R \right) -p_T \right) \right) . \end{aligned}$$
(5)

As with the conventional parton shower approach, the square bracket in Eq. 3 will integrate to unity and hence the total inclusive cross section will be given by the NLO result.

Typically, the counter term, \(C(\Phi _N, \Phi _R)\), can be rewritten as a sum of dipole functions, \({\mathcal {D}}_i \), each of which describes the behaviour of the real-emission matrix element in a singular region of phase space, i.e. when the emitted parton becomes soft or collinear to one of the legs in the Born process. By doing so, the different singular regions can be separated such that Eq. 5 becomes a product of Sudakov form factors

$$\begin{aligned} \Delta ^\mathrm{PO}(p_T) = \prod _i \exp \!\left( - \!\int d\Phi _R \frac{R(\Phi _N, \Phi _R) {\mathcal {D}}_i}{\sum _j {\mathcal {D}}_j B(\Phi _N)} \, \Theta \left( p_T \left( \Phi _R \right) -p_T \right) \right) \!,\nonumber \\ \end{aligned}$$
(6)

each of which describes the non-emission probability in a particular region of phase space specified by the dipole function \({\mathcal {D}}_i\).

When applying the POWHEG method to a parton shower ordered in transverse momentum, the hardest emission is generated first using the POWHEG Sudakov form factor in Eq. 6. Subsequent emissions are generated with the normal parton shower Sudakov given in Eq. 2, with the requirement that no parton shower emission has higher transverse momentum than the emission described by \(R(\Phi _N, \Phi _R)\). However, to allow QCD coherence effects to be included, an angularly ordered parton shower is used in Herwig++. Ordering the parton shower in terms of an angular variable means the first emission in the shower may not be the hardest. The POWHEG approach can be reconciled with angularly ordered parton showers by dividing the shower into several steps [6]. The hardest emission in the shower is generated first using the POWHEG Sudakov form factor and the value of the angular evolution variable corresponding to this emission is determined. An angularly ordered shower, running from the shower starting scale down to the scale of the hardest emission, is then generated. This truncated parton shower simulates coherent soft wide-angle radiation. The hardest emission is then inserted and the shower continues until the IR cut-off of the evolution variable is reached. Finally, in both stages of the parton shower, emissions generated by the shower are discarded if they have higher transverse momentum than the emission generated using the POWHEG Sudakov form factor.

3 Top quark decays

In this section, we describe our implementation of the POWHEG formalism for the example of a top quark decaying to a \(W\) boson and a bottom quark. To implement the full POWHEG correction to this decay, the Born configuration must be generated according to Eq. 4 and the hardest emission in the parton shower simulated using Eq. 6. However, in this work we consider only the effect of the POWHEG correction on the simulation of the hardest emission in the shower and hence generate the Born configuration using only \(B (\Phi _N)\), the leading order contribution in Eq. 4. As such, we use the existing LO Herwig++ implementation of top quark decay and modify the shower such that the hardest emission is generated according to Eq. 6. Justification for excluding the normalisation factor of the POWHEG correction will be given in Sect. 4.

Application of the POWHEG correction to top quark decays, along with top quark pair production in \(e^+ e^-\) collisions, has been previously studied in [28] for massless bottom quarks. In this work, we retain the mass of the bottom quark throughout.

3.1 Implementation in Herwig++

In Herwig++, the decays of fundamental particles are performed in the rest frame of the decaying particle. In this frame, we are free to choose the orientation of the \(W\) boson to be along the negative \(z\)-direction and so, at LO, the bottom quark is orientated along positive \(z\)-direction. The squared, spin and colour averaged matrix element for the LO process is given by

$$\begin{aligned} |{\mathcal {M}}_B|^2 = \frac{g^2}{4m_w^2} \left( m_t^4+m_b^4-2m_w^4+m_t^2m_w^2+m_b^2m_w^2-2m_t^2m_b^2 \right) \!,\nonumber \\ \end{aligned}$$
(7)

where \(m_t\), \(m_b\) and \(m_W\) are the masses of the top quark, bottom quark and \(W\) boson respectively and \(g\) is the weak interaction coupling constant. The relevant CKM factor has been set equal to 1.

The squared, spin and colour averaged matrix element for the \({\mathcal {O}} (\alpha _s)\) real-emission correction to the decay \(t \rightarrow W b\) is

$$\begin{aligned} |{\mathcal {M}}_R|^2&= g^2 g_s^2 C_F \left\{ - \frac{|{\mathcal {M}}_B|^2}{g^2} \left( \frac{p_b}{p_b.p_g} - \frac{p_t}{p_t.p_g} \right) ^2 \right. \nonumber \\&\quad + \left( \frac{p_g.p_t}{p_b.p_g} + \frac{p_b.p_g}{p_t.p_g} \right) \left( 1 + \frac{m_t^2}{2m_w^2} + \frac{m_b^2}{2m_w^2} \right) \nonumber \\&\quad \left. - \frac{1}{m_w^2} \left( m_t^2+m_b^2 \right) \right\} , \end{aligned}$$
(8)

where \(g_s\) is the strong coupling constant, \(C_F= \frac{4}{3}\) and \(p_t\), \(p_b\), \(p_W\) and \(p_g\) are the four-momenta of the top quark, bottom quark, \(W\) boson and gluon. In general, the orientation of the decay products in the three-body final state is such that the emitting parton absorbs the transverse recoil coming from the emission of the gluon and the spectator particle continues to lie along the negative \(z\)-direction. When the radiation originates from the top quark, the bottom quark effectively acts as the emitting particle so that we remain in the rest frame of the top quark. Therefore, for emission from both the top and the bottom quarks, the momenta of the decay products are

$$\begin{aligned} p_W = \left( E_W,0,0,-\sqrt{E_W^2 - m_W^2} \right) , \end{aligned}$$
(9)
$$\begin{aligned} p_b = \left( E_b, -p_T \cos \left( \phi \right) , -p_T \sin \left( \phi \right) , \sqrt{E_b^2 \!-\! p_T^2\! \!-\!m_b^2} \right) ,\nonumber \\ \end{aligned}$$
(10)
$$\begin{aligned} p_g = \left( E_g, p_T \cos \left( \phi \right) , p_T \sin \left( \phi \right) , \sqrt{E_g^2 - p_T^2 } \right) , \end{aligned}$$
(11)

where \(E_x\) is the energy of particle \(x\), and \(p_T\) and \(\phi \) are the transverse momentum and azimuthal angle of the gluon.

The Lorentz invariant phase space element, \({\mathrm {d}}\Phi _R\), describing the emission of the additional gluon is obtained from the relation

$$\begin{aligned} {\mathrm {d}}\Phi _3 = {\mathrm {d}}\Phi _2 {\mathrm {d}}\Phi _R , \end{aligned}$$
(12)

where

$$\begin{aligned} {\mathrm {d}}\Phi _N = \left( 2 \pi \right) ^4 \delta ^4 \left( p_t - \sum _{i=1}^N p_i \right) \prod _{i=1}^{N} {\frac{{\mathrm {d}}^3{\mathbf {p_i}}}{2E_i(2\pi )^3}} \end{aligned}$$
(13)

and \({\mathbf {p_i}}\) is the three-momentum of particle \(i\). We choose to parameterize the radiative phase space in terms of the transverse momentum, \(p_T\), rapidity, \(y\), and azimuthal angle, \(\phi \), of the gluon and so find

$$\begin{aligned} {\mathrm {d}}\Phi _R = J {\mathrm {d}}p_T dy {\mathrm {d}}\phi , \end{aligned}$$
(14)

where the Jacobian factor, \(J\), isFootnote 2

$$\begin{aligned} J = \frac{1}{8\pi ^2} \frac{m^2_t p_T |{\mathbf {p_W}}|^2}{\lambda (m_t^2, m_W^2, m_b^2) [|{\mathbf {p_W}}|(m_t-p_T\cosh y)-E_W p_T\sinh y]}.\nonumber \\ \end{aligned}$$
(15)

This parametrization has the advantage of simplifying the Heaviside function in the POWHEG Sudakov form factor to a lower limit in the integration over \(p_T\).

The final components required for the implementation of the POWHEG Sudakov form factor in Eq. 6 are the dipole functions, \({\mathcal {D}}_i\), which describe the singular behaviour of the real-emission matrix element. We use the dipole functions defined in the Catani–Seymour subtraction scheme, details of which can be found in [39, 40], to describe the singular behaviour resulting from emissions from the decay products. The dipole used to describe radiation from the top quark is as follows

$$\begin{aligned} {\mathcal {D}}_i = \frac{-4 \pi C_F \alpha _s }{E^2_g} |{\mathcal {M}}_B|^2. \end{aligned}$$
(16)

It contains only soft enhancements since, in the top quark rest frame, collinear enhancements are suppressed.

Using the above information, the hardest emission in the shower can then be generated according to Eq. 6 using the veto algorithmFootnote 3, which proceeds as follows:

  1. 1.

    Trial values of the radiative phase space variables are generated. The transverse momentum of the emission is generated by solving

    $$\begin{aligned} \Delta ^\mathrm{over}(p_T) = \exp \left( - \int ^{p^\mathrm{{max}}_T}_{p_T} \frac{C\left( y_\mathrm{{max}}-y_\mathrm{{min}}\right) }{p_T\left( \Phi _R\right) } {\mathrm {d}}p_T\left( \Phi _R\right) \right) = {\mathcal {R}},\nonumber \\ \end{aligned}$$
    (17)

    where \(p^\mathrm{{max}}_T = \frac{(m_t -m_W)^2 -m_b^2}{2 (m_t - m_W )}\) is the maximum possible \(p_T\) of the gluon. \(y_\mathrm{max}\) and \(y_\mathrm{min}\) are the upper and lower bounds on the gluon rapidity, chosen to overestimate the true rapidity range. \(C\) is a constant chosen such that the integrand in Eq. 17 always exceeds the integrand in Eq. 6 and \({\mathcal {R}}\) is a random number distributed uniformly in the range \([0,1]\). Values of \(y\) and \(\phi \) are generated uniformly in the ranges \([y_\mathrm{min}, y_\mathrm{max}]\) and \([0, 2 \pi ]\) respectively;

  2. 2.

    If \(p_T < p^\mathrm{min}_T\), no radiation is generated and the event is hadronized directly. We set \(p^\mathrm{min}_T=1 \,{\mathrm {GeV}}\) throughout this work;

  3. 3.

    If \(p_T \ge p^\mathrm{min}_T\), the momenta of the \(W\) boson, bottom quark and gluon are calculated using the generated values of the radiative variables. Doing so, yields two possible values of \(E_W\) that must both be retained and used in the remainder of the calculation. If the resulting momenta do not lie within the physically allowed region of phase space, we veto this configuration, set \(p^\mathrm{max}_T = p_T\) and return to step 1;

  4. 4.

    Events within the physical phase space are accepted with a probability given by the ratio of the true to overestimated integrands in Eqs. 6 and 17 respectively. If the event is rejected, we set \(p^\mathrm{{max}}_T = p_T\) and return to step 1;

Using this procedure, a trial emission is generated for each dipole, \({\mathcal {D}}_i\), in Eq. 6. The configuration which gives the highest \(p_T\) emission is selected. The existing Herwig++ framework, detailed in [9], is then used to generate the remainder of the parton shower.

3.2 Parton level results

To validate our implementation of the algorithm described in Sect. 3.1, Dalitz style plots were generated for the decay \(t \rightarrow W b\) and are shown in Fig. 1. The Dalitz variables, \(x_W\) and \(x_g\), were defined by the relation, \(x_i = \frac{2E_i}{m_t}\), where \(E_i\) is the energy of particle \(i\) in the rest frame of the top quark. The left-hand plot in Fig. 1 shows the distribution obtained using the POWHEG style correction. In this case, \(x_g\) is the energy fraction of the gluon generated using the full real-emission matrix element. The distribution on the right-hand side of Fig. 1 was generated using the conventional parton shower, limited to one emission in the final state, and so here \(x_g\) is the energy fraction of a gluon produced using the parton shower splitting kernels. On both distributions, the black outline indicates the physical phase space boundaries. The enclosed area is divided into a section populated by radiation from the bottom quark (above the green dashed line), sections populated by radiation from the top quark (below the blue dotted lines) and a dead region (between the blue dotted and green dashed lines) that corresponds to hard gluon radiation and is not populated by the conventional parton shower. These boundaries correspond to the theoretical limits of the Herwig++ parton shower with symmetric phase space partitioning, described in [42], in which the starting values of the shower evolution variables for the top and bottom quarks are chosen such that the volumes of phase space accessible to emissions from each quark are approximately equal.

Fig. 1
figure 1

Dalitz distributions for the decay \(t \rightarrow W b\) with (left) and without (right) the POWHEG style correction. The black outline indicates the physically allowed region of phase space. In the conventional parton shower approach, the region above the green dashed line is populated with radiation from the bottom quark and the regions below the blue dotted lines with radiation from the top quark. These boundaries correspond to the limits of the parton shower with symmetric phase space partitioning

Full size image

As expected, in both plots we see a high density of points in the limit \(x_g \rightarrow 0\), corresponding to soft gluon emission. The POWHEG corrected distribution also has a concentration of points along the upper physical phase space boundary where \(x_W\) is maximal and emissions are collinear to the bottom quark. The density of points along the upper boundary is reduced in the parton shower distribution and points are instead concentrated along the lower boundary of the bottom quark emission region. As discussed in [42], the parton shower approximation agrees with the exact matrix element in the case of collinear radiation from the bottom quark but overestimates it elsewhere in the bottom quark emission region. The factor by which the parton shower approximation exceeds the exact matrix element, increases towards the lower boundary of the region and therefore we see an excess of points near the boundary. The parton shower distribution also has a high density of points in the top quark emission region for \(x_g \lesssim 0.53\). This enhancement is again the result of the parton shower approximation overestimating the exact matrix element in this area [42]. In general, we see that the parton shower in Herwig++ produces areas of high emission density which do not correspond to physically enhanced areas of phase space and therefore has a tendency to overpopulate hard regions of phase space. On the other hand, the POWHEG emission is distributed according to the exact real-emission matrix element and so correctly populates the physically enhanced regions of phase space with no additional spurious high density regions. Finally, we also see that the POWHEG corrected distribution fills the dead region of phase space that is not populated by the standalone parton shower.

To study the impact of the POWHEG style correction on top quark decays, parton level \(e^+ e^- \rightarrow t \bar{t}\) events were simulated and analysed as in [43]. Events were generated at a centre-of-mass energy close to the \(t \bar{t}\) threshold, \(\sqrt{s} =360 \,{\mathrm {GeV}}\), to minimize the effects of radiation from the initial-state shower. Unless otherwise stated, in this study we use the default set of tuned perturbative and non-perturbative parameters, or event tune, in Herwig++ version 2.6 [3]. Final-state partons were clustered into three jets using the FastJet[44] implementation of the \(k_T\) algorithm. The \(W\) bosons were decayed leptonically and their decay products excluded from the jet clustering. Events were discarded if they contained a jet with \(p_{T}<10 \,{\mathrm {GeV}}\) or the minimum jet separation,Footnote 4 \(\Delta R\), did not satisfy \( \Delta R \ge 0.7\). Using events that passed these selection criteria, differential distributions were plotted of \(\Delta R \) and \(\log (y_{32})\), where \(y_{32}\) is the value of the jet resolution parameterFootnote 5 at which a three jet event is classified as a two jet event. The resulting distributions are shown in the left and right-hand plots in Fig. 2. Distributions generated using the normal parton shower and the parton shower including the POWHEG style correction, are shown by the blue dashed and black solid lines respectively. The red dotted lines in Fig. 2 show the distributions obtained when the existing implementation of hard and soft matrix element corrections (MEC) [43] are applied to the normal Herwig++ parton shower. Hard matrix element corrections use the full \(t \rightarrow W b g\) matrix element to distribute emissions in the dead regions of phase space that are not populated by the parton shower. Soft matrix element corrections use the full real-emission matrix element to correct emissions generated by the parton shower that lie outside the areas of phase space where the parton shower approximation is valid, i.e. away from the soft and collinear limits. Applying these corrections ensures that the hardest emission in the shower is generated according to the exact matrix element, therefore, we expect a high level of agreement between the POWHEG and matrix element corrected distributions. The bottom panel in each plot shows the ratio of the parton shower and matrix element corrected distributions to the POWHEG corrected distribution. In both plots, we include error bars indicating the statistical uncertainty.

Fig. 2
figure 2

Comparison of distributions generated using the standalone parton shower with those generated using a matrix element or POWHEG style correction to the decay \(t \rightarrow W b\). Parton level \(e^+e^- \rightarrow t \bar{t}\) events were generated at \(\sqrt{s}=360 \,{\mathrm {GeV}}\). The left-hand plot shows the distribution of the minimum jet separation, \(\Delta R\), and the right-hand plot the logarithm of the jet measure, \(y_{32}\)

Full size image

As discussed in [43], applying the matrix element corrections has the effect of softening both the \(\Delta R\) and \(\log (y_{32})\) distributions. This is due to the soft matrix element correction rejecting a portion of the high \(p_T\) emissions generated by the parton shower. The magnitude of the observed effect illustrates the importance of matching the parton shower to the exact matrix element in high \(p_T\) regions. As expected, the distributions generated using the POWHEG style and matrix element corrections are very similar although, for both variables, the POWHEG style correction yields slightly harder distributions. The discrepancies between the distributions are the result of a number of subtle differences between the POWHEG and matrix element correction schemes. Firstly in the matrix element correction approach, events in the dead region are generated using the fixed-order real-emission matrix element only, without any Sudakov suppression, and subsequent showering of the resulting configuration is simulated starting from the \(1 \rightarrow 3\) process. However, in the POWHEG approach the hardest emission in the shower is reinterpreted such that the conventional parton shower instead begins from the Born hard configuration. The scale of the hardest emission is generated, and then the shower proceeds as normal except that the hardest emission is fixed at the generated scale. In addition to this, the soft matrix element correction is applied to all emissions in the parton shower which are the hardest so far. Normally this leads to the correction of both the hardest emission and a number of other emissions with large values of the evolution parameter, but smaller transverse momentum. These differences all contribute to the discrepancies between the POWHEG style and matrix element corrected distributions although it is unclear which would have the largest effect. However, the difference between the POWHEG style and matrix element corrected results is comparatively small. The agreement between the two approaches serves to further validate our implementation of the POWHEG formalism. Finally, we note that the POWHEG style approach is preferable to the original matrix element correction scheme since it is significantly simpler to implement in Herwig++.

4 Decays of BSM particles

As discussed in Sect. 1, it is important that the simulation of QCD radiation in the decays of BSM particles is done in the most accurate way possible. In this work, we present results illustrating the effect of consistently matching the QCD real-emission matrix element with the parton shower in Herwig++ through the POWHEG formalism. This technique has been applied to a range of decays that occur in most of the well studied BSM scenarios. Table 1 shows the combinations of incoming and outgoing spins for which this method is used and each spin structure is implemented for the colour flows given in Table 2. However, models with coloured tensor particles are beyond the scope of this work and therefore decays involving incoming tensor particles were limited to colour flows in which the tensor is a colour singlet.

Table 1 Spin combinations for which the POWHEG correction has been applied
Full size table
Table 2 Colour flows for which the POWHEG correction has been applied
Full size table

The LO and real-emission matrix elements appearing in the POWHEG Sudakov form factor in Eq. 5 are calculated using helicity amplitude methods to correctly incorporate spin correlations [4]. The dipole functions, \({\mathcal {D}}_i \), are defined as in the Catani–Seymour dipole subtraction method [39, 40] when describing radiation from the decay products. In this approach, dipoles describing quasi-collinear radiation from massive vector bosons are not well defined. Therefore, the Fermion-Fermion-Vector, Scalar-Scalar-Vector and Tensor-Vector-Vector decays are limited to the situation where any final-state coloured vector particles are massless. The Vector-Fermion-Fermion and Vector-Scalar-Scalar decays do, however, include radiation from massive incoming vector particles. Decays are performed in the rest frame of the decaying particle [2] and therefore the dipole describing the singular behaviour of this particle will only contain a universal soft contribution. This is a well defined, spin-independent function given, for the example colour flow \(3 \rightarrow 3\,0\), by Eq. 16.

Finally, in this work we focus solely on the effect of the POWHEG correction on the simulation of the hardest emission in the shower and have not implemented the normalisation factor coming from the presence of \(\bar{B}\) rather than \(B\) in Eq. 3. In many cases, the partial widths and branching ratios used in the simulation are calculated by an external program, for example SDECAY [45], and so already include NLO corrections. These values are then passed to Herwig++ by means of a spectrum file in SUSY Les Houches Accord format [46, 47]. In cases where the calculation of the widths and branching ratios is performed in Herwig++, generated distributions can be rescaled by a global normalisation factor to achieve NLO accuracy for suitably inclusive observables when the necessary calculations exist.

5 Results

5.1 Randall–Sundrum graviton

The effect of applying the POWHEG correction to the decay of the lightest RS graviton was investigated using the Herwig++ implementation of the RS model. LHC proton-proton collisions with a centre-of-mass (CM) energy of \(\sqrt{s}=8 \,{\mathrm {TeV}}\) were simulated. The lightest graviton, \(G\), was produced as a resonance and allowed to decay via \(G \rightarrow gg \) and \(G \rightarrow q \bar{q}\) for \(q = u,d,s,c,b \). The mass of the graviton was chosen to be \(m_G = 2.23 \,{\mathrm {TeV}}\) which corresponds to the lower bound on the allowed graviton mass for the coupling \(k/ \bar{M}_{pl} = 0.1\) in [48]. An analysis based on the ATLAS experiment’s search for new phenomena in dijet distributions [49] was then carried out. Jets were constructed using the FastJet   [44] implementation of the anti-\(k_t\) algorithm [50] with the energy recombination scheme and a distance parameter \(R=0.6\). Jets with \(|y| \ge 4.4\) were discarded, where \(y\) is the rapidity of the jet in the \(pp\) CM frame. Events with less than two jets passing this constraint were vetoed. The rapidities of the two highest \(p_T\) jets in the \(pp\) CM frame are given by \(y_{1}\) and \(y_{2}\). In the dijet CM frame formed by the two hardest jets, their corresponding rapidities are \(y_{*}\) and \(-y_{*}\) where \(y_{*}= \frac{1}{2} (y_1 - y_2)\). Events not satisfying \(|y_{*}| < 0.6\) and \(|y_{1,2}| < 2.8\) were discarded. The dijet invariant mass, \(m_{jj}\), was formed from the vector sum of the two hardest jet momenta and events were vetoed if \(m_{jj} \le 1.0 \,{\mathrm {TeV}}\).

The dijet mass distribution after the above selection criteria were applied, is shown in the left-hand plot in Fig. 3. The blue dashed line shows the invariant mass distribution for the LO matrix element combined with the parton shower while the black solid line shows the result including the POWHEG correction to the graviton decay. Both distributions were generated using the optimum set of tuned perturbative and non-perturbative parameters found in [29]. From Fig. 3, we see that including the POWHEG correction causes a decrease of \({\mathcal {O}} (40\,\%)\) in the number of events in the region \(2.1 \,{\mathrm {TeV}}\le m_{jj} \le 2.3 \,{\mathrm {TeV}}\). This effect is highlighted in the right-hand plot in Fig. 3, which shows the dijet mass distribution in this range. In the conventional parton shower approach, the majority of the graviton’s momentum will be carried by the two partonic decay products. When the POWHEG correction is applied, the highest \(p_T\) emission in the shower will typically be quite hard and so a significant fraction of the graviton’s momentum will be missed by considering the invariant mass of only the hardest two jets, therefore shifting the distribution to lower values of \(m_{jj}\).

Fig. 3
figure 3

Dijet invariant mass distribution for the lightest RS graviton decaying to jets. The left-hand plot shows the distribution in the full range \(1.0 \,{\mathrm {TeV}}\le m_{jj} \le 2.5 \,{\mathrm {TeV}}\) while the right-hand plot emphasises the effect on the peak region \(2.1 \,{\mathrm {TeV}}\le m_{jj} \le 2.3 \,{\mathrm {TeV}}\). The mass of the graviton was \(m_G = 2.23 \,{\mathrm {TeV}}\) and the coupling \(k/ \bar{M}_{pl} = 0.1 \). LHC events were simulated with \(\sqrt{s}=8 \,{\mathrm {TeV}}\). The yellow and orange bands were generated by varying the event tune parameters in the POWHEG corrected and conventional parton shower distributions respectively

Full size image

To give an estimate of the uncertainty arising from our choice of event tune, the dijet mass distributions were generated at ten points in the event tune parameter space and error bands were created showing the maximum and minimum values from the resulting set of distributions. A description of the varied parameters can be found in [29] and their values at each of the ten points are given in Tab. 2 of [29]. The error bands are shown in yellow and orange for the distributions with and without the POWHEG correction respectively. The impact of the POWHEG correction is still clearly evident once this uncertainty has been taken into account.

5.2 Constrained minimal supersymmetric standard model

In addition to the results presented in Sect. 5.1, the effect of the POWHEG correction was also studied in the context of the CMSSM model. The high scale parameters of the model were chosen to be \(m_0 = 1,220 \,{\mathrm {GeV}}\), \(m_{1/2} = 630 \,{\mathrm {GeV}}\), \(\tan {\beta } = 10\), \(A_0 = 0\) and \(\mu >0\). This point lies just outside the exclusion limits set by the ATLAS experiment in [51].The corresponding weak scale parameters and decay modes were calculated using ISAJET 7.80 [52] and the resulting masses of the Supersymmetric (SUSY) particles relevant to this study are given in Table 3. The Herwig++ implementation of the MSSM model was used to generate LHC \(pp\) collisions at a centre-of-mass energy of \(\sqrt{s}=8 \,{\mathrm {TeV}}\). Here we focus on the effect of the correction to the parton shower and so hadronization and the underlying event are not simulated. In the following sections, the impact of the POWHEG correction on two archetypal decays is presented. In both cases, the decaying SUSY particle is pair produced in the hard process and the two subsequent decays are then analysed separately in the rest frame of the decaying particle. Dalitz style distributions were produced, as described in Sect. 3.2, for both the POWHEG corrected emission and the normal parton shower limited to one final-state emission. In addition, transverse momentum distributions of the hardest jet not coming from a visible decay product were also studied. To do so, the full parton shower was generated, with and without the POWHEG style correction, and events were analysed by clustering all visible final-state particles into jets using the FastJet implementation of the anti-\(k_T\) algorithm with the energy recombination scheme and \(R=0.4\). Jets with \(p_T \le 20 \,{\mathrm {GeV}}\) or \(|\eta |>4.0 \) were discarded. Events were required to have at least \(n+1\) jets passing the selection criteria, where \(n\) is the number of visible decay products.

Table 3 Masses of the SUSY particles relevant to the decays studied in Sects. 5.2.1 and 5.2.2
Full size table

5.2.1 \(\tilde{u}_L \rightarrow u \, \tilde{\chi }_1^0\)

Events were generated in which \(\tilde{u}_L\) and its associated anti-particle were produced and then decayed via the mode \(\tilde{u}_L \rightarrow u \, \tilde{\chi }_1^0\). Dalitz style distributions with and without the POWHEG correction were produced and are shown in the left and right-hand plots in Fig. 4. The black outline indicates the kinematic limits of phase space and the green dashed and blue dotted lines are the boundaries of the emission regions of the conventional parton shower with the most symmetric choice of shower phase space partitioning. Emissions from the up quark populate the area above the green dashed line, while the regions below the blue dotted lines are filled by emissions from the \(\tilde{u}_L\). The area between the green and blue lines is the dead zone, unpopulated by the normal parton shower. In the POWHEG corrected distribution, points are concentrated in the soft region as \(x_g \rightarrow 0\) and along the upper boundary of physical phase space where the gluon in collinear to the up quark. However, in the normal parton shower distribution fewer points lie along the upper physical phase space boundary and instead there is an concentration of points in the \(\tilde{u}_L\) emission region with \( x_g \lesssim 0.85\) and along the lower boundary of the up quark emission region. In analogy to the case of top quark decay, it is likely that these unphysical high density regions are due to the parton shower kernels overestimating the exact real-emission matrix element. Finally, we see that including the POWHEG correction ensures that the region of phase space inaccessible to the normal parton shower is populated.

Fig. 4
figure 4

Dalitz distributions for the decay \(\tilde{u}_L \rightarrow u \, \tilde{\chi }_1^0\) with (left) and without (right) the POWHEG style correction. The black outline indicates the physically allowed region phase space. In the conventional parton shower approach, the region above the green dashed line is populated with radiation from the up quark and the regions below the blue dotted lines with radiation from the \(\tilde{u}_L\). These boundaries correspond to the limits of the parton shower with symmetric phase space partitioning

Full size image

Differential distributions of the transverse momentum of the subleading jet,Footnote 6 \(p_{T,2}\), in each decay were also generated and are shown in Fig. 5. The blue dashed line corresponds to the distribution generated using the LO matrix element combined with the parton shower while the black solid line shows the result with the POWHEG correction to the decay applied. The bottom panel in Fig. 5 shows the ratio of the parton shower and POWHEG corrected results and in both distributions error bars are included to indicate statistical uncertainty. As demonstrated in Fig. 4, the parton shower has a tendency to over-populate the hard regions of phase space. Hence, including the POWHEG correction reduces the \(p_T\) of the hardest emission in the decay. This phenomenon is reflected in the \(p_{T,2}\) distributions. When the POWHEG correction is applied, the \(p_{T,2}\) distribution is softened such that there is a reduction in the number of events passing the jet \(p_T\) selection criteria of \({\mathcal {O}} (20\,\%)\). The softening is less pronounced at low values of \(p_{T,2}\) where the parton shower splitting kernels give a good approximation to the exact matrix element. Here the standalone parton shower and POWHEG corrected distributions are similar. At larger values of \(p_{T,2}\), the impact of the POWHEG correction is again reduced as, in this region, the subleading jet in the POWHEG corrected distribution typically has a significant contribution from partons generated by the normal parton shower in addition to the hardest emission coming from the POWHEG correction.

Fig. 5
figure 5

Transverse momentum distributions of the second hardest jet in the decay \(\tilde{u}_L \rightarrow u \, \tilde{\chi }_1^0\) in the rest frame of the \(\tilde{u}_L\). Events were generated with and without the POWHEG correction using the CMSSM model with \(m_0 = 1,220 \,{\mathrm {GeV}}\), \(m_{1/2} = 660 \,{\mathrm {GeV}}\), \(\tan {\beta } = 10\), \(A_0 = 0\) and \(\mu >0\) at the LHC with \(\sqrt{s}=8 \,{\mathrm {TeV}}\)

Full size image

5.2.2 \(\tilde{g} \rightarrow \tilde{t}_1 \, \bar{t}\)

Finally, we investigate the impact of the POWHEG style correction on the decay mode \(\tilde{g} \rightarrow \tilde{t}_1 \, \bar{t}\). The left and right-hand plots in Fig. 6 show Dalitz distributions for this decay with and without the POWHEG correction respectively. In both plots, the black outline indicates the kinematically allowed region phase space. The solid coloured lines show the boundaries of the parton shower emission regions in the scenario where the \(\bar{t}\) absorbs the \(p_T\) of the gluon and the \(\tilde{t}_1\) is orientated along the negative \(z\)-axis in the \(\tilde{g}\) rest frame. The region above the pale green line is populated by emissions from the \(\bar{t}\) and the areas below the dark blue lines are filled by emissions from the \(\tilde{g}\). In this scenario, the two emission regions overlap and there is no region of phase space left unpopulated by the parton shower. The dashed coloured lines indicate the emission boundaries of the parton shower when the \(\tilde{t}_1\) absorbs the transverse recoil of the emission and the \(\bar{t}\) is aligned with the negative \(z\)-axis. The pale green dashed line is the upper limit for emissions coming from the \(\tilde{t}_1\) and the dark blue dashed lines are the lower boundaries from emissions from the \(\tilde{g}\). From the left-hand plot of Fig. 6, we see that the majority of points in the POWHEG corrected distribution are concentrated in the soft region of phase space. High density regions corresponding to emissions collinear to the \(\bar{t}\) or \(\tilde{t}_1\) are suppressed due to the large masses of the decay products. In the parton shower distribution, points are concentrated in the soft region and along the lower boundary of the \(\bar{t}\) and dashed \(\tilde{g}\) emission regions. The latter two unphysical regions of over-population again highlight the importance of correcting hard emissions in the parton shower using the exact real-emission matrix element.

Fig. 6
figure 6

Dalitz distributions for the decay \(\tilde{g} \rightarrow \tilde{t}_1 \, \bar{t}\) with (left) and without (right) the POWHEG style correction applied. The solid (dashed) coloured lines indicate the parton shower emission regions when the \(\bar{t}\) \((\tilde{t}_1)\) absorbs the transverse recoil of the emission. The solid (dashed) pale green line shows the lower (upper) boundary for radiation from the \(\bar{t}\) \((\tilde{t}_1)\). The dark blue solid (dashed) lines are the equivalent upper (lower) boundaries for radiation from the \(\tilde{g}\). All boundaries correspond to the case of symmetric phase space partitioning and the black outline shows the kinematically allowed region of phase space

Full size image

The transverse momentum distribution of the third hardest jet in the rest frame of the \(\tilde{g}\) were also plotted and are shown in Fig. 7. To focus on the effect of the POWHEG correction, the decay products, \(\bar{t}\) and \(\tilde{t}_1 \), were not allowed to decay further. The blue dashed and black solid lines in Fig. 7 correspond to the parton shower and POWHEG correction distributions respectively. The bottom panel of the plot shows the ratio of the parton shower and POWHEG corrected results and in both distributions error bars are included to indicate statistical uncertainty. As in Sect. 5.2.1, we find that the POWHEG correction decreases the total number of events passing the jet \(p_T\) selection criterion. The effect is more pronounced in this case, with an \({\mathcal {O}} (40\,\%)\) reduction. The parton shower distribution significantly exceeds the POWHEG corrected distribution at small \(p_{T,3}\), however, at higher values of \(p_{T,3}\) the two distributions are similar. At lower values of \(p_{T,3}\), the main contribution to the third hardest jet in the POWHEG corrected distribution is from the hardest emission in the decay, generated using the real-emission matrix element. Therefore, we expect the uncorrected distribution to exceed the corrected one in this region. However, the maximum possible \(p_T\) of the gluon generated by the POWHEG correction isFootnote 7 \(p^\mathrm{{max}}_T \approx 75\,{\mathrm {GeV}}\). Jets contributing to the POWHEG corrected distribution above this limit include a number of other partons generated by the normal parton shower in addition to the hardest emission. This reduces the effect of the correction at higher values of \(p_{T,3}\). Therefore, we find that applying the POWHEG correction has a more significant impact on the number of events passing selection criteria when the value of the \(p_{T,3}\) selection criterion lies below \(p^\mathrm{max}_T\) of the gluon produced in the POWHEG correction.

Fig. 7
figure 7

Comparison of parton level distributions generated with and without the POWHEG correction for the decay \(\tilde{g} \rightarrow \tilde{t}_1 \, \bar{t}\) with stable decay products. Results are for the CMSSM model with \(m_0 = 1,220 \,{\mathrm {GeV}}\), \(m_{1/2} = 660 \,{\mathrm {GeV}}\), \(\tan {\beta } = 10\), \(A_0 = 0\) and \(\mu >0\) and LHC events with \(\sqrt{s}=8 \,{\mathrm {TeV}}\). Shown are the \(p_T\) distributions of the third hardest jet in the rest frame of the \(\tilde{g}\)

Full size image

6 Conclusions

In this work, we used the real-emission matrix element to generate hard QCD radiation in a range of particle decays in the Herwig++ event generator. This method is particularly relevant to new physics searches based on the decays of heavy new particles. The POWHEG corrections to these decays can change the shapes of certain experimental observables, thus altering the number of signal events passing selection criteria and modifying the exclusion bounds that can be set on the masses of the new particles. This correction will be available in Herwig++ version 2.7.

The algorithm used to implement the POWHEG style correction in Herwig++ was described in detail for the decay \(t \rightarrow W b\). Dalitz style distributions of the first emission in the conventional parton shower and POWHEG corrected approach were produced and showed that, while the POWHEG style correction ensures the majority of emissions lie in the soft and collinear limits, the parton shower has erroneous, unphysical regions of high emission density. This causes the parton shower to overpopulate the high \(p_T\) regions of phase space. Differential distributions of the minimum jet separation and logarithm of the jet measure were also generated with the POWHEG style correction and compared to those generated with the existing Herwig++ implementation of hard and soft matrix element corrections. The two techniques exhibit a high level of agreement therefore demonstrating the validity of our approach. In addition to this, distributions were generated using the normal parton shower. In agreement with the results from the Dalitz plots, these distributions were found to be considerably harder than those generated with the matrix element or POWHEG style corrections.

The impact of applying the POWHEG style correction to a BSM decay was studied by plotting the invariant mass distribution of dijets produced in the decay of the lightest RS graviton, \(G \rightarrow gg\) and \(G \rightarrow q \bar{q}\). Applying the POWHEG correction was found to have a considerable impact on the height of the distribution in the dijet mass peak. The number of events passing selection criteria in the mass range \(2.1 \,{\mathrm {TeV}}\le m_{jj} \le 2.3 \,{\mathrm {TeV}}\) dropped by \({\mathcal {O}} (40\,\%)\) when the correction was applied. This is a consequence of the dijet invariant mass not including the hardest emission in the shower that carries a significant fraction of the graviton’s momentum when it is simulated using the real-emission matrix element. The sizable impact of the correction in this scenario illustrates the importance of including higher order corrections when optimising experimental searches.

The impact of the POWHEG correction was also investigate for two decays in the CMSSM model by studying the transverse momentum distributions of the hardest jet generated by the shower. At values of the transverse momentum less than the upper limit of the POWHEG correction, it was found that the POWHEG corrected distributions were significantly reduced with respect to those generated with the conventional parton shower. Above this cutoff, the normal parton shower and POWHEG corrected distributions were found to be similar.

In this work, we have used the POWHEG formalism to improve the simulation of hard radiation in particle decays and studied the resulting effect on a number of distributions. However, hard radiation in the initial-state parton shower can also have a significant impact on these distributions. Hence, in order to achieve accurate simulation of hard radiation in BSM processes we must also include effects from the initial-state shower. Using the POWHEG formalism to improve the simulation of the hardest initial-state emission in the shower will be the subject of future work.