Searches for additional Higgs bosons and for vector leptoquarks in $\tau\tau$ final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV

Three searches are presented for signatures of physics beyond the standard model (SM) in $\tau\tau$ final states in proton-proton collisions at the LHC, using a data sample collected with the CMS detector at $\sqrt{s}$ = 13 TeV, corresponding to an integrated luminosity of 138 fb$^{-1}$. Upper limits at 95% confidence level (CL) are set on the products of the branching fraction for the decay into $\tau$ leptons and the cross sections for the production of a new boson $\phi$, in addition to the H(125) boson, via gluon fusion (gg$\phi$) or in association with b quarks, ranging from $\mathcal{O}$(10 pb) for a mass of 60 GeV to 0.3 fb for a mass of 3.5 TeV each. The data reveal two excesses for gg$\phi$ production with local $p$-values equivalent to about three standard deviations at $m_\phi$ = 0.1 and 1.2 TeV. In a search for $t$-channel exchange of a vector leptoquark U$_1$, 95% CL upper limits are set on the dimensionless U$_1$ leptoquark coupling to quarks and $\tau$ leptons ranging from 1 for a mass of 1 TeV to 6 for a mass of 5 TeV, depending on the scenario. In the interpretations of the $M_\mathrm{h}^{125}$ and $M_\mathrm{h, EFT}^{125}$ minimal supersymmetric SM benchmark scenarios, additional Higgs bosons with masses below 350 GeV are excluded at 95% CL.


Introduction
The discovery of a Higgs boson with a mass of around 125 GeV,H(125), at the LHC in 2012 [1][2][3] has turned the standard model (SM) of particle physics into a theory that could be valid up to the Planck scale. In the SM, H(125) emerges from the spontaneous breaking of the electroweak SU(2) L symmetry. While the nature of the underlying mechanism leading to this symmetry breaking and the exact form of the required symmetry-breaking potential are still to be explored, the measured couplings of H (125) to fermions and gauge bosons, with 5-20% experimental precision [4][5][6][7], are in good agreement with the expectation for an SM Higgs boson with a mass of 125.38 ± 0.14 GeV [8]. The SM still leaves several fundamental questions related to particle physics unaddressed, including the presence of dark matter and the observed baryon asymmetry in nature. Many extensions of the SM that address these questions require a more complex structure of the part of the theory that is related to SU(2) L breaking, often referred to as the Higgs sector. Such models usually predict additional spin-0 states and modified properties of H (125) with respect to the SM expectation. Models incorporating supersymmetry (SUSY) [9,10] are prominent examples. In the minimal extension of the SM, the minimal supersymmetric SM (MSSM) [11,12], the model predicts three neutral and two charged Higgs bosons.
Searches for additional heavy neutral Higgs bosons in the context of the MSSM were carried out in electron-positron collisions at the LEP collider at CERN [13] and in proton-antiproton collisions at the Fermilab Tevatron [14][15][16][17]. At the LHC such searches have been carried out by the ATLAS and CMS Collaborations in the b quark [18][19][20][21], dimuon [22][23][24][25], and ττ [22, 26-33] final states. The ττ final state has a leading role in these searches, since τ leptons can be identified with higher purity than b quarks and backgrounds from genuine ττ events can be estimated with higher accuracy, while the branching fractions for the decay into τ leptons are typically larger than those for the decay into muons because of the larger τ lepton mass. There are several other examples of extended Higgs sectors, which are summarized in Ref. [34], that could give appreciable resonant ττ production rates in addition to the known SM processes at the LHC. Furthermore, models that include additional coloured states carrying both baryon and lepton quantum numbers, known as leptoquarks [35,36], can lead to an enhancement in the nonresonant production rates of ττ pairs with large invariant masses via the leptoquark t-channel exchange. Searches for resonant and nonresonant ττ signatures are thus complementary in the exploration of physics beyond the SM (BSM) at the LHC. Recent searches for singleand pair-production of third-generation leptoquarks at the LHC are reported in Refs. [37][38][39][40][41][42][43][44][45].
In this paper the results of three searches for both resonant and nonresonant ττ signatures are presented: i) The first search, which is meant to be as model independent as possible, targets the production of a single narrow spin-0 resonance φ, in addition to H (125), via gluon fusion (ggφ) or in association with b quarks (bbφ). Assumptions that have been made for this search are that the width of φ is small compared with the experimental resolution, and that the φ transverse momentum (p T ) spectrum for ggφ production as well as the relative contributions of t-and b-quarks to ggφ production are as expected for an SM Higgs boson at the tested mass value.
ii) The second search targets the t-channel exchange of a vector leptoquark U 1 .
iii) The third search exploits selected benchmark scenarios of the MSSM that rely on the signal from three neutral Higgs bosons, one of which is associated with H (125).
The results are based on the proton-proton (pp) collision data collected at the LHC during the years 2016-2018, at √ s = 13 TeV, by the CMS experiment. The data correspond to an integrated luminosity of 138 fb −1 . The analysis is performed in four ττ final states: eµ, eτ h , µτ h , and τ h τ h , where e, µ, and τ h indicate τ decays into electrons, muons, and hadrons, respectively. For this analysis the most significant backgrounds are estimated from data, which includes all SM processes with two genuine τ leptons in the final state, and processes where quark-or gluon-induced jets are misidentified as τ h , denoted as jet → τ h .
The paper is organized as follows. Section 2 gives an overview of the phenomenology of the BSM physics scenarios under consideration. Section 3 describes the CMS detector, and Section 4 describes the event reconstruction. Section 5 summarizes the event selection and categorization used for the extraction of the signal. The data model and systematic uncertainties are described in Sections 6 and 7. Section 8 contains the results of the analysis. Section 9 briefly summarizes the paper. A complete set of tabulated results of this search for all tested mass hypotheses is available in the HEPData database [46].

Signal models
Neutral (pseudo)scalar bosons φ appear in many extensions of the SM. They may have different couplings to the upper and lower components of the SU(2) L fermion fields (associated with up-and down-type fermions) and gauge bosons. In several models, like the MSSM models discussed in Section 2.2, the φ couplings to down-type fermions are enhanced with respect to the expectation for an SM Higgs boson of the same mass, while the couplings to up-type fermions and vector bosons are suppressed. This makes down-type fermion final states, such as ττ, particularly interesting for searches for neutral Higgs bosons in addition to H (125). An enhancement in the couplings to down-type fermions also increases the bbφ production cross section relative to ggφ, which is another characteristic signature of these models and motivates the search for enhanced production cross sections in this production mode with respect to the SM expectation. In a first interpretation of the data, which is meant to be as model independent as possible, we search for φ production via the ggφ and bbφ processes in a range of 60 ≤ m φ ≤ 3500 GeV, where m φ denotes the hypothesized φ mass. Diagrams for these processes are shown in Fig. 1. In a second, more specific interpretation of the data, we search for nonresonant ττ production in a model with vector leptoquarks. Finally, in a third interpretation of the data, we survey the parameter space of two indicative benchmark scenarios of the MSSM, which predict multiresonance signatures, one of which is associated with H (125). The most important characteristics of the vector leptoquark model and the MSSM are described in the following.  Figure 1: Diagrams for the production of neutral Higgs bosons φ (left) via gluon fusion, labelled as ggφ, and (middle and right) in association with b quarks, labelled as bbφ in the text. In the middle diagram, a pair of b quarks is produced from the fusion of two gluons, one from each proton. In the right diagram, a b quark from one proton scatters from a gluon from the other proton. In both cases φ is radiated off one of the b quarks.
The Lagrangian for the U 1 coupling to SM fermions is given by [72] with the coupling constant g U , where q L and d R (l L and e R ) denote the left-and right-handed quark (lepton) doublets, and β L and β R are left-and right-handed coupling matrices, which are assumed to have the structures: The motivations for the assumed structures of these matrices are given in Ref. [72]. The normalization of g U is chosen to give β bτ L = 1. Two benchmark scenarios are considered, with different assumptions made about the value of β bτ R . In the first benchmark scenario ("VLQ BM 1"), β bτ R is assumed to be zero. In the second benchmark scenario ("VLQ BM 2"), β bτ R is assumed to be −1, which corresponds to a Pati-Salam-like [49,68] U 1 leptoquark. The β sτ L couplings are set to their preferred values from global fits to the low-energy observables presented in Ref. [72], as summarized in Table 1. The β dτ L , β sµ L , and β bµ L couplings are small and have negligible influence on the ττ signature, and therefore have been set to zero.
If the U 1 leptoquark mass (m U ) is sufficiently small, the U 1 particle will contribute to the ττ spectrum via pair production with each U 1 subsequently decaying to a qτ pair. For larger m U , the pair production cross section is suppressed because of the decreasing probability that the initial-state partons possess sufficiently large momentum fractions of the corresponding protons to produce on-shell U 1 pairs. In this case the dominant contribution to the ττ spectrum is via U 1 t-channel exchange in the bb initial-state as illustrated in Fig. 2, with subdominant contributions from the equivalent bs, sb, and ss initiated processes. In our analysis we target the kinematic region of m U 1 TeV, motivated by the experimental exclusion limits on m U by direct searches, e.g. in Ref. [42]. The contribution to the ττ spectrum from U 1 pair production is negligible in this case, and we therefore consider only production through the t-channel exchange.

