Search for Resonant Production of Dark Quarks in the Dijet Final State with the ATLAS Detector

This paper presents a search for a new $Z^\prime$ resonance decaying into a pair of dark quarks which hadronise into dark hadrons before promptly decaying back as Standard Model particles. This analysis is based on proton-proton collision data recorded at $\sqrt{s}=13$ TeV with the ATLAS detector at the Large Hadron Collider between 2015 and 2018, corresponding to an integrated luminosity of 139 fb$^{-1}$. After selecting events containing large-radius jets with high track multiplicity, the invariant mass distribution of the two highest-transverse-momentum jets is scanned to look for an excess above a data-driven estimate of the Standard Model multijet background. No significant excess of events is observed and the results are thus used to set 95 % confidence-level upper limits on the production cross-section times branching ratio of the $Z^\prime$ to dark quarks as a function of the $Z^\prime$ mass for various dark-quark scenarios.


Introduction
In recent years, exploring beyond the Standard Model (BSM) scenarios involving an additional non-Abelian hidden gauge group has been of increasing interest.This sector, which would contain dark quarks and dark gluons, could provide dark matter candidates in the form of stable dark hadrons, leading to very peculiar signatures at colliders which have not yet been thoroughly explored.For a review of such models and signatures see Ref. [1] and references therein.
At the Large Hadron Collider (LHC) dark quarks could be produced in the decay of a heavy mediator.The dark quarks would then hadronise into dark hadrons, whose decay into Standard Model (SM) particles would produce objects which appear as jets with properties that depend on the underlying theory parameters.Limits were set for signatures which include jets aligned with missing transverse momentum (semi-visible jets) using the full Run 2 data samples by the ATLAS and CMS Collaborations [2,3], and for jets with displaced tracks (emerging jets) using partial Run 2 data samples by the CMS Collaboration [4].
In this paper, the resonant production of a pair of dark quarks,   , through a heavy  ′ mediator is considered, as shown in Figure 1.In the benchmark models considered here, following Ref.[5], the dark hadrons are assumed to decay promptly into SM particles and the fraction of invisible components (stable dark hadrons) produced is assumed to be negligible, thus covering a complementary part of the parameter space relative to semi-visible jet and emerging jet searches, whose sensitivity to fully visible jets is reduced by requirements on the presence of missing transverse momentum.The dark jets considered here are typically wider than the SM QCD jets due in part to the double hadronisation process which would take place, first in the dark sector and then in the SM, after the dark mesons decay in part to high-momentum SM quarks.While tracks associated to SM gluon and quark jets mostly come from charged pions and are thus directly linked to fragmentation, the charged particles inside the dark jets come from dark meson decays.As discussed in Ref. [5], the number of tracks associated to dark jets is thus influenced on the one hand by the multiplicity of dark mesons produced in the dark shower, which depends on the dark fragmentation and dark meson masses, and on the other hand by the decay channels of these dark mesons into SM particles.For the dark sector models considered here, this leads to a higher charged-particle multiplicity for dark jets than for SM jets.This analysis hence looks for a resonant excess above a smooth background in the dĳet invariant mass distribution, after preselecting events containing at least two high transverse momentum ( T ) large-radius jets with high track multiplicity, thus complementing the  ′ →  q dĳet resonance search [6].The invariant mass spectrum of the dominant SM QCD dĳet background is determined using a data-driven method.

ATLAS detector
The ATLAS experiment [7] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a nearly 4 coverage in solid angle. 1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer.The inner tracking detector covers the pseudorapidity range || < 2.5.It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors.Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity.A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|| < 1.7).The endcap and forward regions are instrumented with LAr calorimeters for both the electromagnetic and hadronic energy measurements up to || = 4.9.The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each.The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering.
A two-level trigger system is used to select events.The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz.This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average.An extensive software suite [8] is used in 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.

