Search for a new $Z^{\prime}$ gauge boson in $4\mu$ events with the ATLAS
experiment This paper presents a search for a new $Z^{\prime}$ vector gauge
boson with the ATLAS experiment at the Large Hadron Collider using $pp$
collision data collected at $\sqrt{s}=13$ TeV, corresponding to an integrated
luminosity of 139 fb-1. The new gauge boson $Z^{\prime}$ is predicted by
$L_{\mu}-L_{\tau}$ models to address observed phenomena that can not be
explained by the Standard Model. The search examines the four-muon (4$\mu$)
final state, using a deep learning neural network classifier to separate the
$Z^{\prime}$ signal from the Standard Model background events. The di-muon
invariant masses in the $4\mu$ events are used to extract the $Z^{\prime}$
resonance signature. No significant excess of events is observed over the
predicted background. Upper limits at a 95% confidence level on the
$Z^{\prime}$ production cross-section times the decay branching fraction of
$pp\rightarrow Z^{\prime}\mu\mu\rightarrow 4\mu$ are set from 0.31 to 4.3 fb
for the $Z^{\prime}$ mass ranging from 5 to 81 GeV. The corresponding common
coupling strengths, $g_{Z^{\prime}}$, of the $Z^{\prime}$ boson to the second
and third generation leptons above 0.003 – 0.2 have been excluded.
HDBS-2018-57 CERN-EP-2022-251 JHEP

##  1 Introduction

This paper presents a search for a new vector boson $Z^{\prime}$ in the four-
muon (4$\mu$) final state with data recorded in proton-proton ($pp$)
collisions at $\sqrt{s}$ = 13 TeV by the ATLAS detector [1] at the Large
Hadron Collider (LHC), corresponding to an integrated luminosity of 139 fb-1.
The new gauge boson $Z^{\prime}$, which only interacts with the second and
third generation leptons, is predicted by $L_{\mu}-L_{\tau}$ models [2] which
extend the Standard Model (SM) with an additional $U(1)_{L_{\mu}-L_{\tau}}$
symmetry. In recent years searching for such a new gauge boson $Z^{\prime}$
has attracted attention in both theoretical and experimental communities [3,
4, 5] because it could address the observed muon magnetic dipole moment
($g$-$2$) anomaly [6, 7, 8, 9] through the $Z^{\prime}$ loop corrections
without contradicting other existing data. The new $Z^{\prime}$ is also
ideally suited to address the lepton flavour anomalies measured with the ratio
of B-meson decays to muons and electrons, $R_{K}=\mathcal{B}(B^{+}\rightarrow
K^{+}\mu^{+}\mu^{-})/\mathcal{B}(B^{+}\rightarrow K^{+}e^{+}e^{-})$ [10, 11,
12, 13, 14]. Such anomalies suggest possible lepton flavor universality
violation. The $L_{\mu}-L_{\tau}$ models could provide an explanation for the
observed results, since by construction, the $Z^{\prime}$ arising from the
models only couples to muons and taus, and does not couple to electrons. These
models also aim to probe outstanding questions related to dark matter and
neutrino mass [15, 16, 17].

The $Z^{\prime}$ kinematics, mass, and interactions (with the same coupling
strength to the second and third lepton families), are described by the
Lagrangian below:

$L_{Z^{\prime}}=-{1\over 4}F_{\alpha\beta}F^{\alpha\beta}+{1\over
2}M^{2}_{Z^{\prime}}Z^{\prime\alpha}Z^{\prime}_{\alpha}-g_{Z^{\prime}}Z^{\prime}_{\alpha}(\bar{\ell}_{2}\gamma^{\alpha}\ell_{2}+\bar{\mu}\gamma^{\alpha}\mu-\bar{\ell}_{3}\gamma^{\alpha}\ell_{3}-\bar{\tau}\gamma^{\alpha}\tau),$

where
$F_{\alpha\beta}=\partial_{\alpha}Z^{\prime}_{\beta}-\partial_{\beta}Z^{\prime}_{\alpha}$
is the $Z^{\prime}$ field strength tensor; $\ell_{i}=(\nu_{i},e_{i})^{T}$
($i=2,3$, denoting the second and the third generation left-handed lepton
doublets); $\mu$ and $\tau$ represent the right-handed muon and tau singlets;
and $g_{Z^{\prime}}$ (which will be referred to as $g$ in the rest of this
paper) is the interaction coupling constant. The parameter space of
($m_{Z^{\prime}},g$) has not been strongly constrained in experiments since
the $Z^{\prime}$ does not directly couple to the electron nor to any quarks,
hence it could not be directly produced from $e^{+}e^{-}$ and $pp$ colliding
beams.

In $pp$ collisions at the LHC, the $Z^{\prime}$ could be produced from final
state radiation of $\mu$ or $\tau$ pairs of the Drell-Yan (DY) process as
shown in Figure 1(a) with a $4\mu$ final state, which provides the cleanest
signature to search for the $Z^{\prime}$. For relatively low $Z^{\prime}$
mass, the most promising experimental signature would be an excess of 4$\mu$
events with the invariant mass of one $\mu^{+}\mu^{-}$ pair peaking around the
$Z^{\prime}$ mass. The major background comes from the SM 4$\mu$ production
processes shown in Figure 1(b), 1(c) and 1(d). The $Z^{\prime}$ could also be
produced in $W$ production through the DY process, $pp\rightarrow W\rightarrow
Z^{\prime}\mu\nu\rightarrow 3\mu+\nu$. The experimental signature would have a
final state of 3$\mu$ plus large missing transverse energy. This final state
is not included in this analysis.

((a))

((b))

((c))

((d))

Figure 1: Feynman diagrams of $Z^{\prime}$ production through radiation in a
Drell-Yan process (a), and of the corresponding SM background processes (b -
d) with a 4$\mu$ final state.

Both the ATLAS and CMS Collaborations have measured the cross-sections of the
SM $Z\rightarrow 4\mu$ process [18, 19, 20]. The measurement by ATLAS was used
by theorists to set limits in the parameter space of the $L_{\mu}-L_{\tau}$
model [15]. The CMS Collaboration has directly searched for the $Z^{\prime}$
boson in the mass region between 5 to 70 GeV with the 4$\mu$ final state using
77.3 fb-1 of data [21], and set upper limits on the $Z^{\prime}$ to muon
coupling strength, $g$, of 0.004 – 0.3 at 95% confidence level, depending on
the $Z^{\prime}$ mass.

The organization of this paper is as follows. The ATLAS detector is described
in the next section. The dataset and Monte Carlo samples used in this analysis
are detailed in Section 3. Event reconstruction and selection, followed by
background estimation from data, are described in Section 4 and Section 5.
Event classification using a deep learning approach is presented in Section 6.
Systematic uncertainties and the statistical approach to interpret data to
obtain the results are reported in Section 7 and Section 8. The conclusion is
given in the final section.

##  2 ATLAS Detector

The ATLAS detector [1] at the LHC covers nearly the entire solid angle around
the collision point.111 ATLAS uses a right-handed coordinate system with its
origin at the nominal interaction point (IP) in the centre of the detector and
the $z$-axis along the beam pipe. The $x$-axis points from the IP to the
centre of the LHC ring, and the $y$-axis points upwards. Cylindrical
coordinates $(r,\phi)$ are used in the transverse plane, $\phi$ being the
azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of
the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. Angular distance is
measured in units of $\Delta R\equiv\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$.
It consists of an inner tracking detector surrounded by a thin superconducting
solenoid, electromagnetic and hadron calorimeters, and a muon spectrometer
incorporating three large superconducting air-core toroidal magnets.

The inner-detector system (ID) is immersed in a $2\text{\,}\mathrm{T}$ axial
magnetic field and provides charged-particle tracking in the range
$|\eta|<2.5$. The high-granularity silicon pixel detector covers the vertex
region and typically provides four measurements per track, the first hit
normally being in the insertable B-layer (IBL) installed before Run 2 [22,
23]. It is followed by the silicon microstrip tracker (SCT), which usually
provides eight measurements per track. These silicon detectors are
complemented by the transition radiation tracker (TRT), which enables radially
extended track reconstruction up to $|\eta|=2.0$. The TRT also provides
electron identification information based on the fraction of hits (typically
30 in total) above a higher energy-deposit threshold corresponding to
transition radiation.

The calorimeter system covers the pseudorapidity range $|\eta|<4.9$. Within
the region $|\eta|<3.2$, electromagnetic calorimetry is provided by barrel and
endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an
additional thin LAr presampler covering $|\eta|<1.8$ to correct for energy
loss in material upstream of the calorimeters. Hadron calorimetry is provided
by the steel/scintillator-tile calorimeter, segmented into three barrel
structures within $|\eta|<1.7$, and two copper/LAr hadron endcap calorimeters.
The solid angle coverage is completed with forward copper/LAr and tungsten/LAr
calorimeter modules optimised for electromagnetic and hadronic energy
measurements respectively.

The muon spectrometer (MS) comprises separate trigger and high-precision
tracking chambers measuring the deflection of muons in a magnetic field
generated by the superconducting air-core toroidal magnets. The field integral
of the toroids ranges between $2.0$ and
$6.0\text{\,}\mathrm{T}\text{\,}\mathrm{m}$ across most of the detector. Three
layers of precision chambers, each consisting of layers of monitored drift
tubes, covers the region $|\eta|<2.7$, complemented by cathode-strip chambers
in the forward region, where the background is highest. The muon trigger
system covers the range $|\eta|<2.4$ with resistive-plate chambers in the
barrel, and thin-gap chambers in the endcap regions.

