Exploring CP violation in $\varvec{H \rightarrow \tau ^+\,\tau ^-\,\gamma }$

Abstract

We propose a method of measuring the CP-odd part of the Yukawa interaction of Higgs boson and $\tau$ leptons by observing the forward–backward asymmetry in the decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$. The source of such asymmetry is the interference of the CP-even loop-level contribution coming from $H \rightarrow Z \, \gamma \rightarrow \tau ^+ \, \tau ^- \, \gamma$ decay channel with the contribution from tree-level CP-odd Yukawa interaction. We find that the CP violating effect is maximum when the invariant mass of the $\tau ^+\,\tau ^-$ pair is equal to the mass of the Z boson. We propose and utilise various Dalitz plot asymmetries to quantify the maximal size of the asymmetry and perform Monte Carlo simulations to study the feasibility of measuring it in the high luminosity phase of the Large Hadron Collider (HL-LHC).

1 Introduction

In the Standard Model (SM), violation of the CP symmetry is encoded in the CKM matrix. In principle, a Beyond the Standard Model (BSM) physics may have new sources of CP violation. In particular, BSM CP violation in the Yukawa interactions is welcome for electroweak baryogenesis (it is well known that CP violation in the SM is by far too weak for baryogenesis [1,2,3]). The most general expression for CP violating $H\psi \psi$ Yukawa interaction can be written in the following form,

$$\begin{aligned} {\mathscr {L}}_{H\psi \psi }^{} = - \frac{m_\psi }{v} \, \overline{\psi } \left( a_\psi + i\, \gamma ^5 \, b_\psi \right) \psi \, H , \end{aligned}$$

(1.1)

where $v$ is the vacuum expectation value of the Higgs field, $m_\psi$ denotes the mass of the fermion $\psi$, and $a_\psi , b_\psi$ are two real valued parameters. In the SM, $a_\psi ^\text {SM}=1$, $b_\psi ^\text {SM}=0$. If simultaneously both $a_\psi \ne 0$ and $b_\psi \ne 0$, it implies CP violation in $H\psi \psi$ Yukawa interaction. The parameters $b_\psi$ are strongly constrained by the experimental bounds on the electron and neutron Electric Dipole Moments (for a recent analysis see [4,5,6,7] and references therein). In this context, the $\tau$ lepton Yukawa coupling is of interest as it is large and the EDM bound, $b_\tau <0.3$, is weak enough for the $\tau$ Yukawa to play a role in electroweak baryogenesis [4,5,6,7], see however [8]. CP violation in the $\tau$ Yukawa has also been searched for at the LHC. The recent study by CMS [9] probing $H \rightarrow \tau ^+ \, \tau ^-$ gives $\left| b_\tau \right| \lesssim 0.34$ at $68.3\%$ confidence level (for further prospects see [10]). Majority of the experimental studies on this issue concentrate on measurements of the angle between $\tau$ decay planes determined by the directions of particles produced in subsequent $\tau$ lepton decays, such as in $H\rightarrow \tau ^+ \, \tau ^-\rightarrow \pi ^+ \, \pi ^- \, \nu _\tau \, \bar{\nu }_\tau$ or $H\rightarrow \tau ^+ \, \tau ^- \rightarrow \rho ^+ \, \rho ^- \, \nu _\tau \, \bar{\nu }_\tau$ [9, 11,12,13,14,15,16,17,18].

Fig. 1

The dominant Feynman diagrams contributing to $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$. The contributions of possible box diagrams at one-loop level (not shown here) are negligible, see Refs. [19, 20]

In this paper we propose to measure the forward–backward asymmetry of $\tau$ lepton angular distribution in the $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ decay, as a measure of the CP violation in $H\tau \tau$ Yukawa interaction. To study this asymmetry we utilise the Lorentz invariant Dalitz plot distribution of events. The dominant CP-violating effects which contribute to the forward–backward asymmetry in the $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ decay are proportional to the interference of the tree-level and loop-level diagrams¹ shown in Fig. 1. The lower branching ratio than the $H \rightarrow \tau ^+ \, \tau ^-$ is partially compensated by the fact that one only requires to reconstruct the 4-momenta of the $\tau$ leptons and not the full spatial distributions of the final $\tau$ decay products. Our heuristic simulations for the HL-LHC show that one can possibly probe $b_\tau$ using our proposed methodology. A more thorough Monte Carlo study scanning the full 2-dimensional Dalitz plot distribution is beyond our current expertise, and is hence reserved for future exploration.

Our paper is organised as follows. In Sect. 2 we briefly outline the important phenomenological aspects of the 3-body decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$, showing how the forward–backward asymmetry originates and how can it be probed from the Lorentz invariant Dalitz plot distribution. In Sect. 3 we do a numerical study, looking at the distribution pattern inside the Dalitz plot and assess how large the forward–backward asymmetry could be. In Sect. 4 we perform a heuristic Monte Carlo study of the feasibility of observing the asymmetry in context of HL-LHC. Finally we conclude in Sect. 5 summarising our findings and highlighting the salient features of our proposed methodology.

2 Phenomenological study of $\varvec{H \rightarrow \tau ^+ \, \tau ^- \, \gamma }$

Fig. 2

Kinematic configurations related by CP in the center-of-momentum frame of $\tau ^+ \, \tau ^-$

The decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ is its own CP-conjugate process. Let us study the kinematic configuration of the decay in the center-of-momentum frame of $\tau ^+ \, \tau ^-$ (equivalently called the di-tau rest frame). From Fig. 2 it is clear that the CP transformation takes the angle $\theta$ between $\tau ^+$ and photon to $\pi -\theta$. This implies that any difference (or asymmetry) in the angular distribution of events with respect to $\cos \theta \leftrightarrow -\cos \theta$ (‘forward’ $\leftrightarrow$ ‘backward’) exchange would be a clear signature of CP-violation.

