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

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, \tau \to \nu_{\tau}+\pi, 2\pi, and 3\pi, 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.


Introduction
The recent observation of a standard-model-like Higgs boson with mass around 125 GeV at the LHC [1] has increased our interests in the H → τ + τ − decay mode, which has the modest branching fraction yet to be established [2]. On the other hand, tau leptons have been well known as a particle which can explore the Higgs sector [3,4] as well as supersymmetric models [5,6], through their polarization [7,8,9].
Tauola [10,11,12] 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. 1 The standard Tauola Universal Interface [16] includes the longitudinal spin correlations between pair-produced taus [17], while the extended version also takes into account the transverse spin effects [18,19]. 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 → τ + τ − . a e-mail: tong.li@monash.edu b e-mail: kentarou.mawatari@vub.ac.be c e-mail: junnaka@post.kek.jp 1 The event generators, Herwig++ and Pythia8, have their own τ -decay package [13,14]; TauSpinner for studies on spin effect in tau production was recently reported [15].
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, 2 to simulate polarized τ decays in the FeynRules (FR) [20] and MadGraph5 (MG5) [21] framework. In Section 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 scalartau (stau) decays in scenarios with stau being nearly degenerate in mass with neutralino, where the on-shell tau production is kinematically forbidden. Moreover, in Sections 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. Sect. 6 presents our brief summary.

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 [22], and describe how we implemented them via the FeynRules (FR) [20] and MadGraph5 (MG5) [21] packages. The decay modes which we consider in this article are as the 2π and 3π modes are dominated by the ρ and a 1 vector-meson productions, respectively. The τ -ν τ -π vertex can be introduced by the effective interaction Lagrangian: with the physical constants in Table 1, the chiral projection operator P L , and the constant form factor The effective Lagrangian for the ρ mode is The form factor F ρ (Q 2 ) is parametrized by [23,10] with the Breit-Wigner factor where Q = q 1 + q 2 for τ − → νπ − (q 1 )π 0 (q 2 ). The running width is where the ρ meson line shape factor is We take α = −0.145 in (6) [23]. In practice, we introduce only pions, π ± and π 0 , as new particles and implement the above Lagrangians (2) and (4) into FR, 3 which provides the Ufo (Universal FeynRules Output) model file [24], with f 2 as a constant parameter. The Aloha (Automatic Libraries Of Helicity Amplitudes) program [25] in MG5 reads the model file to create the Helas (HELicity Amplitude Subroutines) library [26] 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 1 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. 4 The effective vertex, or the decay amplitude, for the a 1 mode is where the hadronic current J µ is given by [23,10] with Q = q 1 + q 2 + q 3 for τ − → ν τ π − (q 1 )π − (q 2 )π + (q 3 ).
the corresponding helicity amplitudes. 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. The results from the standard τ decay library Tauola [12] are also displayed for comparison. We note that the difference in the large m π − π 0 invariant mass region for the ρ mode comes from the different parameter choice for m ρ , which is taken as 1.37 GeV in Tauola. Moreover, the fractional energy distributions of the polarized τ decays in Fig. 3 agree with the ones in the collinear limit [3] as well as by Tauola.

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 (τ 1 -χ 0 1 ) coannihilation scenario in the constrained minimal supersymmetric standard model (CMSSM). The neutralinoχ 0 1 is the lightest supersymmetric particle (LSP) and stable by assuming the R-parity conservation. The scalar-tau (stau)τ 1 is the next-to-lightest SUSY particle (NLSP), which exclusively decays into a τ -lepton and a LSPχ 0 1 . The scenario is cosmologically preferred to provide the observed dark matter relic density [27] as well as to solve the 7 Li problem in the standard big-bang nucleosynthesis [28]. The recent global fit analysis for the CMSSM also favours the point [29]. If mτ 1 − mχ0 1 > m τ , there is no obstacle to use Tauola for τ decays, following event generation with τ -leptons as an on-shell particle [30]. For the case of 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 [31], which provides the above mass spectrum with mτ 1 = 230.3 GeV and mχ0 1 = 228.7 GeV [32]. 5 When the two-body decay modeτ 1 →χ 0 1 + τ is closed, the three-body or four-body decays via the off-shell τ (see, e.g. Fig. 4) become dominant and theτ 1 could be longlived [35,36,29]. Figure 5 shows the stau partial decay widths for each τ decay mode against theχ 0 1 mass. Theτ 1 mass and the relevant mixing angles in the neutralino and stau sectors are fixed as in CMSSM40.1.2, providing aτ Rlike stau and a bino-like neutralino. The vertical dashed lines denote the threshold for each decay mode, where a 1 and ρ indicate ∆m = m 3π and m 2π , respectively. As theχ 0 1 mass is approaching theτ 1 mass, i.e. ∆m becomes smaller, the decay width becomes smaller very rapidly. Just below the τ threshold, around the CMSSM40.1. 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 multipion and leptonic modes. Only in the ∆m < m π region the electronic mode becomes dominant. Those agree with the recent study [29], although ρ and a 1 are considered as on-shell particles in their calculations.
The collider signature of the long-livedτ 1 significantly depends on theτ 1 width, i.e. the lifetime; see e.g. Refs. [37,38,39], and as shown in Fig. 5 the lifetime strongly depends on ∆m. For instance, theτ 1 predicted at CMSSM40.1.2, whose lifetime is O(10 −8 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 longlived stau to see the left-right mixing ofτ 1 in this situation.
The scalar partners of left-handed and right-handed taus,τ L andτ R , mix to form two mass eigenstates, and the lighter one is where θτ is the mixing angle. Similarly gauginos,B and W 3 , and neutral Higgsinos,H 0 d andH 0 u , mix to form four mass eigenstates, and the lightest one is The chirality of τ in thẽ τ 1 → τχ 0 1 decay depends on these mixings and determined by the interaction Lagrangian with the chiral-projection operators, P R/L = 1 2 (1 ± γ 5 ). The couplings, a L and a R , are where Y τ = −gm τ / √ 2m W cos β, 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 31 ∼ 0), where the chirality is conserved. Figure 6 shows the energy distributions of the oneprong pion in the π mode (left) and in the ρ mode (right) for the left-handedτ (blue) and right-handedτ (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 righthanded can be still seen especially for the π and ρ modes.

Spin correlations in tau pair decays
Another important application of the FR τ -decay model is spin correlations in tau-pair production [3,4,17,18,19,13,14,15]. 6 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 this section, we present correlations in tau-pair production through a Z boson, a CP -even scalar H, or a CP -odd scalar A, and the subsequent hadronic τ decays. We note that, although we consider the above standard τ -pair production processes in this paper, our τ -decay program is independent of the τ production and can cooperate with all kinds of new physics models.

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 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 the full amplitude can be expressed as the product of the tau-pair production amplitude (M X=Z,H,A ) and two τ → πν decay amplitudes (M 1,2 ): with the τ propagator factor D(q 2 ) = (q 2 − m 2 + imΓ ) −1 .
It is straightforward to obtain the squared matrix elements of the full production plus decay amplitudes, in terms of the production density matrix P λ1λ2 λ1λ2 and the decay density matrices D 1,2 λ1,2 λ1,2 ; In the narrow width limit the propagator factor becomes 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.

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, where β =β  (10), Θ is the scattering angle, and we choose p 1 × q 1 direction as the y-axis. The momenta of the τ − decay products are parametrized in the τ − rest frame, , β 1 sin θ 1 cos φ 1 , β 1 sin θ 1 sin φ 1 , β 1 cos θ 1 ), Fig. 7. Schematic view of the coordinate system.