Interesting events are selected by the first-level trigger system implemented
in custom hardware, followed by selections made by algorithms implemented in
software in the high-level trigger [24]. The first-level trigger accepts
events from the $40\text{\,}\mathrm{MHz}$ bunch crossings at a rate below
$100\text{\,}\mathrm{kHz}$, which the high-level trigger further reduces in
order to record events to disk at about $1\text{\,}\mathrm{kHz}$.

An extensive software suite [25] is used in data simulation, in the
reconstruction and analysis of real and simulated data, in detector
operations, and in the trigger and data acquisition systems of the experiment.

##  3 Dataset and Monte Carlo Simulations

The data used for this analysis were recorded using single-muon and multi-muon
triggers, corresponding to an integrated luminosity of 139 fb-1 after the
application of data quality requirements [26]. The transverse momentum
($p_{\text{T}}$) thresholds for the single-muon trigger vary from 20 to 26
GeV, for the di-muon trigger from 10 to 14 GeV, and for the tri-muon trigger
from 4 to 6 GeV, depending on the data-taking periods [27]. The overall
trigger efficiency for the $4\mu$ events selected in this analysis is higher
than 98%.

###  3.1 Simulation of $Z^{\prime}$ Production

The $L_{\mu}-L_{\tau}$ model is used for Monte Carlo (MC) $Z^{\prime}$ signal
sample production, where the $Z^{\prime}$ couples to the left-handed (LH) muon
or tau leptons and their corresponding neutrinos, and to the right-handed (RH)
muon or tau leptons. In the model, the $Z^{\prime}$ couplings to the first
lepton families (electron and its neutrino) and all quarks are set to zero.
The couplings to second and third generation leptons are assumed to be
identical, therefore the branching fractions of the $Z^{\prime}$ decay to a
pair of muons and a pair of muon neutrinos are set to $\frac{1}{3}$ and
$\frac{1}{6}$, respectively. The signal from tau decays in the 4$\mu$ final
state is found to be negligible in this analysis and is not included as
signal.

The signal events are generated with MadGraph5_aMC@NLO 2.7.3 [28] at leading-
order (LO) accuracy in QCD by using the Universal FeynRules Output (UFO)
format [29, 30]. The interactions [31] mediated by a resonance $Z^{\prime}$
which couples to the second and third generation leptons are used for four-
muon signal event generation. Based on theoretical calculations [32, 33] for
related processes, the appropriate NNLO/LO correction factor K of 1.3 is used
to correct the MC LO signal cross-sections. In the signal production
simulations, the NNPDF2.3nlo set [34] is used as parton distribution function
(PDF) for the $pp$ collisions. Photos++ 3.61 [35] is used as the photon
emissions from electroweak vertices and charged leptons via a QED process.

The MC simulated events are generated for a range of masses and coupling
parameters of the $L_{\mu}-L_{\tau}$ model. The $Z^{\prime}$ mass ranges from
5 to 81 GeV, and the value of coupling constant $g$ ranges from 0.008 to
0.316, as summarized in Table 1. The value of the $g$ at each $Z^{\prime}$
mass was chosen to allow a sensitive search for a very small $Z^{\prime}$
signal in the full Run 2 dataset.

With the chosen $g$, the natural width $\Gamma$ of the $Z^{\prime}$ (smaller
than the experimental resolution) and the cross-section, $\sigma$
($pp\rightarrow
Z^{\prime}\mu^{+}\mu^{-}\rightarrow\mu^{+}\mu^{-}\mu^{+}\mu^{-}$), are
calculated as listed in Table 1.

Table 1: Summary of the chosen $Z^{\prime}$ hypotheses and corresponding
coupling, width, and cross-section (calculated at LO accuracy in QCD) at each
mass point.

$m_{Z^{\prime}}$ [GeV] | $g$ |  $\Gamma$ [GeV] |  $\sigma$ [fb] |  $m_{Z^{\prime}}$ [GeV] | $g$ |  $\Gamma$ [GeV] |  $\sigma$ [fb]
---|---|---|---|---|---|---|---
5 | 0.0080 |  $2.45\text{\times}{10}^{-5}\text{\,}$ $\text{\,}$ | 9.96 | 42 | 0.0900 | $2.71\text{\times}{10}^{-2}\text{\,}$ | 13.38
7 | 0.0085 | $3.99\text{\times}{10}^{-5}\text{\,}$ | 7.06 | 45 | 0.1000 | $3.58\text{\times}{10}^{-2}\text{\,}$ | 11.72
9 | 0.0090 | $5.78\text{\times}{10}^{-5}\text{\,}$ | 5.60 | 48 | 0.1100 | $4.62\text{\times}{10}^{-2}\text{\,}$ | 9.96
11 | 0.0095 | $7.89\text{\times}{10}^{-5}\text{\,}$ | 4.65 | 51 | 0.1200 | $5.84\text{\times}{10}^{-2}\text{\,}$ | 8.24
13 | 0.0100 | $1.03\text{\times}{10}^{-4}\text{\,}$ | 3.95 | 54 | 0.1600 | $1.10\text{\times}{10}^{-1}\text{\,}$ | 10.07
15 | 0.0120 | $1.72\text{\times}{10}^{-4}\text{\,}$ | 4.45 | 57 | 0.2000 | $1.81\text{\times}{10}^{-1}\text{\,}$ | 10.73
17 | 0.0140 | $2.65\text{\times}{10}^{-4}\text{\,}$ | 4.80 | 60 | 0.2665 | $3.39\text{\times}{10}^{-1}\text{\,}$ | 12.92
19 | 0.0160 | $3.87\text{\times}{10}^{-4}\text{\,}$ | 5.00 | 63 | 0.2680 | $3.60\text{\times}{10}^{-1}\text{\,}$ | 8.84
23 | 0.0240 | $1.05\text{\times}{10}^{-4}\text{\,}$ | 7.30 | 66 | 0.2780 | $4.06\text{\times}{10}^{-1}\text{\,}$ | 6.50
27 | 0.0320 | $2.20\text{\times}{10}^{-3}\text{\,}$ | 8.50 | 69 | 0.2890 | $4.59\text{\times}{10}^{-1}\text{\,}$ | 4.89
31 | 0.0400 | $3.95\text{\times}{10}^{-3}\text{\,}$ | 8.72 | 72 | 0.3000 | $5.15\text{\times}{10}^{-1}\text{\,}$ | 3.80
35 | 0.0600 | $1.00\text{\times}{10}^{-2}\text{\,}$ | 12.82 | 75 | 0.3000 | $5.37\text{\times}{10}^{-1}\text{\,}$ | 2.88
39 | 0.0800 | $1.99\text{\times}{10}^{-2}\text{\,}$ | 14.77 | 78 | 0.3080 | $5.89\text{\times}{10}^{-1}\text{\,}$ | 2.40
|  |  |  | 81 | 0.3160 | $6.44\text{\times}{10}^{-1}\text{\,}$ | 2.08

The $Z^{\prime}$ signal samples were simulated with the ATLAS fast simulation
framework (Atlfast-II) [36] to produce predictions that can be directly
compared with the data.

###  3.2 Simulation of Background Events

Dominant SM backgrounds in this analysis come from the SM $Z\rightarrow 4\mu$
processes where the four leptons have an invariant mass close to that of the
$Z$ boson. In the higher mass region the $ZZ^{*}$ production contributes a
sizable number of prompt $4\mu$ events. In addition, there are very small
contributions from the Higgs boson, $t\bar{t}V$ ($V=W,Z$), and tri-boson
($VVV$) production processes. These events are estimated with MC simulations.
The non-prompt muon background events, mostly coming from $Z$ \+ jets,
$t\bar{t}$ and single-top-quark production processes, are estimated from data
in this analysis, as described in Section 5. The $Z$ \+ jets and $t\bar{t}$ MC
samples are also produced for background studies.

Samples of diboson final states ($q\bar{q}\rightarrow VV^{(*)}$), including
the processes shown in Figures 1(b) and 1(c), were simulated with the Sherpa
2.2.2 [37] generator, including off-shell effects and Higgs boson
contributions, where appropriate. Fully leptonic final states and semileptonic
final states, where one boson decays leptonically and the other hadronically,
were generated using matrix elements at next-to-leading-order (NLO) accuracy
in QCD for up to one additional parton and at LO accuracy for up to three
additional parton emissions. Samples for the loop-induced processes $gg\to
VV^{(*)}$, shown in Figure 1(d), were generated using LO-accurate matrix
elements for up to one additional parton emission for both the cases of fully
leptonic and semileptonic final states. The matrix element calculations were
matched and merged with the Sherpa parton shower based on Catani–Seymour
dipole factorisation [38, 39] using the MEPS@NLO prescription [40, 41, 42,
43]. The virtual QCD corrections were provided by the OpenLoops library [44,
45, 46]. The NNPDF3.0nnlo set of PDFs was used [47], along with the dedicated
set of parton-shower parameters (tune) developed by the Sherpa authors.

