Main

In the standard model of particle physics (SM)1,2,3,4, three lepton families (flavours) exist. The number of leptons of each family is conserved in weak interactions, and violation of this assumption is known as lepton flavour violation (LFV). No fundamental principles forbid LFV processes in the SM. The phenomenon of neutrino oscillations, where neutrinos (the neutral leptons) of one flavour transform into those of another5,6, indicates that neutrinos have mass and LFV processes do occur in nature. The mechanisms responsible for neutrinos acquiring mass and weak interactions violating lepton flavour conservation remain unknown. More experimental data are needed to constrain and guide possible generalizations of the SM explaining these phenomena.

An observation of LFV in charged-lepton interactions would be an unambiguous sign of new physics. In particular, decays of the Z boson into a light lepton (electron or muon) and a τ lepton at colliders are of experimental interest. The abundance of Z bosons produced at the Large Hadron Collider (LHC) offers the opportunity to strongly constrain potential LFV Z → eτ or Z → μτ interactions, in particular those proportional to the centre-of-mass energy of the decay7. Moreover, Z → eτ, μτ decays are less constrained by low-energy experiments than Z → eμ decays. According to current knowledge, these decays can occur via neutrino mixing but are too rare to be detected. Only 1 in approximately 1054 Z bosons would decay into a muon and a τ lepton8. An observation of such decays would therefore require new theoretical explanations. For example, theories predicting the existence of heavy neutrinos9 provide a fundamental understanding of the observed tiny masses and large mixing of SM neutrinos. In such theories, up to 1 in 105 Z bosons would be expected to undergo an LFV decay involving τ leptons. The ATLAS experiment can test the predictions of such theories by observing or setting ever more stringent constraints on LFV Z-boson decays.

Constraints on the branching fractions (\({\mathcal{B}}\)) of the LFV decays of the Z boson involving a τ lepton have been set by the experiments at the Large Electron–Positron Collider (LEP): \({\mathcal{B}}(Z\to e\tau )<9.8\times 1{0}^{-6}\) (ref. 10) and \({\mathcal{B}}(Z\to \mu \tau )<1.2\times 1{0}^{-5}\) (ref. 11) at the 95% confidence level (CL). The ATLAS experiment12 at the LHC has set constraints \({\mathcal{B}}(Z\to e\tau )<5.8\times 1{0}^{-5}\) at 95% CL using part of the Run 2 data and \({\mathcal{B}}(Z\to \mu \tau )<1.3\times 1{0}^{-5}\) using the Run 1 data and a subset of the Run 2 data13.

This work uses proton–proton (pp) collision data collected by the ATLAS experiment during Run 2 of the LHC, containing about eight billion Z-boson decays. Only events with a τ lepton that decays hadronically are considered. Neural network (NN) classifiers are used in a novel way for optimal discrimination of signal from background, and to achieve improved sensitivity in the search for LFV effects in the data using a binned maximum-likelihood fit. The result for the μτ channel is combined with a previous LHC Run 1 result to further improve the sensitivity. These results set constraints on LFV Z-boson decays involving τ leptons that supersede the most stringent ones set by the LEP experiments more than two decades ago.

The ATLAS experiment and data sample

To record and analyse the LHC pp collisions, the ATLAS experiment uses a multipurpose particle detector with a forward–backward symmetric cylindrical geometry and a near 4π coverage in solid angle12,14,15. It consists of an inner tracking detector surrounded by a superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer.

The search uses the complete dataset of pp collision events at a centre-of-mass energy of \(\sqrt{s}=\,\text{13}\,\text{TeV}\,\) collected by the ATLAS experiment during LHC Run 2. This dataset was recorded using single-electron or single-muon triggers16 and corresponds to an integrated luminosity of 139 fb−1. For the search in the μτ channel, the results are combined with those of a previous similar search using pp collisions at \(\sqrt{s}=\,\text{8}\,\text{TeV}\,\) during LHC Run 1, corresponding to an integrated luminosity of 20.3 fb−1 (ref. 17).

