TauDecay: a library to simulate polarized tau decays via FeynRules and MadGraph5

Abstract

TauDecay is a library of helicity amplitudes to simulate polarized tau decays, constructed in the FeynRules and MadGraph5 framework. Together with the leptonic mode, the decay library includes the main hadronic modes, τ→ν _τ+π, 2π, and 3π, which are introduced as effective vertices by using FeynRules. The model file allows us to simulate tau decays when the on-shell tau production is kinematically forbidden. We also demonstrate that all possible correlations among the decay products of pair-produced taus through a Z boson and a scalar/pseudoscalar Higgs boson are produced automatically. The program has been tested carefully by making use of the standard tau decay library Tauola.

1 Introduction

The recent observation of a standard-model-like Higgs boson with mass around 125 GeV at the LHC [1, 2] has increased our interests in the H→τ ⁺ τ ⁻ decay mode, which has the modest branching fraction yet to be established [3, 4]. On the other hand, tau leptons have been well known as a particle which can explore the Higgs sector [5–8] as well as supersymmetric models [9–11], through their polarization [12–14].

Tauola [15–17] is a well-known program to simulate polarized τ decays, including various hadronic modes, and recent theoretical and experimental event simulations involving τ decays heavily depend on it.¹ The standard Tauola Universal Interface [22, 23] includes the longitudinal spin correlations between pair-produced taus [24], while the extended version also takes into account the transverse spin effects [25–28]. We point here some limitations of the use of Tauola:

• (i) The produced τ has to be stable, i.e. on-shell, in event generators.

• (ii) The interface for the spin correlation effects is limited to the standard processes such as Z/γ ^∗→τ ⁺ τ ⁻ and H/A→τ ⁺ τ ⁻.

• (iii) For the transverse spin effects, a particular simulation is required by accepting or rejecting a pair of produced taus based on a spin weighting factor.

In this article, we present a new implementation of a τ-decay model file, and construct a library of helicity amplitudes, TauDecay,² to simulate polarized τ decays in the FeynRules (FR) [29] and MadGraph5 (MG5) [30] framework. In Sect. 2 we describe how we implement the effective vertices via FR for the main hadronic decay modes, namely τ→ν _τ+π, 2π, and 3π, and present the validation of our TauDecay package. To address the limitations mentioned above, in Sect. 3 we study scalar-tau (stau) decays in scenarios with stau being nearly degenerate in mass with neutralino, where the on-shell tau production is kinematically forbidden. Moreover, in Sects. 4 and 5 we demonstrate that all possible correlations among the decay products of pair-produced taus via a Z boson and a scalar/pseudoscalar Higgs boson can be produced within our full-fledged package. Section 6 presents our brief summary.

2 Effective vertices for hadronic tau decays

In this section, we consider effective vertices for the three hadronic τ decay modes, which together with the leptonic mode account for 90 % of the τ decays [31], and describe how we implemented them via the FeynRules (FR) [29] and MadGraph5 (MG5) [30] packages.

The decay modes which we consider in this article are

$$\begin{aligned} &{\pi\ \mathrm{mode}{:}\quad\tau^-\to\nu_{\tau}\pi^-,} \end{aligned}$$

(1a)

$$\begin{aligned} &{\rho\ \mathrm{mode}{:}\quad\tau^-\to\nu_{\tau}\rho^- \to \nu_{\tau}\pi^-\pi^0,} \end{aligned}$$

(1b)

$$\begin{aligned} &{a_1\ \mathrm{mode}{:}\quad\tau^-\to\nu_{\tau}a_1^- \to\nu_{\tau}\pi^0\rho^-\to\nu_{\tau} \pi^0\pi^0\pi^-,} \end{aligned}$$

(1c)

$$\begin{aligned} &{\phantom{a_1\ \mathrm{mode}{:}\quad}\tau^-\to \nu_{\tau}a_1^- \to\nu_{\tau}\pi^-\rho^0 \to\nu_{\tau}\pi^-\pi^-\pi^+,} \end{aligned}$$

(1d)

as the 2π and 3π modes are dominated by the ρ and a ₁ vector-meson productions, respectively.

The τ–ν _τ–π vertex can be introduced by the effective interaction Lagrangian:

$$\begin{aligned} \mathcal{L}_{\pi}&= \sqrt{2}G_F f_1 \bar{ \tau}\gamma^{\mu} P_L \nu_{\tau} \partial_{\mu}\pi^-+\mbox{h.c.} \end{aligned}$$

(2)

with the physical constants in Table 1, the chiral projection operator P _L, and the constant form factor

$$\begin{aligned} f_1=f_{\pi}\cos\theta_C. \end{aligned}$$

(3)

Table 1 Physical constants [31]

The effective Lagrangian for the ρ mode is

$$\begin{aligned} \mathcal{L}_{\rho}&= \sqrt{2}G_Ff_2 \bar{\tau} \gamma^{\mu} P_L \nu_{\tau} \bigl(\pi^0 \partial_{\mu}\pi^- - \pi^-\partial_{\mu}\pi^0 \bigr)+\mbox{h.c.} \end{aligned}$$

(4)

with

$$\begin{aligned} f_2=\sqrt{2}\cos\theta_CF_{\rho} \bigl(Q^2\bigr). \end{aligned}$$

(5)

The form factor F _ρ(Q ²) is parametrized by [15, 32]

