Search for Z ′ bosons decaying to pairs of heavy Majorana neutrinos in proton-proton collisions at √ s = 13 TeV

: A search for the production of pairs of heavy Majorana neutrinos ( N ℓ ) from the decays of Z ′ bosons is performed using the CMS detector at the LHC. The data were collected in proton-proton collisions at a center-of-mass energy of √ s = 13 TeV, with an integrated luminosity of 138 fb − 1 . The signature for the search is an excess in the invariant mass distribution of the final-state objects, two same-flavor leptons (e or µ ) and at least two jets. No significant excess of events beyond the expected background is observed. Upper limits at 95% confidence level are set on the product of the Z ′ production cross section and its branching fraction to a pair of N ℓ , as functions of N ℓ and Z ′ boson masses ( m N ℓ and m Z ′ , respectively) for m Z ′ from 0.4 to 4.6 TeV and m N ℓ from 0.1 TeV to m Z ′ / 2 . In the theoretical framework of a left-right symmetric model, exclusion bounds in the m N ℓ - m Z ′ plane are presented in both the electron and muon channels. The observed upper limit on m Z ′ reaches up to 4.42 TeV. These are the most restrictive limits to date on the mass of N ℓ as a function of the Z ′ boson mass.


Introduction
The observation of neutrino oscillations [1][2][3][4] established that at least two of the standard model (SM) neutrinos have mass.The nonzero masses of the neutrinos are clear evidence of physics beyond the SM, in which neutrinos are massless.Upper limits on the neutrino masses have been obtained from cosmological observations [5].A direct measurement from tritium beta decays [6] indicates that the electron-type neutrino mass is less than 0.8 eV.The fact that neutrino masses are so much smaller than the other fermion masses and that right-handed neutrinos are not observed suggest that neutrino masses may have an origin other than a Yukawa coupling to the Higgs field.
In addition to its inability to explain neutrino mass, the SM does not provide any clear answer as to what is the source of the parity violation in the weak sector or of the matter dominance of the universe.One of the leading theoretical solutions is to introduce a left-right symmetry model (LRSM) [7].The LRSM is based on a gauge group of SU(3) C SU(2) L SU(2) R U(1) B−L , established by a symmetry between the left-and righthanded SU (2) groups.In this model, there are three additional gauge bosons, W ± R and Z ′ , as well as three right-handed neutrinos (N e , N µ , and N τ ).The spontaneous symmetry breaking of the LRSM, at some energy scale, can generate the gauge group of the SM in a way that naturally includes the observed parity violation in the weak sector, and also provides an additional source of CP violation that can explain the asymmetry between matter and antimatter in the universe.Heavy right-handed neutrinos and light left-handed neutrinos are also naturally generated by the symmetry breaking, through a process referred to as the see-saw mechanism [8][9][10][11][12].
As the LRSM provides explanations for these basic observations of beyond the SM physics, it is compelling to probe signatures of this model.Although there have been previous searches for both new heavy gauge bosons and heavy right-handed neutrinos [13][14][15][16][17], the unique signature of the LRSM is the presence of events that have an extra gauge boson as well as heavy righthanded neutrinos.In this analysis, we search for the pair production of heavy right-handed neutrinos through an extra neutral gauge boson (Z ′ ).The resonance production of a Z ′ boson and its decay into heavy neutrinos are shown in Fig. 1.In the benchmark model used in this study, the heavy right-handed neutrino is a Majorana particle and decays into a lepton plus two quarks.We therefore select events with two electrons or two muons, and at least two reconstructed jets, to search for an excess in the invariant mass distribution of them.We use opposite-and same-sign leptons inclusively, and mixing between heavy neutrinos of electron type and muon type is not considered.Rich scalar sectors are predicted by LRSM models in general.In this search, results are presented in terms of the product of the cross section and branching fraction of the process, thus no assumptions are made concerning the nature of the scalar sector.Detailed underlying assumptions for this LRSM will be discussed in Section 3.