Candidates for electrons18, muons19, jets20,21,22, and visible decay products of hadronic τ-lepton decays (τhad-vis)23,24 are reconstructed from energy deposits in the calorimeters and charged-particle tracks measured in the inner detector and the muon spectrometer.

Electron candidates are required to pass the Medium likelihood-based identification requirement18 and have pseudorapidity \(\left|\eta \right| <1.37\) or \(1.52<\left|\eta \right| <2.47.\) Muon candidates are required to pass the Medium identification requirement19 and have \(\left|\eta \right| <2.5.\) Both the electron and muon candidates must have transverse momentum pT > 30 GeV and satisfy the Tight isolation requirement18,19. The lower bounds on the electron and muon transverse momenta are driven by the acceptance of the trigger selection.

Quark- or gluon-initiated particle showers (jets) are reconstructed using the anti-kt algorithm20,21 with the radius parameter R = 0.4. Jets fulfilling pT > 20 GeV and \(\left|\eta \right| <2.5\) are identified as containing b hadrons if tagged by a dedicated multivariate algorithm25.

The τhad-vis candidates are reconstructed from jets with pT > 10 GeV, \(\left|\eta \right| <1.37\) or \(1.52<\left|\eta \right| <2.5\), and one or three associated tracks, referred to as ‘1-prong’ (1P) and ‘3-prong’ (3P), respectively. The τhad-vis identification is performed by a recurrent NN algorithm23, which uses calorimetric shower shapes and tracking information to discriminate true τhad-vis candidates from fake candidates from quark- or gluon-initiated jets. The τhad-vis candidates are required to pass the Tight identification selection, which has an efficiency of 60% (45%) for true 1P (3P) τhad-vis candidates, constant in the τhad-vis candidates’ transverse momentum, and a misidentification rate of 1 in 70 (700) for fake 1P (3P) candidates in dijet events. Dedicated multivariate algorithms are used to further discriminate between τhad-vis and electrons, and to calibrate the τhad-vis energy24. The τhad-vis candidate with the largest pT in each event is the selected candidate and is required to have pT > 25 GeV. Based on simulation, in Z → ℓτ decays, the τhad-vis candidate is expected to be correctly selected 98% of the time.

The missing transverse momentum \(({E}_{\,\text{T}}^{\text{miss}})\) is calculated as the negative vectorial sum of the pT of all fully reconstructed and calibrated physics objects26,27. The calculation also includes inner detector tracks that originate from the vertex associated with the hard-scattering process but are not associated with any of the reconstructed objects. The missing transverse momentum is the best proxy for the total transverse momentum of undetected particles (in particular neutrinos) in an event.

Search strategy

The Z → ℓτ → ℓτhad-vis + ν (ℓ = light lepton, e or μ) signal events have a number of key features that can be exploited to separate them from the SM background events. The signal events are characterized by their unique final state, which has exactly one ℓ and one τ lepton, with the invariant mass of the pair being compatible with the Z-boson mass. The ℓ and τ leptons carry opposite electric charges and are emitted approximately back to back in the plane transverse to the proton beam direction. Since the τ lepton is typically boosted due to the large difference between its mass and the mass of its parent Z boson, the neutrino from its decay is usually almost collinear with the visible τ-decay products. The neutrino escapes detection and is reconstructed as part of the \({E}_{\,\text{T}}^{\text{miss}\,}\) of the event. In a signal event, this is the only major source of \({E}_{\text{T}}^{\text{miss}}.\)

The major background contributions for this search are as follows: lepton-flavour-conserving Z → ττ → ℓτhad-vis + 3ν decays, where one of the τ leptons decays leptonically and the other hadronically; Z → ℓℓ decays, where one of the light leptons is misidentified as the τhad-vis candidate; events with a quark- or gluon-initiated jet that is misidentified as the τhad-vis candidate. The last of these are hereafter referred to as events with ‘fakes’ and are mostly W(→ ℓν) + jets events and purely hadronic multijet events. Other SM processes with a real ℓτhad-vis final state, such as decays of a top–antitop-quark pair, two gauge bosons or a Higgs boson, and those with a real τhad-vis and a jet misidentified as a light lepton, such as W(→ τν) + jets, are considered, although their contribution to the overall background is minor.