The MSSM
In the MSSM, which is a concrete example of the more general class of two Higgs doublet models (2HDMs) [81,82], the Higgs sector requires two SU(2) doublets, Φ u and Φ d , to provide masses for up-and down-type fermions. In CP-conserving 2HDMs, this leads to the prediction of two charged (H ± ) and three neutral φ bosons (h, H, and A), where h and H (with At tree level in the MSSM, the masses of these five Higgs bosons and α can be expressed in terms of the known gauge boson masses and two additional parameters, which can be chosen as m A and the ratio of the vacuum expectation values of the neutral components of Φ u and Φ d , Dependencies on additional parameters of the soft SUSY breaking mechanism enter via higherorder corrections in perturbation theory. In the exploration of the MSSM Higgs sector these additional parameters are usually set to fixed values in the form of indicative benchmark scenarios to illustrate certain properties of the theory. The most recent set of MSSM benchmark scenarios provided by the LHC Higgs Working Group has been introduced in Refs. [83][84][85] and summarized in Ref. [86]. The corresponding predictions of masses, cross sections, and branching fractions can be obtained from Ref. [87]. With one exception (the M 125 H scenario), in these scenarios h takes the role of H (125), and H and A are nearly degenerate in mass (m H ≈ m A ) in a large fraction of the provided parameter space.
For values of m A much larger than the mass of the Z boson, the coupling of H and A to downtype fermions is enhanced by tan β with respect to the expectation for an SM Higgs boson of the same mass, while the coupling to vector bosons and up-type fermions is suppressed. For increasing values of tan β, bbφ (with φ = H, A) is enhanced relative to ggφ production. The larger contribution of b quarks to the loop in Fig. 1 (left) in addition leads to softer spectra of the H and A transverse momentum. Extra SUSY particles influence the production and decay via higher-order contributions to the interaction vertices that belong to b quark lines. They also contribute directly to the loop in Fig. 1 (left).

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. 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, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors 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 100 kHz within a fixed latency of about 4 µs [88]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [89]. 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. [90].

Event reconstruction
The reconstruction of the pp collision products is based on the particle-flow (PF) algorithm [91], which combines the information from all CMS subdetectors to reconstruct a set of particle candidates (PF candidates), identified as charged and neutral hadrons, electrons, photons, and muons. In the 2016 (2017-2018) data sets the average number of interactions per bunch crossing was 23 (32). The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Ref. [92]. Secondary vertices, which are displaced from the PV, might be associated with decays of longlived particles emerging from the PV. Any other collision vertices in the event are associated with additional, mostly soft, inelastic pp collisions, referred to as pileup (PU).
Electrons are reconstructed using tracks from hits in the tracking system and energy deposits in the ECAL [93,94]. To increase their purity, reconstructed electrons are required to pass a multivariate electron identification discriminant, which combines information on track quality, shower shape, and kinematic quantities. For this analysis, a working point with an identification efficiency of 90% is used, for a rate of jets misidentified as electrons of ≈1%. Muons in the event are reconstructed by combining the information from the tracker and the muon detectors [95]. The presence of hits in the muon detectors already leads to a strong suppression of particles misidentified as muons. Additional identification requirements on the track fit quality and the compatibility of individual track segments with the fitted track can reduce the misidentification rate further. For this analysis, muon identification requirements with an efficiency of ≈99% are chosen, with a misidentification rate below 0.2% for pions.
The contributions from backgrounds to the electron and muon selections are further reduced by requiring the corresponding lepton to be isolated from any hadronic activity in the detector. This property is quantified by an isolation variable where p e(µ) T corresponds to the electron (muon) p T and ∑ p charged T , ∑ E neutral T , and ∑ E γ T to the p T (or transverse energy E T ) sum of all charged particles, neutral hadrons, and photons, in a predefined cone of radius ∆R = (∆η) 2 + (∆ϕ) 2 around the lepton direction at the PV, where ∆η and ∆ϕ (measured in radians) correspond to the angular distances of the particle to the lepton in the η and azimuthal ϕ directions. The chosen cone size is ∆R = 0.3 (0.4) for electrons (muons). The lepton itself is excluded from the calculation. To mitigate any distortions from PU, only those charged particles whose tracks are associated with the PV are included. Since an unambiguous association with the PV is not possible for neutral hadrons and photons, an estimate of the contribution from PU (p PU T ) is subtracted from the sum of ∑ E neutral T and ∑ E γ T . This estimate is obtained from tracks not associated with the PV in the case of I µ rel and the mean energy flow per area unit in the case of I e rel . For negative values, the I e(µ) rel is set to zero. For further characterization of the event, all reconstructed PF candidates are clustered into jets using the anti-k T algorithm with a distance parameter of 0.4, as implemented in the FASTJET software package [96,97]. To identify jets resulting from the hadronization of b quarks (b jets) the DEEPJET algorithm is used, as described in Refs. [98,99]. In this analysis a working point of this algorithm is chosen that corresponds to a b jet identification efficiency of ≈80% for a misidentification rate for jets originating from light-flavour quarks or gluons of O(1%) [100]. Jets with p T > 30 GeV and |η| < 4.7 and b jets with p T > 20 GeV and |η| < 2.4 are used in 2016. From 2017 onwards, after the upgrade of the silicon pixel detector, the b jet η range is extended to |η| < 2.5.
Jets are also used as seeds for the reconstruction of τ h candidates. This is done by exploiting the substructure of the jets using the "hadrons-plus-strips" algorithm, as described in Refs. [101,102]. Decays into one or three charged hadrons with up to two neutral pions with p T > 2.5 GeV are used. Neutral pions are reconstructed as strips with dynamic size in η-ϕ from reconstructed photons and electrons contained in the seeding jet, where the latter originate from photon conversions. The strip size varies as a function of the p T of the electron or photon candidates. The τ h decay mode is then obtained by combining the charged hadrons with the strips. To distinguish τ h candidates from jets originating from the hadronization of quarks or gluons, and from electrons or muons, the DEEPTAU (DT) algorithm is used, as described in Ref. [102]. This algorithm exploits the information of the reconstructed event record (comprising tracking, impact parameter, and calorimeter cluster information), the kinematic and object identification properties of the PF candidates in the vicinity of the τ h candidate and those of the τ h candidate itself, and quantities that estimate the PU density of the event. It results in a multiclassification output y DT α (α = τ, e, µ, jet) equivalent to a Bayesian probability of the τ h candidate to originate from a genuine τ lepton, the hadronization of a quark or gluon, an electron, or a muon. From this output three discriminants are built according to For the analysis presented here, predefined working points of D e , D µ , and D jet [102] are chosen depending on the ττ final state, for which the τ h selection efficiencies and misidentification rates are given in Table 2. Since the jet → τ h misidentification rate strongly depends on the p T and initiating parton type of the misidentified jet, it should be viewed as approximate.
The missing transverse momentum vector p miss T is also used for further categorization of the events. It is calculated as the negative vector p T sum of all PF candidates, weighted by their probability to originate from the PV [103], and exploits the pileup-per-particle identification algorithm [104] to reduce the PU dependence. With p miss T we refer to the magnitude of this quantity.

Selection of ττ candidates
Depending on the final state, the online selection in the HLT step is based either on the presence of a single electron, muon, or τ h candidate, or an eµ, eτ h , µτ h , or τ h τ h pair in an event. In the offline selection further requirements on p T , η, and I e(µ) rel are applied in addition to the object identification requirements described in Section 4, as summarized in Table 3.
In the eµ final state an electron and a muon with p T > 15 GeV and |η| < 2.4 are required. Depending on the trigger path that has led to the online selection of an event, a stricter requirement of p T > 24 GeV is imposed on one of the two leptons to ensure a sufficiently high trigger efficiency of the HLT selection. Both leptons are required to be isolated from any hadronic activity in the detector according to I e(µ) In the eτ h (µτ h ) final state, an electron (muon) with p T > 25 (20) GeV is required if the event was selected by a trigger based on the presence of the eτ h (µτ h ) pair in the event. From 2017 onwards, the threshold on the muon is raised to 21 GeV. If the event was selected by a singleelectron trigger, the p T requirement on the electron is increased to 26, 28, or 33 GeV for the years 2016, 2017, or 2018, respectively. For muons, the p T requirement is increased to 23 (25) GeV for 2016 (2017-2018), if selected by a single-muon trigger. The electron (muon) is required to be contained in the central part of the detector with |η| < 2.1, and to be isolated according to I e(µ) rel < 0.15. The τ h candidate is required to have |η| < 2.3 and p T > 35 (32) GeV if selected by an eτ h (µτ h ) pair trigger, or p T > 30 GeV if selected by a single-electron (single-muon) trigger. In the τ h τ h final state, both τ h candidates are required to have |η| < 2.1 and p T > 40 GeV. For events only selected by a single τ h trigger, the τ h candidate that has been identified with the triggering object is required to have p T > 120 (180) GeV for events recorded in 2016 (2017-2018).
The selected τ lepton decay candidates are required to be of opposite charge and to be separated by more than ∆R = 0.3 in the η-ϕ plane in the eµ final state and by more than 0.5 otherwise. The closest distance of their tracks to the PV is required to be d z < 0.2 cm along the beam axis. For electrons and muons, an additional requirement of d xy < 0.045 cm in the transverse plane is applied. In rare cases, where more than the required number of τ h candidates fulfilling all selection requirements is found, the candidate with the highest D jet score is chosen. For electrons and muons, the most isolated candidate is chosen.
To avoid the assignment of single events to more than one final state, events with additional electrons or muons, fulfilling looser selection requirements than those given for each corresponding ττ final state above, are rejected from the selection. These requirements also help with the suppression of background processes, such as Z/γ * → ee or Z/γ * → µµ. Table 3: Offline selection requirements applied to the electron, muon, and τ h candidates used for the selection of the τ pair. The expressions first and second lepton refer to the label of the final state in the first column. The p T requirements are given in GeV. For the eµ final state two lepton pair trigger paths, with a stronger requirement on the p T of the electron (muon), are used for the online selection of the event. For the eτ h , µτ h , and τ h τ h final states, the values (in parentheses) correspond to the lepton pair (single lepton) trigger paths that have been used in the online selection. A detailed discussion is given in the text.

Standard categories and signal extraction
To increase the sensitivity of the searches, all selected events are further split into categories. Events with at least one b jet, according to the selection requirements given in Section 4, are combined into a global "b tag" category, used to target bbφ production and to control the background from top quark pair (tt) production. All other events are subsumed into a global "no b tag" category. The events in the τ h τ h final state are not further categorized beyond that point. In the eτ h and µτ h final states, more categories are introduced in the global "b tag" and "no b tag" categories, based on the transverse mass of the e (µ)-p miss T system defined as where ∆ϕ refers to the azimuthal angular difference between p i T and p j T . Events are divided into a tight-m T (m e(µ) T < 40 GeV) and a loose-m T (40 < m e(µ) T < 70 GeV) category. The φ signal is expected to be concentrated in the tight-m T category. However, the loose-m T category increases the acceptance for m φ 700 GeV.
In the eµ final state, events are categorized based on the observable D ζ [105] defined as whereζ corresponds to the bisectional direction between p  isotropically distributed leading to nonpeaking distributions in D ζ . For ττ events from resonant decays, p miss T is expected to roughly coincide withζ, and a stronger correlation between p miss ζ and p vis ζ is expected to lead to a peaking distribution about D ζ ≈ 0 GeV. The inputs to the reconstruction of D ζ are illustrated in Fig. 3. Three further categories are introduced as high-D ζ (D ζ > 30 GeV), medium-D ζ (−10 < D ζ < 30 GeV), and low-D ζ (−35 < D ζ < −10 GeV). A φ signal is expected to be concentrated in the medium-D ζ category. However, the low-and high-D ζ categories still contribute to an increase of the sensitivity of the model-independent φ search in the eµ final state by ≈10%. A control category in the eµ final state with at least one b jet and D ζ < −35 GeV is used to constrain the normalization of tt events in the fit used for signal extraction.
In summary, this leads to 17 event categories per data-taking year. Figure 4 shows the D ζ and m µ T distributions in the eµ and µτ h final states, before splitting the events into the categories described above. The category definitions are indicated by the vertical dashed lines in the figures. An overview of the categories described above is given in Fig. 5.
In all cases the signal is extracted from distributions of the total transverse mass [26] defined as where τ 1(2) refers to the first (second) τ final state indicated in the eµ, eτ h , µτ h , and τ h τ h final state labels, and m T between two objects with transverse momenta p τ 1 T and p τ 2 T is defined in Eq. (6). This quantity is expected to provide superior discriminating power between resonant signals with m φ 250 GeV and nonpeaking backgrounds, such as W+ jets or tt production in the high-mass tails of the distribution.
This strategy is used for the model-independent φ search, to extract the expected signal for hypothesized values of m φ ≥ 250 GeV. It is also used for the extraction of the A and H signal (for m A , m H 250 GeV), when interpreting the data in MSSM benchmark scenarios, and for the vector leptoquark search, which is most sensitive to an excess over the background expectation for m tot T 250 GeV as will be discussed in Section 6.4.2.
To increase the sensitivity of the analyses for the model-independent φ search for hypothesized values of m φ < 250 GeV and the low-mass resonance h for the interpretation of the data in MSSM benchmark scenarios, this signal extraction strategy is modified as discussed in the following sections.

Modifications for the low-mass model-independent φ search
For hypothesized values of m φ < 250 GeV, the background from Z/γ * → ττ production, which features a peaking mass distribution in a region close to the signal mass, starts to exceed the nonpeaking backgrounds. The m tot T distribution loses discrimination power and some of the categories that were introduced to increase the acceptance for high-mass signals are not useful anymore. Therefore, the signal extraction strategy is modified in the following way. The low-D ζ and loose-m T categories are removed. The remaining "no b tag" categories are further split by p τ τ T , obtained from the vectorial sum p T of the visible τ decay products and p miss T , according to p τ τ T is used as an estimate for the φ p T (p φ T ) in data. No further splitting based on p τ τ T is applied to the "b tag" categories because of the lower event populations in these categories. In summary, this leads to 26 event categories per data-taking year. An overview of this modified set of categories is given in Fig. 6. Figure 6: Overview of the categories used for the extraction of the signal for the modelindependent φ search for 60 ≤ m φ < 250 GeV.
In these categories, the signal is extracted from a likelihood-based fit of the invariant mass of the ττ system, m τ τ , before the decay of the τ leptons [106]. This estimate combines the measurement of p miss T and its covariance matrix with the measurements of the visible ττ decay products, utilizing the matrix elements for unpolarized τ decays [107] for the decay into leptons and the two-body phase space [108] for the decay into hadrons. On average the resolution of m τ τ amounts to about 10-25%, depending on the kinematic properties of the ττ system and the ττ final states, where the latter is related to the number of neutrinos that escape detection. This approach exploits the better m φ resolution of m τ τ compared to m tot T , together with the usually harder p φ T , compared to the Z/γ * → ττ p T spectrum.

Modifications for the MSSM interpretation
The MSSM predicts three neutral Higgs bosons φ, one of which is identified with H(125). Each benchmark scenario has to match the observed H (125) properties. To exploit the best possible experimental knowledge about H(125) all events in the global "no b tag" category are split by m τ τ . For events with m τ τ > 250 GeV, the categories described in Section 5.2.1 are used. For events with m τ τ < 250 GeV, the neural-network-based (NN) analysis, which was used for the stage-0 simplified template cross section measurements of Ref. [109], is used to obtain the most precise estimates from data for H(125) production via gluon fusion (ggh), vector boson fusion (VBF), and vector boson associated production (Vh). Although the NN is trained specifically to target events with an SM-like φ with m φ = 125 GeV, signal events for the additional Higgs bosons can also enter the NN categories for m φ 250 GeV, and the y l discriminators contribute to the separation of such events from the background.
This modification adds 18 background and 8 signal categories from the NN-analysis per datataking year. We will refer to these as the "NN categories" throughout this paper. In these categories, the H(125) signal is extracted from distributions of the NN output functions y l in each signal and background category l. For the NN-analysis in Ref. [109], m eµ T calculated from p e T + p µ T and p miss T is required to be less than 60 GeV in the eµ final state, to prevent event overlap with analyses of other H(125) decay modes in the SM interpretation. For the analysis presented here, this requirement is replaced by D ζ > −35 GeV.
A summary of the categories and discriminating variables used for signal extraction for each of the analyses presented in this paper is given in Table 4. Table 4: Event categories and discriminants used for the extraction of the signals, for the searches described in this paper. We note that m φ refers to the hypothesized mass of the modelindependent φ search, while m τ τ refers to the reconstructed mass of the ττ system before the decays of the τ leptons, and thus to an estimate of m φ in data. The variable y l refers to the output functions of the NNs used for signal extraction in Ref. [109].  Table 5: Background processes contributing to the event selection, as discussed in Section 5. The symbol corresponds to an electron or muon. The second column refers to the experimental signature in the analysis, the last four columns indicate the estimation methods used to model each corresponding signature, as described in Sections 6.1-6.4. Diboson and single t production are part of the process group iv) discussed in Section 6. QCD(eµ) refers to QCD multijet production with an eµ pair in the final state.
Estimation method Background process Final-state signature τ-emb. F F SS Sim.