Searches for a W ±
R boson and its decay to a heavy Majorana neutrino (N ℓ ) have been performed in the framework of the LRSM by the ATLAS and the CMS Collaborations at the LHC [18][19][20][21][22][23][24][25][26][27][28].The most stringent lower limit on the W ± R mass from the ATLAS Collaboration is about 6.4 TeV [28], which was obtained using a data set of proton-proton (pp) collisions at a centerof-mass energy of 13 TeV, corresponding to an integrated luminosity of 139 fb −1 .A search for an extra gauge boson Z ′ and its decay to a pair of N ℓ has also been performed by the ATLAS Collaboration, which obtained a lower limit on the Z ′ mass (m Z ′ ) of 2.2 TeV [21] using a data sample corresponding to 20.3 fb −1 at √ s = 8 TeV.
All these results are limited by a reduced efficiency for the signals in the kinematic region where the heavy neutrino mass (m N ℓ ) is much smaller than the mass of the extra gauge boson.A lepton from the decay of a highly Lorentz-boosted heavy neutrino overlaps with jets from the same decay, and lepton identification requires the use of isolation criteria.As a result, the signal selection efficiency is reduced in this kinematic region.It is important to recover the efficiency in this kinematic region, because the dominant decay of extra gauge bosons is to heavy neutrinos.Here we present a method that is sensitive to a broad range of N ℓ masses, including the region where m N ℓ is much smaller than m Z ′ .The analysis presented here uses a pp collision data set at √ s = 13 TeV collected with the CMS detector, corresponding to an integrated luminosity of 138 fb −1 .
Tabulated results are provided in the HEPData record for this analysis [29].

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections.Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors.Muons are detected in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid.A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [30].

Simulated samples
Signal and background events for each data-taking year from 2016 to 2018 are produced using Monte Carlo (MC) event generators and a detailed simulation of the CMS detector.The expected background processes include Drell-Yan (DY) dilepton plus jets production, and pro-duction of tt, dibosons, W+jets, single top quarks, and QCD multijet events.The DY events are simulated with MADGRAPH5 aMC@NLO 2.2.2 [31] with up to four jets at leading order (LO) in quantum chromodynamics (QCD) for 2016.For 2017 and 2018, MADGRAPH5 aMC@NLO 2.4.2 is used.The W events associated with jets are obtained with the MADGRAPH5 aMC@NLO 2.3.3 [31] generator with up to four jets at LO.The MLM scheme [32] is used to match partons from matrix element calculations and parton showers for these two samples.The tt production is simulated with the POWHEG v2 [33][34][35][36] generator at next-to-LO (NLO).Diboson (WW, WZ, and ZZ) events are simulated with the PYTHIA 8.212 [37] generator at LO.The t-channel and tW single top quark events are produced with the POWHEG v1 [33-35, 38, 39] generator.Single top quark production in the s-channel is simulated with the MADGRAPH5 aMC@NLO generator at NLO precision.The tt production in association with a W (Z) boson is simulated using MADGRAPH5 aMC@NLO at NLO (LO) precision.QCD multijet events are produced using MADGRAPH5 aMC@NLO at LO precision.
The NNPDF3.0 [40] parton distribution function (PDF) set is used to simulate background events corresponding to the data recorded in 2016.For the 2017 and 2018 backgrounds and all signal events, the NNPDF3.1 [41] PDF set is used.Parton showering, including photon radiation and hadronization processes, with the choice of NNPDF3.0 (or 3.1) are modeled using PYTHIA 8.226 (8.230) [37] with the CUETP8M1 (CP5) [42,43] tune.Finally, the response of the CMS detector is modeled using GEANT4 [44] for all simulated samples.The signal events for the process shown in Fig. 1 are simulated using MADGRAPH5 aMC@NLO 2.6.0 at NLO precision with the LRSM model cards [45].Heavy neutrinos with electron and muon flavors are considered, and no mixing between the two flavors is assumed.Both oppositeand same-sign configurations of dileptons are taken into account in the simulation.The signal samples are simulated for various m Z ′ and m N ℓ hypotheses, with m Z ′ between 400 and 5000 GeV, and m N ℓ between 100 GeV and m Z ′ /2.In the simulation, the W R mass is assumed to be 5 TeV, so that the heavy neutrinos are always lighter than the W R boson.The product of the NLO cross section and the branching fraction (BF) is the same for both dilepton channels and is 27.77 pb for m Z ′ = 400 GeV, and 8.60 × 10 −5 pb for m Z ′ = 5000 GeV with m N ℓ = 100 GeV.
With m N ℓ = 2400 GeV, it is 4.94 × 10 −7 pb for m Z ′ = 5000 GeV.A product of cross section and branching fraction of 1 fb corresponds typically to the Z ′ boson masses in the range 3 to 4 TeV, depending slightly on m N ℓ .
Minimum bias events generated with PYTHIA are superimposed on each simulated hard scattering event to reproduce the effect of extra pp interactions within the same or neighboring bunch crossings (pileup).The simulated events are weighted such that the distribution of the number of additional pileup interactions, estimated from the measured instantaneous luminosity for each bunch crossing, matches that observed in data sets individually for each year of data-taking in 2016-2018.The simulated events are processed with the same reconstruction software used for pp collision data.