As illustrated in Fig. 1 the decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ proceeds via the tree-level $H\tau \tau$ Yukawa interaction, as well as via the effective vertex of $H \rightarrow {\mathcal {V}}\,\gamma \rightarrow \tau ^+ \, \tau ^- \, \gamma$, with ${\mathcal {V}}=Z,\gamma$. The effective Lagrangian for the later interaction can be, to the lowest mass dimension order, written in the form,

$$\begin{aligned} {\mathscr {L}}_{H{\mathcal {V}}\gamma }&= \frac{H}{4\,v} \Big (2 A_2^{Z\gamma } F^{\mu \nu } Z_{\mu \nu } + 2 A_3^{Z\gamma } F^{\mu \nu } \widetilde{Z}_{\mu \nu } \nonumber \\&\qquad + A_2^{\gamma \gamma } F^{\mu \nu } F_{\mu \nu } + A_3^{\gamma \gamma } F^{\mu \nu } \widetilde{F}_{\mu \nu } \Big ), \end{aligned}$$

(2.1)

where ${\mathcal {V}}_{\mu \nu } = \partial _\mu {\mathcal {V}}_\nu - \partial _\nu {\mathcal {V}}_\mu$, $\widetilde{{\mathcal {V}}}_{\mu \nu } = \frac{1}{2} \epsilon _{\mu \nu \rho \sigma } {\mathcal {V}}^{\rho \sigma }$, and $A_{2,3}^{{\mathcal {V}}\gamma }$ are dimensionless form factors. Such form factors receive contributions from the SM loop-level diagrams (see Fig. 1), and from the interaction beyond the SM, the latter in general possibly also containing CP-violating couplings. We take into account only the SM loop contributions, assuming that BSM loop corrections are small compared to the tree level ones. Thus, we put $A_3^{{\mathcal {V}}\gamma }=0$ while doing numerical study,² but for completeness we will keep the $A_3^{{\mathcal {V}}\gamma }$ dependent terms in our analytical expressions. The expressions for $A_2^{Z\gamma }$ and $A_2^{\gamma \gamma }$ in the SM are given in Ref. [23].

Let us denote the decay amplitude for $H\rightarrow \tau ^+\,\tau ^-\,\gamma$ by ${\mathscr {M}}$. As illustrated in Fig. 1, the amplitude can be split into three parts: (1) tree-level contribution ${\mathscr {M}}^{\text {(Yuk)}}$, (2) loop-level $Z\gamma$ contribution ${\mathscr {M}}^{(Z\gamma )}$, and (3) loop-level $\gamma \gamma$ contribution ${\mathscr {M}}^{(\gamma \gamma )}$, i.e. ${\mathscr {M}} = {\mathscr {M}}^{\text {(Yuk)}} + {\mathscr {M}}^{(Z\gamma )} + {\mathscr {M}}^{(\gamma \gamma )}$. Like any other 3-body decay of a spin-0 particle, the full kinematics of $H(p_H) \rightarrow \tau ^+ (p_+) \, \tau ^- (p_-) \, \gamma (p_0)$ can be described by two independent variables. We choose to work with Lorentz invariant mass squares. Defining

$$\begin{aligned} m_{+-}^2&= \left( p_+ + p_-\right) ^2 = \left( p_H - p_0\right) ^2, \end{aligned}$$

(2.2a)

$$\begin{aligned} m_{+0}^2&= \left( p_+ + p_0\right) ^2 = \left( p_H - p_-\right) ^2, \end{aligned}$$

(2.2b)

$$\begin{aligned} m_{-0}^2&= \left( p_- + p_0\right) ^2 = \left( p_H -p_+\right) ^2, \end{aligned}$$

(2.2c)

where

$$\begin{aligned} m_{+-}^2 + m_{+0}^2 + m_{-0}^2 = m_H^2 + 2 \, m_\tau ^2. \end{aligned}$$

(2.3)

We can express $\cos \theta$, defined in the di-tau rest frame, in terms of the Lorentz invariant variables:

$$\begin{aligned} \cos \theta = \left( 1 - \dfrac{4\,m_\tau ^2}{m_{+-}^2} \right) ^{-\frac{1}{2}} \, \frac{m_{-0}^2 - m_{+0}^2}{m_H^2 - m_{+-}^2}. \end{aligned}$$

(2.4)

At the beginning of this section we have argued that the forward–backward asymmetry in $\cos \theta$ distribution can serve as a probe of CP violation. Therefore, we see that the forward–backward asymmetry would be equivalent to an asymmetry in the distribution or number of events in the $m_{+0}^2$ vs. $m_{-0}^2$ plane (usually called a Dalitz plot) under the exchange $m_{+0}^2 \leftrightarrow m_{-0}^2$. Equivalently, one can consider distribution of events in the $m_{+0}$ vs. $m_{-0}$ plane which may be more convenient from experimental perspective. The ‘forward’ (or ‘backward’) region in Dalitz plot is that region where $m_{-0} > m_{+0}$ (or $m_{-0} < m_{+0}$).

In the rest frame of the Higgs boson, the differential decay rate of $H \rightarrow \tau ^+\,\tau ^-\,\gamma$ in terms of $m_{+0}$ and $m_{-0}$ is given by,

$$\begin{aligned} \dfrac{\textrm{d}^2 \Gamma _{\tau \tau \gamma }}{\textrm{d}m_{+0} \, \textrm{d}m_{-0}} = \frac{m_{+0} \, m_{-0}}{64 \, \pi ^3 \, m_H^3} \left| {\mathscr {M}} \right| ^2 \equiv {\mathcal {D}} \big (m_{+0},m_{-0}\big ), \end{aligned}$$

(2.5)

where the squared amplitude $\left| {\mathscr {M}} \right| ^2$ can be split into six constituents,