Background and signal modelling
All SM background sources that are relevant after the event selection described in Section 5 are listed in Table 5. The expected background composition depends on the ττ final state, event category, and the tested signal mass hypothesis. The most abundant source of background in the "b tag" categories is tt production. In the "no b tag" categories Z/γ * → ττ forms the largest fraction of background processes, followed by W+ jets production and events containing purely quantum chromodynamics (QCD) induced gluon and light-flavour quark jets, referred to as QCD multijet production. These backgrounds are grouped according to their experimental signatures into: i) events containing genuine τ lepton pairs (ττ); ii) events with quark-or gluon-induced jets misidentified as τ h candidates (jet → τ h ) or light leptons (jet → ) in the eµ final state; iii) tt events where an intermediate W boson decays into an electron, muon, or τ lepton, which do not fall into the previous groups (labelled as "tt" in later figures); iv) remaining background processes that are of minor importance for the analysis and not yet included in any of the previous event groups (labelled as "others" in later figures).
Event group (i) mostly contains Z/γ * → ττ events, with smaller contributions from tt, diboson, and single t quark production. These events are modelled using the τ-embedding method [110], labelled "τ-emb." in Table 5, and discussed in Section 6.1. Event group (ii) contains events from QCD multijet, W+ jets, Z+ jets, tt, diboson and single t quark production with jet → τ h misidentification, and QCD multijet production with jet → misidentification in the eµ final state. The events with jet → τ h misidentification are estimated from the "fake factor" (F F ) method, labelled "F F " in Table 5, and discussed in Section 6.2. The events with jet → misidentification in the eµ final state are estimated from the "same-sign" (SS) method, labelled "SS" in Table 5, and discussed in Section 6.3. Event group (iv) comprises diboson and single t quark production (labelled as "electroweak" in Fig. 4 left), H(125) production, Z/γ * → µµ and Z/γ * → ee events, and W+ jets events with jet → misidentification. Events from event groups iii) and iv), and all signal processes are estimated from full event simulation, labelled "Sim." in Table 5, and discussed in Section 6.4.