Data and simulated events
The data sample used in this analysis was collected using unprescaled single-jet triggers, where jets were reconstructed exploiting the anti-  clustering algorithm [9,10] with a radius  = 1.0 and a  T threshold of 460 GeV, with varying jet trimming [11] requirements depending on the year of data taking.These triggers are fully efficient after applying the offline object requirements described in the next section.The integrated luminosity of the data sample is 139 fb −1 after applying the stable beam conditions, detector and data quality requirements [12].
The background from multĳet processes is simulated with the Pythia 8.186 [13] event generator.The A14 set of tuned parameters [14] is used together with the NNPDF2.3LOPDF (parton distribution function) set [15].The contributions of other processes is found to be negligible after selections and not considered further in the analysis.
Signal samples of   →  ′ →   q , considering an SU(3) symmetry in the dark sector, are simulated using the Hidden Valley module of Pythia 8.235 [16] with  ′ masses varying from 1.5 to 3.5 TeV, using the same set of tuned parameters and PDF set as the multĳet samples.Between fully simulated points, which are spaced by 250 or 500 GeV, the signal acceptance and shape are interpolated.The production cross-section is a free parameter as it depends on the coupling of the  ′ to SM quarks,   .This parameter is set indirectly by choosing the  ′ width and its relative branching ratio to SM quarks.For the cross-section 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe.The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards.Cylindrical coordinates (, ) are used in the transverse plane,  being the azimuthal 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 values explored in this paper, the scaling   does not affect the dĳet mass distribution because the true  ′ width remains much smaller than the reconstructed width that is dominated by resolution and reconstruction effects.To probe the dark QCD sector, four combinations of parameters are chosen to cover different possibilities for the number of dark quark flavours,   , the confinement scale Λ  , the mass spectra and the prompt decay modes of the dark pions,   .The values of these parameters are chosen following the benchmark models proposed in Ref. [5], and are summarised in Table 1.While the decay of the dark rho mesons is fixed to   →     , the dark pions are either assumed to decay directly to SM quarks (models A and B) or to SM particles through dark photons,  ′ (models C and D).In the latter case, the dark photons decay according to the branching ratios computed in Ref. [17].In model C, the dark photons decay into any pair of charged fermions excluding third-generation quarks (due to the choice of dark photon mass), while in model D the dark photons decay predominantly into charged pions, and occasionally to electron and muon pairs.
Table 1: Benchmark models taken from the reference in the text (Ref.[5]).In the headings,   corresponds to the number of low-mass dark quarks flavors, Λ  is the dark confinement scale and the masses m ′ ,    ,    correspond to the constituent quark mass-the masses of the dark pion and dark rho mesons respectively. Model For all simulated samples, the EvtGen 1.2.0 program [18] is used for the properties of -and -hadron decays.The generation of the simulated event samples also includes the effect of multiple   interactions per bunch crossing, and the effect on the detector response due to interactions from bunch crossings before or after the one containing the hard interaction.Multiple overlaid   collisions are simulated with the soft QCD processes of Pythia 8.186 [16] using the A2 set of tuned parameters [19] and the MSTW2008LO PDF set [20].All the Monte Carlo events are processed through a full simulation of the ATLAS detector geometry and response [21] based on Geant4 [22].