$$\begin{aligned} \left| {\mathscr {M}} \right| ^2&= \big |{\mathscr {M}}^{(\text {Yuk})}\big |^2 + \big |{\mathscr {M}}^{(Z\gamma )}\big |^2 + \big |{\mathscr {M}}^{(\gamma \gamma )}\big |^2 \nonumber \\&\quad + 2 \, \text {Re}\left( {\mathscr {M}}^{(\gamma \gamma )} \, {\mathscr {M}}^{(Z\gamma )*}\right) \nonumber \\&\quad + 2 \, \text {Re}\left( {\mathscr {M}}^{(\text {Yuk})} \, {\mathscr {M}}^{(Z\gamma )*}\right) \nonumber \\&\quad + 2 \, \text {Re}\left( {\mathscr {M}}^{(\text {Yuk})} \, {\mathscr {M}}^{(\gamma \gamma )*}\right) . \end{aligned}$$

(2.6)

In order to clearly point out the terms responsible for the forward–backward asymmetry and see how it is related to CP-asymmetry, we write down the expression for the individual constituents of amplitude square, as shown in Eq. (2.6), in terms of $m_{+-}^2$ and $\theta$. Using Eqs. (2.3) and (2.4) one can easily rewrite all these expressions in terms of $m_{+0}$ and $m_{-0}$. Neglecting the subdominant $m_\tau$ dependent terms in the numerator, we have:

$$\begin{aligned}&\big |{\mathscr {M}}^{(\text {Yuk})}\big |^2\nonumber \\ {}&= \frac{16 e^2 \left( a_{\tau }^2+b_{\tau }^2\right) m_{\tau }^2 \left( m_H^4+m_{+-}^4\right) \, m_{+-}^4 \, \sin ^2\theta }{v^2 \left( m_H^2-m_{+-}^2\right) ^2 \left( \left( m_{+-}^2 - 4\,m_\tau ^2\right) \, \sin ^2\theta + 4\,m_\tau ^2 \right) ^2}, \end{aligned}$$

(2.7a)

$$\begin{aligned}&\big |{\mathscr {M}}^{(Z\gamma )}\big |^2\nonumber \\ {}&= \frac{g_Z^2 \, \left( \left( c_A^\tau \right) ^2+\left( c_V^\tau \right) ^2\right) \, \left( \left( A_2^{Z\gamma } \right) ^2+\left( A_3^{Z\gamma } \right) ^2\right) }{8 v^2 \left( \left( m_{+-}^2-m_Z^2\right) ^2+\Gamma _Z^2 m_Z^2\right) } \nonumber \\&\quad \times m_{+-}^2 \left( m_H^2-m_{+-}^2\right) ^2 \left( 1+\cos ^2\theta \right) ,\end{aligned}$$

(2.7b)

$$\begin{aligned}&\big |{\mathscr {M}}^{(\gamma \gamma )}\big |^2\nonumber \\ {}&= \frac{e^2 \left( \left( A_2^{\gamma \gamma } \right) ^2+\left( A_3^{\gamma \gamma } \right) ^2\right) }{2 \, m_{+-}^2 \, v^2} \nonumber \\&\quad \times \left( m_H^2-m_{+-}^2\right) ^2 \left( 1+\cos ^2\theta \right) ,\end{aligned}$$

(2.7c)

$$\begin{aligned}&\text {Re}\left( {\mathscr {M}}^{(\gamma \gamma )} \, {\mathscr {M}}^{(Z\gamma )*}\right) \nonumber \\ {}&= -\frac{e \, g_Z \, \left( m_H^2-m_{+-}^2\right) ^2}{4 v^2 \left( \left( m_{+-}^2-m_Z^2\right) ^2+\Gamma _Z^2 m_Z^2\right) } \nonumber \\&\quad \times \Bigg (2 \, c_A^\tau \, \left( A_2^{\gamma \gamma } A_3^{Z\gamma }-A_2^{Z\gamma } A_3^{\gamma \gamma } \right) \, m_Z \, \Gamma _Z \, \cos \theta \, \nonumber \\&\quad + c_V^\tau \, \left( A_2^{\gamma \gamma } A_2^{Z\gamma }+A_3^{\gamma \gamma } A_3^{Z\gamma }\right) \,\nonumber \\&\quad \times \left( m_{+-}^2-m_Z^2\right) \left( 1+\cos ^2 \theta \right) \Bigg ),\end{aligned}$$

(2.7d)

$$\begin{aligned}&\text {Re}\left( {\mathscr {M}}^{(\text {Yuk})} \, {\mathscr {M}}^{(Z\gamma )*}\right) \nonumber \\ {}&= \frac{4 \, e \, g_Z \, m_{\tau }^2 \, m_{+-}^4 \, \sin ^2\theta }{v^2 \left( \left( m_{+-}^2-m_Z^2\right) ^2+\Gamma _Z^2 m_Z^2\right) } \nonumber \\&\quad \times \frac{1}{\left( \left( m_{+-}^2 - 4\,m_\tau ^2 \right) \, \sin ^2\theta + 4\,m_\tau ^2 \right) ^2} \nonumber \\&\quad \times \Bigg ( c_A^\tau \, \left( A_3^{Z\gamma } a_{\tau } - A_2^{Z\gamma } b_{\tau } \right) \, m_Z \, \Gamma _Z \, \left( m_H^2-m_{+-}^2\right) \, \cos \theta \nonumber \\&\quad + c_V^\tau \left( m_{+-}^2-m_Z^2\right) \bigg (A_2^{Z\gamma } a_{\tau } \left( m_H^2-m_{+-}^2 \cos ^2 \theta \right) \nonumber \\&\quad + A_3^{Z\gamma } b_{\tau } \left( m_H^2-m_{+-}^2\right) \bigg )\Bigg ),\end{aligned}$$

(2.7e)