The production of $t\bar{t}$ events was modelled using the Powheg Box v2 [48,
49, 50, 51] generator at NLO in QCD with the NNPDF3.0nlo PDF set and the
$h_{\mathrm{damp}}$ parameter222The $h_{\mathrm{damp}}$ parameter is a
resummation damping factor and one of the parameters that controls the
matching of Powheg matrix elements to the parton shower and thus effectively
regulates the high-$p_{\text{T}}$ radiation against which the $t\bar{t}$
system recoils. set to 1.5 $m_{\text{top}}$ [52]. The events were interfaced
to Pythia 8.230 [53] to model the parton shower, hadronisation, and underlying
event, with parameters set according to the A14 tune [54] and using the
NNPDF2.3lo set of PDFs. The associated production of top quarks with $W$
bosons ($tW$) was modelled by the Powheg Box v2 [55, 49, 50, 51] generator at
NLO in QCD using the five-flavour scheme and the NNPDF3.0nlo set of PDFs. The
diagram removal scheme [56] was used to remove interference and overlap with
$t\bar{t}$ production. The events were interfaced to Pythia 8.230 using the
A14 tune and the NNPDF2.3lo set of PDFs.

The production of $V+$jets was simulated with the Sherpa 2.2.1 generator using
NLO matrix elements for up to two partons, and LO matrix elements for up to
four partons, calculated with the Comix and OpenLoops libraries. They were
matched with the Sherpa parton shower using the MEPS@NLO prescription with the
set of tuned parameters developed by the Sherpa authors. The NNPDF3.0nnlo set
of PDFs was used and the samples were normalised to a next-to-next-to-leading-
order (NNLO) prediction [57].

The production of $t\bar{t}V$ events was modelled using the MadGraph5_aMC@NLO
2.3.3 generator at NLO with the NNPDF3.0nlo PDF. The events were interfaced to
Pythia 8.210 using the A14 tune and the NNPDF2.3lo PDF set. The production of
tri-boson ($VVV$) events was simulated with the Sherpa 2.2.2 generator. Matrix
elements accurate to LO in QCD for up to one additional parton emission were
matched and merged with the Sherpa parton shower based on Catani–Seymour
dipole factorisation using the MEPS@NLO prescription. Samples were generated
using the NNPDF3.0nnlo PDF set, along with the dedicated set of tuned parton-
shower parameters developed by the Sherpa authors.

The generated background MC samples were produced through the full ATLAS
detector simulation based on Geant4 [58]. The effect of multiple interactions
in the same and neighbouring bunch crossings (pile-up) was modelled by
overlaying the simulated hard-scattering event with inelastic $pp$ events
generated with Pythia 8.186 [59] using the NNPDF2.3lo set of PDF and the A3
set of tuned parameters [60]. Simulated events were reweighted to match the
pile-up conditions in the data. All simulated events were processed using the
same reconstruction algorithms and triggering requirements as used in data.

##  4 Event Reconstruction and Pre-selection

Proton-proton interaction vertices are reconstructed in events with at least
two tracks, each with $p_{\text{T}}>0.5\leavevmode\nobreak\
\text{Ge\kern-1.00006ptV}$. The primary hard-scatter vertex is defined as the
one with the largest value of the sum of squared track transverse momenta.

Muons are identified by matching tracks reconstructed in the MS to tracks
reconstructed in the ID (referred to as combined muons). To increase the muon
reconstruction efficiency non-combined muon identification algorithms are also
used in the analysis, including using the MS stand-alone tracks in the region
$2.5<|\eta|<2.7$, and matching the ID tracks with calorimeter hit information
within $|\eta|<$0.1, as well as using the ID tracks associated with at least
one local track segment in the MS. In the $4\mu$ event selection at most one
of the selected muons can be a non-combined muon. Each muon is then required
to satisfy the ‘loose’ identification criteria [61]. Muons are required to be
isolated using a particle-flow algorithm [62] and associated with the primary
hard-scatter vertex by satisfying $|{d_{0}\over\sigma_{d_{0}}}|<$ 3 and
$|z_{0}\times\sin\theta|<$ 0.5 mm, where $d_{0}$ is the transverse impact
parameter calculated with respect to the measured beam-line position,
$\sigma_{d_{0}}$ its uncertainty, and $z_{0}$ is the longitudinal distance
between the point at which $d_{0}$ is measured and the primary vertex. The
minimum muon $p_{\text{T}}$ threshold is 3 GeV.

In addition to muons, electrons, jets and missing transverse momentum
($E_{\text{T}}^{\text{miss}}$) are also used to select control samples for
background estimation in this analysis. The reconstructions of these objects
are described below.

Electrons are identified with a likelihood discriminator built from the shower
shapes of electron energy loss in the calorimeter, track-cluster matching, and
some of the track quality distributions. Each electron is required to satisfy
the ‘medium’ likelihood identification criteria [63], as well as similar
vertex and isolation requirements as muons. The reconstructed electrons are
required to have $p_{\textup{T}}>$ 7 GeV and $|\eta|<$ 2.5, excluding the
transition region between the barrel and endcap calorimeters,
$1.37<|\eta|<1.52$.

Jets are reconstructed with the anti-$k_{t}$ algorithm [64, 65] with a radius
parameter of $R=0.4$. The jet clustering input objects are based on particle-
flow [66] in the ID and the calorimeter. Jets are required to have
$p_{\textup{T}}>20\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ and
$|\eta|<2.5$. Jets containing B hadrons, referred to as $b\text{-jets}$, are
identified with a multivariate discriminant [67]. To reduce the effect of
pile-up an additional quality requirement based on the ‘Jet-Vertex-Tagger’
algorithm (JVT) [68] is applied in jet identification.

$E_{\text{T}}^{\text{miss}}$ is determined as the magnitude of the negative
vectorial sum of the transverse momenta of the selected and calibrated physics
objects (including muons, electrons, photons, and jets including hadronically
decaying tau-leptons) and the ID tracks must come from the main vertex and not
associated with any physics object (soft term) [69].

Events containing at least four muons with kinematics consistent with
$Z(Z^{*}/\gamma^{*})\rightarrow 4\mu$ production are then selected as follows.
The four leading $p_{\text{T}}$-ordered muons are required to pass the
$p_{\text{T}}$ thresholds of 20, 15, 8, and 3 GeV, respectively. If a muon is
selected as a non-combined muon, its $p_{\text{T}}$ must be greater than 15
GeV. Any di-muon pair in the event must have an invariant mass $m_{\mu\mu}$
greater than 4 GeV and an angular separation $\Delta R$ larger than 0.2. To
search for the $Z^{\prime}\rightarrow\mu^{+}\mu^{-}$ signature, two muon pairs
are selected based on their invariant mass values. The first pair (referred to
as $Z_{1}$) is selected from all the possible $\mu^{+}\mu^{-}$ pairs to have
the smallest mass difference between the $Z_{1}$ mass and the $Z$ mass,
$|m_{Z}-m_{Z_{1}}|$. The second $\mu^{+}\mu^{-}$ pair is selected from the
remaining muons that has the highest invariant mass (referred to as $Z_{2}$).
The correct signal di-muon pairing fraction varies with the $Z^{\prime}$ mass,
where the selected di-muon pair that forms $m_{Z1}$ or $m_{Z2}$ originates
from the $Z^{\prime}$. For example, for $m_{Z^{\prime}}$ = 5, 42, 63, 72, and
81 GeV, the correct di-muon pairing fractions are about 78%, 50%, 88%, 82%,
and 90%, respectively. Finally, the selected four muons must have an invariant
mass in a range of 80 to 180 GeV, but excluding the Higgs boson resonance mass
region of 110 to 130 GeV.

The $Z^{\prime}$ signal efficiency at various stages of the event selection is
shown in Table 2 for five representative mass points. The event selection
efficiencies vary significantly depending on the $Z^{\prime}$ mass. At
generator level, an MC filter is applied, which requires at least four muons
with $p_{T}>2$ GeV and $|\eta|<$ 3.0. The MC filter efficiencies of these
representative $Z^{\prime}$ signal samples are listed in the table as well.

Table 2:  The $Z^{\prime}$ signal event selection efficiencies compared to the events passing the previous cut level for several representative mass points. The overall signal efficiencies are the products of the 4$\mu$ MC filter and the combined event selection efficiencies.  $m_{Z^{\prime}}$ [GeV] | 5 | 42 | 63 | 72 | 81
---|---|---|---|---|---
MC filter efficiency | 32.8% | 57.7% | 61.0% | 65.3% | 70.0%
Number of identified muons $\geq$ 4 | 47.3% | 74.1% | 70.8% | 72.4% | 75.4%
$p_{\textup{T}}^{i}(i=1,4)>20,\leavevmode\nobreak\ 15,\leavevmode\nobreak\ 8,\leavevmode\nobreak\ 3$ GeV | 60.0% | 82.6%, | 90.3% | 93.6% | 98.2%
$\Delta R(\mu_{i},\mu_{j})>0.2$ & vertex requirement | 87.2% | 95.4% | 96.2% | 96.6% | 97.2%
Isolation | 54.2% | 76.9% | 79.2% | 84.1% | 87.5%
$m_{4\mu}$ within [80, 110] or [130, 180] GeV | 91.9% | 88.8% | 58.9% | 33.5% | 16.8%
Combined event selection efficiency | 12.3% | 39.9% | 28.7% | 18.4% | 10.6%
Overall 4$\mu$ signal efficiency | 4.1% | 23.0% | 17.5% | 11.9% | 7.4%

