HIG-20-006

HIG-20-006

# Analysis of the CP structure of the Yukawa coupling between the Higgs boson
and leptons in proton-proton collisions at $\sqrt{s}=13\TeV$

###### Abstract

The first measurement of the CP structure of the Yukawa coupling between the
Higgs boson and leptons is presented. The measurement is based on data
collected in proton-proton collisions at $\sqrt{s}=13\TeV$ by the CMS detector
at the LHC, corresponding to an integrated luminosity of 137. The analysis
uses the angular correlation between the decay planes of leptons produced in
Higgs boson decays. The effective mixing angle between CP-even and CP-odd
Yukawa couplings is found to be $-1\pm 19^{\circ}$, compared to an expected
value of $0\pm 21^{\circ}$ at the 68.3% confidence level. The data disfavour
the pure CP-odd scenario at 3.0 standard deviations. The results are
compatible with predictions for the standard model Higgs boson.

##  0.1 Introduction

In the standard model (SM), the electroweak symmetry breaking is postulated
via the Brout–Englert–Higgs mechanism [1, 2, 3, 4, 5, 6]. This mechanism
predicts the existence of a scalar particle, the Higgs boson (). A particle
compatible with this boson was discovered by the ATLAS [7] and CMS [8, 9]
Collaborations at the LHC using proton-proton ($\Pp\Pp$) collision data
collected in 2011 and 2012 at centre-of-mass energies of 7 and 8,
respectively. Since 2012, the couplings of the Higgs boson to heavy quarks,
leptons, and gauge bosons have been measured, including the coupling to
leptons [10, 11, 12], and the most recent measurement of its mass value yields
$125.38\pm 0.14\GeV$ [13].