Tau-pair production and decay amplitudes
The helicity amplitude for the tau-pair production via Zboson (21a) is given by where the reduced amplitude iŝ Here, we neglect the initial fermion mass and take σ = σ 1 = −σ 2 . g q ± and g τ ± are the Z-boson couplings to rightand left-handed quarks and tau leptons, respectively.
The amplitude for the S(= H, A)-boson case (21b) and (21c) is The amplitudes are non-zero only when σ = σ 1 = σ 2 and λ = λ 1 = λ 2 . Here, y τ = √ 2m τ /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.
For the Z production, only for λ 1 = λ 2 the production amplitude is non-zero. After integrating out the scattering angle Θ and the azimuthal angles φ 1,2 , the squared matrix element (25) is with κ = (g τ As seen in (39), the polar angle cos θ 1,2 distribution of the pion arising from τ ± L is (1 − cos θ 1,2 ), while from τ ± R it is (1 + cos θ 1,2 ). Hence in Z → τ + L τ − R or τ + R τ − 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 θ 1 and cos θ 2 in the single differential cross section, and one can see a slightly higher density in the cos θ 1 < 0 and cos θ 2 > 0 region in the left panel of Fig. 8. For the Higgs production, on the other hand, only for λ 1 = λ 2 the production amplitude is non-zero, and therefore Apart from the linear dependence of cos θ 1,2 for the Z case, completely the opposite is favoured for H/A → τ + L τ − L or τ + R τ − 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 .

Polarization correlations
The non-trivial polarization correlations are given by the off-diagonal parts of the density matrices, i.e. λ 1 = −λ 1 and λ 2 = −λ 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 Here, and in the following, summation over repeated indices (λ 1 , λ 2 ,λ 1 ,λ 2 ) = ± is implied. The coefficients F (±) i 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 λ 1 = λ 2 . 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 the coefficients F (±) 1−4 are expressed in terms of the production density matrix and two decay density matrices as 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 φ 1 (or φ 2 ) dependence, however this is very small because the distribution of the azimuthal difference between the two τ -decay planes, φ 1 − φ 2 (≡ ∆φ), is flat. An interesting correlation between the two decay planes for the Z case is After integrating out Θ, θ 1 and θ 2 , the azimuthal asymmetry is given by with The distribution is enhanced around Φ = 0 and 2π, while it is suppressed around Φ = π.
For the H/A production, as mentioned above, only the φ 1 − φ 2 term is non-zero and the coefficient is where the −/+ sign is for CP -even/odd scalar. The azimuthal asymmetry is given by [4] 1 Γ 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 agree well with the above analytic formula.  Fig. 10. The fractional energy correlations in the laboratory frame in pp → X → τ − (→ π − ν) τ + (→ A +ν ) at √ s = 8 TeV for X = Z (top) and H/A (bottom), where A + is the leptonic, π, ρ, and a1 decay modes from the left to right. The resonant mass is taken to be the Z-boson mass.

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.
As a result, the z 1 -z 2 correlation behaves opposite to cos θ 1cos θ 2 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.
Second, an useful observable to measure the parity of Higgs boson is the acollinearity angle, namely the angle δ between π + and π − [4]. The azimuthal asymmetry ∆φ  Fig. 11. The π + π − acollinearity distribution in the Higgs rest frame (left) and in the laboratory frame (right) for pp → 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 → H/A → τ − (→ π − ν) τ + (→ π +ν ) at √ 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 [18].

Spin correlations with TauDecay
To validate our FR τ -decay model file, all the results shown in Sections 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 P λ1λ2 λ1λ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 D i λī λi , according to (25). Figure 12 shows the azimuthal angle ∆φ correlation in pp → H/A → τ − (→ π − ν) τ + (→ π +ν ) for the comparison of the three methods, and the same results are reproduced as expected.

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.