Testing CP properties of the Higgs boson coupling to $\tau$ leptons with heterogeneous graphs

Following the discovery of a scalar resonance with a mass near 125 GeV [1, 2], researchers are now focused on measuring the properties and characteristics to assess whether it aligns with the Higgs boson predicted by the Standard Model (SM) [3, 4]. A key aspect of this investigation is examining its spin and CP transformation properties [5, 6, 7]. For example, top-associated production is a promising channel to extract such properties [8, 9, 10]. Testing its CP properties is especially significant in light of the observed baryon asymmetry of the universe. The SM CP violation has been observed initially in Kaon decays [11], and established by the measurements of the direct CP violation in K system [12], and CP violation in neutral B meson decays [13, 14]. However, it is insufficient to account for the Bayron asymmetry of the Universe [15]. Therefore, additional sources of CP violation are essential ingredients of BSM to address the origin of the matter in our Universe. So far, measurements of the Higgs boson properties, for instance, the interactions with gauge bosons [16, 17, 18, 19, 20, 21], performed by the ATLAS and CMS experiments show no deviations from the SM predictions. Still, the possibility of an extended scalar sector that includes CP violation, and thus, the observed scalar resonance being a CP mixing state, has not been ruled out.

Besides the top-associated Higgs production, studying the $\tau$ spin correlation of the $h\tau\tau$ coupling yields valuable information about the CP state of the Higgs boson and has been proposed for both the LHC and lepton colliders [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Here, the primary focus is on reconstructing acoplanarity angle of the decays of the $\tau$ lepton pairs from the angular correlation of the decay products. A significant challenge in this kind of study at LHC arises from the presence of neutrinos in $\tau$ lepton decays, which complicates the accurate determination of the $\tau$ momentum vector and, consequently, the angular distribution between the $\tau$ lepton pairs at the LHC. However, the one-prong decays of the $\tau$ lepton, $\tau^{\pm}\rightarrow\pi^{\pm}\nu$, $\tau^{\pm}\to\rho^{\pm}(\rho^{\pm}\to\pi^{\pm}\pi^{0})\nu_{\tau},\tau^{\pm}\to a_{1}^{\pm}(a_{1}^{\pm}\to\pi^{\pm}\pi^{0}\pi^{0})\nu_{\tau}$, are promising. For the $\pi^{+}$ final state, the CP information can be reconstructed from charged pion momentum and interaction point (IP) information [26]. For the case of $\rho$ and $a_{1}$ final states, the angle between the decay planes of ($\pi^{+},n\pi^{0}$), ($\pi^{-},n\pi^{0}$) retains the spin correlation information of the $\tau$ lepton pair for $\rho$ and $a_{1}$ mode. It is essential to combine as many final states as possible; the branching ratio of $\tau^{+}\to\pi^{+}\nu$ final state is only 10.8% of the total decay width but increases up to 45.6% if one includes up to $\pi^{+}2\pi^{0}$ final state. Analyses of Higgs boson decays into $\tau$ lepton pairs have been conducted by both ATLAS [35, 16, 36] and CMS [37, 38, 39, 40]. They show agreement with the SM prediction of the Higgs boson but do not rule out CP admixture. ATLAS and CMS will accumulate more than 200 fb^-1 for each experiment in Run 3 (2024-2025), in the future, the High Luminosity Large Hadron Collider (HL-LHC) aims to collect 3000 fb^-1, to reveal the nature of the Higgs boson.

The hadronic final states of $\tau\bar{\tau}$ suffers serious background from various SM processes. Cuts to reduce the background have been developed, such as the isolation of the jet, jet mass, change multiplicities, and momentum distribution of pions. Deep Learning Neural Networks (DNNs) can efficiently increase the signal-to-background yield. Recently, DNNs have been widely used in collider analysis for various tasks, see [41] and references therein. In this paper, we utilize different sets of DNNs to suppress background events and study the CP properties of the Higgs boson at the HL-LHC. The first DNN we consider is the Multi-Layer Perceptron (MLP), which analyzes high-level kinematic and CP distribution. Although MLP can achieve high classification performance between signal and background events, it does not provide optimal performance for an analysis of the CP properties of the Higgs boson. This is because CP information is fully mixed with event kinematics; the learned information about Higgs boson CP properties will be diluted and network performance will be degraded. The issue of degraded performance in MLPs due to inputs containing mixed information has been highlighted, for example, in [42].

Alternatively, one can utilize a heterogeneous graph to analyse the $h\rightarrow\tau\bar{\tau}$. A heterogeneous graph consists of nodes and edges with different types of information, which allows the separate encoding of CP and kinematical information. We construct a graph of the nodes of the final state pions and the reconstructed $\tau$ pair with selective connections. Pion nodes are fully connected to recover the kinematic information, while tau nodes are connected with edges weighted by the value of the reconstructed angular distribution between them. This approach allows CP and kinematic information to be separately encoded within a single graph.

To analyze these constructed graphs, we utilize two Graph Neural Networks (GNNs): Graph Convolutional Network (GCN) and Graph Transformer Network (GTN). One key advantage of the GTN is its ability to dynamically capture the complex information within the graph, as GTNs leverage the attention mechanism to weigh the importance of different nodes irrespective of their position, while GCNs rely on localized neighborhood aggregation. Thanks to the attention mechanism, GTNs can dynamically focus on the most relevant parts of the heterogeneous graph. This flexibility, together with the enhanced capability to model intricate patterns, make GTNs particularly powerful for heterogeneous graph classification tasks. We find that GTN can enhance the current sensitivity on CP nature of ATLAS and CMS analysis.

This paper is organized as follows: In Section 2, we describe the parameterization of the effective Lagrangian relevant to this study and the construction of the CP observable in different hadronically decaying channels of the $\tau$ lepton. Section 3 details the analysis methodology and the Monte Carlo tools used for event simulation. In Section 4, we discuss the three deep learning methods, MLP, GCN and GTN. The results of this study are presented in Section 5 and our conclusion is drawn in Section 6.

The general form of the effective Yukawa interaction between the Higgs boson $h$ and $\tau$ leptons can be written as

Here, $\kappa_{\tau}$ is the reduced Yukawa coupling strength, $v=246$ the vacuum expectation value (vev) of $h$ and $\theta_{\tau}$ is the CP mixing angle, with $\theta_{\tau}$ ranging from $-90^{\circ}$ to $90^{\circ}$. This angle parameterized the relative contributions of the CP-even and CP-odd components to the $h\tau\tau$ coupling. Specifically, $\theta_{\tau}=0^{\circ}$ represents a purely CP-even state, while $\theta_{\tau}=90^{\circ}$ represents a purely CP-odd state. Values of $\theta_{\tau}$ between these extremes indicate an admixture of both components, suggesting a CP-violation in the Higgs sector. The CP-mixing angle $\theta_{\tau}$ affects the correlations between the transverse spin components of $\tau$-leptons in $h\to\tau\tau$ decays. These correlations, in turn, influence the directions of the $\tau$-lepton decay products. The acoplanarity angle $\Phi^{\ast}$, defined between the $\tau$ decay planes, is sensitive to these transverse spin correlations and is influenced by the CP-mixing angle of the Yukawa coupling. The $\Phi^{\ast}$ angle is directly connected to $\theta_{\tau}$ in the $h\to\tau\tau$ differential decay rate, with the relationship taking the form of a first-order trigonometric polynomial in $\theta_{\tau}$. The differential decay rate can be obtained as [43, 44]

where $a$ and $a^{\prime}$ are the tau decay products and $b(E_{\pm})$ is the spectral functions defined in [27].

The reconstruction of $\Phi^{\ast}$ requires the reconstruction of the $\tau$ decay planes, which is challenging at the detector level due to the missing energy associated with $\tau$ neutrinos. Various methods have been developed to approximate the acoplanarity angle based on differential techniques [23, 45, 25, 27, 28, 46, 31]. These methods are tailored to analyze specific $\tau$ lepton decay modes and are adjusted according to the number of visible particles within the $\tau_{\rm jet}$. Among them, we consider two methods for reconstructing $\Phi^{\ast}$[36]. For the $\tau^{\pm}\to\pi^{\pm}\nu_{\tau}$ decay, only one visible charged particle is present. In this case, the $\tau$ plane can be reconstructed using the IP method based on the impact parameter of the charged pions from the two $\tau$. On the other hand, if $\tau$ decays to $\rho^{\pm}\nu_{\tau}$ and $a^{\pm}_{1}\nu_{\tau}$, both charged and neutral pions are present. In this case we use the transverse momenta of the visible particles to estimate $\Phi^{\ast}$.

There is a significant advantage in having two visible particles from each $\tau$ decay. Although the momentum and plane of the $\tau$ lepton decays cannot be fully reconstructed due to the missing energy from the associated neutrino, the momenta of the charged and neutral pions can still be utilized. This allows for retaining information about the $\tau$ polarization and reconstructing $\Phi^{\ast}$. The acoplanarity angle can be derived from the charged and neutral pions at the LHC as [31, 47]

where $\hat{p}^{\pm}$ is the unit vector of the charged pion three momentum in the zero momentum frame of the $\rho$-meson pair, and $\hat{p}^{0+}_{\perp},\hat{p}^{0-}_{\perp}$ are normalised three momentum vectors of neutral pions transverse to the charged pion momentum. Another requirement is the discrimination of phase space with different $\tau$ polarization. The sign of the product of the $\tau$ lepton spin analyzing function, $Y=y_{-}^{\rho}\times y_{+}^{\rho}$ , where $y_{\pm}^{\rho}=\frac{E_{\pi^{\pm}}-E_{\pi^{0}}}{E_{\pi^{\pm}}+E_{\pi^{0}}}$, and $E_{\pi}$ is the pion energy in the laboratory frame, appears in the CP-mixing sensitive terms of the squared matrix element. As $Y$ is not positive definite, integration over pion momenta for both $Y>0$ and $Y<0$ would average out the CP mixing sensitive terms in the matrix element. Accordingly, the events from different classes are separated by shifting the events with $Y<0$ by $\pi$. This way the acoplanarity angle is modified for the case of $Y<0$ only and defined as

With this definition $\Phi^{\ast}$ has the range of $-\pi<\Phi^{\ast}<\pi$.

Analogously to the $\rho$ decay mode, the acoplanarity angle of the $\tau$ decay mode $\tau^{\pm}\to a^{\pm}_{1}\nu_{\tau}\to\pi^{\pm}2\pi^{0}\nu_{\tau}$ can be constructed by considering the four momenta sum of the neutral pions as taken in the neutral component of the $\rho$ method.

For the decay mode $\tau^{\pm}\to\pi^{\pm}\nu_{\tau}$, the IP method is employed to reconstruct $\Phi^{\ast}$. The impact parameter is defined as the shortest distance between the primary vertex and the pion momentum vector extended in the direction of the $\tau$ decay point. Since it is practically impossible to reconstruct the $\tau$ lepton momentum due to the presence of $\tau$ neutrinos among the decay products, the $\tau$ lepton decay plane is reconstructed from the track momentum and the impact parameter of the charged pion. First, the normalized impact parameters of the charged pions, $\hat{n}^{\ast}=(0,\vec{n}^{\ast\pm})$, are measured in the lab frame and then boosted to the zero momentum frame of the visible $\pi^{\pm}$ pair. The transverse components of the boosted impact parameters to the direction of the associated charged pion momentum, $\hat{n}^{\ast\pm}_{\perp}$, are used to define the acoplanarity angle as follows:

Considering this setup, we can analyze the CP properties of Higgs boson decays to $\tau$ lepton pairs from three hadronic decay modes of the $\tau$ lepton. In general, the IP method, which requires only one charged pion to reconstruct $\Phi^{\ast}$, is suitable for analyzing other decay modes such as $\tau^{\pm}\to\rho^{\pm}\nu_{\tau}$ and $\tau^{\pm}\to a^{\pm}_{1}\nu_{\tau}$. However, this method has low efficiency due to the significant uncertainty associated with IP reconstruction. The IP of the $\tau_{\rm jet}$ is relatively small compared to the tracking resolution, limiting the precision of its measurement despite the excellent resolution of the detector tracker. An advantage of the neutral pion method is that it does not rely on the reconstruction of the IP. Instead, it requires determining the direction of the neutral pion. The relatively large distance between the primary interaction point and the electro-magnetic calorimeter (ECAL) ($\mathcal{O}(1)$ m), coupled with the fine ECAL granularity, allows the direction of neutral pions to be reconstructed with smaller relative uncertainties compared to the IP.

With the theoretical framework established, we now proceed with a phenomenological investigation into the CP properties of the Higgs boson through its decay into a pair of $\tau$ leptons. Our analysis centers on the Higgs boson, with a branching ratio BR$(h\to\tau\tau)\sim 6.23\%$, produced via gluon-gluon fusion with a production cross section of $51$ pb at $\sqrt{s}=14$ TeV and $46$ pb at $\sqrt{s}=13.5$ TeV. We consider the hadronic decays of the $\tau$ lepton through three distinct modes: 1) both $\tau$’s decay $\tau^{\pm}\to\pi^{\pm}\nu_{\tau}$, 2) both $\tau$’s decay to $\tau^{\pm}\to\rho^{\pm}\nu_{\tau}$, and 3) both $\tau$’s decay to $\tau^{\pm}(\to a_{1}^{\pm}\nu_{\tau})\rightarrow\pi^{+}2\pi^{0}\nu_{\tau}$, with corresponding branching fractions of a single $\tau$ lepton is to approximately ${\rm Br}(\tau^{-}\rightarrow\pi^{-}\nu)=10.8\%$, ${\rm Br}(\tau^{-}\rightarrow\rho^{-}\nu)=25.49\%$, and ${\rm Br}(\tau^{-}\rightarrow\pi^{-}2\pi^{0}\nu)=9.26\%$ according to Ref. [48]. Our networks analyse these three modes as a single input without changing network structures, and extending the analysis into mixed final states, such as $\pi\rho$, $\pi a_{1}$, $\rho a_{1}$, is straightforward.