Object and event selections
Events are required to have at least one   interaction vertex with at least two tracks with  T greater than 500 MeV.The primary vertex is defined as the one with the largest Σ 2 T of its associated tracks.Events are rejected if they contain jets compatible with noise bursts, beam-induced background or cosmic rays as defined in Ref. [23].
Large- jets are built from three-dimensional topological clusters of energy deposits in the calorimeter, calibrated to the hadronic energy scale with the local cluster weighting (LCW) [24] procedure using the anti-  algorithm [9,10] with a radius parameter  = 1.0.To reduce contributions from pile-up, initial-state radiation and multiple parton interactions that would contaminate the measurement of the hard scatter process, the jets are trimmed [11].The trimming procedure reclusters jet constituents into subjets using the   algorithm [25][26][27] with  = 0.2 and discards the resulting subjets if their  T is less than 5% of the  T of the parent jet [28], before recomputing the four-momenta of the large- jets using the remaining subjets.This analysis selects events containing at least two large- jets satisfying the following requirements: 500 <  T < 3000 GeV for the leading and 400 <  T < 3000 GeV for the subleading jet, || < 2.0 and 50 <  J < 600 GeV, where  J is the invariant mass of the jet.The upper bounds on  T and  J have a small impact on the signal acceptance while ensuring that jet energy and mass scale calibrations and their uncertainties remain valid [29,30].Finally, along with the  T > 500 GeV selection on the leading jet,  JJ is required to be greater than 1.3 TeV to ensure full trigger efficiency.
The distributions of the dĳet invariant mass of the multĳet background and benchmark models with  ′  = 2.5 TeV are shown in Figure 2 after these preselections.The bin width in  JJ is determined by the experimental dĳet mass resolution (evaluated on background simulation) of approximately 100 GeV for the background in this mass range.The signal distributions are relatively broad and shifted to lower values relative to the true  ′ mass.This is mainly due to the removal of some of the jet constituents by the trimming algorithm, which also reduces the background under the bulk of the signal peak [11].The number of tracks matched to the untrimmed jets via ghost association2 [31],  track , is helpful in separating the signal from the background.These tracks are required to have  T > 500 MeV, fulfil hit requirements corresponding to the Loose working point [32] and to be matched to the primary vertex. 3he distributions of  track for the multĳet background and representative signals are shown in Figure 3 after the preselections for the two leading jets are applied.In the high end of the spectrum where the selection is made, the data agree with the background prediction within the modelling uncertainties that are described in Section 6.It was checked that the small fluctuation found at very low values of  track has a negligible impact on the analysis.To enhance the signal-to-background ratio, it is beneficial to require that the jets have a minimal value of  track .However, such a requirement would sculpt the  JJ background distribution since the number of tracks associated with a jet is correlated to its  T .To avoid this undesirable effect, which would impact the identification of a clear resonance over the background spectrum, a new variable,   track , is defined as follows: • First, a target efficiency, , for a background jet to pass the requirement on  track is defined.In this analysis, a value of  = 1% is chosen to reduce the background significantly while maintaining good acceptance across the various signals probed.
• Second, for each bin in  JJ and for the leading jet of the simulated multĳet background, the minimal value  J 1 for which  track,1 >  J 1 leads to the target efficiency, , is derived (where the subscript "1" refers to the leading jet).
• Third, the procedure is repeated with the subleading jet to determine  J 2 , after applying  track,1 >  J 1 on the leading jet.This is done to account for possible correlation effects between the two leading jets.Applying an  track requirement on the leading jet selects preferentially events which also contain higher- track subleading jets, raising the value of  J 2 needed to reach the target efficiency.
• Fourth, the  J i values obtained for the two leading jets ( = 1, 2) are fitted as a function of  JJ to have a smooth distribution.The values obtained for  J 1 and  J 2 are shown in Figure 4.
• Finally, the variable   track, is defined for the two leading jets as   track, =  track, −  J i .Requiring   track,1 > 0 removes 1 −  of the multĳet background; requiring   track,2 > 0 further removes 1 −  of the remaining background.
The signal region (SR) is defined by requiring   track,1 > 0 and   track,2 > 0 after the preselections are applied, which reduces the multĳet background by approximately 99.99%.This selection also reduces the already-negligible backgrounds from top and electroweak processes (< 1% compared to the multĳet Figure 4: Distributions of  J 1 (bottom curve) and  J 2 (top curve) as a function of  JJ for the simulated multĳet background.These variables are used to define a selection on  track which evolves with  JJ in order to achieve a fixed background selection efficiency as a function of  JJ .As mentioned in the text, the value of  J 2 is determined after applying  track,1 <  J 1 on the leading jet, thus preferentially selecting events which also contain higher- track subleading jets.background).For  ′  = 1.75 TeV the signal acceptance times efficiency after all selections lies between 0.3% and 9%, while at  ′  = 2.5 (3.5) TeV it is between 1% (2.5%) and 28% (37%) depending on the signal models.

Background estimation
The requirement on   track can also be reversed by requiring   track < 0 for the two leading jets, in order to construct a control region (CR) which is largely enriched in multĳet background while maintaining the same  JJ shape, as per definition   track is decorrelated from  JJ .This allows a data-driven estimate of the background in the signal region.The background contribution to the SR is determined by taking the  JJ distribution directly from the CR, while the overall background normalisation is left as a free parameter to be extracted from a likelihood fit in the SR.This CR method is adopted for the background estimate because of the relatively large reconstructed width of the signal models considered.Signal injection tests, in which signal models are added to the data-driven expected background distribution in the SR, show that the extracted signal yields obtained when using the CR method are consistent with the injected values within uncertainties.The background shape template from the CR can create a bias in the likelihood fit, a "spurious signal,"  spur , which corresponds to the extracted number of signal events when a signal-plus-background fit is conducted on a background-only sample.A spurious signal test is thus conducted and is shown to satisfy the criterion motivated in Ref. [33] for all models considered.The requirement that  spur / fit < 0.5 is used where  fit is the uncertainty in the fitted number of signal events.In practice it is found to be less than 0.35 for all the models and masses considered.
A validation region (VR), between the CR and the SR, is also built by requiring   track,1 > 0 and   track,2 < 0. This region allows the verification that the decorrelation of   track and  JJ , derived from the simulated multĳet processes, operates well in data: in the VR, the scaled CR template agrees with the data with a chi-squared per degree of freedom of 1.2.This good agreement is shown in Figure 5 in which the  JJ data distribution in the CR is compared to that in the VR after normalising the total number of events in each region to the one expected in the SR.The general  JJ shape agrees well between the CR and the VR, which confirms that the decorrelation method works within the statistical uncertainty expected in the SR.Both the CR and the VR have significantly more background events than the SR: the expected signal-to-background ratio in these regions is thus found to be negligible with respect to that in the SR, for all the models considered.Therefore, any bias from the presence of signal events in these regions will be negligible, in either the construction of the background model or in the above cross-checks.The statistical analysis of the data follows a frequentist approach based on a profile-likelihood ratio test statistic [34].The parameter of interest-the signal strength -is extracted from a binned maximumlikelihood fit of a signal-plus-background model to the  JJ distribution in the SR.The signal strength then acts as a scale factor on the predicted number of events for each model assumption and is allowed to take negative values in the fit.The bin width is fixed at 100 GeV, as motivated by fitting the experimental dĳet mass resolution of the background simulation, while the  JJ range probed extends from 1.3 to 5 TeV.The likelihood is defined from the Poisson probability to observe  data events for a given expected signal   and background   , where  is the background normalisation scale factor that is extracted from the fit, as: where the nuisance parameters, ì , representing the systematic uncertainties , are included with Gaussian constraint terms .The systematic uncertainties are detailed in Section 6.The expected background shape,  , is taken from the CR.