The signal and background events are separated by using a set of event selection criteria that help to define a signal-enhanced sample, referred to as the signal region (SR). The main selection criteria are listed in Table 1 and will be explained in the following. They are primarily based on the multiplicity of reconstructed particle candidates and the event topology, in particular the transverse masses (mT), which are defined as

$${m}_{\text{T}}(X,\,{E}_{\text{T}}^{\text{miss}})\equiv \sqrt{2 {p}_{\text{T}}(X) {E}_{\text{T}}^{\text{miss}} \left(1-\cos ({\phi }_{X}-{\phi }_{{E}_{\text{T}}^{\text{miss}}})\right)}$$
(1)

where X is either a light lepton or a τhad-vis candidate and ɸ denotes the azimuthal angle. A schematic of the expected signal and background topologies is described in Extended Data Figs. 1 and 2.

Table 1 Main selection criteria for events in the signal region
Full size table

NN binary classifiers are used to distinguish signal events from W + jets , Z → ττ and Z → ℓℓ background events. The NNs are trained on simulated events (‘Signal and background predictions’). Each individual NN is optimized to discriminate against a particular background process in a given decay channel. The input to these NNs is a mixture of low-level and high-level kinematic variables, as detailed in Methods. The low-level variables are the momentum components of the reconstructed ℓ, τhad-vis candidate and \({E}_{\,\text{T}}^{\text{miss}\,}\). The high-level variables are kinematic properties of the ℓ–τhad-vis–\({E}_{\,\text{T}}^{\text{miss}\,}\) system, such as the collinear mass mcoll(ℓ, τ), defined as the invariant mass of the ℓ–τhad-vis–ν system, where the ν is assumed to have a momentum that is equal in pT and ϕ to the measured \({E}_{\,\text{T}}^{\text{miss}\,}\) and equal in η to the τhad-vis momentum. Given the finite training-sample size, the high-level variables help the NNs to converge faster, while the NNs exploit any residual correlations between the low-level variables.

The outputs from the individual NNs are numbers between 0 and 1 that reflect the probability for an event to be a signal event; they are combined into a final discriminant, hereafter referred to as the ‘combined NN output’. The combination is parameterized by weights associated with each individual NN and optimized for discrimination among various background processes distributed differently along the range of combined NN output values, as detailed in Methods. This allows the maximum-likelihood fit to determine the background contributions more precisely, which ultimately improves the sensitivity.

Events classified by the NNs as being background-like are excluded from the SR, as indicated in Table 1. The signal acceptance times selection efficiency in the SR is 2.7% for the eτ channel and 3.0% for the μτ channel, as determined from simulated signal samples.

Signal and background predictions

Predictions for signal and background contributions to the event yield and kinematic distributions in the SR are based partly on Monte Carlo (MC) simulations and partly on the use of data in regions that are enriched in background events and do not overlap with the SR.

The signal events were simulated using PYTHIA 828 with matrix elements calculated at leading order (LO) in the strong coupling constant (αs). Parameter values for initial-state radiation, multiparton interactions and beam remnants were set according to the A14 set of tuned parameters (tune)29 with the NNPDF 2.3 LO parton distribution function (PDF) set30. Nominal signal samples were generated with a parity-conserving Zℓτ vertex and unpolarized τ leptons. Scenarios where the decays are maximally parity-violating were considered by reweighting the simulated events using TAUSPINNER31. The event weight was computed as the probability of occurrence of each generated signal event, based on its kinematics, when assuming a specific τ-polarization state (left-handed or right-handed).