In this section, we discuss the construction of the combined signal and background events across different $\tau$ decay modes. Furthermore, we describe simulation tools we used for the simulation preserving $\tau$ spin correlation in the final states for this analysis.

In this section, we explore the application of various deep learning techniques, including MLP, GCN, and GTN, to investigate the CP properties of the Higgs boson. Each network is designed to handle specific types of input data according to its structure. The MLP is adept at analyzing high-level kinematic distributions, while the GCN and GTN are suited for analyzing heterogeneous graphs constructed from final state particles.

For the MLP study, the kinematical variables of the final state $\tau_{\rm jet}$ pairs are reconstructed, and $\Phi^{\ast}$ is calculated for the final state $\tau_{\rm jet}$ pair. The kinematical variables and $\Phi^{\ast}$ are fed into the MLP. Of course, the kinematical data is normalized between the standard 0 and 1 range to ensure effective processing by the neural network. The network is trained to distinguish the background and signal processes. Another method for training a network to distinguish between different CP states is to use a conditional DNN network, as described in [10].

GNNs, on the other hand, analyze graph-like structures. The standard way is that the nodes of the graph represent the final state particles, and all nodes are fully connected. The graph nodes are weighted with the four-momenta of the final state particles and edges are weighted with the angular distance between each node pair. With this approach it is not easy to incorporate the CP properties of the Higgs boson into a fully connected graph. Instead, we utilize a heterogeneous graph in this paper. A heterogeneous graph comprises multiple types of nodes and edges, each representing different entities and interactions within the experimental setup. Nodes in these graphs can still represent final state particles, but also the reconstructed $\tau_{\rm jet}$ and Higgs boson. Each node type has distinct attributes and properties that define its role within the graph. Edges represent the interactions or relationships between these entities. This method has enhanced flexibility representing the physics we focus on.