The $Z^{\prime}$ production signature is searched for in the $Z_{1}$ or
$Z_{2}$ mass spectrum depending on the $Z^{\prime}$ mass. The relatively high-
mass $Z^{\prime}$ signals mostly appear as a peak in the $Z_{1}$ spectrum
while the relatively low-mass signals mostly appear as a peak in the $Z_{2}$
spectrum. Representative examples of the predicted signal over background,
after further selection with a deep learning approach which will be described
in Section 6, are shown in Figure 4. In the analysis the $Z_{1}$ and $Z_{2}$
mass spectra are scanned to search for a $Z^{\prime}$ with mass greater or
smaller than 42 GeV, respectively. The boundary value of 42 GeV is chosen
based on the studies to optimize the search sensitivity.

The numbers of 4$\mu$ events in data and the estimated background yields are
given in Table 3. More details about the estimation of the reducible
backgrounds containing non-prompt muons can be found in Section 5. The total
uncertainties of simulated backgrounds are also listed in the table. The
evaluations of systematic uncertainties will be described in Section 7.

Table 3: The selected 4$\mu$ events in data and the estimated backgrounds and their combined statistical and systematic uncertainties. Data | Total | $qq\rightarrow ZZ^{*}$ | $gg\rightarrow ZZ^{*}$ | $ttV+VVV+H$ | Reducible background
---|---|---|---|---|---
| background | from MC | from MC | from MC | from data
1131 | $1148\pm 70$ | $1065^{+70}_{-69}$ | $15.6\pm 2.5$ | $6.2\pm 2.9$ | $61.1^{+8.3}_{-9.1}$

##  5 Estimation of Reducible Background from Data

A data-driven technique is used to estimate the reducible background from $Z$
\+ jets and $t\bar{t}$ production. Events from these processes may contain two
prompt leptons from $W$ or $Z$ boson decays, together with additional activity
such as heavy-flavor jets or misidentified components of jets yielding
reconstructed muons, collectively referred to as non-prompt muons. Such
backgrounds are estimated from data using a control sample of
$\mu^{+}\mu^{-}\mu_{j}\mu_{j}$ events, selected with the standard signal
requirements except that non-prompt muons, $\mu_{j}$, are selected in place of
two of the signal muons. Non-prompt muons are defined as muon candidates using
the standard selection but fail the requirements on isolation and impact
parameter. However, they are required to pass a much looser impact parameter
significance requirement, $|{d_{0}\over\sigma_{d_{0}}}|<10$. The background in
the signal sample is estimated by scaling each event in the selected
$\mu^{+}\mu^{-}\mu_{j}\mu_{j}$ control sample by $f_{1}\times f_{2}$, where
the $f_{i}(i=1,2)$ is referred to as a "fake-factor" for each of the two non-
prompt muons and is computed as a function of their $p_{\textup{T}}$ and
$|\eta|$.

The fake-factor $f$ is derived from two independent non-prompt-muon enriched
data control samples dominated by $Z$ \+ jets or $t\bar{t}$ events. The $Z$ \+
jets ($t\bar{t}$) sample contains non-prompt muons dominantly from light jets
(heavy-flavor jets). The $Z$ \+ jets sample is tagged by two isolated prompt
leptons from $Z$ decays ($Z\rightarrow\mu^{+}\mu^{-},\leavevmode\nobreak\
e^{+}e^{-}$). The $p_{\text{T}}$-threshold for the two prompt leptons are 27
and 25 GeV and their invariant mass must be within 10 GeV of the $Z$ mass.

The $E_{\text{T}}^{\text{miss}}$ of the event must be less than 25 GeV. The
$t\bar{t}$ sample is tagged by two prompt isolated high- $p_{\text{T}}$
$e^{\pm}\mu^{\mp}$ leptons with $p_{\text{T}}$ thresholds of 28 and 25 GeV and
associated with at least one $b$-jet. The event is required to have
$E_{\text{T}}^{\text{miss}}>$ 50 GeV. If there is an additional isolated
high-$p_{\text{T}}$ lepton in the event, the transverse mass calculated with
the $E_{\text{T}}^{\text{miss}}$ and the third lepton is required to be less
than 60 GeV to reduce the contribution from the $t\bar{t}W$ events. Events
with four muons that satisfy the signal selection criteria are removed from
both $Z$ \+ jets and $t\bar{t}$ control samples. Additional reconstructed muon
objects, passing either the signal-muon or the non-prompt-muon selection
criteria are counted in the control samples. The fake-factor $f$ is calculated
as the ratio of the probability for a non-prompt muon to satisfy the signal
muon selection criteria to the probability for a non-prompt muon to fail the
signal muon selection. The factor $f$ is derived as a function of $|\eta|$ and
$p_{\text{T}}$ of the non-prompt muons, which varies from 0.11 to 0.76 (0.01
to 0.27) in the $Z$ \+ jets ($t\bar{t}$) control sample. Then, the fake-
factors derived from two data control samples are combined using the $m_{Z1}$
spectra obtained from simulated $4\mu$ background events from $Z$ \+ jets and
$t\bar{t}$ processes. A simultaneous fit of the $m_{Z1}$ spectra of the MC
events to data selected in the $\mu^{+}\mu^{-}\mu_{j}\mu_{j}$ control sample
is performed to determine the fractions of the $Z$ \+ jets and $t\bar{t}$
events in each bin of the mass spectrum.

The overall systematic uncertainty of $f$ is about 14%. It is determined with
alternative non-prompt muon selections, such as changing the muon isolation
criteria in data control samples (7.6%), and the uncertainties of the $Z$ \+
jets and $t\bar{t}$ event fractions when combining the fake-factors derived
from two control samples (8.6%). The statistical uncertainties of the control
samples are also accounted as part of the systematic uncertainties (6.1%).

Additional reducible background contributions come from the $WZ$ process.
These events contain three prompt muons from $W$ and $Z$ decays and one non-
prompt muon. This background is estimated by scaling a $3\mu+\mu_{j}$ control
sample by the $f$ derived from the $Z$ \+ jets sample. The total estimated
number of reducible background events, $61.1^{+8.3}_{-9.1}$, is listed in
Table 3.

##  6 Event Classification with Deep Learning Approach

The selected 4$\mu$ events are classified with a "Deep Learning" approach to
further separate the $Z^{\prime}$ signal from the SM background. A
parametrized deep neural network (pDNN) architecture [70] is used in the
analysis. This algorithm allows the training of a single classifier for
multiple signal mass hypotheses in the search range. In practice, the
kinematic inputs together with the $Z^{\prime}$ mass parameters of signal and
background are used for training. The MC generated $Z^{\prime}$ masses (listed
in Table 1) are used as the multiple signal mass parameters, while the
distribution of the mass parameter for background events is randomly drawn
from the same distribution as for the signal class. The algorithm was
implemented in the Keras [71] framework with the TensorFlow [72] backend. Two
classifiers are trained for low (high) $Z^{\prime}$ mass searches using MC
mass parameters lower (higher) than 40 GeV. During the training of the
classifier, the training samples were composed of simulated signal and
background 4$\mu$ events using the pre-selection described in Section 4. A set
of kinematic distributions were used for pDNN training input features. They
are the $p_{\text{T}}$ and $\eta$ of each muon, the $p_{\text{T}}$ of the
$Z_{1}$ and $Z_{2}$, the mass difference of the $Z_{1}$ and $Z_{2}$, $\Delta
R$ and $\Delta\eta$ of each muon pair that forms the $Z_{1}$ and $Z_{2}$, and
the $p_{\text{T}}$ and mass of the $4\mu$ system. Examples of the input
variable distributions, both for representative signals ($m_{Z^{\prime}}$ =
15, and 51 GeV) and background, are shown in Figure 2 to compare the predicted
and observed transverse momenta ($p_{\textup{T}}^{Z_{1}}$ and
$p_{\textup{T}}^{Z_{2}}$) distributions and the mass difference of the $Z_{1}$
and $Z_{2}$. In addition to the major background from the SM
$Z(Z^{*})\rightarrow 4\mu$ production, other backgrounds, including 4$\mu$
events containing non-prompt muons estimated from data, and from $ttV$, $VVV$,
and Higgs boson production processes, are included in the plots.

((a))  ((b))  ((c))

Figure 2: Distributions of $p_{\text{T}}^{Z_{1}}$ (a), and
$p_{\text{T}}^{Z_{2}}$ (b), and the mass difference of the $Z_{1}$ and $Z_{2}$
candidates (c). Small background contributions are denoted as "other
backgrounds", including 4$\mu$ events containing non-prompt muons estimated
from data and from $ttV$, $VVV$, and Higgs boson production processes.

The training results are a set of weights of the pDNN model which are applied
to real data or MC samples to obtain the pDNN output discriminating variable.
The $m_{Z1}$ or $m_{Z2}$ are used as the mass parameter when applying the
model to data.

((a)) low mass: $m_{Z^{\prime}}\leq 42\leavevmode\nobreak\
\text{Ge\kern-1.00006ptV}$

((b)) high mass: $m_{Z^{\prime}}\geq 42\leavevmode\nobreak\
\text{Ge\kern-1.00006ptV}$

Figure 3:  The pDNN output discriminant variable distributions for low mass
(a) and high mass (b) with a signal sample at 35 GeV and 51 GeV, respectively.

To improve the power of the pDNN classifier the model is optimized through an
automatic process developed with the package Tune [73] and HEBO Search
algorithm [74]. The neural-network structure was chosen to have two fully
connected hidden layers, each with 256 (32) nodes for the high (low) mass
search in the analysis. Other training hyper-parameters, such as the learning
rate with decay and class weight, for both high and low mass searches were
also selected. Two pDNN output discriminant classifiers for high and low mass
searches were obtained and are shown in Figure 3 for data and simulated signal
and background events.