Background Z → ττ events were simulated with the SHERPA 2.2.132 generator using the NNPDF 3.0 NNLO PDF set33 and next-to-leading-order (NLO) matrix elements for up to two partons, and LO matrix elements for up to four partons, calculated with the COMIX34 and OPENLOOPS35,36,37 libraries. They were matched with the SHERPA parton shower38 using the MEPS@NLO prescription39,40,41,42 with the default SHERPA tune. This set-up follows the recommendations of the SHERPA authors. Background Z → ℓℓ events were simulated using the POWHEG-BOX43 generator with NLO matrix elements and interfaced to PYTHIA 8 to model the parton showers, hadronization and underlying events. All MC samples include a detailed simulation of the ATLAS detector with GEANT 444, to produce predictions that can be compared with the data. Furthermore, simulated inelastic pp collisions, generated with PYTHIA 8 using the NNPDF 2.3 LO PDF set and the A3 tune45 were overlaid on the hard-scattering events to model the additional pp collisions occurring in the same proton bunch crossing. All simulated events were processed using the same reconstruction algorithms as used for data.

The simulation of Z-boson production is improved with a correction derived from measurements in data. The simulated pT spectra of the Z boson are reweighed to match the unfolded distribution measured by ATLAS in ref. 46. This improves the predictions of signal, Z → ττ and Z → ℓℓ events, which are simulated at different orders in αs using different generators. It also reduces the uncertainties related to missing higher orders in αs.

The predicted overall yields of signal and Z → ττ events are determined by a binned maximum-likelihood fit to data (see next section) in the SR and in a control region enhanced in Z → ττ → ℓτhad-vis + 3ν events (CRZττ), using an unconstrained fit parameter, which accounts for theoretical uncertainties in the total Z-boson production cross-section (σZ), as well as the experimental uncertainties related to the acceptance of the common ℓτhad-vis final state. The selection criteria for events in the CRZττ are the same as those for events in the SR, except that events are required to have \({m}_{\text{T}}({\tau }_{{\rm{had}}{\rm{-}}{\rm{vis}}},{E}_{\text{T}}^{\text{miss}})>\text{35}\,\text{GeV}\), \({m}_{\text{T}}(\ell ,{E}_{\text{T}}^{\text{miss}})<\text{40}\,\text{GeV}\) and 70 GeV < mcoll(ℓ, τ) < 110 GeV.

A much smaller contribution to the total background originates from Z → ℓℓ events. Their predicted overall yield is based on the measured value of σZ (ref. 47) times the measured integrated luminosity. The uncertainty in the measurement is taken into account. The predicted rates of misidentifying electrons and muons in Z → ℓℓ events as 1P τhad-vis candidates are corrected using data in a region enriched in Z → ℓℓ events and orthogonal to the SR, where the last selection criterion in Table 1 is inverted and the outputs of the NN classifiers optimized to reject Z → ττ, and W + jets events are required to be greater than 0.8. The corrections are derived as functions of pT and \(\left|\eta \right|\) of the τhad-vis candidate. Statistical uncertainties in the correction are considered.

Events with fakes are one of the dominant contributions to the background, and are estimated from data using the ‘fake-factor method’, which is described in ref. 13. A fake factor is defined as the ratio of the number of events with a fake τhad-vis candidate passing the Tight τhad-vis identification requirement to those failing it. Four fake factors, one for each of the most important backgrounds with fakes (W(→ ℓν) + jets, multijet, Z(→ ℓℓ) + jets and \(t\bar{t}\) events), are measured in data in four corresponding fakes-enriched regions. Each of these regions has a dominant contribution from one of the four targeted backgrounds with fakes. These regions do not overlap with any of the regions used in the final maximum-likelihood fit. The purity of the multijet-enriched region is improved by introducing two additional selection criteria: events must have a same-sign charged ℓ–τhad-vis pair and \({m}_{\text{T}}(\ell ,{E}_{\text{T}}^{\text{miss}})<\text{40}\,\text{GeV}\). The fake factors are measured as functions of the transverse momentum of the τhad-vis candidate, separately for eτ and μτ events and for events with 1P or 3P τhad-vis candidates.

The number of events with a fake 1P or 3P τhad-vis candidate in a given pT range in the SR or CRZττ is estimated by the number of events with a τhad-vis candidate failing the Tight identification requirement, but otherwise satisfying all other selection criteria for that region, multiplied by an average of the fake factors. To calculate this average, the fake factors are summed with weights equal to the expected relative contribution of the corresponding background to the total yield of events in the region with the inverted identification requirement. This approach is used to model the kinematic properties of the events with fakes. The total predicted yields of these events in the SR and CRZττ are instead determined by a maximum-likelihood fit to data (see next section), separately for events with 1P and 3P τhad-vis candidates. This approach avoids the uncertainties associated with the simulation of events with fakes, and makes full use of the large amount of data collected.