Systematic uncertainties
Systematic uncertainties in the background estimate and in the expected signal yields and shapes arise from detector effects and Monte Carlo (MC) modelling; they are treated as nuisance parameters in the statistical analysis.Given that the shape of the background is taken from the CR and that its normalisation is a free parameter determined by a fit in the SR, for all the uncertainties listed below, only their effects on the shape agreement between the SR and the CR are considered as uncertainties in the background.
The experimental uncertainties include uncertainties on the jets, tracks and luminosity.The uncertainties associated with the large- jets arise from the jet energy and mass scales, mass response and energy resolution [29,35].For the model B, C and D simulated signal samples, an additional jet energy scale uncertainty of 5% is added to cover for a small deviation from unity in the jet energy response with respect to the true energy scale (non-closure) observed after calibration following the SR selections.For model A and for the SM dĳet background, closure was obtained within the nominal uncertainties.Uncertainties in the track reconstruction performance are related to their reconstruction efficiency and to the rate at which fake tracks are reconstructed [32].These also include additional uncertainties related to tracking in dense environments, given the large number of tracks requested in the core of high- T jets in this analysis [36].The limited tracker coverage in combination with the use of large- jets does not affect the results, since its effect is found to be small and well modelled within the CR and the VR.The uncertainty in the combined 2015-2018 integrated luminosity, only affecting the signal samples, is 1.7% [37], obtained using the LUCID-2 detector [38] for the primary luminosity measurements.
The MC generator uncertainties are applied to signal and include uncertainties in the PDFs [39], renormalisation and factorisation scales, and strong coupling constant, as well as uncertainties related to the parton shower set of tuned parameters in initial-or final-state radiation, such as variations of the renormalisation scale for QCD emission ( , ) or the inclusion of non-singular terms [40].
Following the prescription described in Ref. [33], an uncertainty related to the spurious signal test is added.It is estimated by performing a signal-plus-background fit on the VR data which is scaled down to the expected number of events in the SR.
The statistical uncertainty dominates over the systematic uncertainties over the full mass range probed.The impact of the largest relative systematic uncertainties on  are shown in Table 2 for the four models and an example  ′ mass of 2.5 TeV.They are determined by re-running the fit while fixing, in turn, each of these nuisance parameters to their post-fit value plus or minus one standard deviation.Models B and C are most affected by uncertainties which can shift signal events in and out of acceptance and, due to a lower number of tracks and a more displaced  JJ peak, these models have the lowest acceptance.