$$\begin{aligned}&\text {Re}\left( {\mathscr {M}}^{(\text {Yuk})} \, {\mathscr {M}}^{(\gamma \gamma )*}\right) \nonumber \\ {}&= - \frac{8 \, e^2 \, m_{\tau }^2 \, m_{+-}^2 \, \sin ^2\theta }{v^2 \, \left( \left( m_{+-}^2 - 4\, m_\tau ^2 \right) \, \sin ^2\theta + 4\, m_\tau ^2 \right) ^2} \nonumber \\&\quad \times \left( A_2^{\gamma \gamma } a_{\tau } \left( m_H^2 - m_{+-}^2 \, \cos ^2\theta \right) \right. \nonumber \\ {}&\quad \left. + A_3^{\gamma \gamma } b_{\tau } \left( m_H^2-m_{+-}^2\right) \right) , \end{aligned}$$

(2.7f)

where $c_V^\tau = -1/2 + 2 \, \sin ^2\theta _W$, $c_A^\tau = -1/2$, and $g_Z = e/(\sin \theta _W \, \cos \theta _W)$, with $\theta _W$ being the weak mixing angle. Note that we have kept the total width of the Z boson, $\Gamma _Z$, because the Z boson can be on-shell in our case.

We are interested in terms that are odd (linear) in $\cos \theta$ (or, using Lorentz invariant variables, odd in the difference $m_{+0}^2 - m_{-0}^2$). Such terms are found to be proportional to $m_Z \, \Gamma _Z$ as well as the product of CP-even and CP-odd couplings. If we use the narrow-width approximation for the Z boson propagator,

$$\begin{aligned} \frac{1}{\left( m_{+-}^2-m_Z^2\right) ^2 + \Gamma _Z^2 \, m_Z^2} \approx \frac{\pi }{m_Z \, \Gamma _Z} \delta \big (m_{+-}^2-m_Z^2\big ), \end{aligned}$$

(2.8)

the $m_Z \, \Gamma _Z$ factor in the terms linear in $\cos \theta$ cancels out, and it is obvious that maximum CP-violation occurs for $m_{+-}^2 = m_Z^2$. Thus, the dominant contribution to the forward–backward asymmetry comes from the events for which invariant mass of the $\tau$ pair is close to the Z boson mass.

It is clear from Eq. (2.7e) that to a good approximation the asymmetry in the $\cos \theta$ distribution probes the combination $(A_3^{Z\gamma }a_\tau -A_2^{Z\gamma }b_\tau )$. In our numerical study in Sect. 3 we put $A_3^{Z\gamma }=0$.

In the following section we illustrate how the distribution pattern in the ‘forward’ and ‘backward’ regions of the Dalitz plot differ due to CP violation (i.e. $b_\tau \ne 0$) by studying the following distribution asymmetry,

$$\begin{aligned} {\mathcal {A}}\left( m_{+0},m_{-0}\right) = \frac{\big |{\mathcal {D}}\left( m_{+0},m_{-0}\right) - {\mathcal {D}}\left( m_{-0},m_{+0} \right) \big |}{{\mathcal {D}}\left( m_{+0},m_{-0}\right) + {\mathcal {D}}\left( m_{-0},m_{+0}\right) }. \end{aligned}$$

(2.9)

Additionally, we also study the asymmetry integrated over the region where the invariant mass of the $\tau ^+\,\tau ^-$ pair is close to the Z boson mass,

$$\begin{aligned} A(n) = \frac{\displaystyle \Big | \iint \Big ({\mathcal {D}}\big (m_{+0}< m_{-0}\big ) - {\mathcal {D}}\big (m_{+0}>m_{-0} \big ) \Big ) \, \Pi (m_{+-},n) \, \textrm{d}m_{+0} \, \textrm{d}m_{-0} \Big |}{\displaystyle \iint {\mathcal {D}}\big (m_{+0},m_{-0}\big ) \, \Pi (m_{+-},n) \, \textrm{d}m_{+0} \, \textrm{d}m_{-0}},\nonumber \\ \end{aligned}$$

(2.10)

where the function $\Pi \left( m_{+-}, n\right)$ defines the cut on the invariant mass of the $\tau ^+\,\tau ^-$ pair