The SM is even under charge-parity (CP) inversion. A sizeable deviation from a
pure CP-even interaction of the with any of the SM particles would be a direct
indication of physics beyond the SM. Therefore, the CP structure of the
couplings of the is of paramount interest. The CMS and ATLAS Collaborations
have studied the couplings of the to vector gauge bosons, including tests of
CP violation [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. These
studies excluded pure pseudoscalar (CP-odd) interactions of the with the and
bosons (referred to collectively as bosons).

There are strong theoretical motivations to search for CP-violating effects in
couplings of the to fermions rather than bosons. In couplings to bosons, CP-
odd contributions enter either via higher-order operators that are suppressed
by powers of $1/\Lambda^{2}$ [28], where $\Lambda$ is the scale of the physics
beyond the SM in an effective field theory, or via nonrenormalisable
interaction terms [29, 30]. Therefore, these are expected to only yield a
small contribution to the coupling. A renormalisable CP-violating Higgs-to-
fermion coupling may occur at tree level. The lepton and top quark Yukawa
couplings, $\PH\PGt\PGt$ and $\PH\mathrm{tt}$, respectively, are therefore the
optimal couplings for CP studies in $\Pp\Pp$ collisions [31], and measurements
of these two couplings are complementary. Recently, both the CMS [32] and
ATLAS [33] Collaborations presented first measurements of the CP structure of
the coupling to top quarks. The CMS result rejects the purely CP-odd
hypothesis with a significance of 3.2 standard deviations, $\sigma$, while the
ATLAS analysis rejects it with a significance of $3.9\sigma$. The CMS and
ATLAS Collaborations have also studied the CP-nature of the interaction with
gluons [21, 34] which was found to be consistent with the SM expectation,
albeit with limited sensitivity. Such studies may also be interpreted in terms
of the coupling to top quarks, under the assumption that the interaction is
mediated predominantly via top quark loops.

The CP-properties of the $\PH\to\PGt\PGt$ process is commonly described in
terms of an effective mixing angle $\alpha^{\PH\tau\tau}$, which is virtually
equal to $0^{\circ}$ in the SM. The measurement of a nonzero
$\alpha^{\PH\tau\tau}$ would therefore directly contradict the SM predictions,
and have implications for beyond the SM physics models, such as two-Higgs-
doublet models [35], including supersymmetry. For example, in the minimal
supersymmetric model CP violation in the Higgs-to-fermion couplings is
expected to be small and therefore the measurement of a sizeable mixing angle
would disfavour such scenarios. In contrast, in the next-to-minimal
supersymmetric model, $\alpha^{\PH\tau\tau}$ can be as large as $27^{\circ}$
[36]. Feasibility studies have indicated that the LHC experiments can measure
$\alpha^{\PH\tau\tau}$ to a precision of about 5–$10^{\circ}$ with 3of data
[29, 37].

In this paper we present the first measurement of the CP structure of the
coupling to leptons. This analysis uses the $\Pp\Pp$ data sets collected by
the CMS detector at $\sqrt{s}=13\TeV$ in 2016, 2017, and 2018. These
correspond to integrated luminosities of 35.9, 41.5, and 59.7, respectively.
This analysis targets the most sensitive $\tauh\tauh$, $\tau_{\Pgm}\tauh$ and
$\tau_{\Pe}\tauh$ decay channels, where a lepton decaying to hadrons is
denoted as , and a lepton decaying to a muon or an electron as $\tau_{\Pgm}$
or $\tau_{\Pe}$ (or collectively as $\tau_{\ell}$), respectively. In total,
this covers about 70% of all possible lepton pair decay modes. Table 0.1
summarises the lepton decay modes used in this analysis, their branching
fractions, and the shorthand symbols that we use to denote them in the rest of
this paper. The decays into light leptons are accompanied by two neutrinos,
while the hadronic modes involve one neutrino. The charged hadrons are denoted
by the symbol $\Ph^{\pm}$, which consist mainly of charged pions but include a
smaller contribution from charged kaons. Throughout this paper we will assume
that all $\Ph^{\pm}$ are charged pions () since the CMS detector is not able
to distinguish between different types of $\Ph^{\pm}$.

Decay modes of leptons used in this analysis and their branching fractions
$\mathcal{B}$ [38]. Where appropriate, we indicate the known intermediate
resonances. The last row gives the shorthand notation for the decays used
throughout this paper.  Mode $\Pepm\nu\nu$ $\PGmpm\nu\nu$ $\Ph^{\pm}\nu$
$\Ph^{\pm}\Pgpz\nu$ $\Ph^{\pm}\Pgpz\Pgpz\nu$ $\Ph^{\pm}\Ph^{\mp}\Ph^{\pm}\nu$
Type $\tau_{\Pe}$ $\tau_{\Pgm}$ $\mathcal{B}(\%)$ 17.8 17.4 11.5 25.9 9.5 9.8
Resonance Symbol $\mathrm{a_{1}^{1pr}}$ $\mathrm{a_{1}^{3pr}}$

This paper is organised as follows. The parameterisation of the CP properties
of the Yukawa coupling is discussed in Section 0.2. In Section 0.3 the
experimental setup is outlined, and this is followed by a discussion of the
data sets and simulated samples in Section 0.4. Subsequently, the event
reconstruction is presented in Section 0.5. Thereafter, in Section 0.6 the CP-
sensitive observables used to extract the results are outlined. In Section 0.7
the event selection is presented. The estimation of the backgrounds is
discussed in Section 0.8. The techniques used to distinguish the signal from
the background events are outlined in Section 0.9. In Section 0.10 various
distributions that are used to extract the results are displayed and
discussed. In Section 0.11 the systematic uncertainties are presented. The
results are discussed in Section 0.12, and a summary of the analysis is given
in Section 0.13.

##  0.2 Parametrisation of the CP properties of the Yukawa coupling

We parameterise the Lagrangian for the Yukawa coupling in terms of the
coupling strength modifiers $\kappa_{\tau}$ and $\widetilde{\kappa}_{\tau}$
that parameterise the CP-even and CP-odd contributions, respectively [31]:

$\mathcal{L}_{\mathrm{Y}}=-\frac{m_{\PGt}}{v}\PH(\kappa_{\tau}\PAGt\PGt+\widetilde{\kappa}_{\tau}\PAGt
i\gamma_{5}\PGt).$ (1)

In this equation, $m_{\PGt}$ is the mass of the lepton, denotes the Dirac
spinor of lepton fields, and the vacuum expectation value of the Higgs field,
$v$, has a value of 246. The effective mixing angle $\alpha^{\PH\tau\tau}$ for
the $\PH\PGt\PGt$ coupling is defined in terms of the coupling strengths as

$\tan(\alpha^{\PH\tau\tau})=\frac{\widetilde{\kappa}_{\tau}}{\kappa_{\tau}},$
(2)

while the fractional contribution of the CP-odd coupling
$f_{{\textit{CP}}}^{\PH\PGt\PGt}$ is obtained from the mixing angle as
$f_{{\textit{CP}}}^{\PH\PGt\PGt}=\sin^{2}(\alpha^{\PH\tau\tau})$. A mixing
angle of $\alpha^{\PH\tau\tau}=0\,(90)^{\circ}$ corresponds to a pure scalar
(pseudoscalar) coupling. For any other value of $\alpha^{\PH\tau\tau}$, the
has a mixed coupling with CP-even and CP-odd components, with maximal mixing
at a value of $\pm 45^{\circ}$.

The angle $\phi_{{\textit{CP}}}$ denotes the angle between the lepton decay
planes in the rest frame. An illustration of the decay planes in the single
pion channel is depicted in Fig. 1. The relation between
$\alpha^{\PH\tau\tau}$ and $\phi_{{\textit{CP}}}$ may be inferred from the
decay of a via leptons to two outgoing charged particles [39] as

$\dd{\Gamma}{\phi_{{\textit{CP}}}}(\PH\to\PGtp\PGtm)\sim
1-b(E^{+})b(E^{-})\frac{\pi^{2}}{16}\cos(\phi_{{\textit{CP}}}-2\alpha^{\PH\tau\tau}).$
(3)

In this equation, the outgoing charged particles have an energy $E^{\pm}$ in
their respective rest frames. The functions $b$ are spectral functions [40]
that encapsulate the correlation between the spin and the momentum of the
outgoing charged particle. We note that the spectral functions for the
leptonic and various hadronic decays are different.

Figure 1: The decay planes of two leptons decaying to a single charged pion.
The angle $\phi_{{\textit{CP}}}$ is the angle between the decay planes. The
illustration is in the rest frame.

Figure 2 shows the normalised distribution of $\phi_{{\textit{CP}}}$ at the
generator level, calculated in the rest frame of the , for the scalar,
pseudoscalar, and maximally mixed values of $\alpha^{\PH\tau\tau}$, as well as
the $\phi_{{\textit{CP}}}$ distribution from Drell–Yan processes. The
simulated event samples that are used to generate these distributions are
discussed in Section 0.4. These distributions are for the scenario where both
leptons decay to a charged pion and a neutrino.

Figure 2: The normalised distribution of $\phi_{{\textit{CP}}}$ between the
lepton decay planes in the rest frame at the generator level, for both leptons
decaying to a charged pion and a neutrino. The distributions are for a
decaying scalar (CP-even, solid red), pseudoscalar (CP-odd, dash blue), a
maximal mixing angle of $45^{\circ}$ (CP-mix, dash-dot-dot green), and a
vector boson (black dash-dot). The transverse momentum of the visible decay
products $\pt^{\PGt}$ was required to be larger than 33during the event
generation.

There is a phase shift between different mixing scenarios such that the
difference in $\phi_{{\textit{CP}}}$ equals two times the difference in
$\alpha^{\PH\tau\tau}$, as given by Eq. (3). It is important to note that the
distribution of $\phi_{{\textit{CP}}}$ for the Drell–Yan background is
constant; we will exploit this symmetry to reduce statistical fluctuations in
the background estimates, as explained in Section 0.9.

The observable $\phi_{{\textit{CP}}}$ was originally introduced in the context
of collisions [41, 42] where the lepton momenta can be reconstructed and thus
$\phi_{{\textit{CP}}}$ can be calculated in the rest frame. In hadronic
collisions the momenta of the neutrinos cannot be well constrained, except for
the configuration in which both leptons decay via the $\mathrm{a_{1}^{3pr}}$
mode to three charged pions. Therefore, the methods for estimating
$\phi_{{\textit{CP}}}$ have been extended and optimised for hadronic
collisions [37], as discussed in Section 0.6. Throughout this document, we
will denote the angle between the decay planes as $\phi_{{\textit{CP}}}$,
irrespective of the frame in which it is calculated.

##  0.3 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6m
internal diameter, providing a magnetic field of 3.8T. Within the solenoid
volume are a silicon pixel and strip tracker, a lead tungstate crystal
electromagnetic calorimeter (ECAL), and a brass and scintillator hadron
calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward
calorimeters extend the pseudorapidity ($\eta$) coverage provided by the
barrel and endcap detectors. Muons are detected in gas-ionisation chambers
embedded in the steel flux-return yoke outside the solenoid. Events of
interest are selected using a two-tiered trigger system. The first level (L1),
composed of custom hardware processors, uses information from the calorimeters
and muon detectors to select events at a rate of around 100kHz within a fixed
latency of about 4 [43]. The second level, known as the high-level trigger,
consists of a farm of processors running a version of the full event
reconstruction software optimised for fast processing, and reduces the event
rate to around 1kHz before data storage [44]. A more detailed description of
the CMS detector, together with a definition of the coordinate system used and
the relevant kinematic variables, can be found in Ref. [45].

##  0.4 Simulated samples

The signal and relevant background processes are modelled with samples of
Monte Carlo simulated events. The signal samples with a produced through
gluon-gluon fusion ($\Pg\Pg\PH$), vector boson fusion (VBF), or in association
with a or vector boson (denoted as $\PW\PH$ or $\PZ\PH$, or $\PV\PH$ when
combined) are generated at next-to-leading order (NLO) in perturbative quantum
chromodynamics (QCD) with the 2.0 [46, 47, 48, 49, 50, 51, 52] event
generator. The production mechanism is configured to only produce scalar Higgs
bosons, as opposed to pseudoscalar Higgs bosons or mixed couplings; the latter
scenarios would also affect various properties of the production, for example
the production rate and the topology of associated jets when present. The
distributions of the Higgs boson’s transverse momentum () and of the jet
multiplicity are reweighted to match the predictions at next-to-NLO (NNLO)
accuracy obtained from full phase space calculations with the nnlops (version
1) generator [53, 54]. The decay of the does not depend on its production. The
description of the decay of the to leptons is obtained using the generator
version 8.230 [55]. These samples are simulated without accounting for the
spin correlations. After the samples have been generated, the tauspinner
package [56] is used to calculate event weights that can be applied to the
simulated signal samples to model polarisation effects for a boson with CP-
mixing angles of 0, 45, and $90^{\circ}$. There is no normalisation effect
from the reweighting procedure, the integrated $\PH\to\PGt\PGt$ cross section
of the signal samples is invariant under rotations in $\alpha^{\PH\tau\tau}$.
All 2016 samples are generated with the NNPDF3.0 [57] NLO parton distribution
functions (PDFs), while the NNPDF3.1 [58] NNLO distributions are used for
2017–2018.

The [59] generator (version 2.6.0) is used for processes involving a or boson
and up to four outgoing partons generated with the matrix element, and these
processes are denoted $\PZ+\text{jets}$ and $\PW+\text{jets}$, respectively.
Processes involving bosons originating from top quark decays are not
considered in these samples. They are simulated at leading order (LO) with the
MLM jet matching and merging approach [60]. The same generator is used at NLO
for diboson production, whereas 2.0 (1.0) is used for top quark-antiquark pair
production [61] and single top quark production (associated with a boson) [62,
63]. The generators are interfaced with to model the parton showering and
fragmentation, as well as the decay of the leptons. The parameters that affect
the description of the underlying event are set to the CUETP8M1 tune [64] in
2016, and CP5 tune [65] in 2017–2018.

Monte Carlo generated events are processed through a simulation of the CMS
detector that is based on [66], and are reconstructed with the same algorithms
as the ones used for data. Additional $\Pp\Pp$ interactions per bunch crossing
(”pileup”) are included. The effect of pileup is taken into account by
generating concurrent minimum bias collision events with . The pileup
distribution in simulation is weighted to match the pileup in data.

##  0.5 Event reconstruction

The reconstruction algorithms for both observed and simulated events are based
on the particle-flow (PF) algorithm [67], which relies on the information from
the different CMS subdetectors to reconstruct muons, electrons, photons, and
charged and neutral hadrons. These objects are combined to form more complex
ones, such as candidates or missing transverse momentum ().

###  0.5.1 Primary vertex reconstruction

The positions of all $\Pp\Pp$ interactions (vertices) in the event, including
the hard scatter (primary) and soft (pileup) vertices, are reconstructed in a
two-step procedure [68]. The steps consist in clustering the tracks that
appear to originate from the same interaction using the deterministic
annealing algorithm [69], and subsequently fitting the position of each vertex
using tracks associated to its cluster with the adaptive vertex fitter (AVF)
algorithm [70]. The candidate vertex with the largest value of the sum of the
$\pt^{2}$ of all associated physics objects is considered to be the primary
$\Pp\Pp$ interaction vertex (PV). The physics objects included in this sum are
jets, clustered using the anti-jet finding algorithm [71] with the tracks
assigned to candidate vertices as inputs, and the associated , taken as the
negative vector sum of the of those jets.

###  0.5.2 Muon reconstruction

Muons are identified and reconstructed with requirements on the quality of the
track reconstruction and on the number of hits in the tracker and muon systems
[72], and selected within $\abs{\eta}<2.4$. In order to reject muons that
originate from nonprompt interactions, or are misidentified, a relative
isolation is defined as

$I^{\Pgm}\equiv\frac{\sum_{\text{charged}}\pt+\max\left(0,\sum_{\text{neutral}}\pt-\frac{1}{2}\sum_{\text{charged,
PU}}\pt\right)}{\pt^{\Pgm}}.$ (4)

In this equation, $\sum_{\text{charged}}\pt$ is the scalar sum of the charged
particles originating from the PV and located in a cone of size
$\DR=\sqrt{\smash[b]{(\Delta\eta)^{2}+(\Delta\phi)^{2}}}=0.4$ (where $\phi$ is
azimuthal angle in radians) centred on the muon direction. The sum
$\sum_{\text{neutral}}\pt$ is a similar quantity for neutral particles. The
$\sum_{\text{charged, PU}}\pt$ term sums over charged particles originating
from pileup vertices in order to estimate and subtract the contribution of
pileup to the neutral particle sum, which is scaled by $1/2$ to account for
the fraction of neutral to charged energy in pileup interactions. The of the
muon is denoted by $\pt^{\Pgm}$. In the $\tau_{\Pgm}\tauh$ channel, it is
required that $I^{\Pgm}<0.15$.

###  0.5.3 Electron reconstruction

Electrons are reconstructed using tracks from the tracking system and
calorimeter deposits in the ECAL, with a veto on objects with a large HCAL to
ECAL energy ratio. Electrons are identified using a multivariate analysis
(MVA) discriminant combining several quantities that describe the shape of the
energy deposits in the ECAL, the quality of tracks, and the compatibility of
the measurements from the tracker and the ECAL [73]. The energy scale of
electrons is adjusted in data and simulation using the mass peak, while its
resolution in simulation is adjusted to data.

For the electrons, an isolation criterion $I^{\Pe}$ is defined for a cone size
of $R<0.3$ centred on the electron direction. Its definition is analogous to
Eq. (4) for the charged tracks, but the pileup contribution of neutral
particles is estimated via an effective-area method as

$I^{\Pe}=\frac{\sum_{\text{charged}}\pt+\max\left(0,\sum_{\text{neutral}}\pt-\rho\,\text{EA}\right)}{\pt^{\Pe}}.$
(5)

In this equation, the pileup contribution is estimated as $\rho\,\text{EA}$,
where $\rho$ is the event-specific average pileup energy density per unit area
in the $\phi$-$\eta$ plane and EA, which depends on the electron $\eta$, is
the effective area specific to the given type of isolation. In the
$\tau_{\Pe}\tauh$ channel, it is required that $I^{\Pe}<0.15$.

###  0.5.4 Jet and reconstruction

Jets are reconstructed using the anti-algorithm [71] with distance parameter
$R=0.4$ as implemented in the package [74]. The anti-algorithm functions by
taking PF objects and grouping them together based on inverse powers of the of
the objects [71, 75]. Jet momentum is determined as the vectorial sum of all
particle momenta in the jet, and is found from simulation to be, on average,
within 5 to 10% of the true momentum over the whole spectrum and detector
acceptance. Pileup interactions can contribute additional tracks and
calorimetric energy depositions to the jet momentum. To mitigate this effect,
charged particles identified to be originating from pileup vertices are
discarded and an offset correction is applied to correct for remaining
contributions. Jet energy corrections are derived from simulation to bring the
measured response of jets to that of particle level jets on average. In situ
measurements of the momentum balance in dijet, $\text{photon}+\text{jet}$,
$\PZ+\text{jet}$, and multijet events are used to account for any residual
differences in the jet energy scale between data and simulation [76]. The jet
energy resolution amounts typically to 15–20% at 30, 10% at 100, and 5% at 1
[76]. Additional selection criteria are applied to each jet to remove jets
potentially dominated by anomalous contributions from various subdetector
components or reconstruction failures. Data collected in the ECAL endcaps were
affected by large amounts of noise during the 2017 data-taking period, which
led to disagreements between simulation and data. To mitigate this issue, jets
used in the analysis of the 2017 data are discarded if they have $\pt<50\GeV$
and $2.65<\abs\eta<3.10$. Hadronic jets that contain -quarks (-jets) are
tagged using a deep neural network (DNN), called DeepCSV algorithm [77]. The
medium working point used for the DeepCSV algorithm corresponds to a -jet
identification efficiency of about 70% for a misidentification rate for jets
originating from light quarks and gluons of around 1%.

The pileup per particle identification algorithm [78] is applied to reduce the
pileup dependence of the observable. The and its magnitude () are computed
from the PF candidates weighted by their probability to originate from the PV
[79]. The is adjusted for the effect of jet energy corrections.

###  0.5.5 Tau lepton reconstruction

The lepton reconstruction is performed with the Hadron-Plus-Strips (HPS)
algorithm [80]. Starting from the constituents of reconstructed jets, the
algorithm works by combining charged hadrons with the signature of neutral
pions—one or more electron/photon candidates falling within a certain
$\Delta\eta{\times}\Delta\phi$ region (referred to as a “strip”). The
combination of these signatures provides the four-vector of the visible decay
products of the parent . The identification of candidates makes use of
isolation discriminators to reject quark and gluon jets that could be
misidentified as . For this analysis, a DNN called DeepTau [81] is used on the
HPS candidates to provide further discrimination. In order to achieve an
optimal identification performance, the DNN combines information from the
high-level reconstructed features together with the low-level information from
the inner tracker, calorimeters and muon sub-detectors, using PF candidates
reconstructed within the isolation cone. The working point on the output
discriminant is chosen to provide a identification efficiency of about 60% at
a jet misidentification rate of approximately $5\times 10^{-3}$. Two other
DNNs are used to reject electrons and muons misidentified as candidates using
dedicated criteria based on the consistency between the measurements in the
tracker, calorimeters, and muon detectors.

The mass of the $\PGt\PGt$ system $m_{\PGt\PGt}$ is calculated using a
simplified matrix-element algorithm, Svfit [82], which combines the and its
uncertainty matrix with the four-vectors of both candidates to calculate the
parent boson’s mass. The resolution of $m_{\PGt\PGt}$ is 15–20% depending on
the $\PGt\PGt$ final state and the boost of the $\PGt\PGt$ system.

##  0.6 Reconstruction of CP-sensitive observables

In this section we outline the methods used to construct CP-sensitive
observables, collectively referred to as $\phi_{{\textit{CP}}}$ angles.
Various techniques can be used to define $\phi_{{\textit{CP}}}$ depending on
the decay topology of the leptons. In total, four methods are employed in the
analysis: the ”impact parameter” [83, 84], ”neutral-pion” [85, 84], ”combined”
[84], and ”polarimetric vector” [86] methods. We provide a detailed
description of these methods below. We then summarise for which di-final
states each method is utilised, and outline the procedures used to optimise
the resolving power of the $\phi_{{\textit{CP}}}$ observables.

###  0.6.1 Impact parameter method

This method exploits the finite lifetime of the leptons and can be applied to
all events where both leptons decay to a single charged particle. We define
the impact parameter $\vec{j}^{\pm}$ of a track (where $\pm$ refers to the
charge of the track) as the vector between the PV and the point on the track
where distance to the PV is minimal.

For each lepton we define a plane using the impact parameter vector and the
charged-particle momentum vector. This plane, which is constructed in the
laboratory frame, only represents the genuine plane of the decay into a single
charged pion and neutrino when the laboratory frame coincides with the rest
frame of the . This means that this method does not reconstruct the genuine
lepton decay plane, but rather a plane that is correlated with it. In order to
approximate the rest frame of the we use the charged decay products of the
leptons of the to define a zero-momentum frame (ZMF) into which the decay
planes are boosted. The ZMF is used to define $\phi_{{\textit{CP}}}$ for all
channels in this analysis, except the
$\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$ channel, where both leptons decay
to three charged pions and the rest frame can be reconstructed.

We then construct four-component vectors in the laboratory frame as
$\lambda^{\pm}=(0,\vec{j}^{\pm})$. The $\lambda^{\pm}$ four-vectors are
boosted into the ZMF and denoted ${\lambda^{ZMF}}^{\pm}$. We also boost the
respective charged-pion four-vectors to the ZMF, denoted $q^{ZMF\pm}$.
Subsequently, we take the transverse components of ${\lambda^{ZMF}}^{\pm}$
with respect to $q^{ZMF\pm}$. We normalise the vectors to obtain unit vectors
$\hat{\lambda}^{ZMF+}_{\bot}$ and $\hat{\lambda}^{ZMF-}_{\bot}$.

To reconstruct $\phi_{{\textit{CP}}}$, we first define the angle $\phi^{ZMF}$
and $O^{ZMF}$ as

$\displaystyle\phi^{ZMF}$
$\displaystyle=\arccos(\hat{\lambda}^{ZMF+}_{\bot}\cdot\hat{\lambda}^{ZMF-}_{\bot}),\,\text{and}$
(6) $\displaystyle O^{ZMF}$
$\displaystyle=\hat{q}^{ZMF-}\cdot(\hat{\lambda}^{ZMF+}_{\bot}\times\hat{\lambda}^{ZMF-}_{\bot}).$

From $\phi^{ZMF}$ and $O^{ZMF}$ we reconstruct $\phi_{{\textit{CP}}}$ in a
range $[0,360^{\circ}]$ as

$\phi_{{\textit{CP}}}=\left\\{\begin{aligned} &\phi^{ZMF}&&\text{if
}O^{ZMF}\geq 0\\\ &360^{\circ}-\phi^{ZMF}&&\text{if
}O^{ZMF}<0\end{aligned}\right..$ (7)

The lepton spectral functions have opposite signs for single-pion decays and
leptonic decays in the kinematic regions considered in this analysis. This
causes a phase flip between the $\phi_{{\textit{CP}}}$ distributions for
single pion decays and leptonic decays when the impact parameter method is
used [40]. An illustration of the definition of the $\phi_{{\textit{CP}}}$
observable using the impact parameter method is shown in Fig. 3 (left).

Figure 3: Illustration of the lepton decay planes and the angle
$\phi_{{\textit{CP}}}$ for various decay configurations. The decay planes are
illustrated with the shaded regions, and either the vector $\hat{\lambda}$ or
the momentum vector of the neutral pion is in the decay plane. The
illustrations are in the frame in which the sum of the momenta of the charged
particles is zero. Left: the decay plane for the decays $\PGtm\to\Pgpm+\PGn$
and $\PGtp\to\Pgpp+\PAGn$. Middle: the decay plane as reconstructed from the
neutral and charged pion momenta. Right: $\phi_{{\textit{CP}}}$ for the mixed
scenario, in which one lepton decays to a pion while the other decays via an
intermediate meson.

###  0.6.2 Neutral-pion method

This method can be applied to hadronic decay channels in which both leptons
undergo decays involving more than one outgoing hadron. We describe the method
applied to the intermediate meson decay, and the intermediate meson to 1- and
3-prong decay modes.

For the meson decays, the vector $\lambda$ is replaced by the four-momentum
vector of the , which means we use the planes spanned by the decay products
(the and in the case of the $\PGr^{\pm}\to\Pgppm\Pgpz$ decay) to define the
$\phi_{{\textit{CP}}}$ observable. The four-momentum vector of the is obtained
as follows: to estimate the energy, we the sum of the energies of all
electron/photon candidates collected by the HPS algorithm. The direction of
the is then taken as the direction of the leading electron/photon candidate.
In most cases the leading candidate is a photon and the direction is
determined by pointing its associated ECAL clusters back towards the PV.
Finally, the mass is set to the known mass.

The same method is applied to $\mathrm{a_{1}^{1pr}}$ decays involving two
neutral pions by summing the neutral constituents in the decay, as they cannot
be easily resolved experimentally. The angle $\phi_{{\textit{CP}}}$ is then
calculated in an analogous method to that used in the impact parameter method
except that to avoid destructive interference from differently polarised
states of the mesons, the following observables need to be defined:

$y^{\PGtpm}=\frac{E_{\Pgppm}-E_{\Pgpz}}{E_{\Pgppm}+E_{E_{\Pgpz}}},\,y^{\PGt}=y^{\PGtm}y^{\PGtp}.$
(8)

In this equation, $E_{\Pgp}$ is the energy of the pion in the laboratory
frame. If $y^{\PGt}$ is negative, $\phi_{{\textit{CP}}}$ is obtained via the
shift $360^{\circ}-\phi_{{\textit{CP}}}$. The neutral-pion method can also be
successfully adapted to the $\mathrm{a_{1}^{3pr}}$ decay mode. In these decays
we select the oppositely charged pion pair with an invariant mass closest to
the intermediate $\PGr^{0}$, an illustration is depicted in Fig. 4. Of this
pair we treat the pion with the charge opposite of that of the lepton as
though it was a , and the momentum of the pion with the same sign as the is
used for the calculation of the ZMF. After these assignments the neutral-pion
method is applied as described for 1-prong decays.

Figure 4: The decay of $\mathrm{a_{1}^{3pr}}$ via an intermediate $\PGr^{0}$
to three charged pions.

An illustration of the definition of the $\phi_{{\textit{CP}}}$ observable
using the neutral-pion method is shown in Fig. 3 (middle).

###  0.6.3 Combined method

This method combines the impact parameter and neutral-pion methods outlined in
the two previous sections, which is appropriate for events where only one of
the two leptons decay into multiple hadrons. For the lepton decaying into a ,
$\mathrm{a_{1}^{1pr}}$, or $\mathrm{a_{1}^{3pr}}$ mesons, the vector $\lambda$
in Eq. (6) is replaced by four-momentum vectors as described in Section 0.6.2,
and the angle $\phi_{{\textit{CP}}}$ is then calculated using the same
formulae.

Analogously to the neutral-pion method we avoid destructive interference from
differently polarised states of the mesons by applying the shift
$360^{\circ}-\phi_{{\textit{CP}}}$ for events with $y^{\PGtpm}$, where
$y^{\PGtpm}$ is computed for the lepton that decays to the intermediate
resonance.

An illustration of the definition of the $\phi_{{\textit{CP}}}$ observable
using the combined method is shown in Fig. 3 (right).

###  0.6.4 Polarimetric vector method

This method can, in principle, be applied to any lepton decay mode in which
both lepton momenta can be well reconstructed. When leptons decay via the
$\mathrm{a_{1}^{3pr}}$ mode, the lepton rest frames can be reconstructed using
the secondary vertices (SVs), that are extracted by fitting the three tracks
originating from the $\mathrm{a_{1}^{3pr}}$ decays. Therefore, we only apply
the polarimetric vector method to the
$\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$ decay configuration.

The polarimetric vector $\vec{h}$ can be considered as an estimate of the most
likely direction of the spin vector $\vec{s}$ of the lepton in the lepton rest
frame [86]. We start by outlining the reconstruction of the lepton momenta,
which are required to compute the polarimetric vectors. Subsequently, we
describe how $\phi_{{\textit{CP}}}$ is reconstructed in the rest frame from
both the lepton momenta and the polarimetric vectors.

To reconstruct the lepton momentum in the $\mathrm{a_{1}^{3pr}}$ channel we
assume that the reconstructed lepton candidate has a mass $m_{\PGt}$ and
undergoes a two-body decay to a massless neutrino and an intermediate
$\mathrm{a_{1}}$ meson with mass $m_{\mathrm{a_{1}}}$. Furthermore, we define
the Gottfried–Jackson angle $\theta_{\mathrm{GJ}}$ as the angle between the
$\mathrm{a_{1}}$ momentum and the lepton momentum [87]. The latter is
reconstructed from the positions of the lepton production and decay vertex.
The magnitude of the lepton momentum is then given by [87]

$\abs{\vec{p}_{\PGt}}=\frac{(m_{\mathrm{a_{1}}}^{2}+m_{\PGt}^{2})\abs{\vec{p}_{\mathrm{a_{1}}}}\cos\theta_{\mathrm{GJ}}\pm\sqrt{\smash[b]{(m_{\mathrm{a_{1}}}^{2}+\abs{\vec{p}_{\mathrm{a_{1}}}}^{2})((m_{\mathrm{a_{1}}}^{2}-m_{\PGt}^{2})^{2}-4m_{\PGt}^{2}\abs{\vec{p}_{\mathrm{a_{1}}}}^{2}\sin^{2}\theta_{\mathrm{GJ}})}}}{2(m_{\mathrm{a_{1}}}^{2}+\abs{\vec{p}_{\mathrm{a_{1}}}}^{2}\sin^{2}\theta_{\mathrm{GJ}})}.$
(9)

The maximal allowed value $\theta_{\mathrm{GJ}}^{\text{max}}$ of the
Gottfried–Jackson angle is defined as

$\theta_{\mathrm{GJ}}^{\text{max}}=\arcsin\left(\frac{m_{\PGt}^{2}-m_{\mathrm{a_{1}}}^{2}}{2m_{\PGt}\abs{\vec{p}_{\mathrm{a_{1}}}}}\right).$
(10)

For decays in which the reconstructed $\theta_{\mathrm{GJ}}$ exceeds
$\theta_{\mathrm{GJ}}^{\text{max}}$, due to the finite angular resolution of
the charged pions and direction measurements, the value of
$\theta_{\mathrm{GJ}}$ is set to $\theta_{\mathrm{GJ}}^{\text{max}}$.

As can be seen in Eq. (9), there can be two solutions for the lepton momentum.
This can be understood by considering the decay in the lepton rest frame. In
this frame, the $\mathrm{a_{1}}$ meson may be emitted in either the same or
the opposite direction to that of the lepton momentum in the lab frame. When
the $\mathrm{a_{1}}$ meson is emitted in the direction orthogonal to the
lepton, we obtain the uniquely determined solution when the square root in the
numerator of Eq. (9) vanishes. Thus, we may obtain up to four pairs of
solutions for the momenta of the two leptons. This ambiguity is resolved by
selecting the pair of solutions with the mass closest to that of the . The
direction of the lepton in the lab frame is determined by the vector SV-PV.

Once the leptons and $\mathrm{a_{1}^{3pr}}$ momenta have been determined, the
polarimetric vectors $\vec{h}_{1,2}$ may be retrieved using the
$\mathrm{a_{1}^{3pr}}$ resonance model as implemented in the [88, 89, 90]
program, which uses the parameters as measured by the CLEO Collaboration [91].
To reconstruct $\phi_{{\textit{CP}}}$ from the polarimetric vectors and the
lepton momenta vectors, we introduce a vector $\vec{k}$ that is defined as

$\vec{k}_{1,2}=\frac{\vec{h}_{1,2}\times\vec{h}_{1,2}}{\abs{\vec{h}_{1,2}\times\vec{h}_{1,2}}}.$
(11)

In this definition, $\vec{h}_{1,2}$ are the two lepton momentum unit vectors
in the rest frame. We then reconstruct $\phi^{*}$ and $O^{*}$ (in the rest
frame) as

$\displaystyle\phi^{*}=\arccos(\vec{k}_{1}\cdot\vec{k}_{2}),\,\text{and}$ (12)
$\displaystyle O^{*}=-(\vec{h}_{1}\times\vec{h}_{2})\cdot\vec{h}_{1}.$

From $\phi^{*}$ and $O^{*}$ we reconstruct $\phi_{{\textit{CP}}}$ via the
assignments defined in Eq. (7).

In summary, for the configuration involving two $\mathrm{a_{1}^{3pr}}$ decays,
the secondary decay vertices are exploited to reconstruct the momenta in the
rest frame of the . Together with the $\mathrm{a_{1}}$ resonance model, this
allows for the extraction of the polarimetric vectors. Studies on simulated
signal events revealed that fits to $\phi_{{\textit{CP}}}$ measured using the
polarimetric vector method have approximately twice the resolving power
between the CP-even and CP-odd states as compared to applying the neutral-pion
method.

###  0.6.5 Strategy for selecting the CP-sensitive observables

The impact parameter is relatively small compared to the tracking resolution
and therefore the precision to which it can be measured is limited despite the
excellent resolution of the CMS tracker. An advantage of the neutral-pion
method is that it does not rely on the reconstruction of the impact parameter;
instead, the direction of the neutral pion needs to be determined. Due to the
relatively large distance between the primary interaction point and the ECAL
$\left(\mathcal{O}(1\unit{m})\right)$, coupled with the fine ECAL granularity,
the direction of neutral pions can be reconstructed with smaller relative
uncertainties compared to the impact parameter direction.

Studies were performed on signal events to review the CP sensitivity of the
neutral pion and impact parameter methods in regions of phase space where the
latter is expected to perform optimally. The sensitivity normalised to the
number of events was comparable while the selections (explained below) that
are needed for the impact parameter method discard a significant number of
events. The cuts imposed on the impact parameter significance mean that the
neutral-pion method can be applied to about twice as many events as the impact
parameter method. Therefore, although the impact parameter method can in
principle be applied to every lepton decay mode, in this analysis we only use
this method for the $\Pgp\Pgp$, $\Pgm\Pgp$, and $\Pe\Pgp$ final states. For
the $\PGr\PGr$, $\PGr\mathrm{a_{1}^{1pr}}$,
$\mathrm{a_{1}^{1pr}}\mathrm{a_{1}^{1pr}}$,
$\mathrm{a_{1}^{1pr}}\mathrm{a_{1}^{3pr}}$, and $\PGr\mathrm{a_{1}^{3pr}}$
final states, the neutral-pion method is deployed, and the polarimetric vector
method is used exclusively for the $\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$
channel. In other configurations where one lepton decays to a single charged
hadron or lepton and the other to multiple hadrons, the combined methods is
used.

###  0.6.6 Extraction of $\phi_{{\textit{CP}}}$ optimisation

In this section we outline the experimental techniques that are developed for
this analysis to optimise the experimental extraction of
$\phi_{{\textit{CP}}}$. A dedicated MVA discriminant is deployed to improve
the identification of the decay modes. To improve the estimate of the impact
parameters, the reconstruction of the PV coordinates is improved, and a
helical extrapolation of the track to the PV was implemented. These methods
are discussed in detail below.

#### Multivariate discriminant for decay mode identification

In order to optimally discriminate between the different decay modes, a
boosted decision tree (BDT) [92] is deployed. It is trained using the XGBoost
[93] framework, and is applied on top of the selection. The algorithm was
trained to distinguish between the 1- and 3-prong lepton decays: , ,
$\mathrm{a_{1}^{1pr}}$, $\mathrm{a_{1}^{3pr}}$, and $\Pgppm\PGpmp\Pgppm\Pgpz$.
The $\Pgppm\PGpmp\Pgppm\Pgpz$ decay is not used in the extraction of the CP
angle but must be separated from $\mathrm{a_{1}^{3pr}}$ to avoid
contamination.

The inputs to the BDT are the kinematic features of the reconstructed by the
HPS method and its constituents. The BDT exploits angular correlations between
the decay products, invariant mass quantities, and kinematic properties of the
photons.

#### Primary vertex refitting

The finite lifetime of the lepton means that tracks emanating from its decay
do not originate from the PV. These tracks are removed and the PV is refitted
using the remaining tracks as input to the AVF algorithm. The LHC beamspot
represents a three-dimensional (3-D) profile of the luminous region, where the
LHC beams collide in the CMS detector. The parameters of the beamspot are
determined from an average over many events [68]. The uncertainties in the
beamspot parameters are relatively small and are incorporated into the AVF
algorithm to provide an additional constraint on the PV position. The
inclusion of the beamspot information leads to an improvement of the PV
resolution in the transverse plane of a factor of about 3 for signal events
and of about 4 for Drell–Yan events, while the $z$ coordinate of the PV is
largely unaffected. This refitted PV is used when estimating the impact
parameters and the polarimetric vectors.

#### Impact parameter estimate and significance

A dedicated algorithm is deployed to derive the impact parameter of the
charged track from the lepton decay using an analytic extrapolation of track
trajectory towards the PV position. The extrapolation depends on the magnetic
field and the helical parameters of the track. The distance between the
extrapolated track and the PV position is then minimised numerically to
determine the impact parameter.

This procedure has two advantages. Firstly, with this extrapolation, the
minimisation of the impact parameter is performed in three dimensions. For
tracks with large $\eta$ values, the procedure leads to a better estimation of
the $z$ coordinate of the impact parameter than when the minimisation is done
exclusively in the transverse plane. Secondly, the helical extrapolation
allows for the propagation of both the track and PV uncertainties into an
overall impact parameter significance $S_{\text{IP}}$ (defined as the ratio of
the magnitude of the impact parameter divided by its uncertainty). In this
analysis, selections are made on the impact parameter significance, as further
explained in Section 0.9.

##  0.7 Event selection

Events are selected online by the CMS trigger system. For the
$\tau_{\ell}\tauh$ channels, events are triggered by either a paired
$\ell+\tauh$ cross trigger or a single-lepton trigger with a higher threshold
for the lepton compared to the cross trigger. For the $\tauh\tauh$ channel, a
di-trigger is used.

Offline, a pair of oppositely charged leptons separated by $\DR>0.5$ is
required. The offline-reconstructed objects must match the required trigger
objects (the object as reconstructed by the trigger system) within $\DR<0.5$.
The offline-reconstructed light lepton is required to have a value that is at
least 1higher than the online threshold. If an offline candidate is matched to
a trigger object (including the leg of the $\ell+\tauh$ cross trigger for the
semileptonic channels), the must have a at least 5above the trigger threshold.
The offline thresholds are higher than the online thresholds due to the turn-
on curve in the trigger efficiencies.

Table 0.7 summarises the online trigger and offline thresholds for 2016–2018.
The offline requirements apply only to objects that are matched to a trigger
object. If, in the $\tau_{\ell}\tauh$ channels, the event is selected online
by the single-lepton trigger, the offline is required instead to have a of at
least 20.

Kinematic trigger and offline requirements applied to the $\tau_{\Pe}\tauh$,
$\tau_{\Pgm}\tauh$, and $\tauh\tauh$ channels. The trigger requirement is
indicated in parentheses (in ). The thresholds indicated for the apply only
for the object matched to the hadronic trigger or to the hadronic leg from the
cross trigger.  Channel Year Trigger requirement Offline () $\tauh\tauh$ All
years $\tauh(35)\,\&\,\tauh(35)$ $\pt^{\tauh}>40$ $\tau_{\Pgm}\tauh$ 2016
$\Pgm(22)$, $\Pgm(19)\,\&\,\tauh(20)$ $\pt^{\Pgm}>20$, $\pt^{\tauh}>25$ 2017,
2018 $\Pgm(24)$, $\Pgm(20)\,\&\,\tauh(27)$ $\pt^{\Pgm}>21$, $\pt^{\tauh}>32$
$\tau_{\Pe}\tauh$ 2016 $\Pe(25)$ $\pt^{\Pe}>26$ 2017 $\Pe(27)$,
$\Pe(24)\,\&\,\tauh(30)$ $\pt^{\Pe}>25$, $\pt^{\tauh}>35$ 2018 $\Pe(32)$,
$\Pe(24)\,\&\,\tauh(30)$ $\pt^{\Pe}>25$, $\pt^{\tauh}>35$

For the $\tau_{\ell}\tauh$ channels, the large $\PW+\text{jets}$ background is
reduced by rejecting events based on the transverse mass of the light lepton
and system,

$\mT\equiv\sqrt{{2\pt^{\ell}\ptmiss[1-\cos(\Delta\phi)]}}<50\GeV,$ (13)

where $\Delta\phi$ is the azimuthal angle between the direction of the light
lepton and .

The longitudinal and transverse impact parameters $d_{\mathrm{z}}$ and
$d_{\mathrm{xy}}$ of the muon and electron are required to satisfy
$\abs{d_{\mathrm{z}}}<0.2\cm$ and $\abs{d_{\mathrm{xy}}}<0.045\cm$. These
impact parameters originate from a minimisation of the magnitude of the impact
parameters in the transverse plane only, in contrast to the impact parameters
used for calculating $\phi_{{\textit{CP}}}$, which are derived using a 3-D
minimisation. For the purpose of the event selection this factorised approach
is sufficiently precise. For the leading track, only the requirement
$\abs{d_{\mathrm{z}}}<0.2\cm$ is imposed to avoid loss of selection
efficiency. Further, a veto on events containing loosely identified additional
electrons or muons is imposed. For the $\tau_{\ell}\tauh$ channels, a veto on
jets passing -tagging requirements is also applied. When multiple lepton pairs
are present, the pairs are ranked based on the output scores of the DeepTau
algorithm for the candidates, and the relative isolation for the $\tau_{\ell}$
candidates. The highest ranked pair is selected.

##  0.8 Background estimation

The processes that contribute to the background in this analysis are
$\PZ+\text{jets}$, $\PW+\text{jets}$, top quark-antiquark pair production (),
single top quark, and diboson production. Additionally, events comprised
uniquely of jets produced through the strong interaction, referred to as QCD
multijet events, form a significant background. These processes contribute to
the production of genuine leptons, jets and leptons that are misidentified as
, as well as prompt leptons and jets that are misidentified as $\tau_{\ell}$
in the semileptonic channels. All background processes resulting in two
genuine leptons constitute a major background, and are estimated from data
using a -embedding technique [94]. The majority of the backgrounds due to jets
misidentified as candidates are estimated using the ”fake factor”
($F_{\text{F}}$) method, as described in Ref. [95]. The remaining minor
backgrounds are determined using simulated events. In the remainder of this
section we outline the -embedding and $F_{\text{F}}$ methods, and describe the
corrections applied to simulated events in order to improve their description
of the data.

###  0.8.1 The embedding method

In order to obtain the genuine $\PGt\PGt$ background we exploit lepton
universality, and replace oppositely charged muon pairs in data events with
simulated oppositely charged lepton pairs. The dominant process for this
background is $\PZ\to\PGt\PGt$, but there are also small contributions from
and diboson processes.

For all data-taking periods, events containing an oppositely charged dimuon
pair were collected using a dedicated di-trigger. The detector hits belonging
to the muon tracks are removed from these events. A boson is simulated in an
empty detector, which is forced to decay to a pair of oppositely charged
leptons with identical kinematics to the muon pair that was removed. The
leptons are forced to decay fully hadronically or semileptonically in order to
simulate either the $\tauh\tauh$ or $\tau_{\ell}\tauh$ channels. The detector
response to the pair is then simulated and added to the data event.

In order to model the background processes in data well, various corrections
need to be applied to the embedded event samples. Muons and electrons are
corrected for mismodelling of their trigger, tracking and identification, and
isolation requirements efficiencies. The candidates are corrected for
mismodelling of their trigger, reconstruction and identification efficiencies.
The ”tag-and-probe” method [96] is used to derive these corrections. The
energy scales are corrected per decay mode to match the corresponding scales
in data. The electron, muon, and pion impact parameters are corrected using a
samples of $\PZ\to\Pgm\Pgm$ and $\PZ\to\Pe\Pe$ events and quantile-mapping
techniques.

###  0.8.2 The $F_{\text{F}}$ method

This method is designed to provide an estimate of the shape and normalisation
of events in which at least one quark or gluon jet is misidentified as a
lepton based on control samples in data. We refer to such a jet as a
$\text{jet}\to\tauh$.

We define a determination region that is orthogonal to the signal region and
dominated by a background process resulting in $\text{jet}\to\tauh$
misidentifications; the construction of these regions is outlined below. We
define a nominal ID as a object that passes nominal ID requirements as
outlined in Section 0.5, and a relaxed ID as objects that pass a looser
requirement on the DNN output but fail the nominal ID. In this determination
region we obtain the ratio between the nominal ID rate and the relaxed ID
rate. The ratio in the determination region is the $F_{\text{F}}$. To obtain
the rate of misidentified jets in the signal region, an application region is
defined by selecting events that fulfil all event selection criteria except
that they contain a lepton that passes the relaxed instead of the nominal
requirement (for the $\tauh\tauh$ channel it must be the leading ). The rate
of misidentified jet events in the signal region is obtained by applying the
$F_{\text{F}}$ values from the determination region on an event-by-event basis
as an event weight to the events in the application region.

The $\text{jet}\to\tauh$ background in the $\tauh\tauh$ channel originates
almost entirely from QCD multijet events. The determination region is thus
defined by inverting the opposite-sign requirement on the lepton pair to a
same-sign requirement, which effectively selects a control region pure in QCD
jets. The $F_{\text{F}}$ are parameterised for the leading lepton as a
function of the of the , and binned in the reconstructed decay mode, jet
multiplicity, and impact parameter significance. Correction factors are
derived using control regions in data to correct for residual differences in
the spectrum, and to account for the sign inversion used to define the
determination region. The final $F_{\text{F}}$ value for the channel is
obtained by applying the raw $F_{\text{F}}$ and the two corrections
multiplicatively. This $F_{\text{F}}$ also accounts for other processes with a
jet misidentified as the leading lepton, such as $\PW+\text{jets}$ production.
The events in which the subleading is a misidentified jet and the leading
candidate is a genuine lepton are modelled via simulation; these events
constitute only a small fraction $(\mathcal{O}(2\%))$ of the total
misidentified jet background in the $\tauh\tauh$ channel.

In the $\tau_{\ell}\tauh$ channels, the $\PW+\text{jets}$ process, and to a
lesser extent the process, contribute to jet misidentification as well as
events originating from QCD multijet production. Therefore, separate
$F_{\text{F}}$ are derived for these processes, and these individual
$F_{\text{F}}$ values are subsequently weighted into an overall $F_{\text{F}}$
to account for their different contributions in the application region.
Simulated events are used to determine the expected relative contributions of
$\PW+\text{jets}$ and events, and the QCD contribution is estimated by
subtracting all simulated non-QCD processes from the data in the application
region. The overall $F_{\text{F}}$ accounts for the jet misidentification in
all background processes. The procedure for the QCD $F_{\text{F}}$ is similar
to the method described for the $\tauh\tauh$ channel, except that the light
lepton isolation parameter must be larger than 0.05 to reduce processes
resulting in genuine leptons. Correction factors are derived to correct for
residual differences in the lepton and spectra, and to account for the sign
inversion and additional lepton isolation requirement used to define the
determination region. A determination region sufficiently pure in
$\PW+\text{jets}$ is defined by selecting events with $\mT>70\GeV$. Correction
factors are derived to correct for residual differences in the lepton and
spectra, and to account for the inverted selection used to define the
determination region. For the process, it is difficult to define a
sufficiently pure region in data, and thus the $F_{\text{F}}$ values are
estimated from a simulated sample. Correction factors are derived to correct
for residual differences in the spectrum, and to account for differences in
the $F_{\text{F}}$ in data and simulation. The latter is derived by comparing
the $F_{\text{F}}$ values measured for $\PW+\text{jets}$ in data and
simulation.

###  0.8.3 Estimation of minor backgrounds

The -embedding and $F_{\text{F}}$ methods combined describe around 90% of the
backgrounds in this analysis. All events containing a genuine lepton pair are
taken from the embedded samples, while events in which the (leading) is a
misidentified hadronic jet in the ($\tauh\tauh$) $\tau_{\ell}\tauh$ channels
are obtained from data using the $F_{\text{F}}$ method. All other background
events are obtained from simulation.

In addition to the genuine and $\text{jet}\to\tauh$ contributions to the
selected pairs, there are additional sources of misidentifications that may
occur. This includes prompt leptons that may be either misidentified as a
$\tau_{\ell}$ or as a , $\tau_{\ell}$ leptons being misidentified as , and
jets misidentified as $\tau_{\ell}$ candidates. In Tables 0.8.3 and 0.8.3 we
summarise the different background composition configurations and their
modelling for the $\tauh\tauh$ and $\tau_{\ell}\tauh$ channels, respectively.
To avoid double-counting events with a genuine lepton pair, such events are
subtracted from all simulated samples, as well as events in which the is a
misidentified hadronic jet (for the $\tauh\tauh$ channel this must be the
leading ).

The different sources of backgrounds in the $\tauh\tauh$ channel are shown in
the rows and columns. The entries in the table represent the possible lepton
pair background contribution from different processes and misidentifications
and encapsulate the different experimental techniques that are deployed to
estimate the background contributions.  Leading  Subleading  Genuine  Jet
$\to\tauh$ (Prompt lepton/$\tau_{\ell}$) $\to\tauh$ Genuine  -Embedding
Simulation Simulation Jet $\to\tauh$ $F_{\text{F}}$ $F_{\text{F}}$
$F_{\text{F}}$ (Prompt lepton/$\tau_{\ell}$) $\to\tauh$ Simulation Simulation
Simulation

The different sources of backgrounds in the $\tau_{\ell}\tauh$ channel are
shown in the rows and columns. The entries in the table represent the possible
lepton pair background contribution from different processes and
misidentifications and encapsulate the different experimental techniques that
are deployed to estimate the background contributions.  $\tau_{\ell}$ Genuine
Jet $\to\tauh$ (Prompt lepton/$\tau_{\ell}$) $\to\tauh$ Genuine $\tau_{\ell}$
-Embedding $F_{\text{F}}$ Simulation Jet $\to\tau_{\ell}$ Simulation
$F_{\text{F}}$ Simulation Prompt lepton $\to\tau_{\ell}$ Simulation
$F_{\text{F}}$ Simulation

In order to model the background processes in data well, various corrections
need to be applied to the simulated samples. All corrections to the lepton
decay products applied to the embedded samples (described in Section 0.8.1)
are also applied to the simulated samples. The corrections are known to differ
between the embedded samples and the corresponding simulated samples, and
therefore dedicated correction factors are derived in each case.

Jet energy scale corrections are applied to both data and simulated events as
described in Section 0.5. Recoil corrections to the are applied to reduce the
mismodelling of the simulated $\PZ+\text{jets}$, $\PW+\text{jets}$, and Higgs
boson samples. The corrections are applied to the simulated events based on
the vectorial difference of the measured and total of the neutrinos
originating from the decay of the , , or . Their average effect is the
reduction of the obtained from simulation by a few . Recoil corrections to are
measured in $\PZ\to\Pgm\Pgm$ events. The corrections are subsequently applied
to Drell–Yan plus jets events, $\PW+\text{jets}$, and signal event samples.
The $\ell\to\tauh$ misidentification rates are corrected in simulation by
applying the tag-and-probe method to $\PZ\to\ell\ell$ events, and the energy
scales are corrected in simulation to match the scale in data.

The boson mass and spectra in simulation are corrected to better match the
data. To this purpose the mass and are measured in data and simulation in di-
muon events, and the simulated events are corrected to match the spectra in
data. A correction is also applied to the top quark spectrum in the sample,
using a dedicated control region. The procedure used to derive this correction
is detailed in Ref. [97].

After applying all corrections, we obtain a satisfactory description of the
observables that we use to categorise events, which are described in Section
0.9.

###  0.8.4 Validating the modelling of the $\phi_{{\textit{CP}}}$ observables

To validate the modelling of the $\phi_{{\textit{CP}}}$ spectrum, the
$\Pgm\to\tauh$ events in data and the background estimates are divided into
distributions in which the charged is ”nearly coplanar” or ”nearly
perpendicular” to the production plane of the beam axis and the momentum in
the laboratory frame, as described in Ref. [37]. To this purpose we introduce
the observable $\alpha^{\pi}_{-}$ that is defined as

$\cos{\alpha^{\pi}_{-}}=\left|\frac{\hat{z}\times\hat{p}}{\abs{\hat{z}\times\hat{p}}}\cdot\frac{\hat{z}\times\hat{j}}{\abs{\hat{z}\times\hat{j}}}\right|.$
(14)

In this equation, $\hat{z}$ is the unit vector pointing along the beam axis,
$\hat{p}$ is the unit momentum vector of the charged , and $\hat{j}$ is the
unit impact parameter vector. We can define a subset of events in which the
charged is nearly perpendicular or coplanar by requiring
$\alpha^{\pi}_{-}>\pi/4$ or $\alpha^{\pi}_{-}<\pi/4$, respectively. We also
perform equivalent checks for decays into mesons, where we substitute the unit
momentum vector for $\hat{j}$ in Eq. (14) to define an equivalent observable,
$\alpha^{\rho}_{-}$. In Fig. 5 we display the coplanar and perpendicular
distributions in the $\Pgm\Pgp$ and $\Pgm\PGr$ channels.

Figure 5: The angle $\phi_{{\textit{CP}}}$ for $\Pgm\to\tauh$ events in which
the decays to a charged (upper) or a charged meson (lower). The distributions
are decomposed in a subset in which the charged is ”nearly coplanar” (left) or
”nearly perpendicular” (right) to the production plane.

##  0.9 Event categorisation

In order to enhance the sensitivity of this analysis we apply MVA
discriminants to separate signal from background events.

This event categorisation is formulated as a multi-class problem. The output
is a set of scores for each event (one per class) that, by construction, sum
to unity. Each score can therefore be interpreted loosely as the probability
that an event belongs to a given class. We then assign each event to a
category depending on the class that received the highest score. Since both
the $\tau_{\ell}\tauh$ and $\tauh\tauh$ channels are dominated by backgrounds
containing contributions from genuine and $\text{jet}\to\tauh$ production, the
discriminant is trained to categorise events in three classes:

  * •

The “Higgs” category is trained to distinguish events from the $\Pg\Pg\PH$,
VBF, and $\PV\PH$ samples from background events, which are reweighed by their
cross sections before merging them into one sample. Events in this category
are used to infer the CP quantum number of the boson.

  * •

The “Genuine” category includes all background processes involving two genuine
leptons.

  * •

The “Mis-ID” category includes all background processes in which minimally one
hadronic jet is misidentified as a lepton. This category also contains
$\ell\to\tauh$ misidentified events for all channels and prompt light leptons
in the $\tau_{\ell}\tauh$ channels.

The three categories are mutually exclusive and, by definition, the lower
bound for the highest MVA score is 1/3. Subsequently, the three training
categories are normalised to account for the different number of events in
each data set. All event classes are then chosen to contribute to the training
with the same weight, with uniform prevalence. For the semileptonic channels,
the backgrounds for the training are provided from simulated samples, except
for QCD events, which are obtained using same-sign lepton pair candidates in
data. For the hadronic channel, the embedded samples and the $F_{\text{F}}$
method are used in the training. For the latter, the events from the
application region are used and reweighted by their $F_{\text{F}}$ values. A
separate training was performed for each year to account for differences in
the performance of the CMS detector in different data-taking periods. In the
$\tau_{\ell}\tauh$ channels, the event categorisation is performed with a
multiclass neural network. In the $\tauh\tauh$ channel, the event
categorisation is performed using a multi-class BDT algorithm combined with
the XGBoost package. The input variables used in the categorisation of the
$\tau_{\ell}\tauh$ and $\tauh\tauh$ channels are displayed in Table 0.9. The
training is performed inclusively for all the lepton decay modes.

Events are sorted into the three categories depending on which of the three
output scores is closest to unity. The maximum output score is also retained
and used for the purpose of signal extraction. These maximum scores will be
referred to as the ”MVA scores” henceforth.

After the categorisation, a cutoff of $S_{\text{IP}}>1.5$ is applied to the
impact parameter significances of the electron and muon, as well as to the
single pions from a lepton for events that are classified as signal events.
Events with a lower $S_{\text{IP}}$ would dilute the sensitivity of the
analysis. In the background categories, a cutoff on the impact parameter
significance is only applied to the single-pion decays.

Input variables to the MVA discriminants for the $\tau_{\ell}\tauh$ and
$\tauh\tauh$ channels. The Svfit algorithm is used to estimate the di-mass.
Observable $\tau_{\ell}\tauh$ $\tauh\tauh$ of leading  $\checkmark$
$\checkmark$ of trailing  $\checkmark$ of $\tau_{\ell}$ $\checkmark$ of
visible di- $\checkmark$ $\checkmark$ of di-+  $\checkmark$ of
$\tau_{\ell}\tauh$ \+  $\checkmark$ Visible di-mass $\checkmark$ $\checkmark$
Di-mass (using Svfit) $\checkmark$ $\checkmark$ Leading jet  $\checkmark$
$\checkmark$ Trailing jet  $\checkmark$ Jet multiplicity $\checkmark$
$\checkmark$ Dijet invariant mass $\checkmark$ $\checkmark$ Dijet
$\checkmark$ Dijet $\abs{\Delta\eta}$ $\checkmark$ $\checkmark$ $\checkmark$

In Fig. 6, the post-fit MVA score distributions of the Genuine and Mis-ID
categories are displayed for the $\tau_{\Pgm}\tauh$ and $\tau_{\Pe}\tauh$
channels. The best fit signal contributions are overlaid. The fitting
procedure is outlined in Section 0.12. The genuine di-and $\text{jet}\to\tauh$
background contributions are displayed separately as indicated in the legends.
The remaining background contributions are collated and indicated by the
”Others” label. The BDT scores for the $\tauh\tauh$ channel are analogously
displayed in Fig. 7.

Figure 6: The post-fit MVA score distributions for the Genuine (left) and Mis-
ID categories (right) in the $\tau_{\Pgm}\tauh$ (upper) and $\tau_{\Pe}\tauh$
(lower) channels. The distributions are inclusive in decay mode. The best fit
signal distributions are overlaid. In the lower panels, the data minus the
background template divided by the uncertainty in the background template is
displayed, as well as the signal distribution divided by the uncertainty in
the background template. The uncertainty band accounts for all sources of
systematic uncertainty in the background prediction, after the fit to data.

Figure 7: The post-fit MVA score distributions for the Genuine (left) and Mis-
ID categories (right) in the $\tauh\tauh$ channel. The distributions are
inclusive in decay mode. The best fit signal distributions are overlaid. In
the lower panels, the data minus the background template divided by the
uncertainty in the background template is displayed, as well as the signal
distribution divided by the uncertainty in the background template. The
uncertainty band accounts for all sources of systematic uncertainty in the
background prediction, after the fit to data.

##  0.10 The $\phi_{{\textit{CP}}}$ distributions in windows of the MVA
discriminant score

The MVA score distributions described in Section 0.9 allow for a partial
separation of signal from background events. The $\phi_{{\textit{CP}}}$
distributions of the events in the signal categories are then analysed in
windows of increasing MVA score, corresponding to progressively higher signal-
to-background ratios. The result is a set of 2-D distributions built from the
MVA score and $\phi_{{\textit{CP}}}$ variables. These distributions are used
in the fit to data to extract the results.

The statistical fluctuations in the estimates of the background contributions
(denoted as background templates) in the signal and background categories are
sizeable. It has been underlined that backgrounds involving two genuine
leptons are flat in $\phi_{{\textit{CP}}}$ at the generator level [37].
Experimental smearing effects do not modulate this flat shape for decay modes
in which we apply the neutral pion method for at least one lepton. Therefore,
for this background process and these decay modes we flatten the background
templates by merging the bins. The $\phi_{{\textit{CP}}}$ distribution is not
flat for the $\text{jet}\to\tauh$ background for all decay modes due to the
kinematic properties of the events, but the distributions are still symmetric
around $\phi_{{\textit{CP}}}=180^{\circ}$, and so this background is
symmetrised—meaning the symmetry around $\phi_{{\textit{CP}}}=180^{\circ}$ is
enforced. For other background templates, for example the $\Pgm\to\tauh$
contribution, the distributions are found to be flat within the statistical
uncertainties, and therefore these background templates are also flattened.

The backgrounds are not expected to be flat in decay modes in which the impact
parameter method is used or when the polarimetric vector method is applied
when both leptons decay via the $\mathrm{a_{1}^{3pr}}$ mode. This can be
understood from the fact that smearing effects in the PV are correlated for
the decay planes. The smearing of the PV results in a depletion in the region
$\phi_{{\textit{CP}}}=180^{\circ}$ [37], such that the shape of the background
distributions in the $\Pgm\Pgp$, $\Pe\Pgp$, and
$\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$ ($\Pgp\Pgp$) channels tends to
resemble the CP-even (CP-odd) signal rather than the CP-odd (CP-even).
However, for such channels the backgrounds are symmetric around
$\phi_{{\textit{CP}}}=180^{\circ}$, and therefore the background templates are
symmetrised.

For certain decay modes, the statistical fluctuations in the signal templates
are also sizeable. Therefore, the templates for the scalar and pseudoscalar
cases are symmetrised around $\phi_{{\textit{CP}}}=180^{\circ}$ as well. The
maximally mixed signal template, which is not displayed in the plots, is used
in the fitting procedure described in Section 0.12. In order to symmetrise
this template, we reweight the signal sample to another sample with
$\alpha^{\PH\tau\tau}=-45^{\circ}$. The $\phi_{{\textit{CP}}}$ distribution is
shifted by $180^{\circ}$, and the average is taken between the sample with
$\alpha^{\PH\tau\tau}=45^{\circ}$ and $\alpha^{\PH\tau\tau}=-45^{\circ}$.

In Figs. 8–10 we display the post-fit data and background template
distributions, after the bin smearing and symmetrisation, with the best fit
and pseudoscalar signal templates overlaid. The uncertainties have been
adjusted to their value after the fit described in Section 0.12. The most
sensitive decay modes of the analysis are displayed, which are the $\Pgm\PGr$
and $\Pgm\Pgp$ mode in the $\tau_{\Pgm}\tauh$ channel displayed in Fig. 8, the
$\PGr\PGr$ and $\Pgp\PGr$ mode in the $\tauh\tauh$ channel displayed in Fig.
9, and the $\Pe\PGr$ and $\Pe\Pgp$ channels in the $\tau_{\Pe}\tauh$ channel
displayed in Fig. 10. The distributions highlight the effectiveness of the MVA
discriminant in optimising the signal over background ratio, as well as the
CP-sensitivity of the measurement that follows from the visibly different
phases of the best fit signal and CP-odd signal distributions. The
$180^{\circ}$ phase shift between the $\tauh\tauh$ and $\tau_{\ell}\tauh$
channels is manifest in the figures. The correlated effect of the PV smearing
is also visible in the $\Pgm\Pgp$ and $\Pe\Pgp$ modes via the non-flat shapes
of the background distributions.

Figure 8: Distributions of $\phi_{{\textit{CP}}}$ in the $\Pgm\PGr$ (upper)
and $\Pgm\Pgp$ (lower) channels in windows of increasing MVA score, shown on
top of each window. The best fit and pseudoscalar (PS) signal distributions
are overlaid. The $x$-axis represents the cyclic bins in
$\phi_{{\textit{CP}}}$ in the range of $(0,360^{\circ})$. In the lower panels,
the data minus the background template divided by the uncertainty in the
background template is displayed, as well as the signal distributions divided
by the uncertainty in the background template. The uncertainty band accounts
for all sources of systematic uncertainty in the background prediction, after
the fit to data.

Figure 9: Distributions of $\phi_{{\textit{CP}}}$ in the $\PGr\PGr$ (upper)
and $\Pgp\PGr$ (lower) channels in windows of increasing MVA score, shown on
top of each window. The best fit and pseudoscalar (PS) signal distributions
are overlaid. The $x$-axis represents the cyclic bins in
$\phi_{{\textit{CP}}}$ in the range of $(0,360^{\circ})$. In the lower panels,
the data minus the background template divided by the uncertainty in the
background template is displayed, as well as the signal distributions divided
by the uncertainty in the background template. The uncertainty band accounts
for all sources of systematic uncertainty in the background prediction, after
the fit to data.

Figure 10: Distributions of $\phi_{{\textit{CP}}}$ in the $\Pe\PGr$ (upper)
and $\Pe\Pgp$ (lower) channels in windows of increasing MVA score, shown on
top of each window. The best fit and pseudoscalar (PS) signal distributions
are overlaid. The $x$-axis represents the cyclic bins in
$\phi_{{\textit{CP}}}$ in the range of $(0,360^{\circ})$. In the lower panels,
the data minus the background template divided by the uncertainty in the
background template is displayed, as well as the signal distributions divided
by the uncertainty in the background template. The uncertainty band accounts
for all sources of systematic uncertainty in the background prediction, after
the fit to data.

##  0.11 Systematic uncertainties

The uncertainties considered in this analysis can be categorised into
normalisation and shape uncertainties. The former modify only the
normalisation of a distribution while leaving its shape unchanged, whereas the
latter allow for correlated changes across bins that also alter the shapes of
the distributions. The uncertainties are accounted for as nuisance parameters
in the fit to data. The normalisation and shape uncertainties are summarised
in Table 0.11, in which we also state their correlations between the three
different years of data-taking considered in this analysis.

Overview of the systematic uncertainties. The third column indicates if the
source of uncertainty was treated as being correlated between the years in the
fit described in Section 0.12. The forth column indicates if the uncertainty
affects the shapes of the distributions.  Uncertainty Magnitude Correlation
Shape Muon reconstruction 1% Yes No Electron reconstruction 2% Yes No Muon
trigger 2% No No Electron trigger 2% No No $\Pe\to\tauh$ rate in $\tauh\tauh$
10 (2)% 2016 (2017,2018) No No $\Pe\to\tauh$ rate in $\tau_{\Pe}\tauh$ 10% No
No $\Pgm\to\tauh$ rate in $\tau_{\Pgm}\tauh$ up to 40% No No -jet veto 1–9% No
No Luminosity 2.3–2.5% Partial No Embedded yield 4% No No cross section 4.2%
Yes No Diboson cross section 5% Yes No Single top quark cross section 5% Yes
No $\PW+\text{jets}$ cross section 4% Yes No Drell–Yan cross section 2% Yes No
cross sections 2–5% [98] Yes No $\PH\to\PGt\PGt$ branching fraction 2% [98]
Yes No SV reco. eff. in $\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$ 2% No No ID
efficiency 3% No No $S_{\text{IP}}$ In , , and decays Decay-mode dependent,
1–5% No No Muon energy scale 0.4–2.7% Yes Yes Electron energy scale $<$1% Yes
Yes Trigger /Decay-mode dependent No Yes Reconstruction /Decay-mode dependent
(2–3%) Partial Yes Top quark reweighing $\pt^{\text{top}}$-Dependent Yes Yes
and mass reweighing $\pt^{\PZ}/m_{\PZ}$-Dependent Partial Yes Energy scale
/Decay-mode dependent (0.2–1.1%) No Yes $\Pe\to\tauh$ Energy scale 0.5–6.5% No
Yes $\Pgm\to\tauh$ Energy scale 1% No Yes Jet energy scale Event-dependent
Partial Yes Jet energy resolution Event-dependent No Yes Unclustered scale
Event-dependent No Yes Recoil corrections Event-dependent No Yes
$F_{\text{F}}$ uncertainties Described in text Partial Yes /diboson in
embedded 10% Yes Yes L1 trigger timing (2016–2017) Event-dependent (0–4%) Yes
Yes Renorm./Fact. scales Event-dependent Yes Yes Parton showering Event-
dependent Yes Yes

###  0.11.1 Normalisation uncertainties

The integrated luminosity uncertainty amounts to 2.5, 2.3, and 2.5% for 2016,
2017, and 2018 respectively [99, 100, 101], and is applied to all simulated
samples discussed in Section 0.4.

The uncertainty in the muon reconstruction efficiency including the tracking,
identification, and isolation requirements is 1%, while for electrons it is
2%. The uncertainty in the muon and electron trigger efficiencies, which
affect both the single-lepton and cross-triggers, is 2%. An additional
normalisation uncertainty of 4% is applied to the embedded event samples,
originating from the uncertainty in the measurement of the muon trigger and
identification efficiencies used to scale the embedded samples.

For the $\tau_{\ell}\tauh$ channels, which contain a veto on events containing
-jets, an uncertainty in the propagation of the -quark tagging scale factors
of 1–9% is applied on the and diboson event yields (the uncertainties on the
event yields for other simulated processes are found to be negligible).

The 3.1 program [102] was used to calculate the $\PW+\text{jets}$ and
$\PZ+\text{jets}$ cross sections. Uncertainties in the factorisation and
renormalisation scales, the PDF, and the running coupling were propagated and
added in quadrature. The Top++v2.0 program [103] was used to calculate the
cross section and its uncertainty. The extracted uncertainties for the
simulated $\PZ+\text{jets}$, $\PW+\text{jets}$, and cross sections amount to
2, 4, and 4%, respectively. For the diboson and single top quark production
processes, a combined systematic uncertainty in the background yield is
estimated to be 5% using CMS measurements [104, 105]. The uncertainties in the
signal $\Pg\Pg\PH$, VBF, and $\PV\PH$ production cross sections, as well as
the uncertainty in the $\PH\to\PGt\PGt$ branching fraction, are applied as
recommended in Ref. [98].

The uncertainty in the $\Pgm\to\tauh$ misidentification rate in the
$\tau_{\Pgm}\tauh$ channel is split into four independent uncertainties
depending on the MVA decay mode of the $\Pgm\to\tauh$ candidate. The sizes of
the uncertainties are 20% for and , 30% for $\mathrm{a_{1}^{1pr}}$, and 40%
for $\mathrm{a_{1}^{3pr}}$, respectively. An uncertainty of 10 (2)% to the
$\Pe\to\tauh$ misidentification rate is applied for 2016 (2017, 2018) in the
$\tauh\tauh$ channel. In the $\tau_{\Pe}\tauh$ channel, the $\Pe\to\tauh$
misidentification rate is split per decay mode and is, at most, 10%.

For the $\tauh\tauh$ channel, the uncertainty in the $\text{jet}\to\tauh$
background normalisation due to the extrapolation of the $F_{\text{F}}$ from
same-sign to opposite-sign regions ranges between 4 and 7%.

For the decay of the lepton to or a single-pion, an uncertainty in the
correction of $S_{\text{IP}}$ is applied by varying the size of the correction
by $\pm 25\%$, while for the decay to an electron the correction is varied by
40%. The uncertainty is converted into a normalisation uncertainty per decay
mode and ranges 1–5%. For the $\mathrm{a_{1}^{3pr}}\mathrm{a_{1}^{3pr}}$ mode,
the uncertainty in the SV reconstruction efficiency is 2%.

Finally, a 3% uncertainty in the efficiency of the candidates to pass the DNN
discrimination against muon and electron misidentifications is applied.

###  0.11.2 Shape uncertainties

The uncertainty in the reconstruction and identification efficiency is
typically of the order of 3%, and split into several uncertainties in each and
MVA decay mode bin. The uncertainties in these corrections originates from
uncertainties in the fits to the scale factors for these corrections and are
statistically dominated. The uncertainty in the trigger depends on the and
decay mode, and originates from the statistical uncertainty in parameterising
the turn-on curve of the triggers. The energy scale uncertainty is 0.8–1.1
(0.2–0.5)% for simulated (embedded) events, and is decay mode dependent. The
uncertainty in the momentum scale varies as a function of the $\eta$ of the
muon and ranges 0.4–2.7%. The uncertainty in the electron energy scale is less
than 1% and depends on the and $\eta$. The $\Pe\to\tauh$ energy scale
uncertainty ranges 0.5–6.5%, while the $\Pgm\to\tauh$ energy scale uncertainty
is 1%.

Uncertainties in the jet energy scale originate from different sources with
limited correlations. The uncertainties depend on the jet kinematics and are
typically larger in the forward regions. Uncertainties in the jet energy
resolution are also incorporated; these uncertainties are typically smaller
than the jet energy scale uncertainties. Uncertainties related to the hadronic
recoil response and resolution as derived from the $\PZ+\text{jets}$,
$\PW+\text{jets}$ and signal samples, are propagated to and observables
dependent on in the simulated samples that are subject to hadronic recoil
corrections. For the samples in which no hadronic recoil is applied (diboson,
single top quark, and ), the jet energy scale and resolution uncertainties as
well as the uncertainty in the unclustered energy are propagated to and
observables dependent on in the simulated samples instead.

The embedded samples contain small fractions of and diboson events. A shape
uncertainty is therefore applied by adding and subtracting 10% of the
simulated and diboson contributions. The top quark and Drell–Yan and mass
spectra are reweighed. For the top samples, the size of the correction is
taken as the uncertainty, while for the Drell–Yan samples the correction is
varied by 10%.

The $F_{\text{F}}$ values are parameterised with continuous functions, and the
statistical uncertainties in the fitted parameters are treated as nuisance
parameters. The uncertainties are parameterised in a manner that allows for
asymmetric variations above and below the value where the uncertainty is
minimal; the procedure is similar to the method described in detail in Ref.
[11]. The size of the correction in is taken as an uncertainty for all
$F_{\text{F}}$ values. For the $\tauh\tauh$ channel, the shape uncertainty in
the QCD same-sign to opposite-sign region correction is determined as the
difference between a correction binned in the distance between the two leptons
and the jet multiplicity, and the unbinned correction. For the
$\tau_{\ell}\tauh$ channels, the equivalent shape uncertainty is taken as the
size of the same-sign to opposite-sign correction. In addition, a systematic
uncertainty due to the light-lepton correction is taken as the size of the
correction. For the $\PW+\text{jets}$ $F_{\text{F}}$ values, the uncertainty
due to the extrapolation from the high-to the low-region is taken as the size
of residual differences observed when applying $F_{\text{F}}$ values derived
for high-simulated events to low-simulated events. For the $F_{\text{F}}$, a
systematic uncertainty is applied to account for potential differences between
data and simulation. To this purpose, the difference between $F_{\text{F}}$
values derived via data and simulated $\PW+\text{jets}$ samples is applied as
the uncertainty.

For uncertainties that are common to simulated and embedded samples we treat
the lepton and identification uncertainties and the lepton and energy scale
uncertainties as being 50% correlated. All other common uncertainties are
treated as being uncorrelated.

During the 2016 and 2017 data-taking periods, a gradual shift in the timing of
the inputs of the ECAL L1 trigger in the forward endcap region
($2.5<\abs{\eta}<3.0$) led to a specific inefficiency. Additional correction
factors and corresponding uncertainties are applied to the simulation to
account for this inefficiency. The magnitude of the uncertainties ranges
between 0–4% depending on the process, category, and channel.

For the signal samples, renormalisation and factorisation scales and parton
showering uncertainties were incorporated [98] .

The limited number of events in the signal and background templates is
accounted for using the ”Barlow–Beeston” method [106, 107], which assigns a
single nuisance parameter per bin per process. For background templates that
have been flattened as described in Section 0.9 the bin-by-bin uncertainties
are fully correlated such that there is only one independent nuisance
parameter for all $\phi_{{\textit{CP}}}$ bins. For background templates that
are symmetric in $\phi_{{\textit{CP}}}=180^{\circ}$ one nuisance parameter per
pair of symmetrised bins is utilised. It should be noted that for flattened
background templates multiple nuisance parameters are still needed per process
since multiple windows of increasing MVA score are used.

##  0.12 Results

In order to extract the CP-mixing angle $\alpha^{\PH\tau\tau}$, a simultaneous
fit to the data is performed using the likelihood function
$L(\lumi,\vec{\mu},\alpha^{\PH\tau\tau},\vec{\theta})$ that depends on
$\vec{\mu}=(\mu_{\Pg\Pg\PH},\mu_{\Pq\Pq\PH})$, which are the Higgs boson
signal strength modifiers (defined as the cross section times $\PH\to\PGt\PGt$
branching fraction with respect to the SM value), the CP-mixing angle
$\alpha^{\PH\tau\tau}$, and the nuisance parameters $\vec{\theta}$ that
account for the systematic uncertainties. In the fit, all $\PH\to\PGt\PGt$
production processes involving boson couplings, namely VBF and $\PV\PH$, are
scaled by $\mu_{\Pq\Pq\PH}$, while the $\Pg\Pg\PH$ process is scaled by
$\mu_{\Pg\Pg\PH}$. The fit is able to differentiate between these production
modes because the shapes of the MVA score distributions are different; the VBF
signal tends to peak more sharply towards larger MVA scores, whereas the
$\Pg\Pg\PH$ distribution is broader.

The likelihood function is defined as a product of conditional probabilities
$P$ over binned distributions of the discriminating observables in each event
category:

$L(\lumi,\vec{\mu},\alpha^{\PH\tau\tau},\vec{\theta})=\prod^{N_{\text{categories}}}_{\mathrm{j}}\prod^{N_{\text{bin}}}_{\mathrm{i}}P(n_{\text{i,j}}\mid\lumi\,\vec{\mu}\,\vec{A}_{\text{i,j}}(\vec{\theta},\alpha^{\PH\tau\tau})+B_{\text{i,j}}(\vec{\theta}))\,\prod^{N_{\text{nuisance}}}_{\mathrm{m}}C_{\mathrm{m}}(\vec{\theta}).$
(15)

In this equation, the Poisson distributions $P$ correspond to the observation
of $n_{\text{i,j}}$ events in each bin of the discriminating observable given
the expectation for the background $B_{\text{i,j}}(\vec{\theta})$ and the
signal
$S_{\text{i,j}}(\lumi,\alpha^{\PH\tau\tau},\vec{\mu},\vec{\theta})=\lumi\,\vec{\mu}\,\vec{A}_{\text{i,j}}(\vec{\theta},\alpha^{\PH\tau\tau})$,
in which is the integrated luminosity and
$\vec{A}_{\text{i,j}}(\vec{\theta},\alpha^{\PH\tau\tau})$ is the signal
acceptance in each production bin. Constraints on the nuisance parameters
corresponding to the systematic uncertainties described in Section 0.11 are
represented by the functions $C_{\mathrm{m}}(\vec{\theta})$. A more detailed
discussion on the formulation of the statistical inference may be found in
Refs. [108, 107]. The systematic uncertainties affecting the normalisation of
the signal and background templates are incorporated in the fit via nuisance
parameters with a log-normal prior probability density function. The shape-
altering systematic uncertainties are represented by nuisance parameters whose
variations cause continuous morphing of the signal or background template
shapes, and are assigned a Gaussian prior probability density function. The
bin-by-bin statistical uncertainties in the background samples are also
assigned a Gaussian prior probability density function.

Using the negative log-likelihood, which is defined as

$-2\Delta\ln
L=-2\,\left(\ln(L\alpha^{\PH\tau\tau})-\ln(L\alpha^{\PH\tau\tau}_{\text{best
fit}})\right),$ (16)

we find the 68.3, 95.5, and 99.7% confidence intervals when $-2\Delta\ln L$
equals 1.00, 4.02, and 8.81, respectively. A detailed discussion may be found
in Section 3.2 of Ref. [109].

The inputs to the likelihood fits differ for the signal and background
categories. For the signal categories, the $\phi_{{\textit{CP}}}$
distributions in bins of the MVA score are used (a subset of these are
displayed in Figs. 8–10). For the background categories, the MVA score
distributions are used. This allows for the background contributions and
systematic uncertainties to be further constrained, and helps to improve the
fit convergence.

###  0.12.1 $\alpha^{\PH\tau\tau}$ mixing angle results

We present the observed and expected negative log-likelihood scan for the
combination of the $\tau_{\Pe}\tauh$, $\tau_{\Pgm}\tauh$, and $\tauh\tauh$
channels in Fig. 11. The two rate parameters that scale the $\Pg\Pg\PH$ and
$\Pq\Pq\PH$ production signal strength were left to float freely in the fit.

Figure 11:  Negative log-likelihood scan for the combination of the
$\tau_{\Pe}\tauh$, $\tau_{\Pgm}\tauh$, and $\tauh\tauh$ channels. The observed
(expected) sensitivity to distinguish between the scalar and pseudoscalar
hypotheses, defined at $\alpha^{\PH\tau\tau}=0$ and $\pm 90^{\circ}$,
respectively, is $3.0\sigma$ ($2.6\sigma$). The observed (expected) value for
$\alpha^{\PH\tau\tau}$ is $-1\pm 19^{\circ}$ ($0\pm 21^{\circ}$) at the 68.3%
. At 95.5% the range is $\pm 41^{\circ}$ ($\pm 49^{\circ}$), and at the 99.7%
the observed range is $\pm 84^{\circ}$.

The data disfavour the pure CP-odd scenario at $3.0\sigma$. The expected
exclusion assuming the SM is $2.6\sigma$. The observed (expected) value of
$\alpha^{\PH\tau\tau}$ is found to be $-1\pm 19^{\circ}$ ($0\pm 21^{\circ}$)
at the 68.3% , and $\pm 41^{\circ}$ ($\pm 49^{\circ}$) at the 95.5% .
Furthermore, we obtain an observed $\pm 84^{\circ}$ at the 99.7% . The
uncertainty can be decomposed into: statistical; bin-by-bin fluctuations in
the background templates; experimental systematic uncertainties; and
theoretical uncertainties. In this decomposition we obtain

$\alpha^{\PH\tau\tau}=-1\pm 19\stat\pm 1\syst\pm 2\,\text{(bin-by-bin)}\pm
1\thy^{\circ}.$

This result is compatible with the SM predictions within the experimental
uncertainties.

The expected sensitivities of the $\tau_{\Pe}\tauh$, $\tau_{\Pgm}\tauh$, and
$\tauh\tauh$ channels are 1.0, 1.5, and $1.8\sigma$, respectively. The
$\Pgm\PGr$ mode yields the most sensitive expected contribution of
$1.2\sigma$, followed by the $\PGr\PGr$ and $\Pgp\PGr$ modes that contribute
1.1 and $1.0\sigma$, respectively. All other modes have sensitivities below
$1\sigma$.

The statistical uncertainties in the background templates are the subleading
source of systematic uncertainty in this analysis. As the dominant
contributions to the backgrounds are determined themselves from control
samples in data, the amount of data is the limiting factor in this
uncertainty. The next most dominant sources of uncertainty are the hadronic
trigger efficiency, theory uncertainties, the energy scale, and uncertainties
related to the implementation of the $F_{\text{F}}$ method.

It was shown in Ref. [36] that in the next-to-minimal supersymmetric model
mixing angles as large as ${\approx}27^{\circ}$ can be accommodated by the
latest electric dipole moment and Higgs boson measurements. This measurement
is thus sensitive to the larger allowed mixing angles in this model.

A fit to the data is also performed assuming
$\mu_{\Pg\Pg\PH}=\mu_{\Pq\Pq\PH}=\mu$. In this case $\mu$ is the combined
signal strength modifier that scales the total production cross section times
$\PH\to\PGt\PGt$ branching fraction relative to the SM value. In Fig. 12 we
display a scan of $\mu$ versus $\alpha^{\PH\tau\tau}$. We observe that there
is no strong correlation between these parameters.

Figure 12:  The 2-D scan of the signal strength modifier $\mu$ versus
$\alpha^{\PH\tau\tau}$. The 68.3, 95.5, and 99.7% confidence regions are
overlaid.

In order to make a 2-D scan of $\kappa_{\tau}$ and
$\widetilde{\kappa}_{\tau}$, as defined in Eq. (2), we parameterise the
likelihood from Eq. (15) in terms of $\kappa_{\tau}$ and
$\widetilde{\kappa}_{\tau}$. All other couplings are fixed to their expected
SM values.

In the case of a 2-D negative log-likelihood, the 68.3, 95.5, and 99.7%
confidence regions are found when $-2\Delta\ln L_{\mathrm{2D}}$ equals 2.30,
6.20, and 11.62 [109], respectively, defined analogously to Eq. (16) with the
likelihood now a function of both $\kappa_{\tau}$ and
$\widetilde{\kappa}_{\tau}$. All other CP-even (CP-odd) couplings affecting
the production cross sections and/or the $\PH\to\PGt\PGt$ branching fraction
are fixed to their SM values, $\kappa_{i}=1$ ($\widetilde{\kappa}_{i}=0$). The
observed result of the scan is shown in Fig. 13. It should be noted that the
fit is only sensitive to the relative sign between $\kappa_{\tau}$ and
$\widetilde{\kappa}_{\tau}$ and thus the scan has two best fit points for
positive and negative values of $\kappa_{\tau}$.

Figure 13:  The 2-D scan of the (reduced) CP-even ($\kappa_{\tau}$) and CP-odd
($\widetilde{\kappa}_{\tau}$) Yukawa couplings. The 68.3, 95.5, and 99.7%
confidence regions are overlaid.

In Fig. 14 we display the data of the $\PGr\PGr$, $\Pgp\PGr$, $\Pgm\PGr$, and
$\Pe\PGr$ channel together with CP-even and CP-odd predictions. These channels
are chosen as the same number of $\phi_{{\textit{CP}}}$ bins are used in the
fit to data, and collectively they account for most of the sensitivity to
$\alpha^{\PH\tau\tau}$. Events are included from all MVA score bins in these
signal categories. Each MVA score bin is weighed by $A\>S/(S+B)$, where $S$
and $B$ are the signal and background rates, respectively, and $A$ is a
measure for the average asymmetry between the scalar and pseudoscalar
distributions. The definition of the value of $A$ per bin is
$\abs{{\textit{CP}}^{\text{even}}-{\textit{CP}}^{\text{odd}}}/({\textit{CP}}^{\text{even}}+{\textit{CP}}^{\text{odd}})$,
and $A$ is normalised to the total number of bins. In this equation,
${\textit{CP}}^{\text{even}}$ and ${\textit{CP}}^{\text{odd}}$ are the scalar
and pseudoscalar contributions per bin. This distribution shows that the data
favour the CP-even scenario.

Figure 14:  The $\phi_{{\textit{CP}}}$ distributions for the $\PGr\PGr$,
$\Pgp\PGr$, $\Pgm\PGr$, and $\Pe\PGr$ channels weighed by $A\>S/(S+B)$ are
combined. Events are included from all MVA score bins in the four signal
categories. The background is subtracted from the data. The scalar
distribution is depicted in red, while the pseudoscalar is displayed in blue.
In the predictions, the rate parameters are taken from their best fit values.
The grey uncertainty band indicates the uncertainty in the subtracted
background component. In combining the channels, a phase-shift of
$180^{\circ}$ was applied to the channels involving a lepton since this
channel has a phase difference of $180^{\circ}$ with respect to the two
hadronic channels due to a sign-flip in the spectral function of the light
lepton.

##  0.13 Summary

The first measurement of the effective mixing angle $\alpha^{\PH\tau\tau}$
between scalar and pseudoscalar $\PH\PGt\PGt$ couplings has been presented for
a data set of proton-proton collisions at $\sqrt{s}=13\TeV$ corresponding to
an integrated luminosity of 137. The data were collected with the CMS
experiment at the LHC in the period 2016–2018. The following lepton decay
modes were included: , , , $\PGr^{\pm}\to\Pgppm\Pgpz$,
$\mathrm{a_{1}}^{\pm}\to\Pgppm\Pgpz\Pgpz$, and
$\mathrm{a_{1}}^{\pm}\to\Pgppm\PGpmp\Pgppm$. Dedicated strategies were adopted
to reconstruct the angle $\phi_{{\textit{CP}}}$ between the decay planes for
the various decay modes. The data disfavour the pure CP-odd scenario at 3.0
standard deviations. The observed effective mixing angle is found to be $-1\pm
19^{\circ}$, while the expected value is $0\pm 21^{\circ}$ at the 68.3%
confidence level (). The observed and expected uncertainties are found to be
$\pm 41^{\circ}$ and $\pm 49^{\circ}$ at the 95.5% , respectively, and the
observed sensitivity at the 99.7% is $\pm 84^{\circ}$. The leading uncertainty
in the measurement is statistical, implying that the precision of the
measurement will increase with the accumulation of more collision data. The
measurement is consistent with the standard model expectation, and reduces the
allowed parameter space for its extensions. Tabulated results are provided in
HEPDATA [110].

###### Acknowledgements.

We congratulate our colleagues in the CERN accelerator departments for the
excellent performance of the LHC and thank the technical and administrative
staffs at CERN and at other CMS institutes for their contributions to the
success of the CMS effort. In addition, we gratefully acknowledge the
computing centres and personnel of the Worldwide LHC Computing Grid and other
centres for delivering so effectively the computing infrastructure essential
to our analyses. Finally, we acknowledge the enduring support for the
construction and operation of the LHC, the CMS detector, and the supporting
computing infrastructure provided by the following funding agencies: BMBWF and
FWF (Austria); FNRS and FWO (Belgium); CNPq, CAPES, FAPERJ, FAPERGS, and
FAPESP (Brazil); MES and BNSF (Bulgaria); CERN; CAS, MoST, and NSFC (China);
MINCIENCIAS (Colombia); MSES and CSF (Croatia); RIF (Cyprus); SENESCYT
(Ecuador); MoER, ERC PUT and ERDF (Estonia); Academy of Finland, MEC, and HIP
(Finland); CEA and CNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany); GSRI
(Greece); NKFIA (Hungary); DAE and DST (India); IPM (Iran); SFI (Ireland);
INFN (Italy); MSIP and NRF (Republic of Korea); MES (Latvia); LAS (Lithuania);
MOE and UM (Malaysia); BUAP, CINVESTAV, CONACYT, LNS, SEP, and UASLP-FAI
(Mexico); MOS (Montenegro); MBIE (New Zealand); PAEC (Pakistan); MSHE and NSC
(Poland); FCT (Portugal); JINR (Dubna); MON, RosAtom, RAS, RFBR, and NRC KI
(Russia); MESTD (Serbia); SEIDI, CPAN, PCTI, and FEDER (Spain); MOSTR (Sri
Lanka); Swiss Funding Agencies (Switzerland); MST (Taipei); ThEPCenter, IPST,
STAR, and NSTDA (Thailand); TUBITAK and TAEK (Turkey); NASU (Ukraine); STFC
(United Kingdom); DOE and NSF (USA). Individuals have received support from
the Marie-Curie programme and the European Research Council and Horizon 2020
Grant, contract Nos. 675440, 724704, 752730, 758316, 765710, 824093, 884104,
and COST Action CA16108 (European Union); the Leventis Foundation; the Alfred
P. Sloan Foundation; the Alexander von Humboldt Foundation; the Belgian
Federal Science Policy Office; the Fonds pour la Formation à la Recherche dans
l’Industrie et dans l’Agriculture (FRIA-Belgium); the Agentschap voor
Innovatie door Wetenschap en Technologie (IWT-Belgium); the F.R.S.-FNRS and
FWO (Belgium) under the “Excellence of Science – EOS” – be.h project n.
30820817; the Beijing Municipal Science & Technology Commission, No.
Z191100007219010; the Ministry of Education, Youth and Sports (MEYS) of the
Czech Republic; the Deutsche Forschungsgemeinschaft (DFG), under Germany’s
Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306, and under
project number 400140256 - GRK2497; the Lendület (“Momentum”) Programme and
the János Bolyai Research Scholarship of the Hungarian Academy of Sciences,
the New National Excellence Program ÚNKP, the NKFIA research grants 123842,
123959, 124845, 124850, 125105, 128713, 128786, and 129058 (Hungary); the
Council of Science and Industrial Research, India; the Latvian Council of
Science; the Ministry of Science and Higher Education and the National Science
Center, contracts Opus 2014/15/B/ST2/03998 and 2015/19/B/ST2/02861 (Poland);
the Fundação para a Ciência e a Tecnologia, grant CEECIND/01334/2018
(Portugal); the National Priorities Research Program by Qatar National
Research Fund; the Ministry of Science and Higher Education, projects no.
14.W03.31.0026 and no. FSWW-2020-0008, and the Russian Foundation for Basic
Research, project No.19-42-703014 (Russia); the Programa Estatal de Fomento de
la Investigación Científica y Técnica de Excelencia María de Maeztu, grant
MDM-2015-0509 and the Programa Severo Ochoa del Principado de Asturias; the
Stavros Niarchos Foundation (Greece); the Rachadapisek Sompot Fund for
Postdoctoral Fellowship, Chulalongkorn University and the Chulalongkorn
Academic into Its 2nd Century Project Advancement Project (Thailand); the
Kavli Foundation; the Nvidia Corporation; the SuperMicro Corporation; the
Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).