Backgrounds with genuine τ lepton pairs (ττ)
For all events where the decay of a Z boson results in two genuine τ leptons, the τ-embedding method, as described in Ref. [110], is used. For this purpose, µµ events are selected in data. All energy deposits of the muons are removed from the event record and replaced by simulated τ lepton decays with the same kinematic properties as the selected muons. In this way the method relies only on the simulation of the well-understood τ lepton decay and its energy deposits in the detector, while all other parts of the event, such as the identification and reconstruction of jets, including b jets, or the non-τ related parts of p miss T are obtained from data. This results in an improved modelling of the data compared with the simulation of the full process. In turn, several simulation-to-data corrections, as detailed in Section 6.5, are not needed. The selected muons predominantly originate from Z/γ * → µµ events. However, contributions from other processes resulting in two genuine τ leptons, like tt or diboson production, are also covered by this model. For a selection with no (at least one) b jet in the event, as described in Section 5, 97% (84%) of the µµ events selected for the τ-embedding method are expected to originate from Z/γ * → µµ and <1% (14%) from tt production. A detailed discussion of the selection of the original µµ events, the exact procedure itself, its range of validity, and related uncertainties is reported in Ref. [110].

Backgrounds with jets misidentified as hadronically decaying τ leptons (jet → τ h )
The main processes contributing to jet → τ h events in the eτ h , µτ h , and τ h τ h final states are QCD multijet, W+ jets, and tt production. These events are estimated using the F F method described in Refs. [32,111], and adapted to the analyses described in this paper.
For this purpose, the signal region (SR), defined by the event selection given in Section 5, is complemented by three additional regions: the application region (AR) and two determination regions DR i , where i stands for QCD or W+ jets. For the AR a looser working point for the identification of the τ h candidate is chosen and the events from the SR are excluded, which is the only selection difference with respect to the SR. In this way the AR forms an orthogonal, though still adjacent, sideband to the SR that is enriched in jet → τ h events. The events in the AR are then multiplied with a transfer function, which is obtained from each corresponding DR i or simulation, to estimate the contribution of jet → τ h events in the SR. The background processes in the AR and each corresponding DR i that are not targeted by this method are estimated either from simulation or the τ-embedding method and subtracted from the data.
In the τ h τ h final state, where QCD multijet production contributes 95% of the events in the AR, the transfer function is determined from DR QCD only, for which the charges of the two selected τ h candidates are required to be of same sign. This function is assumed to be applicable also for the small fraction of W+ jets and tt events in the AR. In this final state, both τ h candidates usually originate from jet → τ h misidentification. We require only the τ h candidate with the larger p T to fulfil the AR requirements, which provides an estimate for events where only this τ h candidate is misidentified. Events in which the τ h candidate with the larger p T is a genuine τ lepton and the one with the lower p T is misidentified, which constitute ≈2% of the total jet → τ h background, are modelled from simulation.
In the eτ h (µτ h ) final state, where the sharing of processes contributing to the AR is more equal, separate contributions to the transfer function F i F are used, where the index i runs over the processes of QCD multijet, W+ jets, and tt production. For QCD multijet and W+ jets production each F i F is derived in its corresponding DR i . For DR QCD we require 0.05 < I e(µ) rel < 0.15 and the charges of the selected e(µ) and the τ h candidate to be of same sign. For DR W+ jets we require m e(µ) T > 70 GeV and the absence of b jets. The estimate of F tt F is obtained from simulation. Each F i F is then used to estimate the yield N SR and kinematic properties of the combination of the main contributing backgrounds i in the SR from the number of events N AR in the AR according to Each F i F is combined into a weighted sum, using the simulation-based estimate of the fractions w i of each process in the AR. A template fit to the data in the AR yields a similar result for the w i .
Each F i F is computed on an event-by-event basis. It mainly depends on the p T of the τ h candidate with the larger p T , p τ h T , the ratio p jet T /p τ h T where p jet T corresponds to the p T of the jet seeding the τ h reconstruction, and the jet multiplicity N Jets . Each F i F is further subject to a number of residual corrections derived from both control regions in data and simulation to take subleading dependencies of the F i F into account. Depending on the transfer function F i F and the ττ final state these are dependencies on p T , the invariant mass of the visible decay products of the ττ system, I rel , or p τ h T of the second-leading τ h candidate.

Backgrounds with jets misidentified as electron-muon pairs (QCD(eµ))
The background from QCD multijet production where two quark-or gluon-induced jets are misidentified as an eµ pair is estimated using the SS method. In this case, an AR is distinguished from the SR by requiring the charges of the electron and muon to have the same sign. A sideband region DR is defined requiring the muon to be nonisolated (0.2 < I µ rel < 0.5). From this DR an SS to opposite-sign (OS) transfer function F T is obtained to extrapolate the number N AR of events in the AR to the number N SR of events in the SR according to The function F T primarily depends on the distance ∆R(e, µ) between the e and µ trajectories in η-ϕ and N Jets . Additional dependencies on the electron and muon p T enter via a bias correction, ranging from 0.85-0.9. To validate the method, a second transfer function F T is calculated from a modified DR with an isolated muon (I µ rel < 0.2) and nonisolated electron (0.15 < I e rel < 0.5), which is applied to the SS selection of the DR. The resulting event yield and shapes of the m tot T and m τ τ distributions are compared to the OS selection of the DR. This test reveals a consistent result within the statistical uncertainties of the estimate, for events with N b-jets = 0. For events with N b-jets ≥ 1, a global correction factor r b is required, with a value of 0.71-0.75 depending on the year of data-taking.
A potential bias from requiring the muon to be nonisolated in the definition of DR is checked from a third definition of the transfer function F T , in a DR with a nonisolated muon (0.2 < I µ rel < 0.5) and electron (0.15 < I e rel < 0.5). This test reveals another correction of 0.94-0.95, depending on the year of data-taking, to correct for the fact that r b , with an isolated muon, is systematically smaller by ≈5% than in the case of a nonisolated muon.

Simulated backgrounds and signal
In the τ h τ h final state, the τ-embedding and F F methods cover 97% of all expected background events. The fractions of expected background events described by these two methods are 83% in the eτ h and 90% in the µτ h final states. In the eµ final state, 53% of all events are obtained by either the τ-embedding or SS method. All remaining events originate from processes such as Z boson, tt, or diboson production, where at least one decay of a vector boson into an electron or muon is not covered by any of the previously discussed methods. These backgrounds and the signal processes are modelled using the simulation of the full processes.

Background processes
The W+ jets and Z/γ * → processes are simulated at leading order (LO) accuracy in the strong coupling α S , using the MADGRAPH5 aMC@NLO 2.2.2 (2.4.2) event generator [112,113] for the simulation of the data taken in 2016 (2017-2018). To increase the number of simulated events in regions of high signal purity, supplementary samples are generated with up to four outgoing partons in the hard interaction. For diboson production, MADGRAPH5 aMC@NLO is used at next-to-LO (NLO) precision in α S . In each case, the FxFx [114] (MLM [115]) prescription is used to match the NLO (LO) matrix element calculation with the parton shower model. For tt [116] and (t-channel) single t quark production [117], samples are generated at NLO precision in α S using POWHEG 2.0 [118][119][120][121]. The POWHEG version 1.0 at NLO precision is used for single t quark production in association with a W boson (tW channel) [122].
When compared with data, W+ jets, Z/γ * → , tt, and single t quark events in the tW channel are normalized to their cross sections at next-to-NLO (NNLO) precision in α S [123][124][125]. Single t quark (t-channel) and diboson events are normalized to their cross sections at NLO precision in α S or higher [125][126][127].