Event reconstruction and object identification
Candidate events from decays of heavy neutrino pairs are selected using a two-tiered trigger system [46].The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed time interval of less than 4 µs.The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
For dielectron channel events, a diphoton trigger is used, which requires two electron or photon (e/γ) objects detected in the ECAL with a minimum transverse energy (E T ) of 60 GeV for 2016 data or 70 GeV for 2017 and 2018 data.Dimuon channel events must pass a single-muon trigger that requires the presence of a muon with a minimum p T of 50 GeV.
In the dielectron channel, the trigger efficiency of signal events is about 85%, dropping to about 10% in the kinematic region where m N ℓ is much smaller than m Z ′ .The reduction of efficiency is mainly caused by the requirement on H/E, the ratio between energies deposited in HCAL and ECAL, being smaller than 0.15 (barrel) or 0.10 (endcap).Since the requirement on H/E is essential for distinguishing electrons and photons from hadrons, there are inevitable signal losses in the dielectron channel for highly Lorentz-boosted events.This analysis uses a doublephoton trigger, since it has the loosest H/E requirement among all e/γ HLT paths.
The trigger efficiency for signal events in the dimuon channel depends on m N ℓ and m Z ′ , and ranges between 80 and 95%, except at the lightest Z ′ boson mass points, where it is somewhat lower because of the high momentum threshold of the muon trigger.For example an efficiency of about 60% is observed for the mass points at the 400 GeV and 600 GeV.The cross sections for these mass points are higher than those for heavier Z ′ bosons.This allows good sensitivity to be maintained while keeping the same trigger for all masses.
The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [47].
Particles in selected events are reconstructed and identified using the global event reconstruction (also called particle-flow event reconstruction [48]) with an optimized combination of all subdetector information.The identification of the particle type (photon, electron, muon, charged hadron, neutral hadron) plays an important role in the determination of the particle direction and energy.Photons (e.g., coming from π 0 decays or from electron bremsstrahlung) are identified as ECAL energy clusters not linked to the extrapolation of any charged particle trajectory to the ECAL.Electrons (e.g., from Z boson leptonic decays) are identified as a primary charged particle track and potentially many ECAL energy clusters corresponding to this track extrapolation to the ECAL and to possible bremsstrahlung photons emitted while traversing the tracker material.Muons (e.g., from W boson leptonic decays) are identified as tracks in the central tracker matching with either a track or several hits in the muon system, and associated with calorimeter deposits compatible with the muon hypothesis.Charged hadrons are identified as charged particle tracks neither identified as electrons, nor as muons.Finally, neutral hadrons are identified as HCAL energy clusters not linked to any charged hadron trajectory, or as a combined ECAL and HCAL energy excess with respect to the expected charged hadron energy deposit.
The energy of photons is obtained from the ECAL measurement.The energy of electrons is determined from a combination of the track momentum at the PV, the corresponding ECAL cluster energy, and the energy sum of all bremsstrahlung photons associated with the track.The energy of muons is obtained from the corresponding track momentum.The energy of charged hadrons is determined from a combination of the track momentum and the corresponding ECAL and HCAL energies, corrected for the response function of the calorimeters to hadronic showers.Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.

