Search for events with a pair of displaced vertices from long-lived neutral particles decaying into hadronic jets in the ATLAS muon spectrometer in pp collisions at $\sqrt{s}=13$ TeV

A search for events with two displaced vertices from long-lived particles (LLP) pairs using data collected by the ATLAS detector at the LHC is presented. This analysis uses 139~fb$^{-1}$ of proton-proton collision data at $\sqrt{s}=13$ TeV recorded in 2015-2018. The search employs techniques for reconstructing vertices of LLPs decaying to jets in the muon spectrometer displaced between 3 m and 14 m with respect to the primary interaction vertex. The observed numbers of events are consistent with the expected background and limits for several benchmark signals are determined. For the Higgs boson with a mass of 125 GeV, the paper reports the first exclusion limits for branching fractions into neutral long-lived particles below 0.1%, while branching fractions above 10% are excluded at 95% confidence level for LLP proper lifetimes ranging from 4 cm to 72.4 m. In addition, the paper present the first results for the decay of LLPs into into $t\bar{t}$ in the ATLAS muon spectrometer.


Introduction
The discovery of the Higgs boson at the LHC completed the Standard Model (SM) of elementary particles and focused attention on the many central features of our universe that the SM does not address: dark matter, neutrino mass, matter-antimatter asymmetry, and the hierarchy problem (naturalness). Many beyond-the-SM (BSM) theoretical constructs that address these phenomena predict the existence of long-lived particles (LLPs), with macroscopic decay lengths.
This search focuses on DVs occurring in the ATLAS muon spectrometer (MS) and uses 139 fb −1 of 13 TeV collision data. The event selection criteria and vertex reconstruction algorithms choose candidate events with two MS DVs. The results are interpreted in terms of a scalar portal model, but the analysis may be sensitive to additional models. These results supersede those of a previous search [38] for events with two MS DVs in a smaller 13 TeV dataset. Furthermore, the current analysis uses an improved background estimation methodology, more advanced modeling of the trigger and vertex reconstruction in Monte Carlo (MC) simulation, as well as a refined signal efficiency extrapolation as a function of the LLP lifetime.
The paper describes the ATLAS detector in Section 2, followed by an overview of the analysis strategy in Section 3 and the theoretical models in Section 4. Details of the collision data and MC simulation are provided in Section 5, while those of the specialized trigger and the purpose-built displaced vertex reconstruction algorithm follow in Section 6. The baseline event selection is described in Section 7.1, while Section 7.2 describes the signal selection and Section 8 describes the background estimation. The systematic uncertainties are discussed in Section 9, and the results are presented in Section 10.