To efficiently incorporate CP information into a graph-like structure, we consider a fully connected graph of the final state pions with additional heterogeneous nodes representing the reconstructed $\tau_{\rm jet}$ and the Higgs boson. The edges of the fully connected pion graph are weighted with the angular distance between each node pair, while the edge between the two $\tau_{\rm jet}$ nodes is weighted with the value of the reconstructed $\Phi^{\ast}$. By selectively connecting $\tau_{\rm jet}$ nodes to their decaying pions, we construct a graph that integrates both the kinematic information of the signal and the CP properties of different Higgs boson states.

The discrimination power of each network is measured by the background rejection for a given signal efficiency. The discrimination power is intertwined with various kinematic distributions, including CP information. We utilize three different neural networks: MLP, GCN, and GTN, each trained on eleven benchmark points with $\theta_{\tau}$ ranging from $-90^{\circ}\ \rm to\ 90^{\circ}$. Each network was trained and tested on each benchmark point individually. We use the area under curve (AUC) of the receiver operating characteristic curve (ROC) to assess the networks’ performance. ROC is a curve of the True Positive Rate (TPR) as the function of the False Positive rate (FPR), and AUC is the area surrounded by ${\rm FPR}=1$ and ${\rm TPR}>0$. Figure 6 shows the AUC values for all eleven signal points for MLP (green), GCN (orange), and GTN (blue). GTN demonstrates superior performance with an AUC of approximately $88\%$ for all points, while GCN and MLP achieved AUCs of approximately $86\%$ and $84\%$, respectively. Note that the AUC does not depend on $\theta_{\tau}$ so that the use of DL does not introduce additional bias to the analysis.