Lepton selection
Electrons are selected in the region of |η C | < 2.5, where η C is the pseudorapidity of the ECAL cluster with respect to the nominal center of the CMS detector.The transition region between ECAL barrel and endcaps, 1.44 < |η C | < 1.57, is excluded.The electron E T is required to be greater than 65 (75) GeV for the 2016 (2017 and 2018) data.In this analysis, two electron identification criteria are used: "loose" and "tight".A "loose" electron is an object passing the "HEEP" selection or passing the selection of a loose cut-based electron with no isolation criteria [49,50].The "HEEP" criteria have been developed for the analysis of high-energy electrons.A "tight" electron is an object passing the "HEEP" requirement.
In addition, a "veto" electron is defined using the same criteria as the loose electron but lowering the E T requirement to 10 GeV.
Muons are selected in the region |η| < 2.4 with p T > 65 (75) GeV for 2016 (2017 and 2018).Higher p T requirements for muons compared to the trigger thresholds are applied to reduce QCD background.A dedicated muon identification algorithm, HighPtID [51], adapted from searches for resonances using high-momentum leptons [52] is used for this analysis.Two muon identification criteria are used, "loose" and "tight".A "loose" muon is required to satisfy the HighPtID requirements.A "tight" muon must pass both the HighPtID and isolation requirements.The isolation is defined as the p T sum of tracks within a cone of radius ∆R = √ (∆η) 2 + (∆ϕ) 2 = 0.3 around the muon candidate direction, where ϕ is its azimuthal angle in radians.The momentum of the muon candidate is excluded from the sum.If the sum is less than 10% of the muon candidate's p T , it passes the isolation requirement.In addition, a "veto" muon is defined using the particle-flow muons with p T > 10 GeV.

Jet selection
For each event, hadronic jets are reconstructed by clustering particle-flow objects using the antik T algorithm [53,54] with two distance parameters of 0.4 and 0.8, producing AK4 and AK8 jets, respectively.The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from the simulation to be within 5 to 10% of the true momentum, on average, over the entire p T spectrum and detector acceptance.Pileup interactions can contribute extra tracks and calorimetric energy depositions to the jet momentum.
In order to mitigate the effect of pileup, the pileup-per-particle identification algorithm [55,56] is used at the reconstructed-particle level, making use of local shape information, event pileup properties and tracking information.Charged particles identified to be originating from pileup vertices are discarded.For each neutral particle, a local shape variable is computed using the surrounding charged particles compatible with the PV within the tracker acceptance (|η| < 2.5), and using both charged and neutral particles in the region outside of the tracker coverage.The momenta of the neutral particles are then rescaled according to their probability to originate from the PV deduced from the local shape variables, superseding the need for jet-based pileup corrections [57].For this analysis, AK4 (AK8) jets are required to have p T > 40 (300) GeV and |η| < 2.7.
AK8 jets are groomed using the soft-drop algorithm [58,59].In this algorithm, the constituents of the AK8 jets are reclustered using the Cambridge-Aachen algorithm [60,61].Soft radiation and wide-angle radiation from the jet are removed by setting the angular exponent β to zero, the soft cutoff threshold to 0.1, and the characteristic radius R 0 to 0.8 [62].All AK8 jets are required to have mass greater than 40 GeV to reduce the number of AK8 jets that originate from a single parton.
In order to avoid double counting of AK4 jets with a lepton or an AK8 jet, any AK4 jet located within ∆R < 0.4 of a loose lepton or within ∆R < 1.0 of an AK8 jet is not used in this analysis.