The remaining background processes (summarized as ‘Others’ in the following) have relatively small contributions in the SR and are estimated using simulations. They include events from the production and decays of top quarks, pairs of gauge bosons, the Higgs boson and W(→ τν) + jets. The yields of these events are normalized to their theoretical cross-sections.

The modelling of the estimated background is validated using events in regions where a possible contamination from signal is negligible. Especially important to the search is the modelling of the combined NN output distribution of Z → ττ events and events with fakes. This is validated by comparing the predicted distributions with data in the CRZττ and in a region similar to the SR, but with events that have same-sign charged ℓ–τhad-vis pairs, as shown in Fig. 1.

Fig. 1: Distributions of the combined NN output in control regions and validation regions.
figure 1

a,b, CRZττ for the μτ channel with 1P (a) and 3P (b) τhad-vis candidates. c,d, Same-sign validation region (VRSS) for the eτ channel for events with 1P (c) and 3P (d) τhad-vis candidates. The expected contributions are determined in a fit to data (‘Constraints on\({\mathcal{B}}(Z\to \ell \tau)\)’). The panels below each plot show the ratios of the observed yields to the best-fit background yields. The hatched error bands represent a one standard deviation of the combined statistical and systematic uncertainties. The statistical uncertainties on the data are shown as vertical bars. The last bin in each plot includes overflow events. Similarly good agreement is observed in the same-sign validation region for the μτ channel and CRZττ for the eτ channel, which are not shown here.

Full size image

Constraints on \({\mathcal{B}}(Z\to \ell \tau )\)

A statistical analysis of the selected events is performed to assess the presence of LFV signal events. The statistical analysis method is detailed in Methods. A simultaneous binned maximum-likelihood fit to the combined NN output in the SR and mcoll(ℓ, τ) in the CRZττ is used to constrain uncertainties in the models and extract evidence of a possible signal. The fit is performed independently for the eτ and μτ channels. Events with 1P and 3P τhad-vis candidates are considered separately. Hypothesis tests, in which a log-likelihood ratio is used as the test statistic, are used to assess the compatibility between the background and signal models and the data.

There are four unconstrained parameters in the fits: two of them determine the overall yields of events with fake 1P τhad-vis or 3P τhad-vis candidates, one determines σZ times the overall acceptance and reconstruction efficiency of the ℓτhad-vis final state in Z → ττ and signal events, and the last one, the parameter of interest, determines the LFV branching fraction \({\mathcal{B}}(Z\to \ell \tau )\) by modifying an arbitrary pre-fit signal yield.

Constrained parameters are also introduced to account for systematic uncertainties in the signal and background predictions. In the case of no significant deviations from the SM background, exclusion limits are set using the CLS method48.

Systematic uncertainties in this search include uncertainties in simulated events in the modelling of trigger, reconstruction, identification and isolation efficiencies, as well as energy calibrations and resolutions of reconstructed objects. Conservative theory uncertainties ranging between 4% and 20% are also assigned to the predicted cross-sections used for the estimation of minor background processes. These uncertainties are not assigned to events with fakes or Z-boson decays, whose yields are determined from data. These events constitute only a small fraction of the background events in the SR. The dominant uncertainties in this search are those in the overall yields of events with fakes, which are predominantly of statistical nature, and those in the τhad-vis energy calibration, which are independent between 1P and 3P τhad-vis candidates and constrained by the fit of the collinear mass spectrum to the data in the CRZττ. A summary of the uncertainties and their impact on the best-fit LFV branching fraction is provided in Table 2, which shows that the sensitivity of the search is primarily limited by the available amount of data.

Table 2 Summary of the uncertainties and their impacts on the measured branching fraction \({\mathcal{B}}(Z\to \ell \tau)\)
Full size table