Signal processes
The kinematic properties of single h production are simulated at NLO precision in α S using POWHEG 2.0 separately for the production via ggh [128], VBF [129], or in association with a Z (Zh) or W (Wh) boson [130,131]. For ggh production, the distributions of the h boson p T and the jet multiplicity in the simulation are tuned to match the NNLO accuracy obtained from full phase space calculations with the NNLOPS event generator [132,133]. For this purpose, h is assumed to behave as expected from the SM. This applies to the modelling of H(125) as part of the background for the model-independent φ search, as well as for the SM and the MSSM hypotheses for the interpretation of the data in MSSM benchmark scenarios, where h is associated with H(125) with properties as expected from the SM.
The production of φ, H, and A bosons via gluon fusion is simulated at NLO precision in α S using the 2HDM implementation of POWHEG 2.0 [128]. To account for the multiscale nature of the process in the NLO plus parton shower prediction, the p T spectra corresponding to the contributions from the t quark only, b quark only, and tb-interference are each calculated separately. The POWHEG damping factor h damp , which controls the matching between the matrix element calculation and the parton shower, is set specifically for each contribution as proposed in Refs. [134][135][136].
For the model-independent φ search, the individual distributions are combined according to their contribution to the total cross section as expected for an SM-like Higgs boson with given mass. For the tests of MSSM benchmark scenarios, where the contributions of the individual distributions also depend on the model parameters, these distributions are scaled using the effective Yukawa couplings as predicted by the corresponding benchmark model [87], before combining them into one single prediction. In this context, the tan β-enhanced SUSY corrections to the φbb couplings are also included via the corresponding effective Yukawa couplings, where appropriate. Other SUSY contributions have been checked to amount to less than a few percent and are neglected. An example of the A boson p T spectrum for m A = 1.6 TeV and tan β = 30 is shown in Fig. 7 (left). The bbφ production is simulated at NLO precision in α S using the corresponding POWHEG 2.0 implementation [137] in the four-flavour scheme (4FS).
The signal process of the U 1 t-channel exchange is simulated in the five-flavour scheme (5FS) at LO precision in α S using the MADGRAPH5 aMC@NLO event generator, v2.6.5 [138]. Events are generated with up to one additional outgoing parton from the matrix element calculation and matched following the MLM prescription, with the matching scale Q match set to 40 GeV. The contribution from on-shell U 1 → qτ production and decay is excluded during the event generation. Samples are produced with g U = 1, for several values of m U between 1 and 5 TeV. We observe no large dependence, neither of the templates used for signal extraction nor of the overall cross section, on the assumed U 1 decay width Γ, even after variations of factors of 0.5 and 2 and therefore, for each considered value of m U , we choose Γ to approximately match the value predicted for U 1 production with couplings as obtained from the global fit presented in Ref. [72].
We expect a sizeable effect of destructive interference between the U 1 signal and Z/γ * → ττ production, where the relative sizes of the interference and noninterference contributions depend on g U . To include this dependence we generate separate samples for each contribution to form signal templates, which are negative in case of the interference contribution. These are scaled by g 4 U (for the noninterference contribution) and g 2 U (for the interference contribution), respectively, and combined to form the overall signal distributions for any value of g U . Finally, the resulting signal event yields are normalized to the cross sections for the inclusive U 1 mediated pp → ττ process, computed at LO precision in α S . The contribution of the U 1 t-channel exchange to the m tot T distribution in the τ h τ h final state for m U = 1 TeV and g U = 1.5 for the VLQ BM 1 scenario is shown in Fig. 7 (right). As visible from the figure, a complex contribution of the signal to the overall ττ event yield in m tot T is expected, with a reduction for m tot T 250 GeV and an enhancement otherwise. Both features may contribute to the signal inference, while the sensitivity of the analysis relies on the yield enhancement for m tot T 250 GeV, as will be discussed in more detail in Section 8.2. We note that for the φ searches presented in this paper interference effects with ττ backgrounds, e.g., from Z/γ * → ττ production are not an issue due to the different spin configurations of the ττ final states.