Event selections
The signature targeted in this search is two same-flavor leptons and four jets from heavy neutrino decays, as shown in Fig. 1.The kinematic distributions of the final state are strongly dependent on the ratio of twice the heavy neutrino mass and the Z ′ boson mass (2m N ℓ /m Z ′ ).If the ratio is much smaller than one, the lepton and two jets from a heavy neutrino decay are merged into one jet because of the high Lorentz boost provided by the heavy neutrino.As a consequence, leptons are mostly reconstructed as loose leptons since they seldom satisfy the isolation criteria.In contrast, leptons and jets are well separated if the ratio is close to unity.A dedicated strategy is developed to achieve a good signal selection efficiency in the wide kinematic region covered by this search.
The selected events are required to have exactly two same-flavor (ee or µµ) loose leptons using both the opposite-and same-sign configurations.Events with additional veto leptons are removed.To suppress the background contribution from DY events, the mass of the dilepton pair (m ℓℓ ) is required to be greater than 150 GeV.
No b-tagging veto criteria are considered in the analysis, since high mass heavy neutrinos have a significant branching fraction into bottom quark jets, via N ℓ → ℓ − tb and N ℓ → ℓ + tb decays.Furthermore, for lower heavy neutrino masses, final states typically contain jets with very high energies, for which b-tagging performance is degraded and, as a result, systematic uncertainties prevent an improvement in sensitivity.
We define three different signal regions (SRs), depending on the number of AK8 jets found in the event, as shown in Table 1.The region SR1 contains no AK8 jets, SR2 contains exactly one AK8 jet, and SR3 contains at least two AK8 jets.
In SR1, we require two tight leptons and at least four AK4 jets, but no AK8 jets.To reconstruct two heavy neutrinos for an event, the two tight leptons and the leading four AK4 jets are used.Each heavy neutrino candidate is reconstructed using one lepton and two AK4 jets.From the 12 possible lepton plus two-jet combination combinations, we select the one that minimizes the difference between the masses of the two heavy neutrino candidates.
In SR2, we require exactly one AK8 jet, at least one tight lepton and at least two AK4 jets.
In events with only one tight lepton, this lepton is paired with the two leading AK4 jets to reconstruct one heavy neutrino.The AK8 jet is assigned as a proxy to the other heavy neutrino.If an event has two tight leptons, the AK8 jet is paired with the closer lepton, and the other lepton is paired with the two leading AK4 jets.
In SR3, we require at least two AK8 jets, and loose leptons are ignored.In cases where an event has no tight leptons, the two leading AK8 jets are assigned as proxies to the two heavy neutrinos.If there is only one tight lepton, the lepton is paired with the closer AK8 jet out of the two leading AK8 jets to reconstruct one heavy neutrino.The remaining AK8 jet is assigned as a proxy to the other heavy neutrino.If there are two tight leptons, the leading AK8 jet is paired to the closer lepton while the secondary AK8 jet and the other lepton are used to reconstruct the other heavy neutrino.
In all three SRs, the invariant mass of the Z ′ boson is obtained using the two reconstructed heavy neutrinos.If a tight lepton is included in the energy of the AK8 jet, the four-momentum of the AK8 jet is assigned to the four-momentum of the heavy neutrino.For m N ℓ = 100 GeV and m Z ′ = 4 TeV, about 2% of the signal events that enter the SR3 are of this type.Once the Z ′ boson and heavy neutrinos are reconstructed, additional requirements (reconstructed m N ℓ > 80 GeV and reconstructed m Z ′ > 300 GeV) are applied to reduce backgrounds, and we start search range from m Z ′ = 400 GeV and m N ℓ = 100 GeV.
The products of acceptance and efficiency for the simulated signal samples are 20-25 (25-33)% and 1.2-4.7 (7-16)% in the dielectron (dimuon) channel for the most resolved and boosted mass points, respectively, depending on m Z ′ .

Background estimation
The backgrounds are dominated by SM processes containing prompt leptons.These are estimated using the simulated samples described in Section 3. The leading background contribution is from tt events.The simulated tt sample is normalized using the QCD next-to-NLO cross section, which includes resummation of next-to-next-to-leading logarithmic soft gluon terms for the top quark mass of 172.5 GeV [63].A control region (CR), CR1, is defined to validate the tt background using an eµ sideband passing the same single-muon triggers used in the dimuon SR.The region CR1 has exactly the same event selection as the SR, except the flavors of the two leptons are required to be different.The events of CR1 are classified in three orthogonal regions according to the same criteria used for the SR, which are summarized in Table 1.
Since no mixing between generations of heavy neutrinos is assumed and the heavy neutrino always decays into a lepton and two jets in the signal process, the signal in CR1 is negligible.This expectation is supported by the observed consistency between simulated background predictions in CR1 and the corresponding data-driven estimates.Additional scale factors are applied to the tt simulated samples, depending on the AK8 jet multiplicity.The scale factors are extracted from a simultaneous fit to all CRs and SRs and range from 0.83 to 1.06.The uncertainties in these scale factors, determined from the fit, are assumed to cover the theoretical uncertainties in the simulated tt distributions.
Because of the large production cross section, DY events constitute the second largest background, although most of these are removed by the requirement, m ℓℓ > 150 GeV.A second CR (CR2) is chosen to verify that the simulated sample of DY events is well modeled.The selection used to define CR2 is similar to that of the SR except the requirement on m ℓℓ is changed to |m ℓℓ − m Z | < 10 GeV, where m Z is the mass of Z boson.The events of CR2 are classified in three orthogonal regions according to the same criteria used for the SR, which are summarized in Table 1.To correct for the mismodeling of the dilepton p T in simulated events, a dedicated correction estimated using inclusive dimuon events is applied as a function of generator level dilepton p T and mass [64].The overall normalization factors as functions of the AK8 jet multiplicity are fit using all the CRs and SRs.They range from 0.93 to 1.14 for events with fewer than two AK8 jets.For events with two or more AK8 jets, normalization factors have a wider range from 0.87 to 1.28 owing to the small numbers of events.Figure 3 shows the agreement between .The fitting procedure is described in Section 7. The last bin of each plot includes the overflow events.The contribution from DY events is too small to be visible on these plots.The lower panel of each plot shows the ratio of observed events to expected background.the observed and simulated distributions of the reconstructed Z ′ candidate mass in CR2.
Additional prompt backgrounds are considered from di-and triboson production, tt pair production in association with a vector boson, and single top production.The contribution of these backgrounds is found to be almost negligible.We assign a conservative 50% uncertainty in the estimated size of these backgrounds, to account for limited MC samples in the extreme phase space regions of interest.The background contribution with nonprompt leptons, which are leptons not originating from the hard scattering process, is estimated using the simulated samples.We assign an uncertainty of 100% in the estimated size of the nonprompt-lepton background since this depends on detector effects, which are often more difficult to estimate.