ATLAS detector
The ATLAS detector [44], which has nearly 4 steradian coverage, 1 is a multipurpose detector consisting of an inner tracking detector (ID) surrounded by a superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer (MS) based on three large air-core toroidal superconducting magnets, each with eight coils. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.
The ID covers the pseudorapidity range | | < 2.5. It consists of a silicon pixel detector, a silicon microstrip detector, and a straw-tube transition-radiation tracker.
The calorimeter system covers the pseudorapidity range | | < 4.9. It consists of a high-granularity electromagnetic calorimeter (ECal) surrounded by a hadronic calorimeter (HCal). Within the region | | < 3.2, the ECal comprises liquid-argon (LAr) barrel and endcap electromagnetic calorimeters with lead absorbers. An additional thin LAr presampler covering | | < 1.8 is used to correct for energy loss in material upstream of the calorimeters. The ECal extends from 1.5 m to 2.0 m in in the barrel and from 3.6 m to 4.25 m in | | in the endcaps. The HCal is a steel/scintillator-tile calorimeter that is segmented into three barrel structures within | | < 1.7, and two copper/LAr hadronic calorimeters in the endcaps (1.5 < | | < 3.2). The HCal covers the region from 2.25 m to 4.25 m in in the barrel (although the HCal active material extends only up to 3.9 m) and from 4.3 m to 6.05 m in | | in the endcaps. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimized for electromagnetic and hadronic measurements, respectively. Together the ECal and HCal have a thickness of 9.7 interaction lengths at = 0. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal IP in the center of the detector and the -axis along the beam pipe. The -axis points from the IP to the center of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, where is the azimuth angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of Δ ≡ √︁ (Δ ) 2 + (Δ ) 2 .
The MS comprises three stations of separate trigger and tracking chambers that measure the deflection of muons in a magnetic field generated by the air-core toroid magnets. The barrel and endcap chamber systems are subdivided into 16 sectors: 8 large sectors and 8 small sectors. For the barrel, the small sectors are inside each of the 8 magnet coils and the large sectors are between the coils. The endcap TGC chambers are arranged into 24 or 48 sectors depending on their and radial position. Three stations of resistive-plate chambers (RPC) and thin-gap chambers (TGC) are used for triggering in the MS barrel and endcaps, respectively. The first two RPC stations, which are radially separated by 0.5 m, begin at a radius of either 7 m (large sectors) or 8 m (small sectors). The third station is located at a radius of either 9 m (large sectors) or 10 m (small sectors). In the endcaps, the first TGC station is located at | | = 13 m. The other two stations start at | | = 14 m and | | = 14.5 m, respectively. The muon trigger system covers the range | | < 2.4. The muon tracking chamber system covers the region | | < 2.7 with three stations of monitored drift tubes (MDT), complemented by cathode-strip chambers (CSC) in the forward region. The MDT chambers consist of two multilayers separated by a distance ranging from 6.5 mm to 317 mm. Each multilayer consists of three or four layers of drift tubes. The individual drift tubes are 30 mm in diameter and have lengths of 2-5 m (barrel) and 2-6.5 m (endcaps) depending on the location of the chamber in the spectrometer. In each multilayer, charged-particle track-segment reconstruction entails finding the line that is tangent to the drift circles. These single-multilayer segments are local measurements of the position and direction of the charged particle. Because of its design, the MDT measurement provides only a very coarse position of the track hit. In order to reconstruct the position and direction, the MDT measurements are combined with the -coordinate measurements from the trigger chambers.
The ATLAS trigger and data acquisition system [45] consists of a hardware-based first-level (L1) trigger followed by a software-based high-level trigger (HLT) that reduces the rate of events selected for offline storage to 1 kHz.
The implementation of the L1 muon trigger logic is similar for the RPC and TGC systems. Each of the three stations of the RPC system and the two outermost stations of the TGC system consist of a doublet of independent detector layers. The first TGC station contains a triplet of detector layers. The transverse momentum ( T ) of the muon candidate is measured by the L1 muon trigger, using different algorithms for low-T and high-T triggers. In the barrel, a low-T (<10 GeV) muon region-of-interest (RoI) is generated by requiring a coincidence of hits in at least three of the four layers of the two inner RPC stations. In the endcaps, the trigger requires hits in the two outer TGC stations. A high-T muon RoI in the barrel requires additional hits in at least one of the two layers of the outer RPC station, while for the endcaps, additional hits in two of the three layers of the innermost TGC station are required. The muon RoIs have an angular extent of 0.2 × 0.2 in Δ × Δ in the MS barrel and 0.1 × 0.1 in Δ × Δ in the MS endcaps.
An extensive software suite [46] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.
The main background to LLPs decaying into hadronic jets in the MS originates from hadronic or electromagnetic showers not fully contained in the calorimeter volume (punch-through jets) that result in tracks reconstructed in the MS. Multĳet events that result in vertices in the MS will often have ID tracks and calorimeter jets that point towards the MS DV. To reduce the acceptance of such fake vertices from multĳet events, vertices are required to be isolated in and from ID tracks and calorimeter jets. Multĳet events with mismeasured jets are another source of background as they will not be removed by the isolation criteria.
Additional background, referred to as noncollision background, can be generated by electronic noise in the MDT and RPC/TGC chambers, by cosmic-ray muons, and by the LHC beam (beam-induced background [49]). This last contribution is composed of hadronic and electromagnetic showers caused by beam protons interacting with collimators or residual gas molecules inside the vacuum pipe.
To avoid unintended biasing of the results during the development of the analysis, the number of events observed in the signal region was not measured until the entire set of selection criteria were optimized and all the systematic uncertainties completely evaluated.