We also list the number of misclassified events for $y>y_{\rm cut}$ where signal efficiency is 80% in table 2 for a benchmark point with $\theta_{\tau}=90^{\circ}$ for integrated luminosity of $100$ fb^-1 at $\sqrt{s}=14$ TeV¹¹1We take the integrated luminosity close to the current LHC analysis for comparison. . The acceptance of the background events is 1.7% for MLP, 1.0% for GCN and 0.7% for GTN. We also compute the signal significance, following [77, 78]

where $S,B$ is the number of the signal and background events, and $\delta_{B}$ represents the systematic uncertainty of the SM background events and is set to $20\%$ [36]. These results demonstrate that the $h\to\tau\tau$ process can be effectively identified at the HL-LHC by using the proposed networks.

To better understand the outputs of different networks, particularly the feature regions each network focuses on to achieve its classification performance, we use Shapley Additive Explanations (SHAP) [79]. SHAP is a method to estimate an importance value for each feature using the output of deep learning models. For a given prediction $f(x)$, the SHAP value for a feature $i$ is calculated as:

where $N$ is the set of all features, $S$ is a subset of $N$ excluding feature $i$, $f(S)$ is the prediction based on the features in subset $S$, and $f(S\cup\{i\})$ is the prediction with feature $i$ added to $S$. The SHAP value $\phi_{i}$, thus, represents the average contribution of feature $i$ to the prediction over all possible subsets $S$. This method ensures a fair and consistent allocation of feature importance considering the correlation between the input distributions. Since the networks used have a static structure with a fixed input dataset size, feature $i$ is randomly sampled to eliminate its impact on the network output.