Results
Figure 6 shows the observed  JJ distribution in the SR.The highest mass event in data is observed at 4.8 TeV.A BumpHunter test [41] is performed on the distribution to test its compatibility with the background-only hypothesis: the vertical lines on the Figure indicate the most discrepant interval identified which has a global p 0 -value of 0.63 in the interval 1500 GeV <  JJ < 1700 GeV.This p 0 -value and the bin-by-bin significance of the difference between the data and the fit, shown in the lower panel of the plot, consider only statistical uncertainties.Given no significant deviation from the expected background is observed in Figure 6, limits can be set on the various models considered as shown in Figure 7: upper limits on the signal production cross-section times branching ratio of the  ′ to dark quarks are extracted at 95% confidence level (CL) using the CL s method [42], with the asymptotic formulae 4 for the profile likelihood ratio test statistic distribution [34].The  ′ -mass lower bound in the limit plots is chosen such that the mean of the signal  JJ distribution is above the lowest  JJ value probed of 1.3 TeV.Given the  JJ shape of the data in the SR, the signal-plus-background fit is able to accommodate some large-width, low-  ′ signal by lowering the background normalisation at low- JJ values, while still being compatible with the data in the higher- JJ region.The converse is true for higher-  ′ values, resulting in the observed limit shape.
In Figure 7, the limits are also compared to theory predictions for given values of   and    ; these are obtained at leading-order in QCD with MadGraph5_aMC@NLO [43], using a model based upon that used in the  ′ →  q ATLAS dĳet resonance search [6], in which the coupling to dark matter particles is replaced by a coupling to dark quarks.In that model, a universal coupling of the  ′ to each SM quark   is assumed and the dark quarks are all represented by one species with an effective coupling    to the  ′ .For   = 0.05 and    = 0.2, a value of   which evades the dĳet resonance search constraints, the observed (expected) limit on the  ′ mass is at 3.0 (2.9) TeV for model A and reaches up to 2.9 (2.8) TeV for model D (for which the region between 2.1 and 2.2 TeV is not excluded), while for models B and C, due to a lower acceptance, these couplings do not lead to any exclusion.Pushing the couplings up to   = 0.15 and    = 0.5, a regime in which the  ′ width is still below the experimental resolution of about 100 GeV,  ′ masses up to 2.3 TeV and between 2.8 and 3.5 TeV can be excluded for model C, a sensitivity which is close to the one obtained by the dĳet resonance analysis which sets a limit for   = 0.15 at around 3.5 TeV.   Figure 7: Observed and expected limits at 95% CL on the cross-section times branching ratio for the production of dark quarks through a  ′ as a function of the  ′ mass for models A, B, C and D (see the text for details).The darker and lighter shaded bands around the expected limits represent the ±1 and ±2 uncertainty range, respectively.The predicted cross-section times branching ratio is also shown as a thin solid line for   = 0.05 and    = 0.2 for models A and D, while values of   = 0.15 and    = 0.5 are used for models B and C.

Conclusion
A search for a new resonance decaying into a pair of dark quarks was performed using 139 fb −1 of √  = 13 TeV proton-proton collision data recorded with the ATLAS detector at the Large Hadron Collider between 2015 and 2018.The invariant mass spectra of the two highest-transverse-momentum large- jets is probed after selecting jets with a high track multiplicity using a variable which decorrelates this selection from the  JJ distribution.No significant excess is observed above the Standard Model multĳet background whose shape is obtained using a data-driven method from a low-track multiplicity control region.Constraints on the production cross-section times branching ratio to dark quarks are presented for four different dark QCD scenarios as a function of the  ′ mass, excluding masses of up to 3.0 TeV for   = 0.05 and    = 0.2, depending on the model.

Figure 1 :
Figure 1: Resonant dark quark pair production through a  ′ mediator.

Figure 2 :
Figure 2: Distribution of the dĳet invariant mass,  JJ for the data, the simulated multĳet background and of some representative signals (models A, B, C and D with   ′ = 2.5 TeV), shown after applying the preselections described in the text.The simulated background is normalised to the data and the signals are normalised to a production cross-section of 10 fb.

Figure 3 :
Figure 3: Distributions of the number of tracks associated with the two leading jets,  track,1 (a) and  track,2 (b), for the data, the simulated multĳet background and of some representative signals (models A, B, C and D with   ′ = 2.5 TeV), shown after applying the preselections described in the text.All distributions are normalised to unity.The uncertainty band around the background prediction corresponds to the modelling uncertainty described in Section 6.A red arrow indicates an out-of-range data point for a given bin.

Figure 5 :
Figure5: Distribution of  JJ in the CR (black dots) compared with that in the VR (purple dots), with their respective statistical uncertainties, after normalising the total number of events in each region to the one expected in the SR.The uncertainties in the data points correspond to the uncertainties on the number of events in the CR and the VR, and the shaded grey band shows the statistical uncertainty expected in the SR.

Figure 6 :
Figure6: Distribution of  JJ in the SR (black points) compared to the background-only expectation, where the shape is derived from the CR and the normalisation is determined from a fit in the SR.The most discrepant interval found by the BumpHunter algorithm, indicated by the vertical lines, is between 1500 GeV <  JJ < 1700 GeV, with a p 0 -value of 0.63.This p 0 -value and the bin-by-bin significance shown in the lower panel are computed considering statistical uncertainties only.

Table 2 :
Impact of the largest relative systematic uncertainties, Δ/ (in percentage), for the various models and for  ′  = 2.5 TeV.