To maximize the search sensitivity, a scan of the pDNN output scores was
performed to find the optimal cut values for each $Z^{\prime}$ mass
hypothesis. These cut values vary from 0.42 to 0.74 (0.12 to 0.16) for the low
(high) mass region. These cuts keep high signal efficiencies between 98% and
95% (90% and 50%), while the corresponding background reductions range from
10% to 50% (50% to 96%) for the low (high) mass region.

After the 4$\mu$ event selection with the pDNN classifier, the final
discriminant to search for the $Z^{\prime}$ resonance signature is the $Z_{1}$
(for $m_{Z^{\prime}}\geq 42$ GeV) or $Z_{2}$ (for $m_{Z^{\prime}}\leq 42$ GeV)
mass spectrum as shown in Figure 4. Data are compared to the estimated
background together with two representative signals with masses of 15 and 51
GeV, shown in Figure 4(a) and 4(b), respectively. The values of the gauge
coupling strengths ($g$) for the two mass points are chosen for the purpose of
illustration.

One should note that the dimuon masses ($m_{Z1},\leavevmode\nobreak\ m_{Z2}$)
are not used as the pDNN training variables, therefore cutting on the pDNN
output scores does not sculpt the dimuon invariant mass distributions of the
backgrounds (see Figure 4). To interpret the observed mass spectra systematic
uncertainties in the predictions are estimated in Section 7.

((a))

((b))

Figure 4:  Mass spectra of $m_{Z2}$ (left) and $m_{Z1}$ (right) for the pDNN-
selected events with a signal sample at 15 GeV and 51 GeV, respectively.

##  7 Systematic Uncertainties

Systematic uncertainties in the simulated event yields and shapes, for both
signal and background processes, may arise from the calibration of the physics
objects and from the theoretical modelling used in the predictions.

The major experimental uncertainties come from the muon reconstruction,
identification, and isolation requirement efficiencies. These efficiencies are
corrected based on studies performed in data control regions. The energy and
momentum scales and resolutions of the simulated objects are calibrated to
reproduce data from $Z\rightarrow\mu^{+}\mu^{-}$ and
$J/\psi\rightarrow\mu^{+}\mu^{-}$ decays [61]. The uncertainties on the 4$\mu$
detection efficiency are determined by varying the nominal calibrations in the
MC samples by one standard deviation, including muon momentum resolutions and
scales, and the trigger, reconstruction, identification, and isolation
requirement efficiencies. The overall relative experimental uncertainties in
the 4$\mu$ event selection efficiency is about 3.9%, dominated by the
uncertainty of the isolation efficiency (2.9%) and the low $p_{\textup{T}}$
calibration uncertainty (2.0%). The signal event selection efficiency
uncertainties vary from 8.3% to 3.9%, depending on the $Z^{\prime}$ mass. In
addition, the uncertainty in the combined 2015-2018 integrated luminosity is
1.7% [75], obtained using the LUCID-2 detector [76] for the primary luminosity
measurements.

Sources of theoretical uncertainties come from the choice of QCD scales
(renormalization $\mu_{R}$ and factorization $\mu_{F}$), strong coupling
constant $\alpha_{S}$, and PDFs, as well as the parton shower models. These
uncertainties affect the signal and background event selection efficiencies,
normalization, and the shape of their kinematic distributions. The scales are
varied independently from 0.5 to 2.0 times the nominal values and the largest
deviation is chosen as the systematic uncertainty. The PDF uncertainty is
estimated by comparing events generated with different PDF sets, as well as
the uncertainties from the nominal PDF set itself. The maximal variation
(envelope) is accounted for as the systematic uncertainty, following the
PDF4LHC [77] recommendations. The $\alpha_{S}$ uncertainty is estimated by
comparing events generated with different $\alpha_{S}$ values using the
nominal PDF set. The parton shower uncertainty is estimated by comparing
events with different parton shower parameters in the Sherpa MC samples. For
the $Z(Z^{*})\rightarrow 4\mu$ process, the relative uncertainties of event
yields for scale, PDF, $\alpha_{s}$, and parton shower, are 4.6%, 1.8%, 1.0%,
and 1.9%, respectively. The $t\bar{t}V$, $VVV$, and Higgs boson, processes
contribute 5.2, 0.50, and 0.53 background events, respectively. The total
relative theoretical uncertainty for these background is estimated about 16%.
The interference effect between the signal and background is estimated (using
MC samples generated by MadGraph5_aMC@NLO) by evaluating the cross-section
ratio, ${\Delta\sigma}/{\sigma_{Z^{\prime}}}$, where $\Delta\sigma$ is the
difference of the cross-section from the inclusive MC sample and the sum of
the cross-sections of the signal and background samples in the $Z^{\prime}$
detection phase space. The effect varies from 1.8 to 7.5% for the $Z^{\prime}$
mass from 5 - 81 GeV, which is accounted for as additional signal yield
systematic uncertainties. In addition, the MC $Z^{\prime}$ signal filter
acceptance uncertainty is estimated to be about 2%, which is calculated by
varying the QCD scales, the PDF-sets, and the strong coupling constant using
MC events at the generator level for different $Z^{\prime}$ mass points.

All the uncertainties (from both experiment and theory) on the final
discriminant di-muon mass spectra are included as nuisance parameters in the
signal extraction fitting process described in Section 8.1.

##  8 Data Interpretation and Results

The statistical analysis is performed by comparing the data to the sum of the
background prediction and the signal to search for the $Z^{\prime}$ signature
and information about the $pp\rightarrow\mu^{+}\mu^{-}Z^{\prime}\rightarrow
4\mu$ signal production cross-section and the associated coupling strength. In
case of no significant data excess over the background prediction, upper
limits on the signal production cross-section times the decay branching
fraction for different $Z^{\prime}$ masses are set at 95% confidence level
(CL).

###  8.1 Statistical Procedure and Results

The statistical procedure used to interpret the data is described in Ref.
[78]. The likelihood function is constructed from the product of Poisson
probabilities:

$L(\text{data}|\leavevmode\nobreak\
\mu,\theta)=\prod_{i=1}^{N}\text{Poisson}(\text{data}_{i}|\leavevmode\nobreak\
\mu\cdot s_{i}(\theta)+b_{i}(\theta))\times{\text{G}}(\tilde{\theta}|\theta)$

where $i$ is the bin-index of the fitted variable distribution; $\mu$ is the
signal strength, and $\theta$ denotes the nuisance parameters, which represent
the uncertainties of the measurements; ${\text{G}}(\tilde{\theta}|\theta)$ is
the Gaussian scaling function of the nuisance parameters constructed as
deviations from the nominal model of the systematic uncertainties, where
$\tilde{\theta}$ provides a maximum likelihood estimate for $\theta$. The
parameter of interest in the statistical analysis is the global signal
strength factor $\mu$, which acts as a scale factor on the total number of
events predicted by the signal model. This factor is defined such that $\mu=0$
corresponds to the background-only hypothesis and $\mu>0$ corresponds to a
$Z^{\prime}$ signal in addition to the background. Hypothesised values of
$\mu$ are tested based on the profile likelihood ratio [79], which compares
data with background-only ($b$) and signal+background ($s+b$) models using the
following test statistic:

$q_{\mu}=\begin{cases}-2\text{ln}\frac{L(\text{data}|\leavevmode\nobreak\
\mu,\hat{\hat{\theta}}_{\mu})}{L(\text{data}|\leavevmode\nobreak\
\hat{\mu},\hat{\theta})}&\hat{\mu}\leq\mu\\\ 0&\hat{\mu}>\mu\\\ \end{cases}$

where $\hat{\mu}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\hat{\theta}$ are the values for $\mu$ and $\theta$ when maximizing the $L$
with all parameters floating, which are named as the unconditional maximum-
likelihood (ML) estimators; $\hat{\hat{\theta}}_{\mu}$ is the conditional ML
estimator of $\theta$ for a fixed value of $\mu$. This test statistic extracts
the information on the signal strength from a full likelihood fit to the data.
The likelihood function includes all the parameters that describe the
systematic uncertainties and their correlations.

In case of no significant data excess over the background, exclusion limits
are set, based on the CLs prescription [80], which calculates the ratio of the
probabilities based on the $s+b$ background-only models:

$\text{CL}_{s}(\mu)=\frac{\text{CL}_{s+b}}{\text{CL}_{b}}=\frac{\text{P}(q_{\mu}>q_{\mu}^{obs})}{\text{P}(q_{0}>q_{0}^{obs})}$

A value of $\mu$ is regarded as excluded at 95% confident level when CLs is
less than 5%.

The significance of an excess in the data is first quantified with the local
$p_{0}$, the probability that the background can produce a fluctuation greater
than or equal to the excess observed in data. The equivalent formulation in
terms of number of standard deviations is referred to as the local
significance, which is computed from test statistic $q_{0}$:

$p_{0}=\text{P}(q_{0}<q_{0}^{obs}).$

The global probability for the most significant excess to be observed anywhere
in a given search region is estimated with the method described in Ref. [81].