$$\begin{aligned} F_{\rho}\bigl(Q^2\bigr)=\bigl[B_{\rho} \bigl(Q^2\bigr)+\alpha B_{\rho'}\bigl(Q^2\bigr) \bigr]/(1+\alpha) \end{aligned}$$

(6)

with the Breit–Wigner factor

$$\begin{aligned} B_{V}\bigl(Q^2\bigr)=\frac{m_{V}^2}{ m_{V}^2-Q^2-i\sqrt{Q^2}\varGamma_{V}(Q^2)}, \end{aligned}$$

(7)

where Q=q ₁+q ₂ for τ ⁻→νπ ⁻(q ₁)π ⁰(q ₂). The running width is

$$\begin{aligned} \varGamma_{V}\bigl(Q^2\bigr)=\varGamma_{V} \frac{\sqrt{Q^2}}{m_{V}} \frac{g_V(Q^2)}{g_V(m_{V}^2)}, \end{aligned}$$

(8)

where the ρ meson line shape factor is

$$\begin{aligned} g_{\rho}\bigl(Q^2\bigr)=\bar{\beta}\biggl(\frac{m_{\pi^-}^2}{Q^2}, \frac{m_{\pi^0}^2}{Q^2} \biggr)^3 \end{aligned}$$

(9)

with

$$\begin{aligned} \bar{\beta}(a,b)\equiv\bigl(1+a^2+b^2-2a-2b-2ab \bigr)^{1/2}. \end{aligned}$$

(10)

We take α=−0.145 in (6) [32].

In practice, we introduce only pions, π ^± and π ⁰, as new particles and implement the above Lagrangians (2) and (4) into FR,³ which provides the Ufo (Universal FeynRules Output) model file [33], with f ₂ as a constant parameter. The Aloha (Automatic Libraries Of Helicity Amplitudes) program [34] in MG5 reads the model file to create the Helas (HELicity Amplitude Subroutines) library [35, 36] for helicity amplitude computations. At this stage we replace the constant parameter by the momentum-dependent form factor (5).

Unlike the above two cases, the effective vertex for the a ₁ mode cannot be obtained by the Lagrangian since the vertex structure is not symmetric between the two identical pion in the final state. Therefore, instead of introducing the Lagrangian, we implement it by hand in FR and created the Ufo model file.⁴ The effective vertex, or the decay amplitude, for the a ₁ mode is

$$\begin{aligned} \mathcal{M}_{a_1}=\sqrt{2}G_F \bar{\tau} \gamma^{\mu}P_L\nu_{\tau}J_{\mu}, \end{aligned}$$

(11)

where the hadronic current J ^μ is given by [15, 32]

$$\begin{aligned} J^{\mu}=f_3 \bigl[F^{13}\bigl(q_1^{\mu}-q_3^{\mu}-G^{13}Q^{\mu} \bigr) + (1\leftrightarrow2) \bigr] \end{aligned}$$

(12)

with Q=q ₁+q ₂+q ₃ for τ ⁻→ν _τ π ⁻(q ₁)π ⁻(q ₂)π ⁺(q ₃). The form factors are

$$\begin{aligned} &{f_3=\frac{4}{3f_{\pi}}\cos\theta_CB_{a_1} \bigl(Q^2\bigr),} \end{aligned}$$

(13)

$$\begin{aligned} &{F^{i3}=F_{\rho}\bigl(Q_{i3}^2 \bigr),\qquad G^{i3}=\frac{Q\cdot(q_i-q_3)}{Q^2},} \end{aligned}$$

(14)

with Q _i3=q _i+q ₃ (i=1,2). The a ₁ meson line shape factor in (8) is parametrized as [15, 32]

$$\begin{aligned} g_{a_1}\bigl(Q^2\bigr)= \begin{cases} \frac{4.1}{Q^2}(Q^2-9m_{\pi}^2)^3 \\ \quad{} \times [1-3.3(Q^2-9m_{\pi}^2)+5.8(Q^2-9m_{\pi}^2)^2 ] \\ \quad\mathrm{if}\ Q^2<(m_{\rho}+m_{\pi})^2, \\ 1.623+\frac{10.38}{Q^2}-\frac{9.32}{Q^4}+\frac{0.65}{Q^6} \\ \quad\mathrm{if}\ Q^2>(m_{\rho}+m_{\pi})^2. \end{cases} \end{aligned}$$

(15)

The π ⁰ π ⁰ π ⁻ mode is the same, except the mass difference between π ^± and π ⁰. The relevant particle masses and widths are collected in Table 2.

Table 2 Particle masses and widths [31]

Our Ufo model now allows MG5 to generate Feynman diagrams including hadronic τ decays (Fig. 1) and the corresponding helicity amplitudes, while the leptonic decay mode can be simulated within the default Standard Model Ufo via the off-shell W boson. TauDecay is a library to simulate polarized τ decays, based on the τ decay helicity amplitudes created by MG5 with the form factor implementation. Using TauDecay, we produce the partial decay widths in Table 3 as well as the pion invariant mass distributions in Fig. 2 in good agreement with the standard τ decay library Tauola [17].⁵ We note that the slight mismatch between the two programs in the table as well as in the large $m_{\pi^{-}\pi^{0}}$ invariant mass region for the ρ mode comes from the QED correction for the leptonic mode in Tauola and the different parameter choice, e.g. f _π and m _ρ′. Moreover, the fractional energy distributions of the polarized τ decays in Fig. 3 agree with the ones in the collinear limit [5–7] as well as by Tauola.

Fig. 1

Effective vertices generated by MG5

Fig. 2

π ⁻ π ⁰ (left) and π ⁰ π ⁰ π ⁻ (right) invariant mass distributions for the τ decay into two and three pions

Fig. 3

The fractional energy distributions of ρ and a ₁ (top) and 1-prong π ⁻ (bottom) from left-handed τ ⁻ (left) and right-handed τ ⁻ (right) at E _τ=50 GeV, normalized to the respective branching ratios. The leptonic decay mode is also shown. The energy fractions are $z\equiv E_{\rho,a_{1}}/E_{\tau}$ (top) and E _π,e/E _τ (bottom) in the laboratory frame

Table 3 τ decay partial widths

3 Stau decays via an off-shell tau

An important application of our FR τ-decay model file is for the case that on-shell τ production is kinematically forbidden. The model file allows MG5 to treat an intermediate τ as a propagator.

Such a case can happen in the so-called stau-neutralino ($\tilde{\tau}_{1}$-$\tilde{\chi}^{0}_{1}$) coannihilation scenario in the constrained minimal supersymmetric standard model (CMSSM). The neutralino $\tilde{\chi}^{0}_{1}$ is the lightest supersymmetric particle (LSP) and stable by assuming the R-parity conservation. The scalar-tau (stau) $\tilde{\tau}_{1}$ is the next-to-lightest SUSY particle (NLSP), which exclusively decays into a τ-lepton and a LSP $\tilde{\chi}^{0}_{1}$. The scenario is cosmologically preferred to provide the observed dark matter relic density [38] as well as to solve the ⁷Li problem in the standard big-bang nucleosynthesis [39]. The recent global fit analysis for the CMSSM also favors the point [40]. If $m_{\tilde{\tau}_{1}}-m_{\tilde{\chi}^{0}_{1}}>m_{\tau}$, there is no obstacle to use Tauola for τ decays, following event generation with τ-leptons as an on-shell particle [41]. For the case of

$$\begin{aligned} \Delta m=m_{\tilde{\tau}_1}-m_{\tilde{\chi}^0_1}<m_{\tau}, \end{aligned}$$

(16)

however, we cannot generate events unless τ decays are taken into account in event generators. One of the common benchmark points for the SUSY searches at the LHC is e.g. CMSSM40.1.2 [42], which provides the above mass spectrum with $m_{\tilde{\tau}_{1}}=230.3$ GeV and $m_{\tilde{\chi}^{0}_{1}}=228.7$ GeV [43].⁶

When the two-body decay mode $\tilde{\tau}_{1}\to\tilde{\chi}^{0}_{1}+\tau$ is closed, the three-body or four-body decays via the off-shell τ (see, e.g. Fig. 4) become dominant and the $\tilde{\tau}_{1}$ could be long-lived [40, 46, 47]. Figure 5 shows the stau partial decay widths for each τ decay mode against the $\tilde{\chi}^{0}_{1}$ mass. The $\tilde{\tau}_{1}$ mass and the relevant mixing angles in the neutralino and stau sectors are fixed as in CMSSM40.1.2, providing a $\tilde{\tau}_{R}$-like stau and a bino-like neutralino. The vertical dashed lines denote the threshold for each decay mode, where a ₁ and ρ indicate Δm=m _3π and m _2π, respectively. As the $\tilde{\chi}^{0}_{1}$ mass is approaching the $\tilde{\tau}_{1}$ mass, i.e. Δm becomes smaller, the decay width becomes smaller very rapidly. Just below the τ threshold, around the CMSSM40.1.2 point (denoted by a vertical solid line), decay widths of all the modes are comparable with those in the on-shell τ case, while the π mode is dominant in the small Δm region due to the phase space suppression for the multi-pion and leptonic modes. Only in the Δm<m _π region the electronic mode becomes dominant. Those agree with the recent study [40], although ρ and a ₁ are considered as on-shell particles in their calculations.

Fig. 4

Feynman diagram of a four-body $\tilde{\tau}_{1}$ decay generated by MG5

Fig. 5

Stau partial decay widths as a function of the LSP neutralino mass. The sum of the partial widths is also shown by a black dashed line. The stau mass is fixed at 230.3 GeV, based on the CMSSM40.1.2 benchmark point. The vertical dashed lines denote the threshold for each decay mode

The collider signature of the long-lived $\tilde{\tau}_{1}$ significantly depends on the $\tilde{\tau}_{1}$ width, i.e. the lifetime; see e.g. Refs. [48–50], and as shown in Fig. 5 the lifetime strongly depends on Δm. For instance, the $\tilde{\tau}_{1}$ predicted at CMSSM40.1.2, whose lifetime is $\mathcal{O}(10^{-8}~\mathrm{s})$, could leave a charged track with a displaced vertex of its decay inside the detectors. We briefly study the decay distributions for such a long-lived stau to see the left-right mixing of $\tilde{\tau}_{1}$ in this situation.

The scalar partners of left-handed and right-handed taus, $\tilde{\tau}_{L}$ and $\tilde{\tau}_{R}$, mix to form two mass eigenstates, and the lighter one is

$$\begin{aligned} \tilde{\tau}_1 = \cos\theta_{\tilde{\tau}}\tilde{\tau}_L + \sin\theta_{\tilde{\tau}}\tilde{\tau}_R, \end{aligned}$$

(17)

where $\theta_{\tilde{\tau}}$ is the mixing angle. Similarly gauginos, $\tilde{B}$ and $\tilde{W}_{3}$, and neutral higgsinos, $\tilde{H}_{d}^{0}$ and $\tilde{H}_{u}^{0}$, mix to form four mass eigenstates, and the lightest one is

$$\begin{aligned} \tilde{X}_i = U_{i1}\tilde{\chi}^0_1, \end{aligned}$$

(18)

where $\tilde{X}=(\tilde{B},\tilde{W}_{3},\tilde{H}_{d}^{0},\tilde{H}_{u}^{0})$. The chirality of τ in the $\tilde{\tau}_{1}\to\tau\tilde{\chi}^{0}_{1}$ decay depends on these mixings and determined by the interaction Lagrangian

$$\begin{aligned} \mathcal{L} = \overline{\tilde{\chi}_1^0}(a_LP_L+a_RP_R) \tau \tilde{\tau}_1^* + \mbox{h.c.} \end{aligned}$$

(19)

with the chiral-projection operators, $P_{R/L}=\frac{1}{2}(1\pm\gamma^{5})$. The couplings, a _L and a _R, are

$$\begin{aligned} &{a_L = \cos\theta_{\tilde{\tau}}\frac{g}{\sqrt{2}} (U_{21}+U_{11}\tan\theta_W) + \sin \theta_{\tilde{\tau}}U_{31}Y_{\tau},} \end{aligned}$$

(20a)

$$\begin{aligned} &{a_R = -\sin\theta_{\tilde{\tau}}\frac{g}{\sqrt{2}}2 U_{11}^*\tan\theta_W + \cos\theta_{\tilde{\tau}}U_{31}^*Y_{\tau},} \end{aligned}$$

(20b)

where $Y_{\tau}=-gm_{\tau}/\sqrt{2}m_{W}\cos\beta$, g is the SU(2)_L gauge coupling, and tanβ is the ratio of the vacuum expectation values of the two Higgs doublets. For most of SUSY scenarios as well as for CMSSM40.1.2 the lightest neutralino is gaugino-like (U ₃₁∼0), where the chirality is conserved.

Figure 6 shows the energy distributions of the one-prong pion in the π mode (left) and in the ρ mode (right) for the left-handed $\tilde{\tau}$ (blue) and right-handed $\tilde{\tau}$ (red) decay in the stau rest frame. The distributions for the reconstructed ρ is also shown by dashed lines as a reference. The phase space of the decaying off-shell τ is limited and it is no longer energetic. Therefore the energy distributions of the pion are quite different from the ones in Fig. 3, where the collinear limit is safely applied, although the remnant of the difference between left-handed and right-handed can be still seen especially for the π and ρ modes.

Fig. 6

Energy distributions of the one-prong pion in the π mode (left) and in the ρ mode (right) for the left-handed $\tilde{\tau}$ (blue) and right-handed $\tilde{\tau}$ (red) decay in the stau rest frame. The distributions for the intermediate ρ is also shown by dashed lines as a reference

4 Spin correlations in tau pair decays

Another important application of the FR τ-decay model is spin correlations in tau-pair production [5–8, 18–21, 24–28].⁷ As in the previous section, one can generate amplitudes including the τ decays, and hence the spin correlations between the two taus are automatically taken into account. In other words, there is no need for a particular simulation to get the transverse spin effects.

In this section, we present correlations in tau-pair production through a Z boson, a CP-even scalar H, or a CP-odd scalar A,

$$\begin{aligned} &{q\bar{q} \to Z\to\tau^+\tau^-,} \end{aligned}$$

(21a)

$$\begin{aligned} &{gg \to H\to\tau^+\tau^-,} \end{aligned}$$

(21b)

$$\begin{aligned} &{gg \to A\to\tau^+\tau^-,} \end{aligned}$$

(21c)

and the subsequent hadronic τ decays. We note that, although we consider the above standard τ-pair production processes in this paper, a special interface to cooperate with any new physics models is not required when we generate amplitudes with intermediate taus being as propagators.

4.1 Helicity formalism

First, we discuss the case that both taus decay to ν _τ π by using analytic expressions. We assign the four-momentum and the helicity of each particle for the process as

$$\begin{aligned} &{a_1(p_1,\sigma_1)+a_2(p_2, \sigma_2)} \\ &{\quad\to\tau^-(q_1,\lambda_1)+ \tau^+(q_2,\lambda_2) } \\ &{\quad\to\pi^-(k_1)+\pi^+(k_2)+ \nu_{\tau}(k_3,-)+\bar{\nu}_{\tau}(k_4,+),} \end{aligned}$$

(22)

where a _1,2 stand for quarks or gluons (or even electrons or photons in case of e ⁺ e ⁻ or γγ collisions). The helicities take the values σ _i/2 and λ _i/2 for quarks and leptons and σ _i for gluons and photons with σ _i=±1 and λ _i=±1. Although we evaluate the full amplitude for the above processes numerically to present the distributions, we discuss it below in terms of the τ-pair production amplitude and the τ decay amplitudes, which give us better understanding of the distributions. Moreover, we construct our TauDecay package based on the following analytic expressions in the end of this section.

Using the completeness relations

(23a)

(23b)

the full amplitude can be expressed as the product of the tau-pair production amplitude ($\mathcal{M}_{X=Z,H,A}$) and two τ→πν decay amplitudes ($\mathcal{M}_{1,2}$):

$$\begin{aligned} &{\mathcal{M}(p_1,\sigma_1;p_2, \sigma_2;k_i) } \\ &{\quad=D\bigl(q_1^2\bigr)D\bigl(q_2^2 \bigr)\sum_{\lambda_{1,2}} \mathcal{M}_X(p_1, \sigma_1;p_2,\sigma_2;q_1, \lambda_1;q_2,\lambda _2) } \\ &{\qquad{}\times\mathcal{M}_1(q_1, \lambda_1;k_1;k_3) \mathcal{M}_2(q_2, \lambda_2;k_2;k_4)} \end{aligned}$$

(24)

with the τ propagator factor D(q ²)=(q ²−m ²+imΓ)⁻¹.

It is straightforward to obtain the squared matrix elements of the full production plus decay amplitudes,

$$\begin{aligned} \sum|\mathcal{M}|^2 &\equiv \sum _{\sigma_{1,2}} \bigl|\mathcal{M}(p_1,\sigma_1;p_2, \sigma_2;k_i) \bigr|^2 \\ &= \bigl|D\bigl(q_1^2\bigr)D\bigl(q_2^2 \bigr) \bigr|^2 \sum_{\lambda_1,\lambda_2}\sum _{\bar{\lambda}_1,\bar{\lambda}_2} {\mathcal{P}}^{\lambda_1\lambda_2}_{\bar{\lambda}_1\bar{\lambda}_2} { \mathcal{D}_1}^{\lambda_1}_{\bar{\lambda}_1}{\mathcal {D}_2}^{\lambda_2}_{\bar{\lambda}_2} \end{aligned}$$

(25)

in terms of the production density matrix ${\mathcal{P}}^{\lambda_{1}\lambda_{2}}_{\bar{\lambda}_{1}\bar{\lambda}_{2}}$ and the decay density matrices ${\mathcal{D}_{1,2}}^{\lambda_{1,2}}_{\bar{\lambda}_{1,2}}$;

$$\begin{aligned} &{{\mathcal{P}}^{\lambda_1\lambda_2}_{\bar{\lambda}_1\bar{\lambda}_2}=\sum _{\sigma} \mathcal{M}^{\lambda_1\lambda_2}_{\sigma} \bigl( \mathcal{M}^{\bar{\lambda}_1\bar{\lambda}_2}_{\sigma} \bigr)^*, } \end{aligned}$$

(26)

$$\begin{aligned} &{{\mathcal{D}_1}^{\lambda_1}_{\bar{\lambda}_1}= \mathcal{M}_{\lambda_1} ( \mathcal{M}_{\bar{\lambda}_1} )^*, } \end{aligned}$$

(27)

$$\begin{aligned} &{{\mathcal{D}_2}^{\lambda_2}_{\bar{\lambda}_2}= \mathcal{M}_{\lambda_2} ( \mathcal{M}_{\bar{\lambda}_2} )^*.} \end{aligned}$$

(28)

In the narrow width limit the propagator factor becomes

$$\begin{aligned} \bigl \vert D\bigl(q^2\bigr)\bigr \vert ^2\to \frac{\pi}{m\varGamma}\delta\bigl(q^2-m^2\bigr). \end{aligned}$$

(29)

Although we only consider parton-level subprocesses, one can generalize (25) to mixed case and apply it for any processes, including summation over different subprocesses and a product of the relevant parton densities. We can also easily replace the τ→ν _τ π decay amplitude by one for the other hadronic decay mode as well as the leptonic mode.

4.2 Kinematics

Let us define the kinematical variables. In the collision center-of-mass (CM) frame, i.e. in the X rest frame, we choose the τ momentum direction as the z-axis,

$$\begin{aligned} \begin{aligned} &{p_1=\frac{\sqrt{\hat{s}}}{2}(1,-\sin \varTheta,0,\cos\varTheta ),} \\ &{p_2=\frac{\sqrt{\hat{s}}}{2}(1,\sin\varTheta,0,-\cos\varTheta ),} \\ &{q_1=\frac{\sqrt{\hat{s}}}{2} \biggl(1+\frac{q_1^2-q_2^2}{\hat{s}}, 0,0, \beta \biggr),} \\ &{q_2=\frac{\sqrt{\hat{s}}}{2} \biggl(1+\frac{q_2^2-q_1^2}{\hat{s}}, 0,0,- \beta \biggr),} \end{aligned} \end{aligned}$$

(30)

where $\beta=\bar{\beta}(\frac{q_{1}^{2}}{\hat{s}},\frac{q_{2}^{2}}{\hat{s}} )$ with $\bar{\beta}(a,b)$ defined in (10), Θ is the scattering angle, and we choose p ₁×q ₁ direction as the y-axis. The momenta of the τ ⁻ decay products are parametrized in the τ ⁻ rest frame,

$$\begin{aligned} \begin{aligned} &{k_1=\frac{\sqrt{q_1^2}}{2} \biggl(1+\frac{m_{\pi}^2}{q_1^2},\beta_1\sin\theta_1\cos \phi_1, \beta_1\sin\theta_1\sin \phi_1,\beta_1\cos\theta_1\biggr), } \\ &{k_3=\frac{\sqrt{q_1^2}}{2}\beta_1 (1,-\sin \theta_1\cos\phi_1, -\sin\theta_1\sin \phi_1, -\cos\theta_1),} \end{aligned} \end{aligned}$$

(31)

with $\beta_{1}=1-m_{\pi}^{2}/q_{1}^{2}$. Similarly, those of the τ ⁺ decay products are

$$\begin{aligned} \begin{aligned} &{k_2=\frac{\sqrt{q_2^2}}{2} \biggl(1+\frac{m_{\pi}^2}{q_2^2},\beta_2\sin\theta_2\cos \phi_2, \beta_2\sin\theta_2\sin \phi_2,\beta_2\cos\theta_2\biggr), } \\ &{k_4=\frac{\sqrt{q_2^2}}{2}\beta_2 (1,-\sin \theta_2\cos\phi_2, -\sin\theta_2\sin \phi_2, -\cos\theta_2),} \end{aligned} \end{aligned}$$

(32)

with $\beta_{2}=1-m_{\pi}^{2}/q_{2}^{2}$. We note that the polar (z-)axis and the y-axis normal to the scattering plane are chosen common to all the three frames, and the two decay frames differ only by the boost along the τ production axis (see Fig. 7). The τ width is very narrow, $\varGamma\sim\mathcal{O}(10^{-12}~\mathrm{GeV})$, and hence we take the narrow width limit, $q_{1}^{2}=q_{2}^{2}=m^{2}$, in the following analytic amplitudes.

Fig. 7

Schematic view of the coordinate system

4.3 Tau-pair production and decay amplitudes

The helicity amplitude for the tau-pair production via Z-boson (21a) is given by

$$\begin{aligned} \mathcal{M}_Z=\frac{e^2\hat{s}}{\hat{s}-m_Z^2+im_Z\varGamma_Z} \hat{\mathcal{M}}_{\sigma}^{\lambda_1\lambda_2}, \end{aligned}$$

(33)

where the reduced amplitude is

$$\begin{aligned} \hat{\mathcal{M}}_{\sigma}^{\lambda_1\lambda_2} = \begin{cases} -\frac{1}{2} g_{\sigma}^q\{g_+^{\tau}(1+\lambda_1\beta)+g_-^{\tau}(1-\lambda _1\beta)\}\\ \quad{}\times(\sigma\lambda_1+\cos\varTheta) \quad\text{for }\lambda_1\ne\lambda_2, \\ \frac{m}{\sqrt{\hat{s}}} \lambda_1g_{\sigma}^q(g_+^{\tau}+g_-^{\tau})\sin\varTheta \quad\text{for }\lambda_1=\lambda_2. \end{cases} \end{aligned}$$

(34)

Here, we neglect the initial fermion mass and take σ=σ ₁=−σ ₂. $g_{\pm}^{q}$ and $g_{\pm}^{\tau}$ are the Z-boson couplings to right- and left-handed quarks and tau leptons, respectively.

The amplitude for the S(=H,A)-boson case (21b) and (21c) is

$$\begin{aligned} \mathcal{M}_S=\frac{-y_{\tau}g_{S}\hat{s}}{\hat{s}-m_{S}^2+im_{S}\varGamma_{S}} \frac{\sqrt{\hat{s}}}{2} \hat{ \mathcal{M}}_{\sigma}^{\lambda _1\lambda_2}, \end{aligned}$$

(35)

where

$$\begin{aligned} &{\hat{\mathcal{M}}_{\sigma}^{\lambda\lambda} = \beta \lambda\quad \text{for }S=H,} \end{aligned}$$

(36)

$$\begin{aligned} &{\hat{\mathcal{M}}_{\sigma}^{\lambda\lambda} = i \quad \text{for }S=A.} \end{aligned}$$

(37)

The amplitudes are non-zero only when σ=σ ₁=σ ₂ and λ=λ ₁=λ ₂. Here, $y_{\tau}=\sqrt{2}m_{\tau }/v$ is the τ Yukawa coupling, and g _S is the ggS coupling, which is given by g _H=α _s g _Htt/3πv and g _A=α _s g _Att/2πv in the heavy-top limit.

The τ→πν decay amplitudes in the τ rest frame are

$$\begin{aligned} \mathcal{M}_{1,2}=-G_Ff_{\pi}m^2 \sqrt{\beta_{1,2}}\hat{\mathcal{M}}_{\lambda_{1,2}}, \end{aligned}$$

(38)

where

$$\begin{aligned} &{\hat{\mathcal{M}}_{\lambda_1}=\sqrt{1+\lambda_1\cos \theta _1}e^{i\lambda_1\phi_1/2},} \end{aligned}$$

(39a)

$$\begin{aligned} &{\hat{\mathcal{M}}_{\lambda_2}=\sqrt{1+\lambda_2\cos \theta _2}e^{-i\lambda_2\phi_2/2}.} \end{aligned}$$

(39b)

Those amplitudes are invariant with the boost to the collision CM frame.

4.4 Helicity correlations

The trivial helicity correlations are given by the diagonal parts of the density matrices, $\lambda_{1}=\bar{\lambda}_{1}$ and $\lambda _{2}=\bar{\lambda}_{2}$. In the following analytic expressions we take $m_{\tau}/\sqrt{\hat{s}}=0$, i.e. β=1.

For the Z production, only for λ ₁≠λ ₂ the production amplitude is non-zero. After integrating out the scattering angle Θ and the azimuthal angles ϕ _1,2, the squared matrix element (25) is

$$\begin{aligned} &{{\mathcal{P}}^{+-}_{+-}{\mathcal{D}_1}^+_+{ \mathcal{D}_2}^-_- +{\mathcal{P}}^{-+}_{-+}{ \mathcal{D}_1}^-_-{\mathcal{D}_2}^+_+} \\ &{\quad\propto(1-\cos\theta_1\cos\theta_2) + \kappa(\cos\theta_1-\cos\theta_2).} \end{aligned}$$

(40)

with $\kappa=({g_{+}^{\tau}}^{2}-{g_{-}^{\tau}}^{2})/({g_{+}^{\tau}}^{2}+{g_{-}^{\tau }}^{2})\sim -0.15$. As seen in (39a), (39b), the polar angle cosθ _1,2 distribution of the pion arising from $\tau^{\pm}_{L}$ is (1−cosθ _1,2), while from $\tau^{\pm}_{R}$ it is (1+cosθ _1,2). Hence in $Z\to\tau^{+}_{L}\tau^{-}_{R}$ or $\tau ^{+}_{R}\tau^{-}_{L}$ decays the two pions tend to be emitted to the opposite direction along the z-axis. The difference between the left and right coupling, i.e. parity violation, gives a small linear dependence of cosθ ₁ and cosθ ₂ in the single differential cross section, and one can see a slightly higher density in the cosθ ₁<0 and cosθ ₂>0 region in the left panel of Fig. 8.

Fig. 8

cosθ ₁-cosθ ₂ correlation in $pp\to X\to \tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ for X=Z (left) and H/A (right). The pion polar angles cosθ _1,2 are defined in the τ rest frame; see Fig. 7

For the Higgs production, on the other hand, only for λ ₁=λ ₂ the production amplitude is non-zero, and therefore

$$\begin{aligned} {\mathcal{P}}^{++}_{++}{\mathcal{D}_1}^+_+{ \mathcal{D}_2}^+_+ +{\mathcal{P}}^{--}_{--}{ \mathcal{D}_1}^-_-{\mathcal{D}_2}^-_- \propto(1+\cos \theta_1\cos\theta_2). \end{aligned}$$

(41)

Apart from the linear dependence of cosθ _1,2 for the Z case, completely the opposite is favored for $H/A\to\tau^{+}_{L}\tau^{-}_{L}$ or $\tau^{+}_{R}\tau^{-}_{R}$ decays as shown in the right panel of Fig. 8. Note that the helicity correlations for H and A are identical. For comparison the masses of H and A are assumed to be the Z-boson mass, m _H/A=m _Z.

4.5 Polarization correlations

The non-trivial polarization correlations are given by the off-diagonal parts of the density matrices, i.e. $\lambda_{1}=-\bar{\lambda}_{1}$ and $\lambda_{2}=-\bar{\lambda}_{2}$, which produce the azimuthal angle dependence. When we isolate the azimuthal angle dependence in (25), there are nine distributions (including one constant piece) as

$$\begin{aligned} &{{\mathcal{P}}^{\lambda_1\lambda_2}_{\bar{\lambda}_1\bar{\lambda}_2} {\mathcal{D}_1}^{\lambda_1}_{\bar{\lambda}_1}{ \mathcal {D}_2}^{\lambda_2}_{\bar{\lambda}_2} } \\ &{\quad= F_1+ \bigl\{2{\Re}e \bigl[F_2\cos\phi _1+F_3 \cos\phi_2} \\ &{\qquad{}+F_4^{+}\cos(\phi_1+ \phi_2) +F_4^{-}\cos(\phi_1- \phi_2) \bigr]} \\ &{\qquad{}+({\Re}\mathrm{e}\to{\Im}\mathrm{m},\cos\to\sin) \bigr\}.} \end{aligned}$$

(42)

Here, and in the following, summation over repeated indices $(\lambda_{1},\lambda_{2},\bar{\lambda}_{1},\bar{\lambda}_{2})=\pm$ is implied. The coefficients $F_{i}^{(\pm)}$ are the functions of the kinematical variables except the azimuthal angles ϕ _1,2. For the production of a Z boson, they also depend on the production angle Θ. For the spin-0 particle case, only the $F_{4}^{-}$ term in (42) survives due to the helicity selection λ ₁=λ ₂. All the sine terms vanish when CP is conserved and when the absorptive part of the amplitudes are neglected, e.g., in the tree-level approximation.

Because the phase of the product of the two decay density matrices is

$$\begin{aligned} {\mathcal{D}_1}^{\lambda_1}_{\bar{\lambda}_1}{\mathcal {D}_2}^{\lambda_2}_{\bar{\lambda}_2}\propto \exp\bigl[i\bigl\{( \lambda_1-\bar{\lambda}_1)\phi_1-( \lambda_2-\bar{\lambda}_2)\phi_2\bigr\}/2 \bigr], \end{aligned}$$

(43)

the coefficients $F_{1-4}^{(\pm)}$ are expressed in terms of the production density matrix and two decay density matrices as

$$\begin{aligned} \begin{aligned} &{F_1= \mathcal{P}^{\lambda_1\lambda_2}_{\lambda_1\lambda_2} {\mathcal{D}_1}^{\lambda_1}_{\lambda_1}{ \mathcal{D}_2}^{\lambda _2}_{\lambda_2},} \\ &{F_2=\mathcal{P}^{+\lambda_2}_{-\lambda_2} {\mathcal{D}_1}^{+}_{-}{\mathcal{D}_2}^{\lambda_2}_{\lambda_2},} \\ &{F_3=\mathcal{P}^{\lambda_1+}_{\lambda_1-} {\mathcal{D}_1}^{\lambda_1}_{\lambda_1}{\mathcal{D}_2}^{+}_{-},} \\ &{F_4^{\pm}=\mathcal{P}^{+\mp}_{-\pm} {\mathcal{D}_1}^{+}_{-}{\mathcal{D}_2}^{\mp}_{\pm}.} \end{aligned} \end{aligned}$$

(44)

The azimuthal angle correlations are manifestly expressed by quantum interference among different helicity states of the intermediate tau leptons.

The Z production could produce the ϕ ₁ (or ϕ ₂) dependence, however this is very small because F ₂ (F ₃) is proportional to $m/\sqrt{s}$. Similarly, $F_{4}^{-}\propto\mathcal{O}(m^{2}/\hat{s})$, and hence the distribution of the azimuthal difference between the two τ-decay planes, ϕ ₁−ϕ ₂(≡Δϕ), is flat. An interesting correlation between the two decay planes for the Z case is the ϕ ₁+ϕ ₂(≡Φ) correlation, whose coefficient $F_{4}^{+}$ is

$$\begin{aligned} F_4^+ \propto-g_+^{\tau}g_-^{\tau} \sin^2\varTheta\sin\theta_1\sin\theta_2. \end{aligned}$$

(45)

After integrating out Θ, θ ₁ and θ ₂, the azimuthal asymmetry is given by

$$\begin{aligned} \frac{1}{\varGamma}\frac{d\varGamma}{d\varPhi} =\frac{1}{2\pi} \bigl[1+A_4^+ \cos\varPhi \bigr] \end{aligned}$$

(46)

with

$$\begin{aligned} A_4^+ \equiv\frac{2F_4^+}{F_1} =\frac{\pi^2}{16} \frac{-g_+^{\tau}g_-^{\tau}}{{g_+^{\tau}}^2+{g_-^{\tau}}^2} \sim0.30. \end{aligned}$$

(47)

The distribution is enhanced around Φ=0 and 2π, while it is suppressed around Φ=π.

For the H/A production, as mentioned above, only the ϕ ₁−ϕ ₂ term is non-zero and the coefficient is

$$\begin{aligned} F_4^- \propto\mp\sin\theta_1\sin\theta_2, \end{aligned}$$

(48)

where the −/+ sign is for CP-even/odd scalar. The azimuthal asymmetry is given by [8]

$$\begin{aligned} \frac{1}{\varGamma}\frac{d\varGamma}{d\Delta\phi} =\frac{1}{2\pi} \bigl[1+A_4^- \cos\Delta\phi \bigr] \end{aligned}$$

(49)

with

$$\begin{aligned} A_4^- \equiv\frac{2F_4^-}{F_1} =\mp\frac{\pi^2}{16} \sim\mp0.62. \end{aligned}$$

(50)

The CP-even and -odd scalars have opposite modulation, and the distribution is strongly enhanced (suppressed) around Δϕ=π for the CP-even (-odd) scalars.

To examine the validity of the model file, Fig. 9 shows the normalized azimuthal correlations between π ⁻ and π ⁺, and our numerical results (solid histograms) agree well with the above analytic formula (dotted lines).

Fig. 9

Azimuthal angle correlations, ϕ ₁+ϕ ₂ (left) and ϕ ₁−ϕ ₂ (right), in $pp\to X\to \tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ for X=Z (blue), H (black), and A (red). The pion azimuthal angles ϕ _1,2 are defined in the τ rest frame by the pp→τ ⁺ τ ⁻ scattering plane; see Fig. 7

4.6 Spin correlations at the LHC

To present the validation of our program, we have so far discussed the angular distributions, assuming that all the τ-decay products can be reconstructed. Here, we show some distributions of the realistic experimental observables at the LHC.

First, we extend the helicity correlation in τ→ν _τ π discussed above to other tau decay modes. In Fig. 10 we present the helicity correlations in terms of energy fractions in the laboratory frame for $pp\to Z/H/A\to\tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to A^{+}\bar{\nu})$, where A ⁺ is the leptonic, π, ρ, and a ₁ decay modes, respectively, from left to right. The one-prong τ decays are considered and the energy fractions are $z_{1}=E_{\pi^{-}}/E_{\tau^{-}}$ and $z_{2}=E_{\hat{A}^{+}}/E_{\tau^{+}}$ with $\hat{A}=e,\pi,\pi(\rho),\pi(a_{1})$. For comparison, the resonant mass is taken to the Z-boson mass, and hence the both taus are highly boosted, where the collinear approximation can be safely applied. For the A=π mode, in the collinear limit (β=1) we have

$$\begin{aligned} z_1=(1+\cos\theta_1)/2\quad\text{and}\quad z_2=(1-\cos\theta_2)/2. \end{aligned}$$

(51)

As a result, the z ₁-z ₂ correlation behaves opposite to cosθ ₁-cosθ ₂ correlation as shown in Fig. 8. The helicity correlations are different between the spin-1 and spin-0 bosons, while those are identical between the CP-even and -odd scalars.

Fig. 10

The fractional energy correlations in the laboratory frame in $pp\to X\to\tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to A^{+}\bar{\nu})$ at $\sqrt{s}=8$ TeV for X=Z (top) and H/A (bottom), where A ⁺ is the leptonic, π, ρ, and a ₁ decay modes from the left to right. The resonant mass is taken to be the Z-boson mass