Systematic uncertainties and fitting procedure
Binned maximum likelihood fits using the distribution of reconstructed m Z ′ for background and signal as templates are performed to extract the signal contribution for each SR.The choice of bin size is based on the resolution for the reconstructed m Z ′ , which ranges from 150 to 550 GeV.Where necessary, bins are then merged to ensure that the predicted number of MC events is nonzero in all bins of all SRs and CRs for each year.Various sources of uncertainty in both shape and normalization of reconstructed m Z ′ distributions are considered as nuisance parameters and profiled in the fit.Log-normal and Gaussian prior probability densities are used for normalization and shape uncertainties, respectively.The normalizations of tt and DY backgrounds are left freely floating in the fit.The pre-fit and post-fit predictions on the plots of Figs.2-4 correspond to the results obtained before and after likelihood fitting, respectively.No nuisance parameter is pulled beyond one standard deviation from the pre-fit value.
For leptons, the efficiencies for reconstruction, identification, isolation, and trigger are estimated using observed and simulated Z → ee, µµ events that are not included in the SRs.Scale factors (SFs) are applied to simulated events to correct for mismodeling in the simulation.Systematic uncertainties in momentum scale and resolution [49,51] are estimated in simulations by varying the lepton momentum scales within their uncertainties.
The jet energy corrections and resolutions estimated using hadronic jets are used for both hadronic jets and AK8 jets with high p T leptons inside.Such corrections were found to provide a good modeling of jets also in the case of non-resolved lepton(s) within the AK8 jet in CRs.The uncertainty in the jet energy scale and resolution are estimated using the same procedures as those used to determine systematic uncertainty in lepton momenta.An additional uncertainty coming from the soft-drop mass scale is found to be less than 1.3% [65].
The integrated luminosities of the 2016, 2017, and 2018 data-taking periods are individually known with uncertainties in the range 1.2-2.5% [66][67][68], while the total 2016-2018 integrated luminosity has an uncertainty of 1.6%, the improvement in precision reflecting the (uncorrelated) time evolution of some systematic effects.
The simulation of pileup events is based on a cross section of 69.2 mb for inelastic pp collisions with energy deposits inside the CMS acceptance.A systematic uncertainty for the pileup weight is assigned by varying the cross section up and down by 5%.An additional event weight that accounts for the L1 trigger inefficiency observed in the |η| < 2.0 region for the 2016 and 2017 data-taking period is also applied [69].The systematic uncertainty for the L1 trigger inefficiency weight is assigned by varying the weight within its uncertainty.
The uncertainties related to the choice of PDF can affect both the cross section and kinematic distributions of the signal.The effect of the PDF uncertainty and the variation re-  lated to the strong coupling constant (α S ) are considered, following the method introduced by PDF4LHC15 [70].The uncertainties in the renormalization and factorization scales (µ R and µ F , respectively) are determined by varying them by a factor of two on an event-by-event basis using combinations in which the two scales do not move in opposite directions at the same time.The effect of these variations on the shape of the m Z ′ distribution is determined while maintaining the signal normalization.
The normalization uncertainties for the rare SM processes and nonprompt-lepton backgrounds are also incorporated into the likelihood as nuisance parameters.We observe that the fit procedure returns an uncertainty of about 40% in the nuisance parameter of nonprompt-lepton background, which is initially constrained at the 100% level.Since this is consistent with the precision of independent background estimation methods using control samples in data-driven methods [15,16,71], we use the fit to constrain the normalization.
The uncertainty in the dilepton p T shape corrections to the simulated DY events is estimated by varying correction factors within their uncertainties.
A complete list of these uncertainties in the dielectron and dimuon channels is shown in Table 2.The lepton trigger scale factors, jet energy resolution, L1 trigger inefficiency, and nonprompt lepton background normalizations are assumed to be uncorrelated between the 2016, 2017, and 2018 data sets, since their corrections and corresponding uncertainties are affected by differences in operating conditions and detector performance between data-taking years.
Other experimental and all theoretical uncertainties are taken as correlated.