The statistical tests are performed by first searching for data excesses by
means of a $p_{0}$ scan, and in case no such significant data excess is found,
by setting limits on the $Z^{\prime}$ production cross-section times the
branching fraction. In both steps the values and reconstructed mass
distributions of hypothesized $Z^{\prime}$ mass $m_{Z^{\prime}}$ are used. The
steps between the mass points used in these test procedures are set to match
the mass resolutions determined by reconstruction of the di-muon invariant
mass. A linear interpolation of the expected event yields and shapes between
the MC generated signal samples is used in the fitting process. The asymptotic
approximation [79] upon which the results are based has been validated against
the method (detailed in Ref. [78]) of using pseudo-experiments for several
mass points.

In both steps of the statistical tests, data are fit to the $m_{Z1}$ and
$m_{Z2}$ spectrum with background ($b$) only and signal+background ($s+b$)
hypotheses. In the fitting process each mass spectrum is divided in the signal
region (SR) and the background control region (CR). For each $Z^{\prime}$ mass
point the SR is defined in a mass window of $m_{Z^{\prime}}\pm
3\sigma_{m_{\mu\mu}}$ of the di-muon mass spectrum. The sidebands outside of
the SR are defined as the CR. The di-muon mass resolution
$\sigma_{m_{\mu\mu}}$ is determined by the fully simulated $Z^{\prime}$ mass
distribution, which combines the $Z^{\prime}$ natural width and the detector
resolution, ranging from 0.10 to 1.75 GeV. The mass resolution is mostly
dominated by the detector resolution. Finer binning is used in the SR to
enhance the sensitivity. The background CR is used to constrain the overall
normalization for the background in the signal region. The shape of the major
background from $Z(Z^{*})\rightarrow 4\mu$ is fixed with prior uncertainties
included in the fitting process, but the normalization (or strength) floats in
the fit. Other background normalizations and shapes are fixed with prior
uncertainties included in the fitting.

###  8.2 The $p_{0}$ Scan Results

The $p_{0}$-values corresponding to the background-only hypothesis are scanned
in the mass range of this analysis. A binned profile-likelihood fit [79] is
performed simultaneously across the $Z^{\prime}$ signal-region and the
background control region using the predicted and observed mass spectrum as
inputs. Data are fit to the $m_{Z1}$ and $m_{Z2}$ distributions for
$m_{Z^{\prime}}\geq 42\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ and
$m_{Z^{\prime}}\leq 42\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}{}$,
respectively, with a "sliding" mass window as the defined SR changes for
different $Z^{\prime}$ mass points. The chosen bin-size inside the SR is
around 0.3 $\sigma_{m_{\mu\mu}}$, for each mass point. The total number of
bins in the CR is 20. The fit mass range of $m_{Z1}$ ($m_{Z2}$) is [30, 85]
GeV ([0, 45] GeV).

The $p_{0}$-values at different $Z^{\prime}$ mass hypothesis points are
computed and transformed into Gaussian standard deviations to indicate the
significance as shown in Figure 5. The smallest $p_{0}$-value is at 39.6 GeV,
corresponding to a local 2.65$\sigma$ deviation from the background-only
hypothesis, while the global deviation [81] is found to be 0.52$\sigma$,
indicating that no significant data excess over the expected background is
observed.

Figure 5: The $p_{0}$-value scan across the $Z^{\prime}$ mass signal regions.

###  8.3 Upper Limits

The upper limits on the production cross-section times branching fraction of
the $pp\rightarrow\mu^{+}\mu^{-}Z^{\prime}\rightarrow 4\mu$ process are
calculated using a similar fitting procedure described in Section 8.2.
Confidence intervals are computed based on the profile-likelihood-ratio test
statistics [79]. The observed and expected upper limits at 95% CL on the
cross-section times branching fraction, $\sigma(pp\rightarrow
Z^{\prime}\mu\mu\rightarrow 4\mu)$, are shown in Figure 6(a). Assuming the
same coupling strength $g$ of the $Z^{\prime}$ boson to the second and third
lepton families and to the left- and right-handed fermions, the upper limits
on the coupling parameter $g$ are extracted from the limits of the
$Z^{\prime}$ production cross-section times branching fraction using the
$L_{\mu}-L_{\tau}$ model, where the branching fraction of
$\mathcal{B}$($Z^{\prime}\rightarrow\mu^{+}\mu^{-})={1\over 3}$, which is
determined by counting all the possible $Z^{\prime}$ decay modes in this
model. At each generated $Z^{\prime}$ mass point, a limit on the coupling
strength $g$ has been obtained from the cross-section limit. The observed and
expected upper limits on the coupling parameter $g$ are shown in Figure 6(b).
The limits on the coupling $g$ are in the range of 0.003 (for $m_{Z^{\prime}}$
= 5 GeV) to 0.2 (for $m_{Z^{\prime}}$ = 81 GeV) depending on the $Z^{\prime}$
mass ranging from 5 to 81 GeV. This ensures that the ratio of the $Z^{\prime}$
natural width and mass, $\Gamma(Z^{\prime})/m_{Z^{\prime}}$, is well below 1%
in this mass range.

((a))

((b))

Figure 6: 95% CL upper limits (expected and observed) on the cross-sections
times branching fraction (a) and coupling parameter (b). The discontinuity at
42 GeV represents the border of the low/high mass classifiers.

Motivated by theoretical interpretations in Ref. [15], a 2-dimensional
exclusion contour at 95% CL in the parameter-space of $(m_{Z^{\prime}},g)$ of
the $L_{\mu}-L_{\tau}$ model from this analysis is produced and shown in
Figure 7. The parameter space exclusion regions calculated by theorists using
data from the Neutrino Trident experiment [82] and the $B_{s}$ mixing
measurements by a global analysis performed in Ref. [15] are also shown in
Figure 7. This had left a large gap in the parameter space not yet excluded,
that was allowed to explain the LHCb $b\rightarrow s\mu^{+}\mu^{-}$ anomalies
[10, 11]. This gap is now largely excluded by this analysis.

Figure 7:  The coupling parameter $g$ limits from this search as a function of
the $Z^{\prime}$ mass compared to the limits [15] from the Neutrino Trident
(red) and the $B_{s}$ mixing (green) experimental results.

##  9 Conclusion

A search for a new vector gauge boson $Z^{\prime}$ predicted by the
$L_{\mu}-L_{\tau}$ models has been performed with a 4$\mu$ final state in the
invariant mass range of [80, 180] GeV, excluding the Higgs boson mass window
[110, 130] GeV, using 139 fb-1 of $\sqrt{s}=$ 13 TeV proton-proton collision
data collected with the ATLAS detector. No significant excess of events over
the expected SM background is observed. Therefore, upper limits are set on the
$Z^{\prime}$ production cross-section times the decay branching fraction of
the $pp\rightarrow
Z^{\prime}\mu^{+}\mu^{-}\rightarrow\mu^{+}\mu^{-}\mu^{+}\mu^{-}$ process,
varying from 0.31 to 4.3 fb at 95% CL, in a $Z^{\prime}$ mass range of [5, 81]
GeV, from which the coupling strength $g$ of the $Z^{\prime}$ to muons above
0.003 to 0.2 (depending on the $Z^{\prime}$ mass) are excluded in the same
mass range. This search shows significant sensitivity improvements over
previous indirect and direct searches of the $Z^{\prime}$ with 4$\mu$ final
state. An interesting parameter space of the $L_{\mu}-L_{\tau}$ model
prediction that was not excluded before and could explain the $b\rightarrow
s\mu^{+}\mu^{-}$ anomalies is now largely excluded by this search.