The best-fit expected and observed distributions of the combined NN output in the SR are shown in Fig. 2. The best-fit yields of Z → ττ and events with fakes are close to the pre-fit predicted values and are determined with a relative precision of 2–4%. Table 3 shows the best-fit expected background and signal yields and the observed number of events in the SR of the eτ and μτ channels with an additional requirement of a combined NN output > 0.7 to consider the most signal-like events.

Fig. 2: Distributions of the combined NN output in the signal region.
figure 2

a,b, eτ events with 1P (a) and 3P (b) τhad-vis candidates. c,d, μτ events with 1P (c) and 3P (d) τhad-vis candidates. The expected contributions are determined in the fit to data. The expected signal, normalized to \({\mathcal{B}}(Z\to \ell \tau )=5\times 1{0}^{-4}\), is shown as a dashed red histogram in each plot. The panels below each plot show the ratios of the observed yields (dots) and the best-fit background-plus-signal yields (solid red line) to the best-fit background yields. The hatched error bands represent a one standard deviation of the combined statistical and systematic uncertainties. The statistical uncertainties on the data are shown as vertical bars. The last bin in each plot includes overflow events.

Full size image
Table 3 Observed number of events and best-fit expected background and signal yields in the SR
Full size table

The best-fit amount of Z → ℓτ signal corresponds to the branching fractions \({\mathcal{B}}(Z\to e\tau )=(-0.1\pm 3.5(\text{stat})\pm 2.3(\text{syst}))\times 1{0}^{-6}\) and \({\mathcal{B}}(Z\to \mu \tau )=(4.3\pm 2.8(\text{stat})\pm 1.6(\text{syst}))\times 1{0}^{-6}\). The positive best-fit value of \({\mathcal{B}}(Z\to \mu \tau )\) is related to a small excess of observed events relative to the background-only hypothesis. This excess has a significance of 0.9 standard deviations when the events with 1P and 3P τhad-vis candidates are fitted simultaneously.

No statistically significant deviation from the SM prediction is observed, and upper limits on the LFV branching fractions are set. For the μτ channel, a more stringent upper limit is set by combining the likelihood function of the presented measurement and a similar measurement done with ATLAS Run 1 data17. Systematic uncertainties from the two measurements are considered uncorrelated in the combined likelihood function. The upper limits are shown in Table 4 for LFV decays with different assumptions about the τ-polarization state. In the scenario where the τ leptons are unpolarized, the observed upper limits at 95% CL on \({\mathcal{B}}(Z\to e\tau )\) and \({\mathcal{B}}(Z\to \mu \tau )\) are 8.1 × 10−6 and 9.5 × 10−6, respectively.

Table 4 Observed and expected (median) upper limits on the signal branching fraction at 95% CL
Full size table

In conclusion, these results from the ATLAS experiment at the LHC set stringent constraints on LFV Z-boson decays involving τ leptons (using only their hadronic decays), superseding the most stringent ones set by the LEP experiments more than two decades ago. The precision of these results is mainly limited by statistical uncertainties.

Methods

Neural network classifiers

Several binary NN classifiers are trained for both the eτ and μτ channels to discriminate signal from the three major backgrounds: W + jets, Z → ττ and Z → ℓℓ. They are referred to as NNWjets, NNZττ and NNZℓℓ, respectively.

The NNs are trained using simulated events selected with the same criteria as those used in the SR, except that the cuts on mvis(ℓ, τ) and the NN output are omitted, and real τhad-vis candidates from Z → ℓτ and Z → ττ are required to pass less stringent identification criteria so as to increase the training sample size. For the Z → ℓℓ process, only events where the τhad-vis candidate is a misidentified light lepton are used. For the W + jets process, jets misidentified as τhad-vis are modelled by simulations. Different NNs are trained separately for eτ and μτ events as well as for events with 1P or 3P τhad-vis candidates. To increase the signal sample size, the Z → eτ and Z → μτ samples are combined and used for training in both channels, assuming equivalent event topology when exchanging e and μ. Owing to the low expected yield of Z → ℓℓ events with 3P τhad-vis candidates, no classifier is trained to discriminate them from background.