Results and discussion
The background-only post-fit invariant mass distributions of Z ′ candidates in the different SRs are shown in Fig. 4. Representative signal distributions are also shown.
We do not observe any significant excess in data with respect to the SM prediction.The 95% confidence level (CL) upper limits on the product of signal cross section and BF are obtained us-ing the distributions of the likelihood ratio calculated with the asymptotic approximation [72] and the CL s criterion [73,74].
The observed and expected exclusion limits on the signal cross section times BF are shown in Fig. 5.The maximum local significance is 3.32σ, observed in the dielectron channel for the (m Z ′ , m N ℓ ) = (4.6,0.1) TeV.The probability to observe similar or larger excess in the dielectron channel across the full analysis mass range is estimated with pseudo experiments under background-only hypothesis.This probability is 1.12 × 10 −2 , corresponding to a global significance of 2.28σ.For cases where m N ℓ = m Z ′ /4, the observed (expected) lower limits at 95% CL on the mass of the Z ′ boson are 3.59 (3.90) TeV in the dielectron channel and 4.10 (3.86) TeV in the dimuon channel.For the phase space with mostly boosted signals (m N ℓ = 100 GeV), the observed (expected) lower limits are found to be m Z ′ = 2.79 (3.12) TeV and 4.38 (4.22) TeV in the dielectron and dimuon channels, respectively.The use of a dedicated algorithm for boosted signals provides a significant improvement in sensitivity, in particular for the muon channel.We note that the sensitivity for Dirac-type heavy neutrinos is identical as long as the branching fraction of the Z ′ boson to a heavy neutrino pair stays the same.
One of the assumptions in this analysis is the presence of a nonnegligible Z ′ coupling to quarks.A comparison with the sensitivity of the search for high mass dijet resonances performed by CMS using the same data set [75] is therefore possible.For comparable Z ′ branching fractions into quarks and heavy neutrinos, our search is found to be about one order of magnitude more sensitive.
The upper limits across the m N ℓ -m Z ′ plane are shown in Fig. 6.The SR1 is sensitive in the phase space regions m Z ′ = 1 to 2 TeV and m N ℓ = 0.5 to 1 TeV.The SR2 is particularly sensitive for higher m Z ′ with resolved heavy neutrinos.The SR3 is most important where heavy neutrinos are boosted.The dimuon channel shows better sensitivity where heavy neutrinos are mostly boosted because the dielectron channel suffers for the H/E requirements discussed in Section 4. In the dielectron channel, for m N ℓ > 0.9 TeV and m Z ′ > 2.85 TeV, the observed limits obtained in the SR2 are better than those of the combined limits.This indicates that the slight excess observed for the combined limits in this region of phase space is driven by the SR3.
A priori, a further categorization separating same-sign and opposite-sign dilepton pair events might increase the sensitivity.However, with the current data sample size, the evaluation of same-sign dilepton backgrounds using control samples, particularly in a highly boosted regime for heavy neutrinos, is affected by large systematic uncertainties.Thus introducing sign categorization would not improve the present search.