Figure 7 shows the SHAP values for 10000 test events for a signal point with $\theta_{\tau}=90^{\circ}$ for the MLP (left) and GTN (right). SHAP values are computed for the inputs of both networks, including the reconstructed distributions shown in Figure 2 and the graph nodes plus the $\Phi^{\ast}$ edge shown in Figure 3. The results indicate that in both networks, $\Phi^{\ast}$ significantly affects the network output. Interestingly, the MLP focuses primarily on the information from $\not{E}_{T},\Phi^{\ast},p_{T_{(\tau_{h}\tau_{h})}}$ and the transverse momentum of the leading $\tau$ jet, while it is less sensitive to all other input distributions. Conversely, the GTN distributes its attention almost equally across all nodes in the input graph, which may explain the improved performance of the GTN over the MLP.

In this paper, we investigate the CP structure of the $h\tau^{\pm}\tau^{\mp}$ vertex at the HL-LHC with $\sqrt{s}=14$ TeV. We consider three different channels of hadronic $\tau$ lepton decays for the case where both of the $\tau$ lepton decays into the same final state: $\tau^{\pm}\to\rho^{\pm}(\rho^{\pm}\to\pi^{\pm}\pi^{0})\nu_{\tau}$, $\tau^{\pm}\to a_{1}^{\pm}(a_{1}^{\pm}\to\pi^{\pm}\pi^{0}\pi^{0})\nu_{\tau}$, and $\tau^{\pm}\to\pi^{\pm}\nu_{\tau}$. The CP structure of the $h\tau^{\pm}\tau^{\mp}$ vertex can be determined from the angular correlation of the $\tau$ spins. This correlation can be reconstructed from the angular distribution between the charged and neutral pions of the $\tau$ lepton pair decays, even though there are tau neutrinos in the final state.