Signal model
Although this analysis is sensitive to a large variety of models, this paper interprets the results in terms of a scalar portal model [13], where a Higgs boson with a mass of 125 GeV or a lower/higher-mass scalar boson, both indicated by the symbol , decays into two long-lived scalars.
A standard way of introducing long-lived, neutral particles to the SM is through a hidden sector that weakly couples to the SM. Scalar portals, where a 125 GeV Higgs boson or a lower/higher-mass scalar boson mixes weakly with a hidden-sector scalar field, can result in pair production of long-lived hidden-sector scalars or pseudoscalars that carry no SM quantum numbers. An indirect upper limit on the BSM Higgs boson decays of the order of 20% can be derived from combined studies of the Higgs boson production and decay [50,51] under several assumptions, potentially allowing sizable branching fractions for decays into non-SM particles.
Long-lived scalars with a mass from 5 to 30 GeV are also well motivated by naturalness models that are generic extensions of hidden-valley portal models [52].
The mechanism for LLP production is shown in Figure 1. Here, a scalar boson decays with some effective coupling into a pair of long-lived scalars, . The scalars subsequently decay into SM particles. In this model the couplings of the scalar to SM fermions are determined by the Higgs Yukawa coupling. Therefore, each long-lived scalar decays mainly into the heaviest fermion pair that is kinematically accessible. The branching fractions of , , and + − decays depend on the mass of the scalar ( ), and for 25 GeV they are almost independent of its mass and equal to 85%, 5%, and 8%, respectively. While only these decay modes were considered in the previous ATLAS publication [38], this paper also considers the decay into as it is the dominant decay mode when it becomes kinematically available.
The branching fraction for decaying into a pair of long-lived scalars is not constrained in these models. It is therefore interesting to focus on decays into LLPs where is a Higgs boson with a mass of 125 GeV or another scalar with a different mass. The relative masses and lifetimes generated, as well as details about the Monte Carlo (MC) event generation are described in Section 5. Only the gluon-gluon fusion production mode is considered, as it is the dominant production mode.

Data and MC simulation
The analysis presented in this paper uses √ = 13 TeV collision data recorded by the ATLAS detector with stable LHC beams during the 2015-2018 data-taking periods [53]. Events are required to have been taken during stable beam conditions and when all the detector subsystems were operational. After data quality requirements [53], the total integrated luminosity is 139 fb −1 .
An unbiased random trigger (the zero-bias trigger) that fires on the bunch crossing occurring one LHC revolution after a low-threshold calorimeter-based trigger was used to collect data for background estimation with negligible signal contamination. The zero-bias trigger runs throughout the ATLAS data-taking, so the events are acquired with the same beam conditions present in normal physics data-taking and can be used to study the expected background. The zero-bias trigger was prescaled as a function of the instantaneous pileup to have a fixed output event rate throughout the data-taking period. Therefore, to estimate the correct number of background events, the pileup profile from the zero-bias trigger was rescaled to match that from the data used to select candidate events.
Signal MC simulation samples were produced using the Hidden Abelian Higgs Model [54], considering only the gluon-gluon fusion production mechanism. The masses, summarized in Table 1, were chosen to span the accessible parameter space. The mean proper lifetime in each sample was tuned to maximize the number of decays throughout the ATLAS detector volume and covers a range of 0.13-6 m, depending on the sample. All samples were generated at leading order using M G 5_ MC@NLO [55] interfaced to the P 8 [56] parton shower model. The transverse momentum distribution for the samples is reweighted to match that obtained for the corresponding next-to-leading order (NLO) Higgs samples generated using the M G 5_ MC@NLO merging approach [57]. A set of tuned parameters called the A14 tune [58] was used together with the NNPDF3.1 parton distribution function (PDF) set [59]. The E G 1.2.0 program [60] was used to model -and -hadron decays. The generated events were processed through a full simulation of the ATLAS detector geometry and response [61] using the G 4 [62] toolkit. The simulation includes multiple interactions per bunch crossing (pileup), as well as the effect on the detector response due to interactions from bunch crossings before or after the one containing the hard interaction. Pileup was simulated with the soft strong-interaction processes of P 8 using the A3 tune [63] and the NNPDF2.3 [59] PDF set. Per-event weights were applied to the simulated events to match the distribution of the average number of interactions per bunch crossing measured in data.
A multĳet background sample was generated using P 8.230 [64] with leading-order matrix elements for dĳet production which were matched to the parton shower. The NNPDF2.3 PDF set was used in the matrix element generation, the parton shower, and the simulation of the multiparton interactions. The A14 tune was used. Perturbative uncertainties were estimated through event weights [65] that encompass variations of the scales at which the strong coupling constant is evaluated in the initial-and final-state showers, as well as the PDF uncertainty in the shower and the nonsingular part of the splitting functions. Since the analysis is sensitive to a wide range of mean proper lifetimes, and the generation of many samples to cover a broad lifetime range would consume far too much CPU time, a toy MC method was adopted to extrapolate the expected number of events to the range of mean proper lifetimes between 0.01 and 1000 m. For each LLP in the MC sample, a random decay time, sampled from an exponential distribution of a chosen proper lifetime, was generated. The physical decay distance in the detector was then calculated for each simulated LLP using its four-momentum. The overall probability of the event to satisfy the signal selection criteria, parameterized as a function of the LLP decay position and boost, is then evaluated from the LLP trigger and vertex efficiencies described in Sections 6.1 and 7.3, respectively. In order to validate the lifetime-extrapolation procedure described above, a second set of samples with a different lifetime was generated for some scalar masses.