Summary
A search for the pair production of heavy Majorana neutrinos (N ℓ ) via the decay of a Z ′ boson, in a final state with two same-flavor leptons and at least two reconstructed jets has been performed in proton-proton collisions at a center-of-mass energy of 13 TeV, using LHC 2016-2018 data corresponding to an integrated luminosity of 138 fb −1 .No significant excess of events beyond the expected background is observed.Upper limits on the product of signal cross section and branching fraction in the context of a left-right symmetry model scenario are set [45].Exclusion regions in the dielectron and dimuon channels are set at 95% confidence level.These are the first results of a search for this process at 13 TeV and are the most restrictive limits to date on the mass of N ℓ as a function of the Z ′ boson mass. ) = (1000, 100) GeV, (1000, 300) GeV, (4000, 200) GeV, and (4000, 1200) GeV are shown together with a reference cross section times branching fraction of 1 fb.The last bin of each plot includes the overflow events.The lower panel of each plot shows the ratio of observed events to expected background.= m Z ′ /4 (upper row) and m N ℓ = 100 GeV (lower row) for dielectron (left column) and dimuon (right column) channels.The green and yellow bands indicate the 68% and 95% CL regions around the expected limit.The red lines represent values coming from the benchmark LRSM model [45].Both opposite-and same-sign dileptons from decays of heavy neutrino pairs are considered.for dielectron (left) and dimuon (right) channels.One standard deviation (s.d.) limits are also shown.Both opposite-and same-sign dileptons from decays of heavy neutrino pairs are considered.
[28] ATLAS Collaboration, "Search for heavy majorana or dirac neutrinos and right-handed W gauge bosons in final states with charged leptons and jets in pp collisions at √ s = 13 TeV with the ATLAS detector", 2023.arXiv:2304.09553.Submitted to Eur.Phys.J. C.

Figure 1 :
Figure 1: Feynman diagram representing the pair production of N ℓ via Z ′ boson exchange, where ℓ = e or µ.Each N ℓ decays into a lepton and two quarks, and we assume that the W ± R is heavier than the N ℓ .Both opposite-and same-sign dileptons from the decays of the heavy neutrino pair are allowed, because of the Majorana nature of N ℓ .

Figure 2 :
Figure2: Reconstructed mass of the Z ′ candidate in CR1 (eµ), which consists of flavor sidebands of SR1 (upper), SR2 (middle), and SR3 (lower) regions.Pre-fit (post-fit) results are shown on the left (right).The fitting procedure is described in Section 7. The last bin of each plot includes the overflow events.The contribution from DY events is too small to be visible on these plots.The lower panel of each plot shows the ratio of observed events to expected background.

Figure 3 :
Figure3: Reconstructed mass of the Z ′ candidate in CR2 (ee and µµ), which consists of the dilepton mass sidebands of SR1 (upper), SR2 (middle), and SR3 (lower) regions.Post-fit dielectron (dimuon) channel results are shown on the left (right).The fitting procedure is described in Section 7. The last bin of each plot includes the overflow events.The lower panel of each plot shows the ratio of observed events to expected background.

Figure 4 :
Figure 4: Distributions of the reconstructed Z ′ candidate mass in SR1 (upper row), SR2 (middle row) and SR3 (lower row).The left (right) column corresponds to the dielectron (dimuon) channel.Signal samples with (m Z ′ , m N ℓ) = (1000, 100) GeV, (1000, 300) GeV, (4000, 200) GeV, and (4000, 1200) GeV are shown together with a reference cross section times branching fraction of 1 fb.The last bin of each plot includes the overflow events.The lower panel of each plot shows the ratio of observed events to expected background.

Figure 5 :
Figure 5: The observed and expected 95% CL upper limits on the product of Z ′ boson cross section and branching fraction are shown for the case of m N ℓ= m Z ′ /4 (upper row) and m N ℓ = 100 GeV (lower row) for dielectron (left column) and dimuon (right column) channels.The green and yellow bands indicate the 68% and 95% CL regions around the expected limit.The red lines represent values coming from the benchmark LRSM model[45].Both opposite-and same-sign dileptons from decays of heavy neutrino pairs are considered.

Figure 6 :
Figure 6: Observed and expected exclusion regions at 95% CL in the 2D phase space of m Z ′ vs. m N ℓfor dielectron (left) and dimuon (right) channels.One standard deviation (s.d.) limits are also shown.Both opposite-and same-sign dileptons from decays of heavy neutrino pairs are considered.

Table 1 :
Multiplicity requirements for AK8 jets, tight leptons and AK4 jets in the different signal regions considered in the search.

Table 2 :
Systematic uncertainties and their impacts on the total number of events in the signal regions before the binned maximum likelihood fit.