$$\begin{aligned} \Pi \left( m_{+-}, n\right) = {\left\{ \begin{array}{ll} 1 &{}\quad \text {for } \big | m_{+-} - m_Z \big | \, \leqslant n \, \Gamma _Z,\\ 0 &{}\quad \text {otherwise}. \end{array}\right. } \end{aligned}$$

(2.11)

The asymmetry A(n) is directly related to the number of events around the Z pole,

$$\begin{aligned} A(n) = \frac{\left| N_F(n)-N_B(n) \right| }{N_F(n) + N_B(n)}, \end{aligned}$$

(2.12)

where $N_{F/B}(n)$ denote the number of events contained in the forward/backward region which are also contained in the region around Z pole as defined in Eq. (2.11).

3 Numerical study

In this section we do a numerical study of the effect of the CP violating parameter $b_\tau$ on the Dalitz plot distribution in $m_{+0}$ vs. $m_{-0}$ plane. Especially we focus on the size of the asymmetries ${\mathcal {A}}\left( m_{+0},m_{-0}\right)$ and A(n) as defined in Eqs. (2.9) and (2.10). As detailed below, we impose a few kinematic cuts in the Higgs rest frame. In the next section we present the results of a heuristic simple MC simulation as an attempt to be closer to the experimental conditions at the HL-LHC.

We note that by neglecting $m_\tau$ in comparison with Higgs mass $m_H$, one can constrain (to a very good approximation) the sum $a_\tau ^2 + b_\tau ^2$ from the experimentally measured $pp \rightarrow H \rightarrow \tau ^+ \, \tau ^-$ cross-section [24], which yields

$$\begin{aligned} a_\tau ^2 + b_\tau ^2 \approx 0.93^{+0.14}_{-0.12}, \end{aligned}$$

(3.1)

where the experimental errors have been added in quadrature.

To avoid infrared divergence, we impose a cut on the photon energy (i.e. specify a minimum energy for the photon) in the Higgs rest frame,

$$\begin{aligned} E_\gamma ^\text {cut} = 5~\text {GeV}. \end{aligned}$$

(3.2)

As we discuss later, the actual value of this cut has little impact on the decay branching ratios in the range of the di-$\tau$ invariant mass squared $m_{+-}^2$ most sensitive to the CP violation effect. For the sake of reference we note that, with this cut, for the full kinematical range of $m^2_{+-}$ the branching ratio of $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ is $BR_{\tau \tau \gamma }=3.72\times 10^{-3}$. The branching ratio decreases once a cut is imposed on the three relative angles $\theta _X$, with $X \in \{ +-,+0,-0 \}$ (see Fig. 3) among the final particles in the Higgs rest frame. An angular cut $\theta _X^{\text {cut}}$ specifies the minimum angle among the final particles. For $\theta _X^{\text {cut}} = 5^\circ$ we get $BR_{\tau \tau \gamma }=3.24\times 10^{-3}$ which further decreases by approximately 15% for each $5^\circ$ increase in the cut. Both the angular cut $\theta _X^{\text {cut}}$ and photon energy cut $E_\gamma ^\text {cut}$ affect the allowed values of $m_{+0}$ and $m_{-0}$.

Fig. 3

A schematic kinematic configuration of $\tau ^+ \, \tau ^- \, \gamma$ in the Higgs rest frame showing angles subtended by the 3-momenta of the final particles

Fig. 4

The expected differential decay rates ${\mathcal {D}}\big ( m_{+0}, m_{-0} \big )$ as well as the asymmetry ${\mathcal {A}} \big ( m_{+0}, m_{-0} \big )$ inside the Dalitz plot regions for $a_\tau =0.950$, $b_\tau =0.1$. The asymmetry is localised around the Z pole as expected (see discussion surrounding Eq. (2.8)). Note that the colour bar is in logarithmic scale

In Fig. 4 we see that the differential decay distribution have maxima close to the axes when $m_{\pm 0}^2$ approaches $m_\tau ^2$. These peaks are characteristic of the tree-level contribution from Fig. 1. A second peak is also easily discernible in the distributions around $m_{+-}^2 = m_Z^2$ as a slightly darker band, and this corresponds to contribution from the on-shell Z contribution, coming from the one-loop level diagrams of Fig. 1. Furthermore, for $b_\tau \ne 0$ we do find non-zero forward–backward asymmetry. Also as expected, the distribution asymmetry ${\mathcal {A}}\left( m_{+0},m_{-0}\right)$ become significantly large around the Z-pole region. The distribution asymmetry can be as large as $\sim 1\%$ depending on the values of $a_\tau$, $b_\tau$ such as for $a_\tau =0.950$ and $b_\tau =0.20$.

Regarding the asymmetries A(n) around the Z-pole, see Eq. (2.10), we note that the Z-pole cut as encoded in Eq. (2.11) can be rewritten, in terms of the photon energy in the Higgs rest frame, as follows,

$$\begin{aligned} \Pi \big ( m_{+-},n \big )&\equiv \Pi \big ( E_\gamma ,n \big ) \nonumber \\&= {\left\{ \begin{array}{ll} 1 &{} \text { for } \Big |\sqrt{m_H^2 - 2\,m_H\,E_\gamma } - m_Z \Big |\leqslant n\,\Gamma _Z,\\ 0 &{} \text { otherwise}, \end{array}\right. } \end{aligned}$$

(3.3)

From the equation above it is clear that for the invariant mass of the $\tau$ pair close to the Z pole, say $\big | m_{+-} - m_Z \big | \, \leqslant 5 \, \Gamma _Z$, that the photon energy cut $E_\gamma ^{\text {cut}}=5$ GeV has no relevance, since the minimum photon energy required for events around Z-pole corresponds to higher photon energies. Only the angular cuts $\theta _X^{\text {cut}}$ have any bearing in such a case.

In Fig. 5 we show the variation of A(n) for $1 \leqslant n \leqslant 5$ and compare it with with the ratio $\Gamma _{\tau \tau \gamma } \text {(around}$Z$\text {pole)}/\Gamma _{\tau \tau \gamma } \text {(full)}$, where $\Gamma _{\tau \tau \gamma } \text {(around}$Z$\text {pole)}$ is the partial decay rates for the decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ with $m_{+-}$ around the Z pole (imposed using Eq. (2.11)), and $\Gamma _{\tau \tau \gamma } \text {(full)}$ is the full partial decay rate. As expected, the asymmetry decreases with n, as it is strongly localised around the Z-pole, whereas $\Gamma _{\tau \tau \gamma } \text {(around }$Z$\text {pole)}/\Gamma _{\tau \tau \gamma } \text {(full)}$ increases with n. The plot clearly shows the challenge for an experimental analysis to find an optimal balance between the magnitude of the effect and the statistics of the events.

4 Simulation study of $\varvec{H \rightarrow \tau ^+ \, \tau ^- \, \gamma }$ in the context of HL-LHC

To estimate the sensitivity of the proposed Higgs boson decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ to the CP violation effects at the HL-LHC, Monte-Carlo (MC) generators were used to simulate the signal in the actual experimental environment. However, due to the limited computing resources, we have used a simplified MC simulation procedure. The differential cross-sections corresponding to the various $a_\tau$ and $b_\tau$ values are computed using GNU Octave [25] as a function of $m_{+0}$ and $m_{-0}$. The MC signal samples are re-weighted using these cross-sections (which include the kinematic cuts of Sect. 3) to properly model the impact of the interference term, similar to the “interpolation” approach used in [26]. The validity of the approach is verified by the comparison of relevant kinematic distributions with the analytical calculations.

We project the Dalitz plot distribution of events in forward and backward regions onto the $m_{+-}$ axis to do a 1-dimensional binned study of the forward–backward asymmetry. A more detailed and thorough MC study taking the full 2-dimensional Dalitz plot distribution into account and exploring unbinned Dalitz plot analysis techniques such as the Miranda method [27, 28], the method of energy test statistic [29,30,31,32,33] and the earth mover’s distance [34] are reserved for future explorations.

In the following, all additional cuts are defined in the laboratory frame. Reconstruction of the Higgs rest frame, that was used in the previous section would require the knowledge of the Higgs boson three-momentum which is not known experimentally. Besides, the observed distribution of events in $m_{+0}$ vs. $m_{-0}$ Dalitz plot can be obtained in any frame of reference.

4.1 Monte Carlo simulation

Fig. 5

The comparison of forward–backward asymmetry A(n), defined in Eq. (2.12) and the ratio $\Gamma _{\tau \tau \gamma } \text {(around}$Z$\text {pole)}/\Gamma _{\tau \tau \gamma } \text {(full)}$, where $\Gamma _{\tau \tau \gamma } \text {(around}$Z$\text {pole)}$ is the partial decay rates for the decay $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ with $m_{+-}$ around the Z pole satisfying Eq. (2.11), and $\Gamma _{\tau \tau \gamma } \text {(full)}$ is the full partial decay rate

Fig. 6

MC-based predictions of the a,b invariant mass of the di-$\tau$ system as well as the c–f asymmetry sensitivity for different values of $a_{\tau }$ and $b_{\tau }$ as a function of di-$\tau$ invariant mass. In c–e the asymmetry is normalised such that it is always 0 in the SM, while in f the statistical fluctuation of the MC generator are included, with the number of events around 5 times larger than the yields expected at HL-LHC

For the calculation of gluon-fusion production of the Higgs boson, the PowhegBox v2 [35,36,37,38] generator was used with the NNPDF3.0NNLO [39] PDF set. Proton-proton collisions are set to happen at center-of-mass energy of 14 TeV, as is expected for HL-LHC. For the simulation of the decay of the Higgs boson, modelling of the parton showers, and hadronization, the simulated events were processed with the Pythia v8.306 [40] program with the CTEQ6L1 [41] PDF set. DELPHES 3.5 [42] framework is then used to emulate the resolution and reconstruction of physical objects (such as photons, $\tau$ leptons, and jets) by a general-purpose particle detector (such as ATLAS or CMS) using the “HLLHC” card. FastJet 3.3.4 [43] package is used to perform the jet clustering using the anti-$k_t$ algorithm [44].

In the simulation studies photons are required to have $p_T > 10$ GeV and to be isolated with an angular cone defined by the condition³$\Delta R\le 0.3$. The reconstructed $\tau$ leptons are required to have $p_T > 15$ GeV. Their reconstruction is based on seed jets with the radius parameter [44] $\textsf {R} = 0.4$. This $p_T$ selection represents a realistic lower limit of what a general purpose detector can achieve. We assume that hadronically decaying $\tau$ leptons can be identified with $100\%$ efficiency. In reality this efficiency will be heavily dependent on the desired jet rejection power achievable with the conditions of the HL-LHC. The results presented in this section scale trivially with the $\tau$ identification efficiency. This optimisation is left for the future, more realistic, simulations of the performance of $\tau$ identification algorithms at the HL-LHC. All plots in this subsection are based on the MC simulation described above.

Fig. 7

The invariant masses $m_{+-}, m_{+0}$ and $m_{-0}$ using true, visible, and fitted physical objects

The HL-LHC is expected to deliver about 3000 fb$^{-1}$ integrated luminosity of data [45]. This corresponds to over 160 million events with gluon-gluon fusion production of the Higgs boson. With hadronically reconstructed $\tau$s and taking the same kinematic constraints as considered in Sect. 3, we estimate that $2.24 \times 10^{5}$ of these Higgs bosons will eventually decay into the $\gamma \, \tau ^+_\text {had} \, \tau ^-_\text {had}$ final state⁴. Approximately $10\%$ of the events will have the di-$\tau$ system with the invariant mass $m_{+-}$ within 5 GeV of the Z-boson mass peak where the forward–backward asymmetry manifests, see Fig. 6a. For any selected range of $m_{+-}$ we can estimate the number of events in ‘forward’ and ‘backward’ regions, say $N_F$ and $N_B$ respectively. Thus we can easily estimate the following forward–backward asymmetry,

$$\begin{aligned} A = \frac{N_F - N_B}{N_F + N_B}. \end{aligned}$$

(4.1)

The laboratory frame kinematic requirements applied to the reconstructed objects, such as the $\tau$ and photon $p_T$ and isolation requirements, further reduce the number of available events by a factor of 3 in the Z mass peak region, see Fig. 6b. The photon $p_T$ requirement by itself is responsible for a $50\%$ decrease in the selection efficiency.⁵

4.2 Kinematic fit

Although the true invariant mass of the di-$\tau$ system ($m_{+-}$) offers a good way to access the forward–backward asymmetry, see Fig. 6c, it is not accessible experimentally. The short lifetime of the $\tau$ leptons means that they will decay before reaching the detector, with $\nu _{\tau }$ escaping undetected. For hadronically decaying taus that are used in the present study, the particles registered in the detector will be predominantly charged and neutral pions. The detectors have limited acceptance and resolution, meaning that energies and momenta of these particles will be reconstructed with a limited accuracy. The visible invariant mass of the di-$\tau$ system ($m_{+-}^\text {vis}$), constructed from the visible decay products of $\tau$ decays, offers a degraded sensitivity to the forward–backward asymmetry, with almost no visible Z peak, see Fig. 6d. A fit procedure to recover the sensitivity to the asymmetry based on the kinematic constraints of the system is described in the following,⁶

Fig. 8

Gaussian fits to a truth-based asymmetry shapes and b data-like shapes in MC with $b_{\tau }=0.1$ with the $b_{\tau }=0.4$ fit overlaid. c Shows the relationship between true $b_{\tau }$ and the scaling factor c from the fits of analytically computed asymmetry A distributions. Note that the uncertainty bars in b represent the statistical uncertainty expected at HL-LHC

The final state of $H \rightarrow \gamma \, \tau ^+ \, \tau ^-$ is subject to two constraints:

• 1. The true invariant mass of the three final particles must be equal to the mass of the Higgs boson.

• 2. The energy in the transverse plane, perpendicular to the beam line, should be conserved and equal to 0, with any deviations coming from either the missing neutrinos ($\nu _\tau , \overline{\nu }_\tau$) or mismeasurements of the particle’s energies.

A fit procedure using Minuit2 [46] is performed based on these two conditions with the overall energy of the two $\tau$ leptons as free parameters. Since the opening angle between the neutrinos and visible tau decay product has to be of the order of $m_\tau / E_\tau$, both $\nu _{\tau }$ and $\bar{\nu }_\tau$ are predominantly collinear with the visible parts of the hadronically decaying $\tau$’s, for the energies considered here. Therefore, the approach of treating the contributions from the $\tau$ neutrino and $\tau$ energy smearing as one common parameter that only affects the energy of the $\tau$-lepton and not its spacial direction is justified.

This simple fit procedure allows us to restore the true energies of the $\tau$-leptons and the fitted two-body invariant masses match well with the true invariant masses, as demonstrated in Fig. 7. Further the fitted invariant masses $m^\text {fit}_{+0}$ and $m^\text {fit}_{-0}$ are used to identify events in forward and backward regions. This information is then used to estimate the asymmetry $A^\text {fit}$ while selecting $m^\text {fit}_{+-}$ in the region around Z-boson mass where the asymmetry is maximal. The asymmetry defined using the fitted masses, $A^\text {fit}$ behaves similarly as expected for the true asymmetry as a function of the fitted di-$\tau$ mass, see Fig. 6e. Thus $A^\text {fit}$ is a reasonable estimator of the forward–backward asymmetry. In the following we evaluate this asymmetry in a real-data-like environment.

4.3 Asymmetry calculation

Table 1 True and predicted values of $b_{\tau }$ from two different methods. Integrated luminosity of 3000fb$^{-1}$ is assumed

We compare two methods of quantifying the asymmetry and estimating the corresponding values of $b_{\tau }$. The first approach uses a simple selection of events with $m^\text {fit}_{+-}$ close to the Z mass peak, where the asymmetry is maximised. Here, the window of $\pm 9$ GeV around $m_Z$ was chosen, i.e. $\left| m^\text {fit}_{+-} - m_Z \right| \leqslant 9$ GeV. The width of this window was inspired by the range of the di-tau mass $m^\text {fit}_{+-}$ where the asymmetry is enhanced, see Fig. 6e, f. The asymmetry estimate $A^\text {fit}$ is then computed as in the Eq. (4.1) and used to predict $b_\tau$.

Fig. 9

MC-based comparisons of the $H\rightarrow \tau \tau \gamma$ signal and $Z+\gamma$ background samples, showing a, b the angular distance between the $\tau$-leptons and photon, c angular distance between the di-tau system and photon, d the transverse momentum of the leading photon. Both samples are normalized to 1 for ease of comparison

The second approach involves widening the di-$\tau$ mass selection to the range of 72–114 GeV. The asymmetry $A^\text {fit}$ is computed in bins of 3 GeV. A skewed Gaussian $f_\text {skew}(t)$ is then fitted to the shape of the asymmetry distribution,

$$\begin{aligned} f_\text {skew}(t) = c \phi (t) \Phi (\alpha t), \qquad t = \frac{m_{+-}^\text {fit}-a}{b}, \end{aligned}$$

(4.2)

where $\phi (t)$ is the normal probability density function, $\Phi (t)$ is the normal cumulative distribution function, and $a,b,c,\alpha$ are the free parameters in the fit. The parameters can be determined by fitting the true asymmetry distributions for specific values of $b_\tau$ as shown in Fig. 8a. Except for the parameter c, which determines the height of the skewed Gaussian and depends on $b_\tau$, all other parameters agree for different $b_\tau$ values (within statistical uncertainty). Keeping all parameters except the overall scaling factor c fixed at the determined values, we fit the distribution of asymmetry $A^\text {fit}$ in the reconstructed MC, which gives us the fitted value of c. As an example the dataset corresponding to $b_\tau = 0.1$ is shown in Fig. 8b along with the fitted skewed Gaussian distribution. For comparison, the skewed Gaussian corresponding to the case of $b_\tau = 0.4$ fit is overlaid. In Fig. 8c we show that the scaling factor c obtained from fitting is directly proportional to $b_\tau$. Thus knowing c one can directly infer the value of $b_\tau$.

The results of the two approaches are summarised in Table 1. The uncertainties of the measurements include the statistical uncertainty of the expected HL-LHC event yields, which is the dominant one. To estimate the HL-LHC uncertainty contribution we rescale the yields to match those expected at 3000fb$^{-1}$ and recompute the statistical uncertainty accordingly. Both of the approaches produce comparable central values with the fit to $f_\text {skew}$ resulting in lower uncertainties. Note that an offset between the theoretical input values and predicted central values of $b_\tau$ is present, due to a statistical fluctuation in the MC sample. The offset remains constant for all tested $b_\tau$ values and is reproduced by both methods tested, further validating their stability.

4.4 Backgrounds contributions and its separation

The dominant SM background we have to consider is $Z(\tau \tau )+\gamma$ production, similar to the ATLAS [47] and CMS [48] $H \rightarrow Z(\ell \ell ) \gamma$ searches. However, compared to the $\ell \ell$ channel, in the $\tau \tau$ channel the non-resonant Higgs decay contribution is dominating, resulting in the much larger (Higgs)/(non-Higgs backgrounds) fraction in the Z mass peak region of $m_{+-}$. The kinematics of the events also change, allowing for easier discrimination of the background. Developing an algorithm for suppression of the $Z+\gamma$ backgrounds (which most likely would have to be done with the application of machine learning techniques to fully utilize several correlated kinematic values) is beyond the scope of this paper. However to illustrate the possibility of such classifiers we perform a simple MC comparison.

The generation of $Z+\gamma$ MC is done using MadGraph5 aMC@NLO 3.5 [49] with the NNPDF3.0NLO [39] PDF set. Further decay chains and hadronization are handled by Pythia and DELPHES with the same setup as described in Sect. 4.1. We note several key kinematic variables related to $\tau$-leptons and photons that can help discriminate between signal and background, such as $\Delta R(\tau ,\gamma )$, $\Delta R(Z^\text {vis},\gamma )$, photon $p_T$, see Fig. 9. Higgs candidate $p_T$ and invariant mass can also be useful, as has been shown in the ATLAS and CMS searches. Based on this we believe that an efficient classifier can be built to significantly suppress the background contribution. The uncertainties quoted in Table 1 will increase in the presense of background events, but the exact effect depends strongly on how well the signal/background separation can be performed.

5 Conclusions

We have analysed the 3-body decay of the Higgs boson $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ as an additional source of information about the CP violation in the $H\tau \tau$ Yukawa coupling, independent from the existing experimental studies on the 2-body decay $H \rightarrow \tau ^+\,\tau ^-$ [9, 11,12,13,14,15,16,17,18]. The forward–backward asymmetry in the $\tau$ angular distribution in our case arises due to the interference of the tree-level contribution (which includes the CP violating $H\tau \tau$ Yukawa coupling $b_\tau \ne 0$) and the CP-even SM loop-level contributions. We have proposed a novel method of measuring forward–backward asymmetry in the Dalitz plot distribution of events in the plane of $\gamma \,\tau ^\pm$ Lorentz invariant masses ($m_{+0}$ vs. $m_{-0}$ plane). Such a Dalitz plot distribution is frame independent, making the method of extraction of forward–backward asymmetry clean and attractive from the experimental point of view. The asymmetry is directly proportional to the CP-odd $H\tau \tau$ coupling parameter $b_\tau$. In principle, the asymmetry can also appear from the interference of CP-even tree-level contribution and CP violating loop-level contributions. However, for our numerical study, we assume no CP-violation at loop-level and focus only on the effects of non-zero $b_\tau$ and whether this can be experimentally probed at HL-LHC.

The forward–backward asymmetry is predicted to be the largest when the di-$\tau$ invariant mass $m_{+-}$ is close to $m_Z$ (it could reach $\sim 1\%$ for high values of $b_\tau$) and it rapidly diminishes as one moves farther away from the Z pole. To estimate the feasibility of such asymmetry measurements at the HL-LHC we have performed a simplified MC simulation with kinematic cuts meant to mimic the experimental conditions. A kinematic fit was used to constrain the hadronically reconstructed $\tau$-leptons and account for the missing $\nu _{\tau }$ information not available in the detector. We estimated the asymmetry directly in the region with di-$\tau$ mass in the range of $m_Z \pm 9$ GeV for different values of $b_\tau$. We also looked for the asymmetry by performing a shape fit in a wider mass region, $72~\text {GeV} \leqslant m_{+-} \leqslant 114~\text {GeV}$.

From our MC studies we find that the statistical uncertainties we currently expect to get with the HL-LHC dataset are significantly larger than the effect itself. Nevertheless, our simplistic MC study suggests that our proposed methodology is experimentally doable, and our results could be encouraging for more detailed and in-depth explorations in the future. Instead of the one-dimensional binned shape fit used in this study a full two-dimensional unbinned Dalitz plot analysis could instead be envisaged using for example the Miranda method [27, 28], the method of energy test statistic [29,30,31,32,33] and the earth mover’s distance [34]. The asymmetry can also appear from the interference of CP-even tree-level contribution and CP violating loop-level contributions, this effect has not been considered in our numerical studies yet. In our MC simulation, we have only considered final states with both of the $\tau$-leptons decaying hadronically, the dataset can be doubled by also considering one of the $\tau$s to decay leptonically, i.e. adding the $H \rightarrow \tau _\text {had} \tau _\text {lep} \gamma$ decay channel. With the better understanding of the technical capabilities of particle detectors such as ATLAS and CMS after the Phase-2 upgrades, the kinematic selections can be further optimised.

Finally, once the asymmetry can be probed with reduced uncertainty, it would be interesting to compare its prediction for $b_\tau$ with that obtained from the already ongoing experimental study of $H \rightarrow \tau ^+ \, \tau ^- \rightarrow m^+ \, \overline{\nu }_\tau \, m^- \, \nu _\tau$ where $m=\pi ,\rho$ etc. If there is significant deviation in the two $b_\tau$ values, one can assume that there is some significant CP-violation coming from the loop-level contribution, which we have neglected in our numerical study in this paper. It is interesting to note that the same loop-level diagrams also contribute to $H \rightarrow \ell ^+ \, \ell ^- \, \gamma$ for $\ell =e,\mu$, and for these decay modes the tree-level contributions are negligible. Moreover, the same Dalitz plot techniques developed for $H \rightarrow \tau ^+ \, \tau ^- \, \gamma$ can also be applied to probe the asymmetry in the Dalitz plots of $H \rightarrow \ell ^+ \, \ell ^- \, \gamma$ to constrain or discover the CP violation at loop-level. Therefore, our formalism of probing the forward–backward asymmetry inside the Lorentz invariant Dalitz plot distribution of events would certainly help explore CP property of the Higgs boson in a more systematic and unified manner.