Trigger selection and event reconstruction
Hadronic LLP decays in the MS typically produce narrow, high-multiplicity hadronic showers. The track multiplicity and shower width vary with the mass and boost of the decaying LLP and the final states to which the LLP decays. Dedicated trigger [47] and vertex [48] algorithms were developed to select and reconstruct displaced decays in the MS. Due to the amount of material in the calorimeter, only decays occurring in or after the last sampling layer of the hadronic calorimeter produce enough hits in the MS to allow reconstruction of the DV.

Muon RoI Cluster trigger
The Muon RoI Cluster trigger is a signature-driven trigger that selects candidate events for decays of LLPs in the MS [47]. Events must contain at L1 two muon RoIs with a T higher than 10 GeV, and in the HLT a cluster of three (four) muon RoIs within a Δ = 0.4 cone centered on the L1 object in the barrel (endcaps). The isolation criteria for jets and tracks discussed in Ref.
[47], used to reduce background DVs from punch-through jets, are not considered in the analysis presented in this paper. The trigger selects both isolated, signal-like events and nonisolated, background-like events used for data-driven background estimations.
The trigger efficiency, defined as the fraction of LLP decays in the MS fiducial volume that are selected by the trigger as a function of the LLP decay position, is shown in Figures 2(a) and 2(b) for four MC simulated benchmark samples with LLP decays in the MS barrel and endcap regions, respectively. The efficiency is parameterized as a function of the transverse decay radius ( ) in the barrel and the longitudinal decay position (| |) in the endcaps. The events are required to pass the data quality requirements [53] and have a reconstructed primary vertex as described in Section 6.3. The trigger is efficient for hadronic decays of LLPs that occur anywhere from the outer regions of the HCal to the middle stations of the MS. These efficiencies are obtained from a subset of simulated signal events with only a single LLP decay in the muon spectrometer in order to ensure that the result of the trigger is due to a single burst of MS activity. The uncertainties shown are statistical only. The relative differences between the efficiencies of the benchmark samples are a result of the different masses of the LLPs, which in turn affect their momenta and consequently the opening angles of the decay products. The trigger efficiency increases as the LLPs decay closer to the end of the hadronic calorimeter, located at ≈ 4 m for the barrel and | | ≈ 6 m for the endcaps. The efficiency decreases rapidly as the decay occurs closer to the middle stations of the muon spectrometer (barrel: ≈ 7 m; endcaps: | | ≈ 13 m). To quantify the mismodeling of the L1 muon trigger efficiency in MC simulation, the distributions of the number of muon-RoI clusters within a cone of Δ = 0.4 around the axis of a punch-through jet in multĳet MC and data events are compared. High-energy jets in data were selected using jet triggers with a T threshold of 400 GeV or 420 GeV, depending on the data-taking period. MC simulation shows a rate of muon-RoI clusters that is 24% (20%) lower than in data in the barrel (endcaps). This mismodeling does not depend on the , or T of the jet, and it is used as a systematic uncertainty on the signal trigger efficiency.

Reconstruction of MS displaced vertices
A dedicated algorithm [48], capable of reconstructing low-momentum tracks in a busy environment, is used to reconstruct the MS DV. This algorithm was also employed by previous searches for displaced decays in the MS [27, 29, 38]. The algorithm takes advantage of the spatial separation between the two multilayers inside a single MDT chamber. Single-multilayer straight-line segments that contain three or more MDT hits are reconstructed using a 2 fit. Segments from multilayer 1 are then matched with those from multilayer 2. The matched pair of a multilayer-1 segment with a multilayer-2 segment and their corresponding track parameters is called a tracklet. Tracklets that are spatially clustered within |Δ | < 0.7 and |Δ | < /3 in the detector are extrapolated backwards to reconstruct the ( , ) position of the MS DV with a 2 fit. At least three (four) tracklets are required to be used in the fit in the barrel (endcaps). Detectable vertices arise from decays that occur in the region between the outer edge of the HCal and the middle stations of muon chambers. Because their detector technology is different (no spatially separated multilayers), the CSC chambers are not used for the MS DV reconstruction.

Reconstruction of the primary vertex and prompt hadronic jets
Events are required to have a primary vertex (PV) with at least two tracks with T > 500 MeV. If more than one PV candidate is reconstructed, the candidate with the largest sum of the squares of the transverse momenta of all tracks associated with the vertex is selected.
Hadronic jets are constructed by using FastJet [66] to apply the anti-jet algorithm [67] with a radius parameter = 0.4. A collection of three-dimensional topological clusters of neighboring energy deposits in the calorimeter cells containing a significant energy above a noise threshold [68] provide input to the anti-algorithm. After reconstruction, jets are calibrated using the procedure outlined in Refs. [69][70][71][72][73][74][75].