A mixture of low-level and high-level kinematic variables are used as input to the NNs. The low-level variables include the four-momenta of the reconstructed ℓ (refs. 18,19), τhad-vis candidate23,24 and \({E}_{\,\text{T}}^{\text{miss}\,}\) (refs. 26,27). To remove known spatial symmetries for optimal training, the low-level variables are transformed in a way that preserves the Lorentz invariance before they are fed into the NNs. The transformation consists of the following steps: first, the ℓ–τhad-vis–\({E}_{\,\text{T}}^{\text{miss}\,}\) system is boosted in a direction in the plane transverse to the beam line such that the total transverse momentum of the system is zero; the system is then rotated about the z axis such that the direction of \({E}_{\,\text{T}}^{\text{miss}\,}\) is aligned with the x axis; if the τhad-vis candidate’s momentum has a negative z component, the entire system is rotated about the new x axis by 180°. After the transformation, only six independent non-vanishing components are left (the τhad-vis candidate is assumed to have zero rest mass), which are the inputs to the NNs.

The high-level variables include Δα, which is a kinematic discriminant defined7 as

$${{\Delta }}\alpha =\frac{{m}_{Z}^{2}-{m}_{\tau }^{2}}{2p(\ell )\times p({\tau }_{{\rm{had}}{\rm{-}}{\rm{vis}}})}-\frac{{p}_{\text{T}}(\ell )}{{p}_{\text{T}}({\tau }_{{\rm{had}}{\rm{-}}{\rm{vis}}})}$$
(2)

where mZ and mτ are the nominal masses of the Z boson and τ lepton, respectively, and p denotes four-momentum. It is specifically defined to test the assumptions that the missing momentum of the event is collinear with the τhad-vis candidate, and that the τ and light leptons in the event are decay products of an on-shell Z boson. For a signal event, where these assumptions are approximately true, it is expected that Δα ≈ 0. Meanwhile, for an SM background event, the value is expected to deviate from zero in general. The other high-level variables are the invariant mass of the ℓ − τhad-vis system, the collinear mass mcoll(ℓ, τ) and the invariant mass of the light lepton and the track associated with the τhad-vis candidate (only used by the Z → ℓℓ classifier).

The training and optimization of the NN classifiers are performed using the open-source software package KERAS49. All of the NNs used in the analysis share the same architecture. Each NN consists of an input layer, two hidden layers of 20 nodes each, and an output layer with a single node. Each layer is fully connected to the neighbouring layers. Low-level and high-level variables are treated in the same way in the input layer. The hidden-layer nodes use rectified linear activation functions, while the output node uses a sigmoid activation function. The NNs are trained using the Adam algorithm50 to optimize the binary cross entropy. All the NNs are trained with a batch size of 256 and 200 epochs. The number of hidden layers, the number of nodes per layer, the training batch size and the learning rate parameter of the optimizer are simultaneously chosen by maximizing the area under the expected receiver operating characteristic curve. The optimization is done with a grid scan. No regularization or dropout is added, and no sign of overtraining is observed. For other configurations and hyperparameters that have not been mentioned, the default settings in KERAS 1.1.0 are used.

Each NN classifier outputs a score between 0 and 1 for each event, where a higher score indicates that the event is more signal-like. The output scores from the different classifiers are combined into the final discriminant (combined NN output) using the formula

$$\text{Combined NN output}=1-\sqrt{\frac{{\sum }_{b}{w}_{b}\times {\left(1-{\text{NN}}_{b}\,\text{output}\right)}^{2}}{{\sum }_{b}{w}_{b}}}$$
(3)

where b = Wjets, Zττ, Zℓℓ and wb are constant parameters. Output scores for events with 1P τhad-vis candidates and those with 3P τhad-vis candidates are combined separately. The summation is over Wjets, Zττ and Zℓℓ for events with 1P τhad-vis candidates, and only over Wjets and Zττ for events with 3P τhad-vis candidates.