Common processing
The PDF4LHC15 [139] (NNPDF3.1 [140] Parton showering and hadronization, as well as the τ lepton decays, are modelled using the PYTHIA event generator [144], where versions 8.212 and 8.226 are used for the simulation of the data taken in 2016, and version 8.230 is used for the data taken in 2017-2018. For all simulated events, additional inclusive inelastic pp collisions generated with PYTHIA are added according to the expected PU profile in data. All events generated are passed through a GEANT4based [145] simulation of the CMS detector and reconstructed using the same version of the CMS event reconstruction software used for the data.

Corrections to the model
The capability of the model to describe the data is monitored in various control regions orthogonal to the signal and background classes, and corrections and corresponding uncertainties are derived where necessary. All corrections that have been applied to the model are described in the following. Their uncertainties are discussed in Section 7.
The following corrections apply equally to simulated and τ-embedded events, where the τ decay is also simulated. Since the simulation part of τ-embedded events happens under detector conditions that are different from the case of fully simulated events, corrections and related uncertainties may differ, as detailed in Ref. [110]. Corrections are derived for residual differences in the efficiencies of the selected triggers, differences in the electron and muon tracking efficiencies, and in the efficiencies of the identification and isolation requirements for electrons and muons. These corrections are obtained in bins of p T and η of the corresponding lepton, using the "tag-and-probe" method, as described in Ref. [146], with Z/γ * → ee and Z/γ * → µµ events. They usually do not amount to more than a few percent. The electron energy scale is adjusted to the scale measured in data using the Z boson mass peak in Z/γ * → ee events.
In a similar way, corrections are obtained for the efficiency of triggering on the τ h decay signature and for the τ h identification efficiency. The trigger efficiency corrections are obtained from parametric fits to the trigger efficiency, as a function of p T , derived for simulated events and data. The identification efficiency corrections are also derived as a function of the p T of the τ h candidate. For p τ h T > 40 GeV, a correction is moreover derived for each τ h decay mode individually, which is used only in the τ h τ h final state. For each data-taking year and each τ h decay mode, corrections to the energy scale of the τ h candidates and of electrons misidentified as τ h candidates are derived from likelihood scans of discriminating observables, such as the reconstructed τ h candidate mass. For muons misidentified as τ h candidates, the energy scale correction has been checked to be negligible.
Corrections are applied to the magnitude and resolution of p miss T in τ-embedded events to account for rare cases of an incomplete removal of the energy deposits from the muons that are replaced by simulated τ decays during the embedding procedure. These corrections are derived by comparing p miss T in τ-embedded events with fully simulated events.
The following corrections only apply to fully simulated events. During the 2016 and 2017 data taking, a gradual shift in the timing of the inputs of the ECAL L1 trigger in the region at |η| > 2.0 caused a specific trigger inefficiency [88]. For events containing an electron (a jet) with p T larger than ≈50 (≈100) GeV, in the region of 2.5 < |η| < 3.0 the efficiency loss is 10-20%, depending on p T , η, and time. Corresponding corrections have been derived from data and applied to the simulation, where this effect is not present.
The energies of jets are corrected to the expected response of the jet at the stable hadron level, using corrections measured in bins of the jet p T and η. These corrections are usually less than 10-15%. Residual data-to-simulation corrections are applied to the simulated event samples. They usually range from subpercent level at high jet p T in the central part of the detector to a few percent in the forward region. The energy resolution of simulated jets is also adjusted to match the resolution in data. A correction is applied to the direction and magnitude of p miss T based on differences between estimates of the hadronic recoil in Z/γ * → µµ events in data and simulation. This correction is applied to the simulated Z/γ * → , W+ jets, h, and φ signal events, where a hadronic recoil against a single particle is well defined. The efficiencies for genuine and misidentified b jets to pass the working points of the b jet identification discriminant, as given in Section 5, are determined from data, using tt events for genuine b jets and jet-associated Z boson production for jets originating from light-flavour quarks. Data-tosimulation corrections are obtained for these efficiencies and used to correct the number of b jets in the simulation.
Data-to-simulation corrections are further applied to simulated events where an electron (muon) is reconstructed as a τ h candidate, to account for residual differences in the e(µ) → τ h misidentification rate between data and simulation. In a similar way, a correction is applied to account for residual differences in the µ → e misidentification rate between data and simulation.
The dilepton mass and p T spectra in simulated Z/γ * → events are corrected to better match the data. To do this, the dilepton mass and p T are measured in data and simulation in µµ events, and the simulated events are corrected to match the spectra in data. In addition, all simulated tt events are weighted to better match the top quark p T distribution observed in data [147]. The overall normalization of tt events is constrained using the tt control region described in Section 5.

Systematic uncertainties
The uncertainty model used for the analysis comprises theoretical and experimental uncertainties, and uncertainties due to the limited population of template distributions available for the background model. The last group of uncertainties is incorporated for each bin of each corresponding template individually following the approach proposed in Refs. [148,149]. All other uncertainties lead to correlated changes across bins either in the form of normalization changes or as general nontrivial shape-altering variations. Depending on the way they are derived, correlations may also arise across data-taking years, samples, or individual uncertainties.

Uncertainties related to the τ-embedding method or the simulation
The following uncertainties, related to the reconstruction of electrons, muons, and τ h candidates after selection, apply to simulated and τ-embedded events. Unless stated otherwise, they are partially correlated across τ-embedded and simulated events.

Uncertainties common to signal and background events
Uncertainties in the identification efficiency of electrons and muons amount to 2%, correlated across all years. Since no significant dependence on the p T or η of each corresponding lepton is observed, these uncertainties are introduced as normalization uncertainties. A similar reasoning applies to uncertainties in the electron and muon trigger efficiencies, which also amount to 2% each. Because of differences in the online selections they are treated as uncorrelated for single-lepton and dilepton triggers. This may result in shape-altering effects in the overall model, since the two trigger types act on different ranges of lepton p T .
For fully simulated events, an uncertainty in the electron energy scale is derived from the calibration of ECAL crystals, and applied on an event-by-event basis. For τ-embedded events, uncertainties of 0.5-1.25%, determined separately for the ECAL barrel and endcap regions, are derived for the corrections described in Section 6.5. Because of the varying detector conditions, and the different ways the uncertainties are determined, they are treated as uncorrelated across simulated and τ-embedded events. They lead to shape-altering variations and are treated as correlated across data-taking years. The muon momentum is precisely known, and a variation within the expected uncertainties was verified to have no effect on the analysis.
Uncertainties in the τ h identification efficiency are between 3-9% in bins of τ h lepton p T . These are dominated by statistical uncertainties and are, therefore, treated as uncorrelated across decay modes, p T bins, and data-taking years. The same is true for the uncertainties in the τ h energy scale, which amount to 0.2-1.1%, depending on the τ h lepton p T and decay mode. For the energy scale of electrons misidentified as τ h candidates, the uncertainties are 1-6.5%. All τ h energy scale uncertainties are also treated as uncorrelated across data-taking years as they are predominantly statistical in nature. The uncertainty in the energy scale of muons misidentified as τ h is 1%. Uncertainties in the τ h trigger efficiencies are typically O(10%), depending on the τ h lepton p T . They are obtained from parametric fits to data and simulation, and are treated as uncorrelated across triggers and data-taking years. All uncertainties discussed in this paragraph lead to shape-altering variations.
Four further sources of uncertainty are considered for τ-embedded events. A 4% normalization uncertainty arises from the efficiency of the µµ selection in data, which is unfolded during the τ-embedding procedure. Most of this uncertainty originates from the triggers used for selection. Since the trigger configurations changed over time, this uncertainty is treated as uncorrelated across data-taking years. An additional shape uncertainty is introduced to quantify the consistency of the embedding method in a sample of µµ events. For this purpose, dedicated event samples are produced where the muons selected in data are replaced by simulated muons instead of τ lepton decays. These events are compared with the originally selected µµ events in data and residual differences in the µµ mass and p T spectra are used as uncertainties. Another shape-and normalization-altering uncertainty in the yield of tt → µµ + X decays, which are part of the τ-embedded event samples, ranges from subpercent level to 8%, depending on the event composition of the model. For this uncertainty, the number and shape of tt events contained in the τ-embedded event samples are estimated from simulation, where the corresponding decay has been selected at the parton level. This estimate is then varied by ±10% to account for the tt cross section and acceptance uncertainties. Finally, an uncertainty in the p miss T correction for the τ-embedded events described in Section 6.5 is applied. Since this correction is derived from a comparison with fully simulated events, this uncertainty is related to the imperfect p miss T reconstruction in the simulation.
For fully simulated events, the following additional uncertainties arise. Uncertainties in the e(µ) → τ h misidentification rate are 18-40% for electrons and 7-65% for muons, depending on the η of the τ h candidate. These uncertainties apply only to simulated Z/γ * → ee and Z/γ * → µµ events, which are of marginal importance for the analysis. The same is true for the uncertainty in the reweighting in the Z/γ * → dilepton mass and p T , discussed in Section 6.5, which is typically smaller than 1%. A normalization uncertainty due to the timing shift of the inputs of the ECAL L1 trigger described in Section 6.5 amounts to 2-3%.
Uncertainties in the energy calibration and resolution of jets are applied with different correlations depending on their sources, which arise from the statistical limitations of the measurements used for calibration, the time-dependence of the energy measurements in data due to detector ageing, and bias corrections introduced to cover residual differences between simulation and data. They range from subpercent level to O(10%), depending on the kinematic properties of the jets in the event. Similar uncertainties, with similar ranges, are applied for the identification rates for b jets and for the misidentification rates for light-flavour quark or gluon jets.
Depending on the process under consideration, two independent uncertainties in p miss T are applied. For processes that are subject to recoil corrections, i.e. Z/γ * → , W+ jets, h, or φ production, uncertainties in the calibration and resolution of the hadronic recoil are applied; they typically result in changes to the event yields ranging from 0-5%. For all other processes, an uncertainty in p miss T is derived from the amount of energy carried by unclustered particle candidates, which are not contained in jets, in the event [103]. This uncertainty typically results in changes to the event yields ranging from 0-10%.
The integrated luminosities of the 2016, 2017, and 2018 data-taking periods are individually known with uncertainties in the 1.2-2.5% range [150][151][152], while the total integrated luminosity for the years 2016-2018 has an uncertainty of 1.6%; the improvement in precision reflects the (uncorrelated) time evolution of some systematic effects. Uncertainties in the predictions of the normalizations of all simulated processes amount to 4% for Z/γ * → and W+ jets production [123], 6% for tt production [124,125], and 5% for diboson and single t quark production [125][126][127], where used in the analyses. These uncertainties are correlated across datataking years. A shape-altering uncertainty is derived in the reweighting of the top quark p T described in Section 6.5 by applying the correction twice or not applying it at all. This uncertainty has only a very small effect on the final discriminant.

Uncertainties in the signal modelling
Theoretical uncertainties in the acceptance of bbφ signal events are obtained from variations of the renormalization (µ R ) and factorization (µ F ) scales, the h damp factor, and the PDFs. The scale uncertainty is obtained from the envelope of the six variations of µ R and µ F by factors of 0.5 and 2, omitting the variations where one scale is multiplied by 2 and the corresponding other scale by 0.5, as recommended in Ref. [153]. The scale h damp is varied by factors of 1/ √ 2 and √ 2. The uncertainty from the variation of µ R and µ F , and the uncertainty from the variation of h damp are added linearly, following the recommendation in Ref. [153], resulting in an overall uncertainty that ranges from 1-8% (1-5%) for the b tag ("no b tag") categories depending on the tested mass. The uncertainties due to PDF variations and the uncertainty in α S are obtained following the PDF4LHC recommendations, taking the root mean square of the variation of the results when using different replicas of the default PDF4LHC sets as described, e.g., in Ref. [139]. They range from 1-2%.
Uncertainties in the acceptance of the ggφ process are also obtained from variations of µ R , µ F , and h damp . The µ R and µ F scales are varied as described above for the bbφ process, whereas the h damp scale is varied by factors of 0.5 and 2 as suggested in Ref. [135]. The influence of the former (latter) variation on the signal acceptance amounts to 20% (35%) for the smallest m φ values. For larger m φ values, the variation is at the subpercent level. In both cases the uncertainties also result in shape-altering effects in the overall model.
For the parameter scan in the MSSM interpretations, theoretical uncertainties in the ggφ and bbφ cross sections are included, as described in Ref. [86]. This includes uncertainties in the µ R and µ F scales, PDFs, and α S . The uncertainties are evaluated separately for each m A -tan β point under consideration. They are typically 5-20% (10-25%) for ggφ (bbφ) production.
Several sources of theoretical uncertainty in the U 1 signal prediction are included. The uncertainty due to the µ R and µ F scale variations is about 15%. The uncertainties due to the PDFs and α S variations are about 15 and 4%, respectively. The Q match and parton shower uncertainties affect the signal acceptances in the "b tag" categories, with magnitudes of about 11 and 1% respectively, and in the "no b tag" categories, with magnitudes of 5 and 6% respectively. The uncertainty in the β sτ L parameter is estimated by varying the coupling strength by the uncertainties obtained in the fit presented in Ref. [72] and summarized in Table 1. The resulting uncertainty varies the signal yields by 4-12%. The uncertainty in the signal acceptance due to the choice of flavour scheme is estimated by comparing the predictions in the 4FS and 5FS calculations, which mainly affect the N b-jets distribution. The resulting uncertainty has a magnitude of 25% (18%) for the "b tag" ("no b tag") categories.
For all results shown in the following, the expectation for SM Higgs boson production is included in the model used for the statistical inference of the signal. Uncertainties due to different choices of µ R and µ F for the calculation of the production cross section of the SM Higgs boson amount to 3.9% for ggh, 0.4% for VBF, 2.8% for Zh, and 0.5% for Wh production [128-130, 154, 155]; uncertainties due to different choices for the PDFs and α S amount to 3.2%, 2.1%, 1.6%, and 1.9% for these four production modes, respectively.

Uncertainties related to jets misidentified as an electron, muon, or τ h candidate
For the F F method, the following uncertainties apply. The F i F and their corrections are subject to statistical fluctuations in each corresponding DR i and simulation. The corresponding uncertainties are split into a normalization and a shape-altering part and propagated into the final discriminant. They are typically 1-10% and are treated as uncorrelated across the kinematic and topological bins where they are derived. An additional uncertainty is defined by varying the choice of the functional form for the parametric fits.
Uncertainties are also applied to cover residual corrections and extrapolation factors, varying from a few percent to O(10%), depending on the kinematic properties of the τ h candidate and the topology of the event. These are both normalization and shape-altering uncertainties.
An additional source of uncertainty concerns the subtraction of processes other than the enriched process in each corresponding DR i . These are subtracted from the data using simulated or τ-embedded events. The combined shape of the events to be removed is varied by 10%, and the measurements are repeated. The impacts of these variations are then propagated to the final discriminant as shape-altering uncertainties.
An uncertainty in the estimation of the three main background fractions in the AR is estimated from a variation of each individual contribution by 10%, increasing or decreasing the remaining fractions such that the sum of all contributions remains unchanged. The amount of variation is motivated by the uncertainty in the production cross sections and acceptances of the involved processes and the constraint on the process composition that can be clearly obtained from the AR. The effect of this variation is found to be very small, since usually one of the contributions dominates the event composition in the AR.
Since the background from QCD multijet events in the eµ final state is also determined from a DR, uncertainties that account for the statistical uncertainty in the data and the subtracted backgrounds in this DR are applied in a similar way. These uncertainties amount to 2-4%. In addition, this background is subject to uncertainties related to the extrapolations from the DR to the corresponding SRs. These uncertainties are O(10%) depending on p e T , p µ T , and N b-jets . Because of their mostly statistical nature, all uncertainties related to the F F and SS methods are treated as uncorrelated across data-taking years.
In the eµ final state, the subdominant contribution to the jet → and µ → e backgrounds is estimated from simulation. Uncertainties in the simulated jet → e and jet → µ misidentification rates are 10 and 12%, respectively. They are treated as correlated across data-taking years. The uncertainty in the µ → e misidentification rate is 15-45%, and is treated as uncorrelated across data-taking years since it is mostly statistical in nature. A summary of all systematic uncertainties that have been discussed in this section is given in Table 6, in which we also state the correlations between the data-taking years.

Results
The statistical model used to infer the signal from the data is defined by an extended binned likelihood of the form where i labels the bins of the discriminating distributions of all categories, split by ττ final state and data-taking year. The function P (k i | ∑ µ s S si ({θ j }) + ∑ B bi ({θ j })) corresponds to the Poisson probability to observe k i events in bin i for a prediction of ∑ µ s S si signal and ∑ B bi background events. The predictions for S si and B bi are obtained from the signal and background models discussed in Section 6. The parameters µ s act as linear scaling parameters of the corresponding signal yields S s . Systematic uncertainties are incorporated in the form of penalty terms for additional nuisance parameters {θ j } in the likelihood, appearing as a product with O(10%) -predefined probability density functions C( θ j |θ j ), where θ j corresponds to the nominal value for θ j . The predefined uncertainties in the θ j , as discussed in Section 7, may be constrained by the fit to the data.
The test statistic used for the inference of the signal is the profile likelihood ratio, as discussed in Refs. [156,157]: where one or more parameters µ s are the parameters of interest (POIs) andμ s ,θ j,µ s , andθ j,μ s are the values of the given parameters that maximize the corresponding likelihood. The index of q µ s indicates that the test statistic is evaluated for a fixed value of µ s . In the large number limit, the sampling distribution of q µ s can be approximated by analytic functions, from which the expected median and central intervals can be obtained as described in Ref. [158]. The signal is inferred from the data in three different ways: i) the model-independent φ search features a signal model for a single narrow resonance φ; ii) for the search for vector leptoquarks, the data are interpreted in terms of the nonresonant U 1 t-channel exchange; iii) the interpretation of the data in terms of MSSM benchmark scenarios relies on three resonances in the ττ mass spectrum with mass values and rates determined by the parameters of the corresponding scenario.
In all cases the tt control region, as defined in Section 5 and shown in Figs. 5-6 is used to constrain the normalization of tt events and all tt related uncertainties. Detailed descriptions of the specific statistical procedures and the results obtained in each case are given in the following sections.