## References

  * [1] F. Englert and R. Brout, “Broken symmetry and the mass of gauge vector mesons”,  Phys. Rev. Lett. 13 (1964) 321, 10.1103/PhysRevLett.13.321.
  * [2] P. W. Higgs, “Broken symmetries, massless particles and gauge fields”,  Phys. Lett. 12 (1964) 132, 10.1016/0031-9163(64)91136-9.
  * [3] P. W. Higgs, “Broken symmetries and the masses of gauge bosons”,  Phys. Rev. Lett. 13 (1964) 508, 10.1103/PhysRevLett.13.508.
  * [4] G. S. Guralnik, C. R. Hagen, and T. W. B. Kibble, “Global conservation laws and massless particles”,  Phys. Rev. Lett. 13 (1964) 585, 10.1103/PhysRevLett.13.585.
  * [5] P. W. Higgs, “Spontaneous symmetry breakdown without massless bosons”,  Phys. Rev. 145 (1966) 1156, 10.1103/PhysRev.145.1156.
  * [6] T. W. B. Kibble, “Symmetry breaking in non-abelian gauge theories”,  Phys. Rev. 155 (1967) 1554, 10.1103/PhysRev.155.1554.
  * [7] ATLAS Collaboration, “Observation of a new particle in the search for the standard model Higgs boson with the ATLAS detector at the LHC”,  Phys. Lett. B 716 (2012) 1, 10.1016/j.physletb.2012.08.020, arXiv:1207.7214.
  * [8] CMS Collaboration, “Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC”,  Phys. Lett. B 716 (2012) 30, 10.1016/j.physletb.2012.08.021, arXiv:1207.7235.
  * [9] CMS Collaboration, “Observation of a new boson with mass near 125 GeV in $\Pp\Pp$ collisions at $\sqrt{s}$ = 7 and 8 TeV”,  JHEP 06 (2013) 081, 10.1007/JHEP06(2013)081, arXiv:1303.4571.
  * [10] ATLAS and CMS Collaborations, “Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV”,  JHEP 08 (2016) 045, 10.1007/JHEP08(2016)045, arXiv:1606.02266.
  * [11] CMS Collaboration, “Observation of the Higgs boson decay to a pair of $\tau$ leptons with the CMS detector”,  Phys. Lett. B 779 (2018) 283, 10.1016/j.physletb.2018.02.004, arXiv:1708.00373.
  * [12] ATLAS Collaboration, “Cross-section measurements of the Higgs boson decaying into a pair of $\tau$-leptons in proton-proton collisions at a centre-of-mass energy of 13 TeV with the ATLAS detector”,  Phys. Rev. D 99 (2019) 072001, 10.1103/physrevd.99.072001, arXiv:1811.08856.
  * [13] CMS Collaboration, “A measurement of the Higgs boson mass in the diphoton decay channel”,  Phys. Lett. B 805 (2020) 135425, 10.1016/j.physletb.2020.135425, arXiv:2002.06398.
  * [14] CMS Collaboration, “On the mass and spin-parity of the Higgs boson candidate via its decays to Z boson pairs”,  Phys. Rev. Lett. 110 (2013) 081803, 10.1103/PhysRevLett.110.081803, arXiv:1212.6639.
  * [15] CMS Collaboration, “Measurement of the properties of a Higgs boson in the four-lepton final state”,  Phys. Rev. D 89 (2014) 092007, 10.1103/PhysRevD.89.092007, arXiv:1312.5353.
  * [16] CMS Collaboration, “Constraints on the spin-parity and anomalous HVV couplings of the Higgs boson in proton collisions at 7 and 8 TeV”,  Phys. Rev. D 92 (2015) 012004, 10.1103/PhysRevD.92.012004, arXiv:1411.3441.
  * [17] CMS Collaboration, “Limits on the Higgs boson lifetime and width from its decay to four charged leptons”,  Phys. Rev. D 92 (2015) 072010, 10.1103/PhysRevD.92.072010, arXiv:1507.06656.
  * [18] CMS Collaboration, “Combined search for anomalous pseudoscalar $\PH\PV\PV$ couplings in $\PV\PH(\PH\to\bbbar$) production and $\PH\to\PV\PV$ decay”,  Phys. Lett. B 759 (2016) 672, 10.1016/j.physletb.2016.06.004, arXiv:1602.04305.
  * [19] CMS Collaboration, “Constraints on anomalous Higgs boson couplings using production and decay information in the four-lepton final state”,  Phys. Lett. B 775 (2017) 1, 10.1016/j.physletb.2017.10.021, arXiv:1707.00541.
  * [20] CMS Collaboration, “Constraints on anomalous HVV couplings from the production of Higgs bosons decaying to $\tau$ lepton pairs”,  Phys. Rev. D 100 (2019) 112002, 10.1103/physrevd.100.112002, arXiv:1903.06973.
  * [21] CMS Collaboration, “Constraints on anomalous Higgs boson couplings to vector bosons and fermions in its production and decay using the four-lepton final state”,  Phys. Rev. D 104 (2021) 052004, 10.1103/PhysRevD.104.052004, arXiv:2104.12152.
  * [22] ATLAS Collaboration, “Evidence for the spin-0 nature of the Higgs boson using ATLAS data”,  Phys. Lett. B 726 (2013) 120, 10.1016/j.physletb.2013.08.026, arXiv:1307.1432.
  * [23] ATLAS Collaboration, “Study of the spin and parity of the Higgs boson in diboson decays with the ATLAS detector”,  Eur. Phys. J. C 75 (2015) 476, 10.1140/epjc/s10052-015-3685-1, arXiv:1506.05669. [Erratum: 10.1140/epjc/s10052-016-3934-y].
  * [24] ATLAS Collaboration, “Test of $CP$ invariance in vector-boson fusion production of the Higgs boson using the optimal observable method in the ditau decay channel with the ATLAS detector”,  Eur. Phys. J. C 76 (2016) 658, 10.1140/epjc/s10052-016-4499-5, arXiv:1602.04516.
  * [25] ATLAS Collaboration, “Measurement of inclusive and differential cross sections in the $\PH\to\PZ\PZ^{*}\to 4\ell$ decay channel in $\Pp\Pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”,  JHEP 10 (2017) 132, 10.1007/JHEP10(2017)132, arXiv:1708.02810.
  * [26] ATLAS Collaboration, “Measurement of the Higgs boson coupling properties in the $\PH\to\PZ\PZ^{*}\to 4\ell$ decay channel at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  JHEP 03 (2018) 095, 10.1007/JHEP03(2018)095, arXiv:1712.02304.
  * [27] ATLAS Collaboration, “Measurements of Higgs boson properties in the diphoton decay channel with 36 of $\Pp\Pp$ collision data at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  Phys. Rev. D 98 (2018) 052005, 10.1103/PhysRevD.98.052005, arXiv:1802.04146.
  * [28] C. Zhang and S. Willenbrock, “Effective-field-theory approach to top-quark production and decay”,  Phys. Rev. D 83 (2011) 034006, 10.1103/physrevd.83.034006, arXiv:1008.3869.
  * [29] R. Harnik et al., “Measuring $CP$ violation in $\PH\to\tau^{+}\tau^{-}$ at colliders”,  Phys. Rev. D 88 (2013) 076009, 10.1103/PhysRevD.88.076009, arXiv:1308.1094.
  * [30] T. Ghosh, R. Godbole, and X. Tata, “Determining the spacetime structure of bottom-quark couplings to spin-zero particles”,  Phys. Rev. D 100 (2019) 015026, 10.1103/physrevd.100.015026, arXiv:1904.09895.
  * [31] A. V. Gritsan, R. Röntsch, M. Schulze, and M. Xiao, “Constraining anomalous Higgs boson couplings to the heavy-flavor fermions using matrix element techniques”,  Phys. Rev. D 94 (2016) 055023, 10.1103/physrevd.94.055023, arXiv:1606.03107.
  * [32] CMS Collaboration, “Measurements of $\ttbar\PH$ production and the $CP$ structure of the Yukawa interaction between the Higgs boson and top quark in the diphoton decay channel”,  Phys. Rev. Lett. 125 (2020) 061801, 10.1103/PhysRevLett.125.061801, arXiv:2003.10866.
  * [33] ATLAS Collaboration, “$CP$ properties of Higgs boson interactions with top quarks in the $\ttbar\PH$ and $\PQt\PH$ processes using $\PH\to\PGg\PGg$ with the ATLAS detector”,  Phys. Rev. Lett. 125 (2020) 061802, 10.1103/PhysRevLett.125.061802, arXiv:2004.04545.
  * [34] ATLAS Collaboration, “Constraints on Higgs boson properties using $\PW\PW^{*}(\rightarrow\Pe\nu\mu\nu)jj$ production in 36.1 of $\sqrt{s}$=13 TeV pp collisions with the ATLAS detector”, 2021. arXiv:2109.13808. Submitted to Eur. Phys. J. C.
  * [35] D. Fontes, J. C. Rom$\tilde{\text{a}}$o, R. Santos, and J. P. Silva, “Large pseudoscalar Yukawa couplings in the complex 2HDM”,  JHEP 06 (2015) 060, 10.1007/JHEP06(2015)060, arXiv:1502.01720.
  * [36] S. F. King, M. Mühlleitner, R. Nevzorov, and K. Walz, “Exploring the $CP$-violating NMSSM: EDM constraints and phenomenology”,  Nucl. Phys. B 901 (2015) 526, 10.1016/j.nuclphysb.2015.11.003, arXiv:1508.03255.
  * [37] S. Berge, W. Bernreuther, and S. Kirchner, “Determination of the Higgs $CP$-mixing angle in the tau decay channels at the LHC including the Drell–Yan background”,  Eur. Phys. J. C 74 (2014) 3164, 10.1140/epjc/s10052-014-3164-0, arXiv:1408.0798.
  * [38] Particle Data Group, P. A. Zyla et al., “Review of particle physics”,  Prog. Theor. Exp. Phys. 2020 (2020) 083C01, 10.1093/ptep/ptaa104.
  * [39] S. Berge, W. Bernreuther, and S. Kirchner, “Determination of the Higgs $CP$-mixing angle in the tau decay channels”,  Nucl. Part. Phys. Proc. 273-275 (2016) 841, 10.1016/j.nuclphysbps.2015.09.129, arXiv:1410.6362.
  * [40] S. Berge, W. Bernreuther, B. Niepelt, and H. Spiesberger, “How to pin down the $CP$ quantum numbers of a Higgs boson in its tau decays at the LHC”,  Phys. Rev. D 84 (2011) 116003, 10.1103/PhysRevD.84.116003, arXiv:1108.0670.
  * [41] J. R. Dell’Aquila and C. A. Nelson, “Distinguishing a spin-0 technipion and an elementary Higgs boson: ${V}_{1}$${V}_{2}$ modes with decays into $\textit{I}^{-}_{A}\boldsymbol{l}_{B}$ and/or $q^{-}_{A}\boldsymbol{q}_{B}$”,  Phys. Rev. D 33 (1986) 93, 10.1103/PhysRevD.33.93.
  * [42] M. Krämer, J. Kühn, M. L. Stong, and P. M. Zerwas, “Prospects of measuring the parity of Higgs particles”,  Z. Phys. C 64 (1994) 21, 10.1007/bf01557231, arXiv:hep-ph/9404280.
  * [43] CMS Collaboration, “Performance of the CMS Level-1 trigger in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JINST 15 (2020) P10017, 10.1088/1748-0221/15/10/P10017, arXiv:2006.10165.
  * [44] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [45] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [46] P. Nason, “A new method for combining NLO QCD with shower Monte Carlo algorithms”,  JHEP 11 (2004) 040, 10.1088/1126-6708/2004/11/040, arXiv:hep-ph/0409146.
  * [47] S. Frixione, P. Nason, and C. Oleari, “Matching NLO QCD computations with parton shower simulations: the POWHEG method”,  JHEP 11 (2007) 070, 10.1088/1126-6708/2007/11/070, arXiv:0709.2092.
  * [48] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX”,  JHEP 06 (2010) 043, 10.1007/JHEP06(2010)043, arXiv:1002.2581.
  * [49] E. Bagnaschi, G. Degrassi, P. Slavich, and A. Vicini, “Higgs production via gluon fusion in the POWHEG approach in the SM and in the MSSM”,  JHEP 02 (2012) 088, 10.1007/JHEP02(2012)088, arXiv:1111.2854.
  * [50] P. Nason and C. Oleari, “NLO Higgs boson production via vector-boson fusion matched with shower in POWHEG”,  JHEP 02 (2010) 037, 10.1007/JHEP02(2010)037, arXiv:0911.5299.
  * [51] T. Ježo and P. Nason, “On the treatment of resonances in next-to-leading order calculations matched to a parton shower”,  JHEP 12 (2015) 065, 10.1007/JHEP12(2015)065, arXiv:1509.09071.
  * [52] F. Granata, J. M. Lindert, C. Oleari, and S. Pozzorini, “NLO QCD+EW predictions for HV and HV+jet production including parton-shower effects”,  JHEP 09 (2017) 012, 10.1007/JHEP09(2017)012, arXiv:1706.03522.
  * [53] K. Hamilton, P. Nason, E. Re, and G. Zanderighi, “NNLOPS simulation of Higgs boson production”,  JHEP 10 (2013) 222, 10.1007/JHEP10(2013)222, arXiv:1309.0017.
  * [54] K. Hamilton, P. Nason, and G. Zanderighi, “Finite quark-mass effects in the NNLOPS POWHEG+MiNLO Higgs generator”,  JHEP 05 (2015) 140, 10.1007/JHEP05(2015)140, arXiv:1501.04637.
  * [55] T. Sjöstrand et al., “An introduction to PYTHIA 8.2”,  Comput. Phys. Commun. 191 (2015) 159, 10.1016/j.cpc.2015.01.024, arXiv:1410.3012.
  * [56] T. Przedzinski, E. Richter-Was, and Z. Was, “Documentation of TauSpinner algorithms: program for simulating spin effects in $\tau$-lepton production at LHC”,  Eur. Phys. J. C 79 (2019) 91, 10.1140/epjc/s10052-018-6527-0, arXiv:1802.05459.
  * [57] NNPDF Collaboration, “Parton distributions for the LHC Run II”,  JHEP 04 (2015) 040, 10.1007/JHEP04(2015)040, arXiv:1410.8849.
  * [58] NNPDF Collaboration, “Parton distributions from high-precision collider data”,  Eur. Phys. J. C 77 (2017) 663, 10.1140/epjc/s10052-017-5199-5, arXiv:1706.00428.
  * [59] J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”,  JHEP 07 (2014) 079, 10.1007/JHEP07(2014)079, arXiv:1405.0301.
  * [60] J. Alwall et al., “Comparative study of various algorithms for the merging of parton showers and matrix elements in hadronic collisions”,  Eur. Phys. J. C 53 (2008) 473, 10.1140/epjc/s10052-007-0490-5, arXiv:0706.2569.
  * [61] S. Alioli, S.-O. Moch, and P. Uwer, “Hadronic top-quark pair-production with one jet and parton showering”,  JHEP 01 (2012) 137, 10.1007/JHEP01(2012)137, arXiv:1110.5251.
  * [62] E. Re, “Single-top $Wt$-channel production matched with parton showers using the POWHEG method”,  Eur. Phys. J. C 71 (2011) 1547, 10.1140/epjc/s10052-011-1547-z, arXiv:1009.2450.
  * [63] R. Frederix, E. Re, and P. Torrielli, “Single-top $t$-channel hadroproduction in the four-flavour scheme with POWHEG and aMC@NLO”,  JHEP 09 (2012) 130, 10.1007/JHEP09(2012)130, arXiv:1207.5391.
  * [64] CMS Collaboration, “Event generator tunes obtained from underlying event and multiparton scattering measurements”,  Eur. Phys. J. C 76 (2016) 155, 10.1140/epjc/s10052-016-3988-x, arXiv:1512.00815.
  * [65] CMS Collaboration, “Extraction and validation of a new set of CMS PYTHIA8 tunes from underlying-event measurements”,  Eur. Phys. J. C 80 (2020) 4, 10.1140/epjc/s10052-019-7499-4, arXiv:1903.12179.
  * [66] GEANT4 Collaboration, “—a simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [67] CMS Collaboration, “Particle-flow reconstruction and global event description with the CMS detector”,  JINST 12 (2017) P10003, 10.1088/1748-0221/12/10/P10003, arXiv:1706.04965.
  * [68] CMS Collaboration, “Description and performance of track and primary-vertex reconstruction with the CMS tracker”,  JINST 9 (2014) P10009, 10.1088/1748-0221/9/10/P10009, arXiv:1405.6569.
  * [69] K. Rose, “Deterministic annealing for clustering, compression, classification, regression, and related optimization problems”,  Proc. IEEE 86 (1998) 2210, 10.1109/5.726788.
  * [70] W. Waltenberger, R. Frühwirth, and P. Vanlaer, “Adaptive vertex fitting”,  J. Phys. G: Nuc. Part. Phys. 34 (2007) N343, 10.1088/0954-3899/34/12/n01.
  * [71] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-jet clustering algorithm”,  JHEP 04 (2008) 063, 10.1088/1126-6708/2008/04/063, arXiv:0802.1189.
  * [72] CMS Collaboration, “Performance of CMS muon reconstruction in pp collision events at $\sqrt{s}=7$ TeV”,  JINST 7 (2012) P10002, 10.1088/1748-0221/7/10/P10002, arXiv:1206.4071.
  * [73] CMS Collaboration, “Performance of electron reconstruction and selection with the CMS detector in proton-proton collisions at $\sqrt{s}=8$ TeV”,  JINST 10 (2015) P06005, 10.1088/1748-0221/10/06/P06005, arXiv:1502.02701.
  * [74] M. Cacciari, G. P. Salam, and G. Soyez, “FastJet user manual”,  Eur. Phys. J. C 72 (2012) 1896, 10.1140/epjc/s10052-012-1896-2, arXiv:1111.6097.
  * [75] CMS Collaboration, “Jet algorithms performance in 13 TeV data”, CMS Physics Analysis Summary CMS-PAS-JME-16-003, 2017.
  * [76] CMS Collaboration, “Jet energy scale and resolution in the CMS experiment in pp collisions at 8 TeV”,  JINST 12 (2017) P02014, 10.1088/1748-0221/12/02/P02014, arXiv:1607.03663.
  * [77] CMS Collaboration, “Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV”,  JINST 13 (2018) P05011, 10.1088/1748-0221/13/05/p05011, arXiv:1712.07158.
  * [78] D. Bertolini, P. Harris, M. Low, and N. Tran, “Pileup per particle identification”,  JHEP 10 (2014) 059, 10.1007/JHEP10(2014)059, arXiv:1407.6013.
  * [79] CMS Collaboration, “Performance of missing transverse momentum reconstruction in proton-proton collisions at $\sqrt{s}=13\TeV$ using the CMS detector”,  JINST 14 (2019) P07004, 10.1088/1748-0221/14/07/P07004, arXiv:1903.06078.
  * [80] CMS Collaboration, “Performance of reconstruction and identification of $\tau$ leptons decaying to hadrons and $\nu_{\tau}$ in pp collisions at $\sqrt{s}=$ 13 TeV”,  JINST 13 (2018) P10005, 10.1088/1748-0221/13/10/P10005, arXiv:1809.02816.
  * [81] CMS Collaboration, “Performance of the DeepTau algorithm for the discrimination of taus against jets, electron, and muons”, CMS Detector Performance Note CMS-DP-2019-033, 2019.
  * [82] L. Bianchini, J. Conway, E. K. Friis, and C. Veelken, “Reconstruction of the Higgs mass in $\PH\to\PGt\PGt$ events by dynamical likelihood techniques”,  J. Phys. Conf. Ser. 513 (2014) 022035, 10.1088/1742-6596/513/2/022035.
  * [83] S. Berge and W. Bernreuther, “Determining the $CP$ parity of Higgs bosons at the LHC in the to 1-prong decay channels”,  Phys. Lett. B 671 (2009) 470, 10.1016/j.physletb.2008.12.065, arXiv:0812.1910.
  * [84] S. Berge, W. Bernreuther, and S. Kirchner, “Prospects of constraining the Higgs boson’s $CP$ nature in the tau decay channel at the LHC”,  Phys. Rev. D 92 (2015) 096012, 10.1103/PhysRevD.92.096012, arXiv:1510.03850.
  * [85] G. R. Bower, T. Pierzchala, Z. Wa̧s, and M. Worek, “Measuring the Higgs boson’s parity using $\PGt\to\PGr\PGn$”,  Phys. Lett. B 543 (2002) 227, 10.1016/S0370-2693(02)02445-0, arXiv:hep-ph/0204292.
  * [86] V. Cherepanov, E. Richter-Was, and Z. Was, “Monte Carlo, fitting and machine learning for tau leptons”,  SciPost Phys. Proc. 1 (2019) 018, 10.21468/SciPostPhysProc.1.018, arXiv:1811.03969.
  * [87] V. Cherepanov and A. Zotz, “Kinematic reconstruction of $\PZ/\PH\to\tau\tau$ decay in proton-proton collisions”, 2018. arXiv:1805.06988.
  * [88] S. Jadach, J. H. Kuhn, and Z. Wa̧s, “TAUOLA: A library of Monte Carlo programs to simulate decays of polarized leptons”,  Comput. Phys. Commun. 64 (1990) 275, 10.1016/0010-4655(91)90038-M.
  * [89] M. Jezabek, Z. Wa̧s, S. Jadach, and J. H. Kuhn, “The decay library TAUOLA, update with exact O(alpha) QED corrections in $\Pgt\to\PGm(\Pe)\Pgn\Pagn$ decay modes”,  Comput. Phys. Commun. 70 (1992) 69, 10.1016/0010-4655(92)90092-D.
  * [90] S. Jadach, Z. Wa̧s, R. Decker, and J. H. Kuhn, “The decay library TAUOLA: Version 2.4”,  Comput. Phys. Commun. 76 (1993) 361, 10.1016/0010-4655(93)90061-G.
  * [91] CLEO Collaboration, “Hadronic structure in the decay $\Pgt\to\PGnGt\PGpm\PGpz\PGpz$ and the sign of the helicity”,  Phys. Rev. D 61 (2000) 012002, 10.1103/PhysRevD.61.012002, arXiv:hep-ex/9902022.
  * [92] CMS Collaboration, “Identification of hadronic tau decay channels using multivariate analysis (MVA decay mode)”, CMS Detector Performance Note CMS-DP-2020-041, 2020.
  * [93] T. Chen and C. Guestrin, “XGBoost: A scalable tree boosting system”, in  Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, p. 785.  2016\.  arXiv:1603.02754.  10.1145/2939672.2939785.
  * [94] CMS Collaboration, “An embedding technique to determine $\tau\tau$ backgrounds in proton-proton collision data”,  JINST 14 (2019) P06032, 10.1088/1748-0221/14/06/p06032, arXiv:1903.01216.
  * [95] CMS Collaboration, “Measurement of the $\PZ/\PGg^{*}\to\tau\tau$ cross section in pp collisions at $\sqrt{s}=$ 13 TeV and validation of $\tau$ lepton analysis techniques”,  Eur. Phys. J. C 78 (2018) 708, 10.1140/epjc/s10052-018-6146-9, arXiv:1801.03535.
  * [96] CMS Collaboration, “Measurements of inclusive W and Z cross sections in pp collisions at $\sqrt{s}=7$ TeV”,  JHEP 01 (2011) 080, 10.1007/JHEP01(2011)080, arXiv:1012.2466.
  * [97] CMS Collaboration, “Measurement of the differential cross section for top quark pair production in pp collisions at $\sqrt{s}=8\TeV$”,  Eur. Phys. J. C 75 (2015) 2339, 10.1140/epjc/s10052-015-3709-x, arXiv:1505.04480.
  * [98] LHC Higgs Cross Section Working Group, “Handbook of LHC Higgs cross sections: 4. deciphering the nature of the Higgs sector”,  CERN (2016) 10.23731/CYRM-2017-002, arXiv:1610.07922.
  * [99] CMS Collaboration, “CMS luminosity measurements for the 2016 data-taking period”, CMS Physics Analysis Summary CMS-PAS-LUM-17-001, 2017.
  * [100] CMS Collaboration, “CMS luminosity measurement for the 2017 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-17-004, 2018.
  * [101] CMS Collaboration, “CMS luminosity measurement for the 2018 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-18-002, 2019.
  * [102] Y. Li and F. Petriello, “Combining QCD and electroweak corrections to dilepton production in the framework of the FEWZ simulation code”,  Phys. Rev. D 86 (2012) 094034, 10.1103/physrevd.86.094034, arXiv:1208.5967.
  * [103] M. Czakon and A. Mitov, “Top++: A program for the calculation of the top-pair cross-section at hadron colliders”,  Comput. Phys. Commun. 185 (2014) 2930, 10.1016/j.cpc.2014.06.021, arXiv:1112.5675.
  * [104] CMS Collaboration, “Measurement of the WZ production cross section in pp collisions at $\sqrt{s}=13$ TeV”,  Phys. Lett. B 766 (2017) 268, 10.1016/j.physletb.2017.01.011, arXiv:1607.06943.
  * [105] CMS Collaboration, “Cross section measurement of $t$-channel single top quark production in pp collisions at $\sqrt{s}=$ 13 TeV”,  Phys. Lett. B 772 (2017) 752, 10.1016/j.physletb.2017.07.047, arXiv:1610.00678.
  * [106] R. Barlow and C. Beeston, “Fitting using finite Monte Carlo samples”,  Comput. Phys. Commun. 77 (1993) 219, 10.1016/0010-4655(93)90005-W.
  * [107] J. S. Conway, “Incorporating nuisance parameters in likelihoods for multisource spectra”, in  PHYSTAT 2011, p. 115.  2011\.  arXiv:1103.0354.  10.5170/CERN-2011-006.115.
  * [108] ATLAS and CMS Collaborations, and LHC Higgs Combination Group, “Procedure for the LHC Higgs boson search combination in Summer 2011”, Technical Report CMS-NOTE-2011-005, ATL-PHYS-PUB-2011-11, 2011.
  * [109] CMS Collaboration, “Precise determination of the mass of the Higgs boson and tests of compatibility of its couplings with the standard model predictions using proton collisions at 7 and 8 TeV”,  Eur. Phys. J. C 75 (2015) 212, 10.1140/epjc/s10052-015-3351-7, arXiv:1412.8662.
  * [110] “HEPData record for this analysis”, 2021.  10.17182/hepdata.104978.