By construction, the combined NN output ranges between 0 and 1, where 0 represents the most background-like (and 1 the most signal-like) event possible. The choice of values of wb affects the expected sensitivity of the analysis because they change how events from the different background processes are distributed along the range of combined NN output values, and thus impact the ability of the binned maximum-likelihood fit to determine the background contributions. The values of wb are chosen with a grid scan to minimize the expected upper limit on the branching fraction in the absence of a signal. The chosen values have the ratio wZττ:wWjets:wZℓℓ = 1.0:1.5:0.33. As could be expected, the optimized weights loosely reflect the impact of the uncertainties in the corresponding backgrounds on the determination of the signal branching fraction.

Maximum-likelihood fit

Binned maximum-likelihood fits are implemented using the statistical analysis packages ROOFIT51, ROOSTATS52 and HISTFITTER53. The expected binned distributions of the combined NN output in the SR and the collinear mass in the CRZττ are fit to data to extract evidence of signal events. Fitting the data in the CRZττ and in part of the SR with low combined NN output values (where no signal is expected) benefits the overall sensitivity to the signal, because it reduces the uncertainties of the background model in the high combined NN output value region, where most of the signal is expected. Owing to the differences in background composition, acceptance and efficiencies, regions with 1P and 3P τhad-vis candidates are fit separately but simultaneously. The probabilities of compatibility between the data and the background-only or background-plus-signal hypotheses are assessed using the modified frequentist CLS method48, and exclusion upper limits on \({\mathcal{B}}(Z\to \ell \tau )\) are set by the inversion of these hypothesis tests.

The background-plus-signal model has four unconstrained parameters before the fit. Two of the parameters determine the overall yields of events with 1P and 3P fakes separately. A third parameter determines σZ times the overall acceptance and reconstruction efficiency of events with a true ℓτhad-vis final state. It is applied to the normalizations of both the signal and Z → ττ events to ensure that the same σZ times acceptance is estimated for both processes. The last unconstrained parameter is the parameter of interest μsig, which controls the normalization of signal events. Given the similarity between the signal and Z → ττ → ℓτhad-vis + 3ν final states and that both processes are estimated with the same σZ and acceptance and efficiency corrections, this choice of parameterization reduces the impact on the determined \({\mathcal{B}}(Z\to \ell \tau )\) from detector effects and uncertainties in predicting σZ. The parameter of interest represents

$${\mu }_{\text{sig}} =\frac{{\mathcal{B}}(Z\to \ell \tau )}{{{\mathcal{B}}}_{\text{pre-fit}}(Z\to \ell \tau )}$$
(4)

where \({{\mathcal{B}}}_{\text{pre-fit}}(Z\to \ell \tau )\) is an arbitrary branching fraction to which the signal prediction is normalized. Although the physical branching fraction must be positive, the parameter of interest in the fit is not constrained to be positive.

Systematic uncertainties are represented by nuisance parameters (NPs) with Gaussian constraints in the likelihood function. The impact of uncertainties on both the shape and normalization of the predicted distributions are taken into account. Uncertainties in the energy calibration and resolution as well as in the trigger, reconstruction, identification and isolation efficiencies of jets, electrons, muons, τhad-vis and \({E}_{\,\text{T}}^{\text{miss}\,}\) are considered. Theoretical uncertainties in the production cross-sections affect only the predictions of the minor backgrounds, because the Z → ττ and signal yields are determined in the maximum-likelihood fit to data and the Z → ℓℓ yield is determined by the measured value of σZ. Statistical uncertainties in the determination of the fake factors are also considered. They are modelled by one NP per pT bin in which the fake factors are measured. As noted in the section ‘Constraints on\({\mathcal{B}}(Z\to \ell \tau )\)’ the dominant uncertainties in the analysis are the statistical uncertainties in determining how many events have fakes and the systematic uncertainties in the reconstructed τhad-vis energy.

For the μτ channel, the likelihood functions of the presented measurement and of the measurement in ref. 17 are combined. As the two measurements are statistically uncorrelated and the predictions are based on different methods, NPs in the individual likelihood functions are considered uncorrelated in the combination. The method of combination is the same as in ref. 13.