Model-independent φ search
For the model-independent φ search, we investigate ggφ and bbφ production corresponding to two independent POIs µ ggφ and µ bbφ in the likelihood of Eq.       The solid histograms show the stacked background predictions after a signal-plus-background fit to the data for m φ = 100 GeV. The best fit ggφ signal is shown by the red line. The total background prediction as estimated from a background-only fit to the data is shown by the dashed blue line for comparison. For all histograms, the bin contents show the event yields divided by the bin widths. The lower panel shows the ratio of the data to the background expectation after the signal-plus-background fit to the data. The signal-plus-background and background-only fit predictions are shown by the solid red and dashed blue lines, respectively, which are also shown relative to the background expectation obtained from the signal-plusbackground fit to the data. expected m tot T (m τ τ ) distributions for a ggφ or bbφ signal with m φ = 1200 (100) GeV are also shown. Figure 10 shows the expected and observed 95% confidence level (CL) upper limits on the product of the cross sections and branching fraction for the decay into τ leptons for ggφ and bbφ production in a mass range of 60 ≤ m φ ≤ 3500 GeV. These limits have been obtained following the modified frequentist approach described in Refs. [159,160]. When setting the limit in one production mode the POI of the other production mode is profiled. The limits are shown with a separation into the low-mass (m φ < 250 GeV) and high-mass (m φ ≥ 250 GeV) regions of the search.
The expected limits in the absence of a signal span four orders of magnitude between ≈10 pb (at m φ = 60 GeV) and ≈0.3 fb (at m φ = 3.5 TeV) for both production modes, with a falling slope for increasing values of m φ . In general, the observation falls within the central 95% interval of the expectation. For the low-mass search, the largest deviation from the expectation is observed for ggφ production at m φ = 100 GeV with a local (global) p-value equivalent to 3.1 (2.7) standard deviations (s.d.). To turn the local into a global p-value, a number N trial of pseudo-data from the input distributions of the background model to the maximum likelihood fit is created. For each mass hypothesis in consideration, a fit of the signal model to these pseudo-data is performed and the fraction of cases, where the outcome of these fits with the maximal significance exceeds the observed significance, with respect to N trial is determined. Finally, the local p-value is reduced by this fraction. The best fit value of the product of the cross section with the branching fraction for the decay into τ leptons is σ ggφ B(φ → ττ) = (5.8 ± 2.5 2.0 ) pb. The excess at m φ = 100 GeV exhibits a p-value of 50% (58%) for the compatibility across ττ final states (datataking years). Within the resolution of m τ τ this coincides with a similar excess observed in a previous search for low-mass resonances by the CMS Collaboration in the γγ final state, where

Low-mass
High-mass Figure 10: Expected and observed 95% CL upper limits on the product of the cross sections and branching fraction for the decay into τ leptons for (left) ggφ and (right) bbφ production in a mass range of 60 ≤ m φ ≤ 3500 GeV, in addition to H (125). The expected median of the exclusion limit in the absence of signal is shown by the dashed line. The dark green and bright yellow bands indicate the central 68% and 95% intervals for the expected exclusion limit. The black dots correspond to the observed limits. The peak in the expected ggφ limit emerges from the loss of sensitivity around 90 GeV due to the background from Z/γ * → ττ events.         the smallest local p-value corresponds to a significance of 2.8 s.d. for a mass of 95.3 GeV [161]. The local (global) significance for the ττ search evaluated at m φ = 95 GeV is 2.6 (2.3) s.d. and the best fit value of the product of the cross section with the branching fraction for the decay into τ leptons is σ ggφ B(φ → ττ) = (7.8 ± 3.9 3.1 ) pb. For the high-mass search, the largest deviation from the expectation is observed for ggφ production at m φ = 1.2 TeV with a local (global) p-value equivalent to 2.8 (2.2) s.d., where the best fit value of the product of the cross section with the branching fraction for the decay into τ leptons is σ ggφ B(φ → ττ) = (3.1 ± 1.0 1.1 ) fb. The excess at m φ = 1.2 TeV exhibits a p-value of 11% (63%) for the compatibility across ττ final states (data-taking years). For bbφ production, no deviation from the expectation beyond the level of 2 s.d. is observed. Figure 11 shows the same results in the form of maximum likelihood estimates with 68% and 95% CL contours obtained from scans of the signal likelihood along the ggφ and bbφ cross sections, for selected values of m φ between 60 GeV and 3.5 TeV.

Search for vector leptoquarks
The inputs for the search for vector leptoquarks are the binned template distributions of m tot T in the categories shown in Fig. 5 resulting in 51 input distributions for signal extraction, for the years 2016-2018. Based on these inputs a signal is searched for in the range of 1 < m U < 5 TeV.
Due to the destructive interference with the Z/γ * → ττ process discussed in Section 6.4, depending on m tot T , a signal from U 1 t-channel exchange may result in an enhancement or a reduction of the yields as expected from the SM. A typical example of this effect for a signal with m U = 1 TeV, g U = 1.5, for the VLQ BM 1 scenario is shown in Fig. 7 (right). There, a reduction in yield (with respect to the SM) is expected for m tot T 250 GeV and an enhancement for 250 m tot T 1000 GeV. In principle both effects contribute to the sensitivity of the analysis to the signal. However, the region of m tot T 250 GeV, which features the deficit, suffers from large backgrounds, reducing the contribution of this effect to the overall sensitivity. Studies confirm that the sensitivity of the analysis to the signal is driven by the high m tot T region, which consistently features an enhancement of the yields expected from the SM. Since the interference reduces this enhancement it also reduces the sensitivity of the analysis, compared to a signal without interference.
No statistically significant signal is observed and 95% CL upper limits on g U are derived for the VLQ BM 1 and 2 scenarios, as shown in Fig. 12, again following the modified frequentist approach as for the previously discussed search. The expected sensitivity of the analysis drops for increasing values of m U following a linear progression with values from g U = 1.3 (0.8) to 5.6 (3.2) for the VLQ BM 1 (2) scenario. The observed limits fall within the central 95% intervals for the expected limits in the absence of signal. The expected limits are also within the 95% confidence interval of the best fit results reported by Ref. [72], indicating that the search is sensitive to a portion of the parameter space that can explain the b physics anomalies.