Second, an useful observable to measure the parity of Higgs boson is the acollinearity angle, namely the angle δ between π ⁺ and π ⁻ [8]. The azimuthal asymmetry Δϕ can be given as a function of angle δ. To display the measurable difference of Higgs parities at the LHC, in Fig. 11, we show the π ⁺ π ⁻ acollinearity distribution in the Higgs rest frame (left) and in the lab frame (right) for $pp\to H/A\to\tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ at $\sqrt{s}=8$ TeV. Here we take the Higgs mass as 125 GeV. One can see the distinguishable difference between CP-even and CP-odd scalars around δ=π, although in the lab frame the boost effect leads to relatively unclear discrepancy. While the reconstruction of the Higgs rest frame at hadron colliers is difficult, it is possible in the Z-associated Higgs production at e ⁺ e ⁻ colliders and the acollinearity distribution in the reconstructed H/A rest frame has been discussed in [25].

Fig. 11

The π ⁺ π ⁻ acollinearity distribution in the Higgs rest frame (left) and in the laboratory frame (right) for $pp\to H/A\to\tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ with m _H/A=125 GeV at $\sqrt{s}=8$ TeV

5 Spin correlations with TauDecay

To validate our FR τ-decay model file, all the results shown in Sects. 3 and 4 are generated in MG5 with intermediate taus being as propagators, denoted by method A. Unless intermediate taus are off-shell such as the case in Sect. 3, a library which collects all the possible τ decay channels would be much more efficient for practical event generations since τ production and its decay can be simulated independently as we use Tauola.