Event selection
All events considered in this search must pass an event-level selection to distinguish signal from background.

Baseline event selection
A common baseline selection is applied to the events considered in this search and summarized in Table 2.
Events with at least one MS DV are required to pass the data quality requirements [53] and the Muon RoI Cluster trigger, and contain a PV. The PV selection has very little impact on the signal efficiency but it helps to reject background events; in simulation, the selected PV corresponds to the signal interaction in about 95%-99% of the cases, depending on the sample. Even though the LLPs are invisible in the ID, the scalar boson is produced with an average T of 20 GeV, and a number of tracks are left by recoiling particles produced in the same interaction.
An MS DV that arises from a displaced decay typically has many more hits than an MS DV from background.
To take advantage of this difference, a minimum number of MDT ( MDT ) and RPC/TGC ( RPC / TGC ) hits is required. The MDT hits are counted in the MDT chambers that have their center within Δ = 0.6 and Δ = 0.6 of the DV ( , ) direction. The RPC and TGC hits are the sum of hits that are within Δ = 0.6 of the DV position. A requirement on the maximum number of MDT hits is also applied to remove background events caused by coherent noise bursts in the MDT chambers. This selection has a negligible impact on the signal events.
A displaced decay that occurs in the transition region between the MS barrel and endcaps results in hits in both regions. Vertex reconstruction is performed separately in the barrel and endcaps, and only the barrel (endcap) hits are used in the barrel (endcap) vertex reconstruction algorithm. Therefore, MS vertices reconstructed from either of the two algorithms have fewer hits, as they are reconstructed from a subset of the full set of hits. This results in a decrease in the reconstruction efficiency, and occasionally two vertices being reconstructed from a single LLP decay. Therefore, the MS DVs with pseudorapidity | vx | between 0.8 and 1.3 are not considered in the analysis. This has a negligible effect on the total signal efficiency, because the average MS DV reconstruction efficiency in this region is less than 2%. Moreover, in the transition region between the barrel and the endcap hadronic calorimeters, 0.7 < | vx | < 1.2, the probability of having a jet that does not fulfill the minimal selection criteria to be considered for isolation and also punches through into the MS is much higher than in other regions of the detector. This region overlaps the already excluded MS transition region, except for 0.7 < | vx | < 0.8. Therefore, vertices reconstructed with pseudorapidity in the range 0.7 < | vx | < 0.8 are also not considered.
The number of events passing the baseline selection is 1,385,587 and 5,138,794 in the barrel and endcaps, respectively. After these requirements, the main background contribution is from punch-through jets.

Signal displaced vertex selection
To reject background vertices created by punch-through jets, a set of vertex isolation criteria was established. These criteria are based on the angular distance, Δ , between the direction of the tracks or jets and the vertex axis, defined as the line from the IP to the DV. No jets or tracks should be present in a Δ cone around the MS DV axis. For isolation from tracks, two criteria are used. One is for isolation from tracks with T > 5 GeV (high-T tracks). The other takes the vector T sum of all tracks associated with the primary vertex (PV) that have 500 MeV < T < 5 GeV (low-T tracks) and are in a cone of Δ = 0.2 around the MS DV axis . The use of two different isolation criteria stems from the fact that some jets have most of their energy in a few hadrons, while others can consist of multiple low-T tracks. The isolation criteria are summarized in Table 3. An MS DV that satisfies these criteria is considered in the analysis.