95% CL excluded:
Observed 68% expected Expected 95% expected 2 scenarios, in a mass range of 1 < m U < 5 TeV. The expected median of the exclusion limit in the absence of signal is shown by the dashed line. The dark and bright grey bands indicate the central 68% and 95% intervals of the expected exclusion limit. The observed excluded parameter space is indicated by the coloured blue area. For both scenarios, the 95% confidence interval for the preferred region from the global fit presented in Ref. [72] is also shown by the green shaded area.

MSSM interpretation of the data
For the interpretation of the data in MSSM benchmark scenarios, the signal is based on the binned distributions of m tot T in the categories shown in Fig. 5, complemented by distributions of the NN output function used for the stage-0 simplified template cross section measurement of Ref. [109], as discussed in Section 5.2, resulting in 129 input distributions for signal extraction.
In the MSSM, the signal constitutes a multiresonance structure with contributions from h, H, and A bosons. For the scenarios chosen for this paper h is associated with H(125). Any MSSM prediction has to match the observed properties of H(125), in particular its mass, cross sections for various production modes, and branching fraction for the decay into τ leptons. For the benchmark scenarios summarized in Ref. [86], all model parameters have been chosen such that m h is compatible with the observed H(125) mass of 125.38 GeV [8], within an uncertainty of ±3 GeV in most of the provided parameter space. The uncertainty of ±3 GeV in the prediction of m h is supposed to reflect the unknown effect of higher-order corrections, as discussed in Ref. [162]. The value of m h is allowed to vary within these boundaries, according to a flat distribution. For the interpretation this is taken into account by simulating the h signal at the observed H(125) mass. For h production, the modes via ggh, b associated production (bbh), VBF, and Vh production are included, and all cross sections and the branching fraction for the decay into τ leptons are scaled according to the MSSM predictions. To remove any dependencies of these predictions on the exact value of m h , they are scaled to the expectation for m h = 125.38 GeV, following the prescription of Ref. [86]. For A and H production, gluon fusion (ggA, ggH) and b associated production (bbA, bbH) are included.
All kinematic distributions are modelled within the accuracies discussed in Section 6.4. In particular, the H (A) boson p T spectra in ggH (ggA) production are modelled as a function of tan β for each tested value of m A , resulting in a softer progression for increasing values of tan β. In the "no b tag" categories for m τ τ > 250 GeV the h signal is expected to be negligible so it is dropped from the signal templates. A summary of the association of signals to the templates used for signal extraction is given in Table 7. To interpolate the simulated mass points to the exact predicted values of m H , a linear template morphing algorithm, as described in Ref. [163], is used.
The m A -tan β plane is scanned and for each tested point in (m A , tan β), the CL s [160] value is calculated. Those points where CL s falls below 5% define the 95% CL exclusion contour for the benchmark scenario under consideration. The underlying test compares the MSSM hypothesis, with signal contributions for h (S h where for brevity the dependence on the nuisance parameters {θ j } has been omitted. Equation (13) represents a nested likelihood model from which the MSSM hypothesis (with µ = 1) evolves through continuous transformation from the SM hypothesis (with µ = 0). We note that the only physically meaningful hypotheses in Eq. (13) correspond to µ = 0 and 1. On the other hand, in the large number limit this construction allows the application of the asymptotic formulas given in Ref. [158], as analytic estimates of the sampling distributions for the MSSM and SM hypotheses, when using the profile likelihood ratio given in Eq. (12) as the test statistic.
We have verified the validity of the large number limit for masses of m A > 1 TeV with the help of ensemble tests. Since we are using the same template distributions for S SM and S h the transition from µ = 0 to 1 corresponds to a normalization change of the signal contribution related to H (125), only.  [164][165][166][167][168][169][170][171]. Branching fractions for the decay into τ leptons and other final states have been obtained from a combination of the codes FEYNHIGGS (and HDECAY [172,173]) for the M 125 h, EFT (M 125 h ) scenario, as described in Ref. [86] following the prescriptions given in Refs. [153,174,175].
Inclusive cross sections for the production via ggφ have been calculated using the program SUSHI 1.7.0 [176,177], including NLO corrections in α S for the t-and b-quark contributions to the cross section [178,179], NNLO corrections in α S in the heavy t quark limit, for the t quark contribution [180][181][182][183][184], and next-to-NNLO contributions in α S for h production [185][186][187]. Electroweak corrections mediated by light-flavour quarks are included at two-loop accuracy reweighting the SM results of Refs. [188,189]. Contributions from squarks and gluinos are included at NLO precision in α S following Refs. [190][191][192]. The tan β-enhanced SUSY contributions to the Higgs-b couplings have been resummed using the one-loop ∆ b terms from Ref. [193] as provided by FEYNHIGGS. Uncertainties in these ∆ b terms which range ≈10% are not included in the overall uncertainties in the predictions as they are subdominant with respect to the other theoretical uncertainties.
For bbH production, cross sections have been calculated for the SM Higgs boson as a function of its mass, based on soft-collinear effective theory [194,195] which combines the merits of both the 4FS [196,197] and 5FS [198,199] calculations. These cross sections coincide with the results of the so-called "fixed order plus next-to-leading log" approach of Refs. [200,201]. The pure t-and loop-induced tb-interference contributions are separately reweighted with effective  h, EFT scenario, the dashed blue line indicates the threshold at m A = 2m t whereby the A → tt decay starts to influence the A → ττ branching fraction. The H → ττ branching fraction is influenced more gradually close to this threshold since A and H are not completely degenerate in mass.
Higgs couplings, using an effective mixing angle α, and including the resummation of tan βenhanced SUSY contributions as in the ggφ case. The same SM cross sections are also used to obtain the reweighted cross section for bbA production. A more detailed discussion is given in Ref. [86]. All Higgs boson masses, effective mixing angles α, Yukawa couplings, branching fractions, cross sections, and their uncertainties, which are included for the exclusion contours, are obtained from Ref. [87].
In the figure, the exclusion sensitivity, estimated from the expected median in the absence of a signal, is indicated by the dashed black line. We note that the central 68% and 95% intervals, also given for the exclusion sensitivity, should not be misinterpreted as an uncertainty in the analysis, but they rather reflect the variation of the expected signal yield in the probed parameter space of the chosen benchmark scenarios. For the M 125 h, EFT scenario the sensitivity sharply drops at m A = 2 m t , caused by a drop of the branching fractions for the decay of A and H into τ leptons where the A and H decays into two on-shell t quarks become kinematically accessible. The distinct boundary is related to the fact that in FEYNHIGGS, which is used for the calculation of all branching fractions for this benchmark scenario, only the decay into on-shell tt pairs is implemented. The parameter space of each benchmark scenario that is excluded at 95% CL by the data is indicated by the coloured blue area.
Both scenarios are excluded at 95% CL for m A 350 GeV. The local excess observed at 1.2 TeV causes the deviation of the observed exclusion from the expectation. For m A 250 GeV, most of the ggH/ggA events do not enter the "no b tag" categories due to the m τ τ > 250 GeV requirement, although these events still contribute to the signal yields in the NN categories. In this parameter space the sensitivity to the MSSM is driven by the measurements of the H (125) production rates, while the sensitivity to the H and A enters mainly via the bbφ signal in the "b tag" categories, especially for increasing values of tan β.

Summary
Three searches have been presented for signatures of physics beyond the standard model (SM) in ττ final states in proton-proton collisions at the LHC, using a data sample collected with the CMS detector at √ s = 13 TeV, corresponding to an integrated luminosity of 138 fb −1 . Upper limits at 95% confidence level (CL) have been set on the products of the branching fraction for the decay into τ leptons and the cross sections for the production of a resonance φ in addition to the observed Higgs boson via gluon fusion (ggφ) or in association with b quarks, ranging from O(10 pb) for a mass of 60 GeV to 0.3 fb for a mass of 3.5 TeV each. The data reveal two excesses for ggφ production with local p-values equivalent to about three standard deviations at m φ = 0.1 and 1.2 TeV. Within the resolution of the reconstructed invariant mass of the ττ system, the excess at 100 GeV coincides with a similar excess observed in a previous search for low-mass resonances by the CMS Collaboration in the γγ final state at a mass of ≈95 GeV. In a search for t-channel exchange of a vector leptoquark U 1 , 95% CL upper limits are set on the U 1 coupling to quarks and τ leptons ranging from 1 for a mass of 1 TeV to 6 for a mass of 5 TeV, depending on the scenario. The search is sensitive to and excludes a portion of the parameter space that can explain the b physics anomalies. In the interpretation of the  [4] 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 √ s = 7 and 8 TeV", JHEP 08 (2016) 045, doi:10.1007/JHEP08(2016)045, arXiv:1606.02266.