TauDecay is a Fortran library of τ-decay helicity amplitudes, constructed in the FR/MG5 framework as described in Sect. 2. The input of the subroutines for each decay channel is the four momenta of the decaying τ as well as its decay products and their helicities. The output is the amplitude, i.e. a complex number, at a given phase space point with a given helicity configuration. In this section, we present how the TauDecay package can reproduce spin correlation discussed in Sect. 4.

We show two methods by means of TauDecay; without and with spin density matrix, denoted by method B and C, respectively. The both methods are carefully checked by method A.

B::

(without density matrix) The τ production amplitude as well as the decay amplitudes from TauDecay are evaluated for a given process for a given phase space point. The τ helicities in the product of the production and decay amplitudes should be summed over as indicated in (24), and then the τ helicity summed amplitudes should be squared. This is theoretically identical to the propagator method in the narrow width limit, but it is useful for systematically studying different tau decay modes for each production process.

C::

(with density matrix) First, the τ production density matrix is computed for each phase space point. If there is only one tau in the final state, this will be a 2×2 matrix, while it will be a 4×4 double density matrix for two taus as ${\mathcal{P}}^{\lambda_{1}\lambda_{2}}_{\bar{\lambda}_{1}\bar{\lambda}_{2}}$ in (26). Then, the correlated tau decays are simulated by multiplying the production density matrix with the corresponding τ decay helicity amplitudes and their complex conjugates provided by TauDecay, namely the decay density matrices ${\mathcal{D}_{i}}^{\lambda_{i}}_{\bar{\lambda}_{i}}$, according to (25).

Figure 12 shows the azimuthal angle Δϕ correlation in $pp\to H/A\to\tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ for the comparison of the three methods, and the same results are reproduced as expected.

Fig. 12

Azimuthal angle Δϕ correlation in $pp\to H/A\to \tau^{-}(\to\pi^{-}\nu)\tau^{+}(\to\pi^{+}\bar{\nu})$ by three different methods; A. propagator (solid), B. without density matrix (dashed), and C. with density matrix (dotted)

6 Summary

In this paper, we have implemented the main hadronic tau decays, τ→ν _τ+π, 2π, and 3π, by using the effective vertices in the FeynRules and MadGraph5 framework, and constructed TauDecay, a library of helicity amplitudes to simulate polarized tau decays. The model file allows us to simulate tau decays when the on-shell tau production is kinematically forbidden, and the stau decay in the stau-neutralino coannihilation region was discussed as an application. We also demonstrated that all possible correlations among the decay products of pair-produced taus through a Z-boson and a scalar/pseudoscalar Higgs boson can be produced within our full-fledged package. The program has been tested carefully by making use of the standard tau decay library Tauola.