## References

  * [1] ATLAS Collaboration  “The ATLAS Experiment at the CERN Large Hadron Collider”  In _JINST_ 3, 2008, pp. S08003  DOI: 10.1088/1748-0221/3/08/S08003
  * [2] Xiao-Gang He, G.. Joshi, H. Lew and R.. Volkas  “Simplest ${Z}^{{}^{\prime}}$ model”  In _Phys. Rev. D_ 44 American Physical Society, 1991, pp. 2118–2132  DOI: 10.1103/PhysRevD.44.2118
  * [3] Fatemeh Elahi and Adam Martin  “Constraints on $L_{\mu}-L_{\tau}$ interactions at the LHC and beyond”  In _Phys. Rev. D_ 93.1 American Physical Society (APS), 2016  DOI: 10.1103/physrevd.93.015022
  * [4] Fatemeh Elahi and Adam Martin  “Using the modified matrix element method to constrain $L_{\mu}-L_{\tau}$ interactions”  In _Phys. Rev. D_ 96.1 American Physical Society (APS), 2017  DOI: 10.1103/physrevd.96.015021
  * [5] Andreas Crivellin, Giancarlo D’Ambrosio and Julian Heeck  “Explaining $h\rightarrow{\mu}^{\pm{}}{\tau}^{\mp}$, $B\rightarrow{K}^{*}{\mu}^{+}{\mu}^{-}$, and $B\rightarrow K{\mu}^{+}{\mu}^{-}/B\rightarrow K{e}^{+}{e}^{-}$ in a Two-Higgs-Doublet Model with Gauged ${L}_{\mu}\text{$-$}{L}_{\tau}$”  In _Phys. Rev. Lett._ 114 American Physical Society, 2015, pp. 151801  DOI: 10.1103/PhysRevLett.114.151801
  * [6] Muon g-2 Collaboration  “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm”  In _Phys. Rev. Lett._ 126.14, 2021, pp. 141801  DOI: 10.1103/PhysRevLett.126.141801
  * [7] Alexander Keshavarzi  “The Muon g-2 experiment at Fermilab”  In _EPJ Web of Conferences_ 212 EDP Sciences, 2019, pp. 05003  DOI: 10.1051/epjconf/201921205003
  * [8] Graziano Venanzoni  “The Fermilab Muon g-2 Experiment”  In _PoS_ EPS-HEP2015, 2016, pp. 568  DOI: 10.22323/1.234.0568
  * [9] Muon g-2 Collaboration  “Final report of the E821 muon anomalous magnetic moment measurement at BNL”  In _Phys. Rev. D_ 73 American Physical Society, 2006, pp. 072003  DOI: 10.1103/PhysRevD.73.072003
  * [10] LHCb Collaboration  “Test of lepton universality in beauty-quark decays”  In _Nat. Phys._ 18.3, 2022, pp. 277–282  DOI: 10.1038/s41567-021-01478-8
  * [11] LHCb Collaboration  “Measurement of Form-Factor-Independent Observables in the Decay $B^{0}\rightarrow K^{*0}\mu^{+}\mu^{-}$”  In _Phys. Rev. Lett._ 111.19 American Physical Society (APS), 2013, pp. 191801  DOI: 10.1103/physrevlett.111.191801
  * [12] LHCb Collaboration  “Test of lepton universality using $B^{+}\rightarrow K^{+}\ell^{+}\ell^{-}$ decays”  In _Phys. Rev. Lett._ 113.15 American Physical Society (APS), 2014, pp. 151601  DOI: 10.1103/physrevlett.113.151601
  * [13] LHCb Collaboration  “Measurement of $CP$-Averaged Observables in the ${B}^{0}\rightarrow{K}^{*0}{\mu}^{+}{\mu}^{-}$ Decay”  In _Phys. Rev. Lett._ 125 American Physical Society, 2020, pp. 011802  DOI: 10.1103/PhysRevLett.125.011802
  * [14] Thomas Blake, Gaia Lanfranchi and David M. Straub  “Rare $B$ Decays as Tests of the Standard Model”  In _Prog. Part. Nucl. Phys._ 92 Elsevier BV, 2017, pp. 50–91  DOI: 10.1016/j.ppnp.2016.10.001
  * [15] Wolfgang Altmannshofer, Stefania Gori, Stefano Profumo and Farinaldo S. Queiroz  “Explaining dark matter and B decay anomalies with an $L_{\mu}-L_{\tau}$ model”  In _JHEP_ 12, 2016, pp. 106  DOI: 10.1007/JHEP12(2016)106
  * [16] Ernest Ma, D.P. Roy and Sourov Roy  “Gauged $L_{\mu}-L_{\tau}$ with large muon anomalous magnetic moment and the bimaximal mixing of neutrinos”  In _Phys. Lett. B_ 525.1-2 Elsevier BV, 2002, pp. 101–106  DOI: 10.1016/s0370-2693(01)01428-9
  * [17] Paul Langacker  “The physics of heavy Z’ gauge bosons”  In _Rev. Mod. Phys._ 81.3 American Physical Society (APS), 2009, pp. 1199–1228  DOI: 10.1103/revmodphys.81.1199
  * [18] CMS Collaboration  “Observation of Z Decays to Four Leptons with the CMS Detector at the LHC”  In _JHEP_ 12, 2012, pp. 034  DOI: 10.1007/JHEP12(2012)034
  * [19] ATLAS Collaboration  “Measurements of Four-Lepton Production at the Z Resonance in pp Collisions at $\sqrt{s}=$7 and 8 TeV with ATLAS”  In _Phys. Rev. Lett._ 112.23, 2014, pp. 231806  DOI: 10.1103/PhysRevLett.112.231806
  * [20] CMS Collaboration  “Measurements of the $pp\rightarrow ZZ$ production cross section and the $Z\rightarrow 4\ell$ branching fraction, and constraints on anomalous triple gauge couplings at $\sqrt{s}=13\,\text{TeV}$”  In _Eur. Phys. J. C_ 78, 2018, pp. 165  DOI: 10.1140/epjc/s10052-018-5567-9
  * [21] CMS Collaboration  “Search for an $L_{\mu}$–$L_{\tau}$ gauge boson using $Z\to 4\mu$ events in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$”  In _Phys. Lett. B_ 792, 2019, pp. 345  DOI: 10.1016/j.physletb.2019.01.072
  * [22] ATLAS Collaboration  “ATLAS Insertable B-Layer: Technical Design Report”, 2010  URL: https://cds.cern.ch/record/1291633
  * [23] B. Abbott  “Production and integration of the ATLAS Insertable B-Layer”  In _JINST_ 13, 2018, pp. T05008  DOI: 10.1088/1748-0221/13/05/T05008
  * [24] ATLAS Collaboration  “Performance of the ATLAS trigger system in 2015”  In _Eur. Phys. J. C_ 77, 2017, pp. 317  DOI: 10.1140/epjc/s10052-017-4852-3
  * [25] ATLAS Collaboration  “The ATLAS Collaboration Software and Firmware”, ATL-SOFT-PUB-2021-001, 2021  URL: https://cds.cern.ch/record/2767187
  * [26] ATLAS Collaboration  “ATLAS data quality operations and performance for 2015-2018 data-taking”  In _JINST_ 15.04 IOP Publishing, 2020, pp. P04003–P04003  DOI: 10.1088/1748-0221/15/04/p04003
  * [27] ATLAS Collaboration  “Performance of the ATLAS muon triggers in Run 2”  In _JINST_ 15, 2020, pp. P09015  DOI: 10.1088/1748-0221/15/09/p09015
  * [28] J. Alwall et al.  “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”  In _JHEP_ 07, 2014, pp. 079  DOI: 10.1007/JHEP07(2014)079
  * [29] Celine Degrande et al.  “UFO - The Universal FeynRules Output”  In _Comput. Phys. Commun._ 183, 2012, pp. 1201–1214  DOI: 10.1016/j.cpc.2012.01.022
  * [30] Adam Alloul et al.  “FeynRules 2.0 - A complete toolbox for tree-level phenomenology”  In _Comput. Phys. Commun._ 185, 2014, pp. 2250–2300  DOI: 10.1016/j.cpc.2014.04.012
  * [31] Francisco Aguila, Mikael Chala, Jose Santiago and Yasuhiro Yamamoto  “Collider limits on leptophilic interactions”  In _JHEP_ 03, 2015, pp. 059  DOI: 10.1007/JHEP03(2015)059
  * [32] Massimiliano Grazzini, Stefan Kallweit and Dirk Rathlev  “ZZ production at the LHC: Fiducial cross sections and distributions in NNLO QCD”  In _Phys. Lett. B_ 750 Elsevier BV, 2015, pp. 407–410  DOI: 10.1016/j.physletb.2015.09.055
  * [33] A. H. Ajjath et al.  “Resummed Drell-Yan cross-section at N3LL”  In _JHEP_ 10 arXiv, 2020, pp. 153  DOI: 10.48550/ARXIV.2001.11377
  * [34] Richard D. Ball  “Parton distributions with LHC data”  In _Nucl. Phys. B_ 867, 2013, pp. 244  DOI: 10.1016/j.nuclphysb.2012.10.003
  * [35] N. Davidson, T. Przedzinski and Z. Was  “PHOTOS Interface in C++: Technical and physics documentation”  In _Comput. Phys. Commun._ 199, 2016, pp. 86–101  DOI: 10.1016/j.cpc.2015.09.013
  * [36] ATLAS Collaboration  “The ATLAS Simulation Infrastructure”  In _Eur. Phys. J. C_ 70, 2010, pp. 823  DOI: 10.1140/epjc/s10052-010-1429-9
  * [37] Enrico Bothmann  “Event generation with Sherpa 2.2”  In _SciPost Phys._ 7.3, 2019, pp. 034  DOI: 10.21468/SciPostPhys.7.3.034
  * [38] Tanju Gleisberg and Stefan Höche  “Comix, a new matrix element generator”  In _JHEP_ 12, 2008, pp. 039  DOI: 10.1088/1126-6708/2008/12/039
  * [39] Steffen Schumann and Frank Krauss  “A parton shower algorithm based on Catani–Seymour dipole factorisation”  In _JHEP_ 03, 2008, pp. 038  DOI: 10.1088/1126-6708/2008/03/038
  * [40] Stefan Höche, Frank Krauss, Marek Schönherr and Frank Siegert  “A critical appraisal of NLO+PS matching methods”  In _JHEP_ 09, 2012, pp. 049  DOI: 10.1007/JHEP09(2012)049
  * [41] Stefan Höche, Frank Krauss, Marek Schönherr and Frank Siegert  “QCD matrix elements + parton showers. The NLO case”  In _JHEP_ 04, 2013, pp. 027  DOI: 10.1007/JHEP04(2013)027
  * [42] S. Catani, F. Krauss, R. Kuhn and B.. Webber  “QCD Matrix Elements + Parton Showers”  In _JHEP_ 11, 2001, pp. 063  DOI: 10.1088/1126-6708/2001/11/063
  * [43] Stefan Höche, Frank Krauss, Steffen Schumann and Frank Siegert  “QCD matrix elements and truncated showers”  In _JHEP_ 05, 2009, pp. 053  DOI: 10.1088/1126-6708/2009/05/053
  * [44] Federico Buccioni et al.  “OpenLoops 2”  In _Eur. Phys. J. C_ 79.10, 2019, pp. 866  DOI: 10.1140/epjc/s10052-019-7306-2
  * [45] Fabio Cascioli, Philipp Maierhofer and Stefano Pozzorini  “Scattering Amplitudes with Open Loops”  In _Phys. Rev. Lett._ 108, 2012, pp. 111601  DOI: 10.1103/PhysRevLett.108.111601
  * [46] Ansgar Denner, Stefan Dittmaier and Lars Hofer  “Collier: A fortran-based complex one-loop library in extended regularizations”  In _Comput. Phys. Commun._ 212, 2017, pp. 220–238  DOI: 10.1016/j.cpc.2016.10.013
  * [47] Richard D. Ball  “Parton distributions for the LHC run II”  In _JHEP_ 04, 2015, pp. 040  DOI: 10.1007/JHEP04(2015)040
  * [48] Stefano Frixione, Paolo Nason and Giovanni Ridolfi  “A positive-weight next-to-leading-order Monte Carlo for heavy flavour hadroproduction”  In _JHEP_ 09, 2007, pp. 126  DOI: 10.1088/1126-6708/2007/09/126
  * [49] Paolo Nason  “A new method for combining NLO QCD with shower Monte Carlo algorithms”  In _JHEP_ 11, 2004, pp. 040  DOI: 10.1088/1126-6708/2004/11/040
  * [50] Stefano Frixione, Paolo Nason and Carlo Oleari  “Matching NLO QCD computations with parton shower simulations: the POWHEG method”  In _JHEP_ 11, 2007, pp. 070  DOI: 10.1088/1126-6708/2007/11/070
  * [51] Simone Alioli, Paolo Nason, Carlo Oleari and Emanuele Re  “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX”  In _JHEP_ 06, 2010, pp. 043  DOI: 10.1007/JHEP06(2010)043
  * [52] ATLAS Collaboration  “Studies on top-quark Monte Carlo modelling for Top2016”, ATL-PHYS-PUB-2016-020, 2016  URL: https://cds.cern.ch/record/2216168
  * [53] Torbjörn Sjöstrand et al.  “An introduction to PYTHIA 8.2”  In _Comput. Phys. Commun._ 191, 2015, pp. 159  DOI: 10.1016/j.cpc.2015.01.024
  * [54] ATLAS Collaboration  “ATLAS Pythia 8 tunes to $7\;\mbox{TeV}$ data”, ATL-PHYS-PUB-2014-021, 2014  URL: https://cds.cern.ch/record/1966419
  * [55] Emanuele Re  “Single-top $Wt$-channel production matched with parton showers using the POWHEG method”  In _Eur. Phys. J. C_ 71, 2011, pp. 1547  DOI: 10.1140/epjc/s10052-011-1547-z
  * [56] Stefano Frixione et al.  “Single-top hadroproduction in association with a $W$ boson”  In _JHEP_ 07, 2008, pp. 029  DOI: 10.1088/1126-6708/2008/07/029
  * [57] Charalampos Anastasiou, Lance J. Dixon, Kirill Melnikov and Frank Petriello  “High precision QCD at hadron colliders: Electroweak gauge boson rapidity distributions at next-to-next-to leading order”  In _Phys. Rev. D_ 69, 2004, pp. 094008  DOI: 10.1103/PhysRevD.69.094008
  * [58] S. Agostinelli  “GEANT4 – a simulation toolkit”  In _Nucl. Instrum. Meth. A_ 506, 2003, pp. 250  DOI: 10.1016/S0168-9002(03)01368-8
  * [59] Torbjorn Sjöstrand, Stephen Mrenna and Peter Z. Skands  “A brief introduction to PYTHIA 8.1”  In _Comput. Phys. Commun._ 178, 2008, pp. 852  DOI: 10.1016/j.cpc.2008.01.036
  * [60] ATLAS Collaboration  “The Pythia 8 A3 tune description of ATLAS minimum bias and inelastic measurements incorporating the Donnachie–Landshoff diffractive model”, ATL-PHYS-PUB-2016-017, 2016  URL: https://cds.cern.ch/record/2206965
  * [61] ATLAS Collaboration  “Muon reconstruction and identification efficiency in ATLAS using the full Run 2 $pp$ collision data set at $\sqrt{s}=13\,\text{TeV}$”  In _Eur. Phys. J. C_ 81.7, 2021, pp. 578  DOI: 10.1140/epjc/s10052-021-09233-2
  * [62] ATLAS Collaboration  “Jet reconstruction and performance using particle flow with the ATLAS Detector”  In _Eur. Phys. J. C_ 77.7 Springer ScienceBusiness Media LLC, 2017  DOI: 10.1140/epjc/s10052-017-5031-2
  * [63] ATLAS Collaboration  “Electron and photon performance measurements with the ATLAS detector using the 2015–2017 LHC proton-proton collision data”  In _JINST_ 14.12 IOP Publishing, 2019, pp. P12006–P12006  DOI: 10.1088/1748-0221/14/12/p12006
  * [64] ATLAS Collaboration  “Track assisted techniques for jet substructure”, ATL-PHYS-PUB-2018-012, 2018  URL: https://cds.cern.ch/record/2630864
  * [65] Matteo Cacciari, Gavin P. Salam and Gregory Soyez  “FastJet user manual”  In _Eur. Phys. J. C_ 72, 2012, pp. 1896  DOI: 10.1140/epjc/s10052-012-1896-2
  * [66] ATLAS Collaboration  “Jet reconstruction and performance using particle flow with the ATLAS Detector”  In _Eur. Phys. J. C_ 77, 2017, pp. 466  DOI: 10.1140/epjc/s10052-017-5031-2
  * [67] ATLAS Collaboration  “Optimisation of the ATLAS $b$-tagging performance for the 2016 LHC Run”, ATL-PHYS-PUB-2016-012, 2016  URL: https://cds.cern.ch/record/2160731
  * [68] ATLAS Collaboration  “Tagging and suppression of pileup jets with the ATLAS detector”, ATLAS-CONF-2014-018, 2014  URL: https://cds.cern.ch/record/1700870
  * [69] ATLAS Collaboration  “Performance of missing transverse momentum reconstruction with the ATLAS detector using proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$”  In _Eur. Phys. J. C_ 78, 2018, pp. 903  DOI: 10.1140/epjc/s10052-018-6288-9
  * [70] Pierre Baldi et al.  “Parameterized neural networks for high-energy physics”  In _Eur. Phys. J. C_ 76.5 Springer ScienceBusiness Media LLC, 2016  DOI: 10.1140/epjc/s10052-016-4099-4
  * [71] François Chollet  “Keras”, https://keras.io, 2015
  * [72] Mart’ın Abadi et al.  “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems” Software available from tensorflow.org, 2015  URL: https://www.tensorflow.org/
  * [73] Richard Liaw et al.  “Tune: A Research Platform for Distributed Model Selection and Training”, 2018  arXiv:1807.05118 [cs.LG]
  * [74] Alexander I. Cowen-Rivers et al.  “HEBO Pushing The Limits of Sample-Efficient Hyperparameter Optimisation”, 2020  arXiv:2012.03826 [cs.LG]
  * [75] ATLAS Collaboration  “Luminosity determination in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ using the ATLAS detector at the LHC”, ATLAS-CONF-2019-021, 2019  URL: https://cds.cern.ch/record/2677054
  * [76] G. Avoni et al.  “The new LUCID-2 detector for luminosity measurement and monitoring in ATLAS”  In _JINST_ 13.07 IOP Publishing, 2018, pp. P07017–P07017  DOI: 10.1088/1748-0221/13/07/p07017
  * [77] Jon Butterworth  “PDF4LHC recommendations for LHC Run II”  In _J. Phys. G_ 43, 2016, pp. 023001  DOI: 10.1088/0954-3899/43/2/023001
  * [78] ATLAS Collaboration  “Combined search for the Standard Model Higgs boson in $pp$ collisions at $\sqrt{s}=7\,\text{TeV}$ with the ATLAS detector”  In _Phys. Rev. D_ 86, 2012, pp. 032003  DOI: 10.1103/PhysRevD.86.032003
  * [79] Glen Cowan, Kyle Cranmer, Eilam Gross and Ofer Vitells  “Asymptotic formulae for likelihood-based tests of new physics”  In _Eur. Phys. J. C_ 71, 2011, pp. 1554  DOI: 10.1140/epjc/s10052-011-1554-0
  * [80] Alexander L. Read  “Presentation of search results: the $CL_{S}$ technique”  In _J. Phys. G_ 28, 2002, pp. 2693  DOI: 10.1088/0954-3899/28/10/313
  * [81] Eilam Gross and Ofer Vitells  “Trial factors for the look elsewhere effect in high energy physics”  In _Eur. Phys. J. C_ 70.1-2 Springer ScienceBusiness Media LLC, 2010, pp. 525–530  DOI: 10.1140/epjc/s10052-010-1470-8
  * [82] S.. Mishra et al.  “Neutrino tridents and W-Z interference”  In _Phys. Rev. Lett._ 66 American Physical Society, 1991, pp. 3117–3120  DOI: 10.1103/PhysRevLett.66.3117
  * [83] CMS Collaboration  In _Eur. Phys. J. C_ 78.6, 2018, pp. 515  DOI: 10.1140/epjc/s10052-018-5769-1
  * [84] ATLAS Collaboration, 2012  URL: https://cds.cern.ch/record/1451888
  * [85] Glen Cowan, Kyle Cranmer, Eilam Gross and Ofer Vitells  In _Eur. Phys. J. C_ 73, 2013, pp. 2501  DOI: 10.1140/epjc/s10052-013-2501-z