To improve the projected reach of measuring the CP mixing angle at the HL-LHC, we consider utilizing advanced deep-learning networks to enhance the signal-to-background yield. For this purpose, we employ three different networks to analyze different data structures. The MLP is used to analyze the kinematic distribution and reconstruct the CP mixing angle. However, due to the fully connected nature of the MLP, it fully mixes the kinematic and CP information, diluting the learned CP information and hindering the overall classification performance. To overcome this, we consider heterogeneous graphs constructed from the information stored in the final and decayed particles. With a selective connection of the nodes we fix the processing of the kinematic and CP information. For heterogeneous graph analysis, we adopt two networks, GCN and GTN. We find that GTN shows superior performance in background rejection, achieving a signal significance of $8.6\sigma$ at the HL-LHC with $\mathcal{L}=100\,\rm{fb}^{-1}$ for a benchmark point with $\theta_{\tau}=90^{\circ}$. GCN and MLP have lower signal significance for the same benchmark point with $7.19\sigma$ and $5.6\sigma$, respectively.

After improving the background rejection, we perform a shape analysis for the remaining events. We use a binned negative log-likelihood analysis to estimate the probability of seeing $\theta_{\tau}\neq 0$ at the HL-LHC. Keeping the limitation of simplified analysis for $\tau$ jet tagging using truth information of $\tau\rightarrow\rho$ and $a_{1}$ decays and $\pi^{0}$ reconstruction, our results show that the pure CP-odd state is excluded at nearly $3\sigma$ using the MLP, while GCN and GTN exclude the pure CP-odd state at above the $3\sigma$ level showing stronger significance than current LHC measurements [36, 40]. Moreover, GTN shows superior performance excluding $\theta_{\tau}\geq|22^{\circ}|$ at $95\%$ C.L., while GCN excludes $\theta_{\tau}\geq|31^{\circ}|$ at $95\%$ C.L. and MLP excludes $\theta_{\tau}\geq|43^{\circ}|$ at $95\%$ C.L. GTN’s improved performance over GCN is due to the fact that GTN applies an attention mechanism during training. The main advantage of the attention mechanism is that it assigns weights to different elements in the input graph, emphasizing the more relevant parts while downplaying the less relevant ones. Conversely, GCN treats all neighbouring vertices equally during the feature aggregation process, which can lead to suboptimal performance if certain neighbours provide more valuable information than others, as is the case encoded in the considered heterogeneous graphs. To ensure the reproducibility of our results, we have made our codes and files publicly available on https://github.com/wesmail/HiggsCP