Isolation requirements Barrel Endcaps
Isolation from high-T tracks ( T > 5 GeV) The isolation criteria were optimized for the MC benchmark samples by comparing simulated signal events with simulated multĳet events. Figure 3 shows the cumulative vertex efficiency as a function of the isolation requirements in the barrel for data events and simulated multĳet and signal events.
All jets considered for DV isolation must satisfy T > 30 GeV and log 10 ( HAD / EM ) < 0.5, where HAD is the jet energy deposited in the HCal, and EM the energy deposited in the ECal. The latter criterion is commonly used to identify the decay of a neutral LLP in the hadronic calorimeter [41] that results in little or no energy deposited in the electromagnetic calorimeter and consequently in an anomalously large value of the hadronic to electromagnetic energy ratio. It ensures that vertices originating from LLPs that decay near the outer edge of the hadronic calorimeter and also have significant MS activity are not rejected by the isolation requirement. In addition, in order to reduce the probability that pileup jets prevent a signal vertex from meeting the isolation criteria, jets with 20 < T < 60 GeV must be matched to the PV by using a jet vertex tagger (JVT) discriminant [71]. Standard jet-quality criteria [76] are not applied because jets that do not fulfill these requirements can also produce a background MS DV, and therefore need to be considered when computing the isolation.
At least two isolated MS DVs must be present in the event. One MS DV must be matched to the trigger-level muon-RoI cluster by satisfying Δ (cluster, vertex) < 0.4. If there are two distinct clusters, each MS DV must be matched to a different cluster. To ensure that the two MS DVs and/or two muon-RoI clusters do not come from the same background activity, the two vertices are required to be separated by at least Δ = 1.0, which has minimal impact ( 3%) on the overall signal acceptance. After the vertex isolation criteria are  Figure 3: Cumulative efficiency of vertices in the barrel where the vertex is required to be isolated from (a) jets and (b) high-T tracks, as a function of the selected Δ , when the sum of low-T tracks in a Δ = 0.2 cone around the vertex direction is required to be less than a specified cut in T (c), and when the sum of low-T tracks in a Δ cone around the vertex is required to be less than 10 GeV, as a function of the Δ value (d). The vertical lines show the selected Δ value that is used in analysis. The differences between Pythia multĳet and data distributions are attributable to the noncollision background, which is not present in MC simulations. applied, most vertices from punch-through jets are eliminated, and the main background contribution is from noncollision backgrounds.

MS displaced vertex reconstruction efficiency
The efficiency for vertex reconstruction [48] is defined as the fraction of simulated LLP decays in the MS fiducial volume that correspond to a reconstructed vertex passing the baseline event selection reported in Table 2 and satisfying the vertex isolation criteria summarized in Table 3. A reconstructed vertex is considered matched to a displaced decay if the vertex is within Δ = 0.4 of the simulated decay position. The MS DV efficiency is parameterized as a function of the and | | LLP decay position in the barrel and endcaps, respectively. Figure 4(a) shows the efficiency for reconstructing a vertex in the MS barrel for a selection of benchmark samples. Figure 4(b) shows the efficiency for reconstructing a vertex in the MS endcaps. For the MC samples considered in this paper, the MS barrel vertex reconstruction efficiency is O (2 − 15)% near the outer edge of the hadronic calorimeter ( ≈ 4 m) and decreases substantially as the decay occurs closer to the middle stations ( ≈ 7 m). The decrease occurs because the charged hadrons and photons are not spatially separated and overlap when they traverse the middle stations. This results in a reduction of the efficiencies for tracklet reconstruction and, consequently, vertex reconstruction. The efficiencies are also shaped by the mass and boost of the LLP. The efficiency for reconstructing vertices is higher in the MS endcaps due to a more efficient selection and vertex reconstruction, and it reaches 40% for higher-mass benchmark models. Since the magnetic field in the region in which endcap tracklets are reconstructed is very weak, the vertex reconstruction does not have the curvature constraints on the tracklets from the charge and momentum measurements that are present in the barrel. Therefore, in the endcaps, the vertex reconstruction algorithm uses straight-line fits so that low-momentum tracks are not rejected, while in the barrel the curvature plus combinatorics provide better rejection of misreconstructed tracks. Consequently, the vertex reconstruction in the endcaps is more efficient for signal, but also less robust in rejecting background events. More details are provided in Ref. [48].
Potential MC simulation inaccuracies in the DV reconstruction are estimated by comparing the distribution of the number of tracklets in a punch-through jet Δ cone of 0.4 in data and MC events using a strategy similar to the one used for the trigger, which was described in Section 6.1. MC simulation shows a rate of tracklets that is 20% (15%) higher than in data in the barrel (endcaps). This mismodeling does not depend on the , or T of the jet. Its impact on the DV reconstruction efficiency is estimated by randomly dropping tracklets used to reconstruct the vertex in accordance with the measured mismodeling factor in the barrel and endcaps and counting the number of reconstructed vertices. The efficiency variation between the nominal reconstruction and the one considering the mismodeling factors, averaged over the MC benchmark samples, is 27% in the barrel, and 9% in the endcaps. The corrected efficiencies are used to compute the expected number of signal events. The systematic uncertainty associated to this correction is discussed in Section 9.
Over the range of the extrapolated lifetimes, the maximum acceptance times efficiency for vertices passing the signal selection ranges from 0.005% (for the MC sample with Φ = 60 GeV, = 5 GeV at a proper lifetime of 0.5 m) to 0.3% (for the MC sample with Φ = 600 GeV, = 50 GeV at a proper lifetime of 1.1 m).