##  .14 The CMS Collaboration

Yerevan Physics Institute, Yerevan, Armenia
A. Tumasyan Institut für Hochenergiephysik, Vienna, Austria
W. Adam, J.W. Andrejkovic, T. Bergauer, S. Chatterjee, M. Dragicevic, A.
Escalante Del Valle, R. Frühwirth1, M. Jeitler1, N. Krammer, L. Lechner, D.
Liko, I. Mikulec, P. Paulitsch, F.M. Pitters, J. Schieck1, R. Schöfbeck, D.
Schwarz, S. Templ, W. Waltenberger, C.-E. Wulz1 Institute for Nuclear
Problems, Minsk, Belarus
V. Chekhovsky, A. Litomin, V. Makarenko Universiteit Antwerpen, Antwerpen,
Belgium
M.R. Darwish2, E.A. De Wolf, T. Janssen, T. Kello3, A. Lelek, H. Rejeb Sfar,
P. Van Mechelen, S. Van Putte, N. Van Remortel Vrije Universiteit Brussel,
Brussel, Belgium
F. Blekman, E.S. Bols, J. D’Hondt, M. Delcourt, H. El Faham, S. Lowette, S.
Moortgat, A. Morton, D. Müller, A.R. Sahasransu, S. Tavernier, W. Van Doninck,
P. Van Mulders Université Libre de Bruxelles, Bruxelles, Belgium
D. Beghin, B. Bilin, B. Clerbaux, G. De Lentdecker, L. Favart, A. Grebenyuk,
A.K. Kalsi, K. Lee, M. Mahdavikhorrami, I. Makarenko, L. Moureaux, L. Pétré,
A. Popov, N. Postiau, E. Starling, L. Thomas, M. Vanden Bemden, C. Vander
Velde, P. Vanlaer, L. Wezenbeek Ghent University, Ghent, Belgium
T. Cornelis, D. Dobur, J. Knolle, L. Lambrecht, G. Mestdach, M. Niedziela, C.
Roskas, A. Samalan, K. Skovpen, M. Tytgat, B. Vermassen, M. Vit Université
Catholique de Louvain, Louvain-la-Neuve, Belgium
A. Benecke, A. Bethani, G. Bruno, F. Bury, C. Caputo, P. David, C. Delaere,
I.S. Donertas, A. Giammanco, K. Jaffel, Sa. Jain, V. Lemaitre, K. Mondal, J.
Prisciandaro, A. Taliercio, M. Teklishyn, T.T. Tran, P. Vischia, S. Wertz
Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, Brazil
G.A. Alves, C. Hensel, A. Moraes Universidade do Estado do Rio de Janeiro, Rio
de Janeiro, Brazil
W.L. Aldá Júnior, M. Alves Gallo Pereira, M. Barroso Ferreira Filho, H.
Brandao Malbouisson, W. Carvalho, J. Chinellato4, E.M. Da Costa, G.G. Da
Silveira5, D. De Jesus Damiao, S. Fonseca De Souza, D. Matos Figueiredo, C.
Mora Herrera, K. Mota Amarilo, L. Mundim, H. Nogima, P. Rebello Teles, A.
Santoro, S.M. Silva Do Amaral, A. Sznajder, M. Thiel, F. Torres Da Silva De
Araujo, A. Vilela Pereira Universidade Estadual Paulista a, Universidade
Federal do ABC b, São Paulo, Brazil
C.A. Bernardesa,a,5, L. Calligarisa, T.R. Fernandez Perez Tomeia, E.M.
Gregoresa,b, D.S. Lemosa, P.G. Mercadantea,b, S.F. Novaesa, Sandra S. Padulaa
Institute for Nuclear Research and Nuclear Energy, Bulgarian Academy of
Sciences, Sofia, Bulgaria
A. Aleksandrov, G. Antchev, R. Hadjiiska, P. Iaydjiev, M. Misheva, M. Rodozov,
M. Shopova, G. Sultanov University of Sofia, Sofia, Bulgaria
A. Dimitrov, T. Ivanov, L. Litov, B. Pavlov, P. Petkov, A. Petrov Beihang
University, Beijing, China
T. Cheng, T. Javaid6, M. Mittal, L. Yuan Department of Physics, Tsinghua
University
M. Ahmad, G. Bauer, C. Dozen7, Z. Hu, J. Martins8, Y. Wang, K. Yi9,10
Institute of High Energy Physics, Beijing, China
E. Chapon, G.M. Chen6, H.S. Chen6, M. Chen, F. Iemmi, A. Kapoor, D. Leggat, H.
Liao, Z.-A. Liu6, V. Milosevic, F. Monti, R. Sharma, J. Tao, J. Thomas-
wilsker, J. Wang, H. Zhang, J. Zhao State Key Laboratory of Nuclear Physics
and Technology, Peking University, Beijing, China
A. Agapitos, Y. An, Y. Ban, C. Chen, A. Levin, Q. Li, X. Lyu, Y. Mao, S.J.
Qian, D. Wang, Q. Wang, J. Xiao Sun Yat-Sen University, Guangzhou, China
M. Lu, Z. You Institute of Modern Physics and Key Laboratory of Nuclear
Physics and Ion-beam Application (MOE) - Fudan University, Shanghai, China
X. Gao3, H. Okawa Zhejiang University, Hangzhou, China
Z. Lin, M. Xiao Universidad de Los Andes, Bogota, Colombia
C. Avila, A. Cabrera, C. Florez, J. Fraga Universidad de Antioquia, Medellin,
Colombia
J. Mejia Guisao, F. Ramirez, J.D. Ruiz Alvarez, C.A. Salazar González
University of Split, Faculty of Electrical Engineering, Mechanical Engineering
and Naval Architecture, Split, Croatia
D. Giljanovic, N. Godinovic, D. Lelas, I. Puljak University of Split, Faculty
of Science, Split, Croatia
Z. Antunovic, M. Kovac, T. Sculac Institute Rudjer Boskovic, Zagreb, Croatia
V. Brigljevic, D. Ferencek, D. Majumder, M. Roguljic, A. Starodumov11, T. Susa
University of Cyprus, Nicosia, Cyprus
A. Attikis, K. Christoforou, E. Erodotou, A. Ioannou, G. Kole, M. Kolosova, S.
Konstantinou, J. Mousa, C. Nicolaou, F. Ptochos, P.A. Razis, H. Rykaczewski,
H. Saka Charles University, Prague, Czech Republic
M. Finger12, M. Finger Jr.12, A. Kveton Escuela Politecnica Nacional, Quito,
Ecuador
E. Ayala Universidad San Francisco de Quito, Quito, Ecuador
E. Carrera Jarrin Academy of Scientific Research and Technology of the Arab
Republic of Egypt, Egyptian Network of High Energy Physics, Cairo, Egypt
A. Ellithi Kamel13, E. Salama14,15 Center for High Energy Physics (CHEP-FU),
Fayoum University, El-Fayoum, Egypt
A. Lotfy, M.A. Mahmoud National Institute of Chemical Physics and Biophysics,
Tallinn, Estonia
S. Bhowmik, R.K. Dewanjee, K. Ehataht, M. Kadastik, S. Nandan, C. Nielsen, J.
Pata, M. Raidal, L. Tani, C. Veelken Department of Physics, University of
Helsinki, Helsinki, Finland
P. Eerola, L. Forthomme, H. Kirschenmann, K. Osterberg, M. Voutilainen
Helsinki Institute of Physics, Helsinki, Finland
S. Bharthuar, E. Brücken, F. Garcia, J. Havukainen, M.S. Kim, R. Kinnunen, T.
Lampén, K. Lassila-Perini, S. Lehti, T. Lindén, M. Lotti, L. Martikainen, M.
Myllymäki, J. Ott, H. Siikonen, E. Tuominen, J. Tuominiemi Lappeenranta
University of Technology, Lappeenranta, Finland
P. Luukka, H. Petrow, T. Tuuva IRFU, CEA, Université Paris-Saclay, Gif-sur-
Yvette, France
C. Amendola, M. Besancon, F. Couderc, M. Dejardin, D. Denegri, J.L. Faure, F.
Ferri, S. Ganjour, A. Givernaud, P. Gras, G. Hamel de Monchenault, P. Jarry,
B. Lenzi, E. Locci, J. Malcles, J. Rander, A. Rosowsky, M.Ö. Sahin, A. Savoy-
Navarro16, M. Titov, G.B. Yu Laboratoire Leprince-Ringuet, CNRS/IN2P3, Ecole
Polytechnique, Institut Polytechnique de Paris, Palaiseau, France
S. Ahuja, F. Beaudette, M. Bonanomi, A. Buchot Perraguin, P. Busson, A.
Cappati, C. Charlot, O. Davignon, B. Diab, G. Falmagne, S. Ghosh, R. Granier
de Cassagnac, A. Hakimi, I. Kucher, J. Motta, M. Nguyen, C. Ochando, P.
Paganini, J. Rembser, R. Salerno, U. Sarkar, J.B. Sauvan, Y. Sirois, A.
Tarabini, A. Zabi, A. Zghiche Université de Strasbourg, CNRS, IPHC UMR 7178,
Strasbourg, France
J.-L. Agram17, J. Andrea, D. Apparu, D. Bloch, G. Bourgatte, J.-M. Brom, E.C.
Chabert, C. Collard, D. Darej, J.-C. Fontaine17, U. Goerlach, C. Grimault,
A.-C. Le Bihan, E. Nibigira, P. Van Hove Institut de Physique des 2 Infinis de
Lyon (IP2I ), Villeurbanne, France
E. Asilar, S. Beauceron, C. Bernet, G. Boudoul, C. Camen, A. Carle, N. Chanon,
D. Contardo, P. Depasse, H. El Mamouni, J. Fay, S. Gascon, M. Gouzevitch, B.
Ille, I.B. Laktineh, H. Lattaud, A. Lesauvage, M. Lethuillier, L. Mirabito, S.
Perries, K. Shchablo, V. Sordini, L. Torterotot, G. Touquet, M. Vander Donckt,
S. Viret Georgian Technical University, Tbilisi, Georgia
I. Bagaturia18, I. Lomidze, Z. Tsamalaidze12 RWTH Aachen University, I.
Physikalisches Institut, Aachen, Germany
V. Botta, L. Feld, K. Klein, M. Lipinski, D. Meuser, A. Pauls, N. Röwert, J.
Schulz, M. Teroerde RWTH Aachen University, III. Physikalisches Institut A,
Aachen, Germany