Background estimation
The two-MS-vertex strategy is designed to be sensitive to models where two LLPs are produced and decay hadronically between the outer region of the HCal and the middle stations of the MS. Requiring two MS DVs significantly reduces the background. In addition, background from punch-through jets is further reduced using the isolation criteria described in Section 7.2.
Residual background can arise from collision or noncollision processes and cannot be simulated accurately. Thus, a data-driven method is used to estimate the expected background. This has the advantage of avoiding the systematic uncertainties related to simulation accuracy present when using MC-based background estimates.
As described in Section 7.1, the signal selection requires that one vertex is always matched to a RoI cluster, and only in the case of two RoI clusters the second vertex is also required to be matched. Therefore, we can naturally factorizes the background estimation into three terms: The first term in Eq. (1) gives the number of background events that contain one vertex-cluster pair anywhere in the detector fiducial region and one other isolated MS vertex anywhere in the same fiducial region. The second and third terms estimate the number of background events that contain one vertex-cluster pair anywhere in the detector fiducial region and an additional vertex-cluster pair in the barrel or endcap region. The background from events with one muon-RoI cluster (first term in Eq. (1)) is estimated by multiplying the number of events with one muon-RoI cluster ( 1cl ) by the probability of finding an isolated MS DV vertex in events not selected by the Muon RoI Cluster trigger ( Vx noMStrig ). The probability Vx noMStrig is determined using data collected with the zero-bias trigger and it is found by dividing the number of events with an isolated MS DV ( 1Vx ) by the total number of zero-bias events that satisfy standard event-quality requirements ( Events ). The background from events with two muon-RoI clusters (second and third terms in Eq. (1)) is estimated by multiplying the number of events with two muon-RoI clusters where one of the two clusters is not matched to an isolated vertex in the barrel ( 2cl 1UMBcl ) or endcaps ( 2cl 1UMEcl ) by the probability of finding an MS DV vertex given a muon-RoI cluster in the barrel or endcaps respectively ( Vx Bcl and Vx Ecl ). The probabilities Vx Bcl and Vx Ecl are found by dividing the number of events with an isolated MS DV matched to a muon-RoI cluster ( 1BclVx and 1EclVx ) by the number of events with a muon-RoI cluster ( 1Bcl and 1Ecl ).
All the numbers of events and probabilities used in Eq. (1) are reported in Table 4. All events used for the background estimation are orthogonal to the signal region that requires two isolated DVs. Therefore, events that pass the Muon RoI Cluster trigger and have only one isolated MS DV are assumed to be fake vertices not due to displaced decays. Potential signal contamination would result in an overestimation of the background rates, but its contribution has been estimated to be smaller than 0.02% for injected signal with rates of the order of the observed 95% CL upper limits.
The background estimation strategy was validated using an orthogonal signal-free data sample obtained by inverting the isolation criteria used to define the signal region. In this validation region, the expected number of background events was calculated to be 0.99 ± 0.20 and 1.56 ± 0.27 for one-cluster and two-cluster events, respectively, where the uncertainty is purely statistical. One event, which contained two clusters, was observed. Therefore, no systematic error is assigned to the background estimation.
The number of expected background events with two MS DVs is 0.32 ± 0.05, where the uncertainty is statistical only.

Systematic uncertainties
The signal efficiency systematic uncertainties are dominated by the modeling of the signal physics processes, pileup and detector response and the extrapolation of the expected number of signal events as a function of the LLP proper lifetime.
One of the sources of systematic uncertainty associated with the Muon RoI Cluster trigger is the modeling of the minimum-bias interactions used to simulate pileup; another stems from the systematic uncertainty due to the PDF used to generate signal MC events. These were estimated by varying the pileup and PDF weights in accord with the respective ±1 systematic uncertainties and evaluating the resultant change in the trigger efficiency. For the latter, uncertainties in the nominal PDF set were evaluated using 100 replica variations. In each case, the systematic uncertainty was determined to be negligible. The systematic uncertainty from the pileup and PDF contributions to the MS DV reconstruction efficiency for signal was evaluated with a procedure similar to that for the trigger, and again the systematic uncertainties were determined to be negligible.
The mismodeling of the L1 muon trigger efficiency in MC simulation described in Section 6 contributes a systematic uncertainty of 20% (24%) in the barrel (endcaps). Another source of systematic uncertainty is the uncertainty in the corrected efficiencies that is due to the MC mismodeling in the vertex reconstruction described in Section 7.3. The corrected vertex reconstruction efficiencies were also estimated with a second method, in which events are weighted such that the MC distribution of the number of reconstructed tracklets matches the data distribution. The efficiency was re-evaluated after reweighting each vertex according to how many tracklets are associated with it. The average efficiency variation is computed and its difference from the variation calculated using the nominal method described in Section 7.3 is taken as the systematic uncertainty, which is 11% in the barrel and 13% in the endcaps.
A systematic uncertainty associated with the efficiency extrapolation method was estimated by comparing the signal efficiency computed using the fully simulated MC samples with the one extrapolated (using toy MC samples) at the same proper lifetime. The difference between the two efficiencies is used as the systematic uncertainty, and it varies from 1.9% to 30%, depending on the kinematics of the sample. For several of the signal samples, two proper lifetime points were fully simulated: one nominal sample and a secondary sample with longer proper lifetime, as described in Section 5. The secondary sample was used to cross-check the lifetime-extrapolation procedure, and good closure was found.
An uncertainty on the NLO-reweighting of the signal samples is obtained by comparing the 125-GeV mediator NLO M G predictions to next-to-next-to-leading-order (NNLO) accuracy in QCD using P B v2 [77][78][79][80][81]. This results in an additional uncertainty on the signal efficiency ranging from 0.1 to 4%.
The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [82], obtained using the LUCID-2 detector [83] for the primary luminosity measurements.
The systematic uncertainties described above change the production cross-section limits by about 15% on average.

Results
In ATLAS Run 2 data, 0.32 ± 0.05 background events are expected for this analysis. Zero events are observed in the signal region.
Upper limits on the production cross-section times branching fraction are derived using the CL s prescription [84], implemented with the pyhf [85, 86] package using a profile likelihood function [87]. The likelihood includes a Poisson probability term describing the total number of observed events. Background and signal uncertainties are taken into account using Gaussian terms, as commonly done in these cases. Pseudo-experiments which sample the distribution of the profile likelihood ratio are generated to compute the -value and derive the exclusion limits.
For scalar boson benchmark samples with ≠ 125 GeV, upper limits are set on × , where represents the branching fraction for → assuming 100% branching fraction of into the heaviest fermion pairs kinematically accessible. As discussed in Section 4, the long-lived scalar mainly decays into , except when > 2 (i.e. for the sample with = 1000 GeV and = 475 GeV) where the dominant decay is into . For scalar boson benchmark samples with = 125 GeV, upper limits are set on ( / SM ) × , where SM = 48.61 pb [88] is the SM Higgs boson gluon-gluon fusion production cross-section. Figure 5 shows a comparison between the expected and observed 95% CL upper limits for one representative sample as well as the observed 95% CL limits for all the MC benchmark samples. Observed limits are consistent with the expected ones within the uncertainty. For = 125 GeV the limits are stronger for intermediate LLP masses, while they become weaker for very low and very high masses. Moreover, the mean proper lifetime at which the upper limit is strongest increases with the long-lived scalar mass. These patterns of behavior are correlated with changes in the trigger and reconstruction efficiencies that depend mainly on the kinematics of the LLP decay. Table 5 summarizes the lifetime ranges excluded by the analysis for branching fractions ( (125) → ) = 10%, 1% and 0.1% for the scalar boson with = 125 GeV decaying into two long-lived scalars. Table 5: Ranges of mean proper lifetime excluded at 95% CL for scalar boson benchmark models with = 125 GeV, assuming a production cross-section times branching fraction equal to 10%, 1% and 0.1% of the SM Higgs boson production cross-section [88].

Summary
This paper presents a search for events with two displaced vertices from pair-produced long-lived particles decaying into hadronic jets using 139 fb −1 of proton-proton collisions at √ = 13 TeV recorded at the LHC by the ATLAS detector during 2015-2018 data-taking period. The search is performed with the same strategy that ATLAS adopted in Run 1 and in a Run 2 search using 2015-2016 data, and benefits from very low background. A data-driven method is used to estimate the expected background in the signal region. No events with two reconstructed displaced vertices in the muon spectrometer are observed in the signal region, and exclusion limits on the LLP production cross-section as a function of its proper lifetime are computed for a scalar portal model where a Higgs boson with a mass of 125 GeV or another short-lived scalar can decay into two long-lived scalars. For the 125 GeV Higgs boson, the paper reports the first exclusion limits for branching fractions below 0.1%, while branching fractions above 10% are excluded at 95% confidence level for LLP mean proper lifetimes ranging from 4 cm to 72.4 m. In addition, the paper presents the first results for the decay of LLPs into in the ATLAS muon spectrometer.