Search for Higgs boson decays into two new low-mass spin-0 particles in the 4$b$ channel with the ATLAS detector using $pp$ collisions at $\sqrt{s}= 13$ TeV

This paper describes a search for beyond the Standard Model decays of the Higgs boson into a pair of new spin-0 particles subsequently decaying into $b$-quark pairs, $H \rightarrow aa \rightarrow (b\bar{b})(b\bar{b})$, using proton-proton collision data collected by the ATLAS detector at the Large Hadron Collider at center-of-mass energy $\sqrt{s}=13$ TeV. This search focuses on the regime where the decay products are collimated and in the range $15 \leq m_a \leq 30$ GeV and is complementary to a previous search in the same final state targeting the regime where the decay products are well separated and in the range $20 \leq m_a \leq 60$ GeV. A novel strategy for the identification of the $a \rightarrow b\bar{b}$ decays is deployed to enhance the efficiency for topologies with small separation angles. The search is performed with 36 fb$^{-1}$ of integrated luminosity collected in 2015 and 2016 and sets upper limits on the production cross-section of $H \rightarrow aa \rightarrow (b\bar{b})(b\bar{b})$, where the Higgs boson is produced in association with a $Z$ boson.


Introduction
The Higgs boson is a particle with a particularly narrow natural width, and its branching fractions to new light particles can be sizable even if they interact weakly with it.Because of this, several new weakly interacting light particles that would not be visible in inclusive searches can be probed by searching for "beyond the Standard Model" (BSM) Higgs boson decays at the LHC [1].These new light particles are predicted in several BSM theories with extended Higgs sectors [2][3][4][5][6] that address open questions in high-energy physics.Theories with new light particles weakly coupled to the Higgs boson provide an explanation for electroweak baryogenesis [7,8] and contain fields that mediate interactions between Standard Model (SM) particles and dark matter [9][10][11][12][13].This paper presents a search for a new spin-0 singlet a that couples to the SM Higgs boson.
When the mass of the spin-0, m a , is less than half of the mass of the Higgs boson, m H , i.e. 2m a < m H , the decay H → aa is kinematically allowed.The search in this paper is performed with events in which each a-boson decays into a pair of b-quarks, and the Higgs boson is produced in association with a Z boson which decays into electrons or muons.The final state with multiple b-quarks has the highest branching ratio in several BSM theories when it is kinematically accessible.The Z boson with leptonic decay provides a simple strategy for triggering and selecting events, as well as powerful background rejection.Figure 1 depicts the main production mechanism of the events sought in this paper.The Higgs boson has been observed by the ATLAS and CMS collaborations [14,15].A comprehensive program is being pursued to measure its branching ratios to SM particles and to search for decays into exotic or non-SM particles.Current measurements constrain the non-SM branching ratio of the Higgs boson to be less than approximately 21% at 95% confidence level (C.L.) with several assumptions [16], leaving enough room for exotic Higgs boson decays.
ATLAS has previously performed a search where each of the four b-quarks was experimentally identified as an individual jet in the detector [17].The search set upper limits on the production cross-section of Z H, followed by H → aa → (b b)(b b), of approximately 0.5 pb at 95% C.L. for m a 30 GeV.However, when the mass of the a-boson is small, it is produced with large momentum and the jets created in the hadronization of the two b-quarks from a single a → b b decay are reconstructed as a single jet in the calorimeter using the standard ATLAS reconstruction algorithms.Because of this, the previous search that covered the range 20 ≤ m a ≤ 60 GeV rapidly loses efficiency for masses m a 30 GeV.This article extends the previous analysis in the mass regime 15 ≤ m a ≤ 30 GeV by relying on a novel strategy for the reconstruction and identification of a → b b decays.The article is structured as follows.Section 2 describes the relevant features of the ATLAS detector.Section 3 lists the data collected for this search and details the simulated signal and background event samples that were used to describe the composition of the selected events.Section 4 describes the basic reconstruction and identification of leptons and jets using the ATLAS detector.Section 5 presents the dedicated method for the reconstruction and identification of low-mass a → b b decays.Section 6 explains the strategy for event selection and categorization.Section 7 discusses the systematic uncertainties considered in this search, and Section 8 presents the results.Finally, Section 9 presents the conclusion.

ATLAS detector
The ATLAS experiment [18][19][20] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4π coverage in solid angle.1It 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 sampling calorimeter provides hadronic energy measurements in the central pseudorapidity range (|η| < 1.7).The endcap and forward regions are instrumented with LAr calorimeters for EM and hadronic energy measurements up to |η| = 4.9.The muon spectrometer surrounds the calorimeters and is based on three large air-core toroidal superconducting 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 custom hardware and uses a subset of the detector information to keep the accepted rate below 100 kHz.This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions.

Dataset and simulated event samples
Events are selected from proton-proton (pp) collisions collected by the ATLAS detector at the LHC at √ s = 13 TeV in 2015 and 2016.Only collisions recorded when all relevant subsystems were operational are considered in the analysis.The dataset corresponds to an integrated luminosity of 3.2 ± 0.1 fb −1 recorded in 2015 and 32.9 ± 0.7 fb −1 recorded in 2016, for a total of 36.1 ± 0.8 fb −1 [21].The uncertainty is obtained from the primary luminosity measurements using the LUCID-2 detector [22].The data used for this search were collected using the single-electron or single-muon triggers with transverse momentum (p T ) thresholds of 20 (26) GeV for muons and 24 (26) GeV for electrons in 2015 (2016) [23].
Simulated event samples are used to study the characteristics of signal events and to calculate the signal efficiency and acceptance, as well as for most aspects of the background estimation.Monte Carlo (MC) samples were produced using the full ATLAS detector simulation [24] based on G 4 [25].To simulate the effects of simultaneous inelastic collisions (pileup), additional interactions were generated using P 8.186 [26] with the A2 set of tuned parameters [27] and the MSTW2008LO [28] parton distribution function (PDF) set, and overlaid on the simulated hard-scatter event.Simulated events were reweighted to match the pileup conditions observed in the data.All simulated events are processed through the same reconstruction algorithms and analysis chain as the data.Decays of band c-hadrons were performed by E G v1.2.0 [29], except in events simulated with the S event generator [30].
Signal samples of associated Higgs boson production with a Z boson, pp → Z H, were generated with P -B v2 [31][32][33][34] using the CT10 PDF set [35] at next-to-leading order (NLO).The sample include gluon-initiated processes at LO.The Higgs boson decay into two spin-0 a-bosons and the subsequent decay of each a-boson into a pair of b-quarks were simulated with P 8.186.The a-boson decay was performed in the narrow-width approximation and the coupling to the b-quarks is assumed to be that of a pseudoscalar.The information about the parity of the a-boson assumed in the simulation is lost in the hadronization of the b-quarks and, therefore, the results of this search apply equally to scalars and pseudoscalars.P 8.186 was also used for the parton showering, hadronization, and underlying-event simulation with the A14 tune [36].Signal events were generated for several a-boson mass hypothesis: 15, 17.5, 20, 22.5, 25, 27.5, and 30 GeV.
The background samples were generated following exactly the same procedure as described in Ref. [17] and only a summarized description is given here.A sample of top-quark pair events was generated using P -B v2 [37] with the NNPDF3.0NLOPDF set.The parton showers and hadronization were modeled by P 8.210 [38] with the A14 tune.To model the t t + b b background with better precision, the relative contributions of the different heavy-flavor categories in the t t sample are scaled to match the predictions of an NLO t t + b b sample including parton showering and hadronization [39], generated with S +O L [30,40], using the procedure described in Ref. [41].
The production of Z bosons in association with jets was simulated with S 2.2.1 [30,42] using the NNPDF3.0NNLOPDF set [43].The matrix element calculation was performed with C [44] and O L [40] and was matched using the MEPS@NLO prescription [45].
Several subleading backgrounds were also simulated.The diboson + jets samples were generated using S 2.1.1 [46] and the CT10 PDF set.Associated production of t tW and t t Z (t tV) were generated with an NLO matrix element using M 5_aMC@NLO interfaced to P 8.210 and the NNPDF3.0NNLOPDF set.Samples of Wt single-top-quark backgrounds were generated with P -B v1 at NLO accuracy using the CT10 PDF set.The production of four top quarks (t tt t) and t tWW was simulated with M 5_aMC@NLO at LO accuracy and interfaced to P 8.186.Background sources with non-prompt leptons contribute negligibly to this search.
Multijet samples are used to compare the data identification efficiency of the a → b b decays with simulation.These samples were generated using P 8.186, with the LO NNPDF2.3PDF set and the A14 tune.To increase the number of simulated events with semileptonically decaying hadrons used in this analysis, samples of multijet events filtered to have at least one muon with p T above 3 GeV and |η| < 2.8 were produced with P using the same version, PDF set, and underlying-event tunes as the unfiltered multijet samples.Both the filtered and unfiltered multijet samples produced with P were processed through the same ATLAS detector simulation.

Object reconstruction and selection
This search relies on the efficient reconstruction of electrons and muons in order to identify leptonically decaying Z bosons and the reconstruction of jets to identify a → b b decays.
Electrons are reconstructed from energy deposited in clusters of cells in the electromagnetic calorimeter matched to tracks in the inner detector [47] and are required to have p T > 15 GeV and |η| < 2.47.Candidates in the transition region between the barrel and endcap calorimeters, 1.37 < |η| < 1.52, are excluded.Electrons are identified using the "Tight" criterion based on a likelihood discriminant [48].Muons are reconstructed by combining matching tracks in the inner detector and the muon spectrometer, and are required to have p T > 10 GeV and |η| < 2.4.Muon candidates must satisfy the "Medium" identification criterion [49].An isolation requirement based on the momentum of the tracks and the calorimeter energy around each lepton candidate is imposed to distinguish between leptons coming from the decay of a Z boson and those from non-prompt sources [48,49].Additionally, all lepton candidates are required to be consistent with the primary vertex, chosen as the reconstructed vertex with the highest sum of the p 2 T of its associated tracks.Jets are reconstructed from three-dimensional topological energy clusters [50] in the calorimeter using the anti-k t jet algorithm [51] implemented in the FastJet package [52] with a radius parameter of 0.4.Jets are calibrated using energy-and η-dependent corrections [53] and are required to have p T > 20 GeV and |η| < 2.5.Events containing jets arising from non-collision sources or detector noise are removed [54].
Finally, a track-based criterion, the jet vertex tagger (JVT), is used to reduce contributions from jets arising from pileup [55].In the region |η| < 2.5, jets are tagged as containing b-hadrons using a multivariate discriminant (MV2) score [56].The MV2 score is obtained from a boosted decision tree (BDT) that combines several algorithms that identify tracks with large impact parameters, secondary vertices, and the topological structure of weak band c-hadron decays inside jets.The BDT was trained using jets reconstructed with the anti-k t algorithm with a radius parameter R = 0.4 from t t simulated events to discriminate b-jets from c-jets and light-flavor jets [57].In this search, the same BDT is used with a novel strategy described in Section 5.1.

Reconstruction and identification of a → b b decays
For low-mass a-bosons, the b-quarks from a-boson decay tend to have small angular separation ∆R and can be reconstructed either as a single jet or as multiple jets in the calorimeter depending on their angular separation and the clustering algorithm used.In order to include both cases, all calibrated jets reconstructed using the anti-k t jet algorithm with radius parameter R = 0.4 and p T > 20 GeV are clustered again, using an anti-k t algorithm with radius parameter R = 0.8 [58].The radius parameter was chosen to optimize the signal acceptance in the mass range considered.Each R = 0.8 jet is considered as a reconstructed a → b b candidate.The R = 0.8 jet will often contain a single anti-k t constituent jet with radius parameter R = 0.4 when the angular separation ∆R between the b-quarks from the a → b b decay is less than 0.4.The four-momentum of an a → b b candidate is the sum of all four-momenta of the set of constituent R = 0.4 jets.Since the R = 0.4 constituent jets are calibrated, no additional momentum calibration is necessary.
The hadronization of the two b-quarks which come from an a-boson decay is identified using variables sensitive to the number of b-hadrons and the mass of the a-boson.The values of these variables are calculated using tracks with p T > 0.4 GeV matched to the reconstructed a → b b candidate.The matching is performed using the ghost-association method [59], which treats the tracks as four-vectors of infinitesimal magnitude during the jet reconstruction and assigns them to the a → b b candidate with which they are clustered.Tracks from the hadronization of different b-quarks are separated by splitting the set of tracks matched to an a → b b candidate into multiple track-jets.Ideally, the decay of each b-quark should be associated with a different track-jet.In this search, the track-jets are reconstructed by clustering all matched tracks using an exclusive-k t algorithm that produces either two (Exk (2)  t ) or three (Exk (3)  t ) final jets [60].The exclusive-k t algorithm implements a sequential clustering in which the two tracks with the smallest k t -distance, defined as the product of the minimum p T of the two tracks and their distance ∆R, are clustered together if this distance is smaller than the transverse momentum of all tracks.If two tracks are clustered together, their momenta is summed and the two considered as a single object in the next iteration of the sequential clustering.If the transverse momentum of a track is smaller than all k t -distances, the track is discarded.Tracks clustered together are considered a final-state track-jet.The sequential clustering is interrupted after the step in which all the tracks have been clustered in the desired number of final-state track-jets [61].The splitting into three final jets attempts to capture cases where significant additional radiation is present.The strategy presented here to identify the two b-quark flight directions as different track-jets differs from the method documented in Ref. [62], where the inclusive version of the anti-k t algorithm is used.At low a-boson momenta, the exclusive-k t algorithm is able to identify the two b-quark flight directions in separate track-jets more often than the inclusive anti-k t algorithm.For a simulated signal event sample with m a = 20 GeV, the inclusive anti-k t algorithm associates the b-quark flight directions with different track-jets in 46% of cases.In contrast, the flight directions are associated with different exclusive-k t track-jets in nearly 100% of cases.
The variables used for the a → b b identification are calculated using the exclusive-k t track-jets.For the track-jets calculated with the Exk (2)  t algorithm, the variables used are the MV2 scores of the two track-jets, as well as their angular separation ∆R and their p T asymmetry, defined as (p 1 T − p 2 T )/(p 1 T + p 2 T ).For Exk (3)   t track-jets, the same variables are used, but they are calculated with the two track-jets with highest and lowest MV2 scores among the three track-jets.The eight variables are used simultaneously.The MV2 scores identify the presence of a b-hadron in the track-jets.Track-jets with large ∆R separation occur in the decay of a massive state.Track-jets with very large p T asymmetry can arise from final-state radiation.The variables calculated with Exk (2)  t track-jets provide most of the discriminating power between signal and background, while the variables calculated in Exk (3)  t help disentangle cases where Exk (2)  t would fail to identify the flight direction of the a → b b decay products.
A BDT is trained with these variables to obtain an efficient identification criterion that distinguishes a → b b candidates in signal events that have two b-quarks produced in the decay of a low-mass resonance, from those in top-quark pair events that contain a single b-quark decay.A sample of simulated SM t t events is used as a source of a → b b candidates with a single b-quark decay, while a simulated signal event sample with m a = 20 GeV is used as a source of a → b b candidates with two b-quark decays.The transverse momentum and angular distributions are not included as inputs for the BDT training, but the differences in these distributions among signal and background are partially taken into account since they are correlated with other variables.In order to classify the b-quark multiplicity of an a → b b candidate, b-hadrons in the simulation of the b-quark hadronization with p T > 5 GeV are matched to the candidates using the same ghost-association method as described above.Figure 2 shows the predicted efficiency for signal and background events using the trained BDT.The BDT discriminator is also efficient in rejecting events without b-quarks, even if such sample was not explicitly included in the BDT training.Two event categories based on the BDT score are defined for the analysis using a tight and a loose working point (WP).A high-purity category (HPC) for a → b b candidates is selected by requiring a BDT score larger than the tight WP, while a low-purity category (LPC) is selected from candidates with BDT score between the loose and the tight WPs.The tight WP is defined by a BDT score of 0.3 while the loose WP is defined by a BDT score of 0.1.The tight WP is chosen such that it provides a background rejection 1/100 in order to reduce the backgrounds from Z + jets and t t events.The LPC contains a relatively large number of events from processes with zero or one b-quark and is used to select background-enriched samples in the search.Reconstructed In order to measure the efficiency of the identification criterion for both signal and background, a → b b candidates are categorized according to the flavor of the track-jets that are reconstructed with the Exk (2)  t algorithm, while the Exk spectrum in signal samples and on the number of events in the multijet data sample used for the efficiency measurement.The complete procedure described below is applied independently in each transverse momentum range.

Efficiency measurement of a → b b identification
The strategy for the efficiency measurement in data closely follows that used in the measurement of the identification efficiency for boosted 125 GeV Higgs boson decays into a pair of b-quarks [62].A multijet sample is selected from a suite of single-jet triggers that differ by their jet p T threshold.Only a small fraction of the events identified by the triggers with low p T threshold are recorded.The choice of which jet events to keep is random and results in an effective integrated luminosity smaller than the total recorded by the ATLAS experiment, but does not introduce any selection bias.The fraction of events kept is known as the trigger prescale fraction.Triggers with prescale fraction less than one are called prescaled triggers.The lowest jet p T threshold for which all events are kept is 300 GeV.When comparing events recorded with prescaled and unprescaled triggers, each event is weighted by the inverse of the prescale fraction of the corresponding trigger used to record it.
The a → b b reconstruction described in Section 5.1 is applied to the multijet sample.The events recorded by the multijet triggers are dominated by LL candidates.Since muons are often produced in semileptonic decays of b-hadrons, a sample with a larger fraction of BB and BL candidates is selected by requiring exactly one muon matched to one of the Exk (2)  t track-jets.The track-jet matched to the muon is called the muon-matched track-jet, while the other one is called the non-muon-matched track-jet.The selected events are compared with simulated multijet events.In order to account for possible mismodeling of the flavor fractions in simulation relative to those in data, a correction is applied to the simulated event sample.The correction is described in detail in Ref. [62] and only a brief summary is given here.The simulated jet sample is split into subsamples depending of the flavor classification of the a → b b candidate: BB, BL, CC, CL, and LL.The selected BC fraction in the multijet sample is negligible.The fraction of each subsample is corrected by fitting the distribution of the signed transverse impact parameter significance S jet d0 = d 0 /σ(d 0 ) of the two Exk (2)  t track-jets to data.The S jet d0 of a track-jet is defined as the average of the three largest signed transverse impact parameter significances S trk d0 of its constituent tracks, since this observable is used to identify the long lifetime of band c-hadrons.The average is used to minimize the impact of misreconstructed tracks on this observable.The track impact parameter d trk 0 is calculated using the vector from the primary vertex to the point of closest approach of the track.The absolute value of d trk 0 is the norm of the projection of this vector in the transverse plane, while the sign depends on the angle between this vector and the track-jet ì p T .If this angle is less than π/2, d trk 0 is taken as positive.For angles larger than π/2, the track impact parameter is considered negative.Large negative impact parameters are often obtained from interactions with the detector material and not from a long-lived bor c-hadron decay, since the direction of the decay is not correlated with the jet axis.
A total of four flavor correction factors that scale the BB, BL, CC, and CL subsamples are determined from a Poisson likelihood fit to data.The scale parameter for the LL subsample is determined implicitly by requiring that the total number of candidates in simulation is the same as in data.The covariance matrix of these four parameters is determined from the statistical uncertainties and correlations, but also from the limited knowledge of the jet energy scale, and from the uncertainty in the impact parameter resolution following the method described in Ref. [62]. Figure 3 shows the result of this fit to data for the transverse momentum range 30 ≤ p a→b b T < 90 GeV.After the flavor correction is applied, the a → b b identification BDT is used to select events in both the HPC and LPC.Once the identification criteria are used, only the BB and the BL subsamples contribute significantly.Any residual disagreement in these regions is the result of a difference in the a → b b identification efficiency between data and simulation.A scale factor (SF) is defined as the ratio of the two efficiencies, SF = ε DATA /ε MC , for each flavor subsample.Only the BB and BL SFs are measured for both HPC and LPC.All other flavors are subleading after applying the identification criterion, and for these the efficiency in data is considered the same as in simulation.In order to measure the BB and BL SFs, in both the HPC and LPC, a second Poisson likelihood fit of the S jet d0 distribution to data is performed after using the identification BDT to select events in both simulation and data.The four SFs measured in each of the three p a→b b T ranges constitute 12 parameters in total.The complete list of uncertainties is described in Section 5.3.Figure 4 shows the measured efficiencies in both data and simulation, for BB and BL candidates.The bottom panel in the same figure shows the SF as defined above.

Systematic uncertainties in the a → b b identification
Several sources of uncertainty are considered when building the covariance matrix of the 12 SFs.The statistical uncertainties and correlations are interpreted directly from the likelihood fit to data.The impact of systematic uncertainties is considered by varying the simulated event samples within ±1σ for each source separately.The impact of each systematic uncertainty is assessed as the difference in the measured SF when fitting the nominal sample and the one with the corresponding source variation.The covariance matrix from the four flavor-fraction corrections described in Section 5.2 is propagated to the SF covariance matrix.The impact of the limited knowledge of the jet energy scale is, once again, considered in the covariance matrix.The uncertainty arising from the choice of hadronization model is included through its effect on the MV2 scores and propagated to the SF covariance matrix.This uncertainty changes the MV2 score by 5-10% depending on the track-jet p T [63] and has a minor impact on the SFs.Two additional sources of uncertainty are considered.First, there is a possible mismodeling of the efficiency for candidates with flavor other than BB or BL.These components are highly suppressed by the BDT selection.An uncertainty of 50% in the efficiency of other flavor components is propagated to the covariance matrix, with negligible impact.The chosen value of 50% is based on the level of agreement between the S jet d0 distribution in data and simulation, after the BDT criteria are applied, but before the SF fit.Second, there is a possible bias from the selection of a → b b candidates with muons.In order to assess it, the measurement is repeated by selecting a → b b candidates with two muons, one inside each track-jet, and comparing the result with the one obtained above.The sample with two muons contains a negligible number of BL candidates and it is only possible to measure the BB SFs.The difference between the SF measured with the one-muon sample and the one measured with the two-muon sample is taken as an estimate of the systematic uncertainty due to a possible bias in the procedure.The same uncertainty is applied to the BL SFs.This uncertainty varies in each p a→b b T range, but it is approximately 20% and is the dominant uncertainty in the p T ranges with a large number of candidates.Figure 5 shows the S jet d0 distribution in the 30 ≤ p a→b b T < 90 GeV range, for the HPC, after fitting for both the BB and BL SFs and including all the uncertainties described here.These SFs are used when comparing simulated signal and background events with data.For the selected background events, the distributions of the variables used for the a → b b identification are similar to the ones in g → b b events.However, for the signal, due to the nonzero mass of the a-boson, the distributions can be quite different, especially for the variables that are sensitive to the mass of the particle.The method presented here relies on the fact that any residual disagreement accounted for by the SFs is independent of the a-boson mass.In order to test this hypothesis, the efficiency measurement is repeated replacing data with a pseudo-data built using the same multijet simulated sample used to obtain the S jet d0 templates but where gluons were replaced by a spin-0 a boson with mass m a = 20 GeV before the decay to two b quarks.Figure 6 shows the results of using this pseudo-data in each of the categories considered above, which can be interpreted as the ratio between the SF for a particle with mass m a = 20 GeV and the one for a massless gluon.The overall distribution is consistent with unity within the statistical uncertainty of the simulated event sample.

Analysis strategy
The analysis strategy targets events where a Higgs boson is produced in association with a Z boson.The Higgs boson is required to decay into two a-bosons that decay into a pair of b-quarks and the Z boson decays into electrons or muons.Events are selected using triggers that require at least one electron or muon.The event is further required to have two oppositely charged electrons or two oppositely charged muons and two reconstructed a → b b candidates.The leading electron or muon is required to have p T > 27 GeV and be matched to the lepton candidate reconstructed by the trigger algorithms.The lepton momentum requirement and trigger matching are used so that all events have at least one lepton with p T above the trigger thresholds.The dilepton mass must be consistent with the Z boson mass and is required to be in the range 85 which probes the difference between the invariant mass of the two a → b b candidates, m a 1 ,a 2 , and the Higgs boson mass m H = 125 GeV.It should be noted that m a is the mass hypothesis for the a-boson.
The reduced mass is required to satisfy −40 < m red < 20 GeV, ensuring that the selection is highly efficient.The presence of m a in the event selection means that different events are used to search for different mass hypotheses.No conditions on the individual values of m a 1 and m a 2 are imposed.The selected mass window, as a function of mass difference and reduced mass, is shown in Figure 7 for signal events with m a = 20 GeV and top-quark pair events.Two signal-enriched regions are defined for this search.One requires the two reconstructed a → b b candidates to be identified in the HPC, while the other requires one a → b b candidate identified in the HPC and one in the LPC.The two main sources of background for this search are top-quark pair and Z boson events produced in association with additional quarks or gluons.In this search, the normalizations of these two backgrounds are measured in dedicated control regions which are selected to be enriched in the specific background.Three regions dominated by top-quark pair events are selected by requiring the dilepton mass to be outside the Z boson mass window, i.e. m ≤ 85 GeV or m ≥ 100 GeV.These three control regions differ by the identification of the two a → b b candidates, with one requiring both to be in the HPC, a second requiring one a → b b candidate in the HPC and one in the LPC and, finally, the third requiring the two a → b b candidates to be identified in the LPC.The three control regions probe t t events produced in association with different numbers of heavy-flavor jets.A dedicated control region for Z boson events is formed by requiring the dilepton mass to be consistent with the Z boson mass and the two a → b b candidates in the LPC. Figure 8 shows the expected background yield fractions in each of the regions described here.Regions labeled "on-Z" require the dilepton mass to be in the range 85 < m < 100 GeV, while regions "off-Z" require the dilepton mass to be outside this window.For a → b b candidates, the HPC and LPC are defined using ranges of the identification BDT score, as described in Section 5.

Systematic uncertainties
Several sources of systematic uncertainty are considered.The identification efficiency for leptons is measured in Z boson data events using a tag-and-probe method [47,49].Small residual disagreements between efficiencies in simulation and those measured in data are corrected as a function of the lepton p T and η.The uncertainties in these corrections are propagated through this search.Uncertainties in the lepton momentum scale and resolution are similarly considered.
Uncertainties associated with jets arise from their reconstruction and identification efficiencies.These are due to the uncertainty in the jet energy scale (JES), mass scale, energy resolution and the efficiency of the JVT requirement that is meant to remove jets from pileup.The JES and its uncertainty are derived by combining information from test-beam data, LHC collision data and simulation [53].
The identification efficiency for a → b b candidates in simulation is also corrected by using the SFs measured with the methods described in Section 5.The full covariance matrix for the 12 SFs is propagated after diagonalization in order to obtain uncorrelated sources of systematic uncertainty.Only a → b b candidates with BB and BL flavors have their identification efficiency corrected in simulation.Candidates with flavors other than BB and BL represent a subleading fraction of candidates selected in this analysis, mostly from BC candidates.In this case, an uncertainty of 50% per candidate is applied, similarly to the uncertainty used when measuring the identification efficiency.
Several sources of systematic uncertainty affecting the modeling of the relative normalization of the background sources in control and signal regions are considered.Since the Z +jets background normalization is measured in a region with two a → b b candidates in the LPC, where a larger fraction of the candidates do not contain two b-quarks, an uncertainty of 50% in the fraction of events which have two or more associated b-hadrons is applied.This uncertainty is derived by comparing the level of agreement between data and simulation for m red < −40 GeV calculated with m a = 20 GeV.Similarly, for the top-quark pair background, three uncorrelated relative uncertainties of 50% are assigned to events with one associated b-hadron, to events with two or more associated b-hadrons, and to events with associated c-hadrons.The number of associated hadrons in each event is determined following the procedure described in Ref. [64].These uncertainties are derived from a comparison of the t t + heavy-flavor production cross-sections predicted by P +P and by S +O L at NLO [64].
Beyond the uncertainties associated with heavy-flavor fractions, several sources of systematic uncertainty affecting the relative normalization between control and signal regions are considered.The procedure closely follows the description in Ref. [17].For the t t background, it includes systematic uncertainties from variations of the factorization and renormalization scales, the PDF set used for simulation, α S , the value of the top-quark mass, the choice of generator, the choice of parton shower and hadronization models, and the effects of initial-and final-state radiation.For the Z + jets backgrounds, additional relative uncertainties are based on variations of the factorization and renormalization scales and of the parameters used in matching the matrix element to the parton showers in the S simulation.
Uncertainties in secondary background sources are also considered, affecting their normalization in both the signal and control regions.A 50% normalization uncertainty in the diboson background is assumed [65].
The uncertainties in the t tW and t t Z NLO cross-section predictions are 13% and 12%, respectively [66,67], and are treated as uncorrelated between the two processes.An additional modeling uncertainty for t tW and t t Z, related to the choice of event generator, parton shower and hadronization models, is derived from comparisons of the nominal samples with alternative ones generated with S .
Several sources of systematic uncertainty affect the theoretical modeling of the signal.Uncertainties originate from the choice of PDFs, the factorization and renormalization scales, and the parton shower, hadronization and underlying-event models.The combined uncertainty in the expected signal yield from these sources is approximately 8%.Higher-order corrections to the decay of the a-boson are small compared to the Higgs boson production uncertainties and, therefore, no additional uncertainty is included.

Results
The results are obtained from a binned maximum-likelihood fit to the data using the two signal regions and four control regions.The likelihood function is constructed from the product of Poisson probabilities in each region.The parameter of interest (POI) scales the signal H → aa → (b b)(b b) yield.The overall normalizations of the Z + jets and t t backgrounds are modeled as unconstrained nuisance parameters.Simulation is used to determine the relative yields of Z + jets and t t backgrounds in each signal and control region.Systematic uncertainties described in Section 7 are incorporated as nuisance parameters with Gaussian priors with a standard deviation equal to the value of the uncertainty, and these nuisance parameters multiply the product of Poisson probabilities.Uncertainties arising from the finite number of simulated events are modeled using gamma distribution priors [68].Gamma distributions are used as a generalization of the Poisson distribution since the expected yield predicted in simulated event samples may not be an integer number.Figure 9 and Table 1 show a comparison of data and simulation when the nuisance parameters have the values that maximize the likelihood function and only the SM background processes are considered, i.e. the POI is fixed at zero.The data in all regions agrees with the prediction within 1 standard deviation.The data observed in each region is included for comparison.The hashed area represents the total uncertainty in the background.
Table 1: Expected yields and total uncertainty for the different background components in each signal and control region after the profile likelihood fit to data.The expected yield for signal with m a = 20 GeV is calculated before the profile likelihood fit and normalized to the total Z H cross-section (B(H → aa → (b b)(b b)) = 1).The data observed in each region is included for comparison.

Signal Regions
Control Regions on-Z, 1HPC1LPC on-Z, 2HPC on-Z, 2LPC off-Z, 2HPC off-Z, 2LPC off-Z, 1HPC1LPC Limits on the production cross-section of Z H, H → aa → (b b)(b b) events are calculated using the test statistic t µ = −2 ln(L(µ, θµ )/L( μ, θ)), where L is the likelihood described above, µ is the single POI and θ is the vector of nuisance parameters (NPs).In addition, μ and θ are the values which maximize the likelihood function, and θµ are the values of the NPs which maximize the likelihood function for a given value of µ [17].Upper limits at 95% C.L. on the production cross-section as a function of the mass hypothesis are determined using the asymptotic distribution for t µ [69][70][71].
The impact of systematic uncertainties on the upper limits is evaluated by varying the corresponding NP when building the Asimov dataset [69] used to estimate the asymptotic distribution for t µ .The NPs are varied by the value of their uncertainties in the fit performed to obtain L( μ, θ).In order to partially account for the correlation between the fitted values of the NPs, the variations are performed after diagonalizing the correlation matrix obtained in the same fit.The diagonalization is performed in blocks of NPs that share a similar origin and that may have large correlations.The impact is defined as the relative variation of the expected upper limit when the modified asymptotic distribution is used.Variations in each block are summed in quadrature and the results are shown in Table 2.The number of events in each of the four control regions is the main factor in determining the impact from the unconstrained nuisance parameters that model the normalization of the Z + jets and t t backgrounds and, therefore, their values are highly correlated.Since they are individually important for the modeling of the background yields, their impacts are reported separately.A correlation of 44% is observed between the two unconstrained nuisance parameters.The impact of the statistical uncertainty is defined as the 1σ uncertainty in the expected upper limit after removing all nuisance parameters, both constrained and unconstrained, from the profile likelihood.
Figure 10 shows the exclusion limits for the production cross-section times the branching ratio for Z H, H → aa → (b b)(b b) as a function of the a-boson mass hypothesis.For comparison, the SM NNLO cross-section for pp → Z H is σ SM (pp → Z H) = 0.88 pb [66].The figure also includes the expected exclusion limit calculated from an Asimov dataset when all the constrained nuisance parameters are fixed to their expected values and the unconstrained nuisance parameters that scale the Z + jets and t t backgrounds are fixed to one.For m a = 20 GeV, an upper limit of 0.71 pb (0.52 +0.31 −0.14 pb) is observed (expected) at 95% C.L. The reduced sensitivity for heavier a-boson mass hypotheses is due to a lower acceptance caused by the increased separation of the b-jets, while the reduced sensitivity for lighter a-boson mass hypotheses is due to a lower efficiency to identify the two b-jets inside an a → b b candidate.The figure includes the results from a previous analysis targeting the higher range of m a [17].)) from this result and a previous search targeting the higher masses (dash-dotted line) [17].For the results of this search, the observed limits are shown, together with the expected limits (dashed line).In the case of the expected limits, one-and two-standard-deviation uncertainty bands are also displayed.The SM NNLO cross-section for pp → Z H of 0.88 pb [66] is also shown.for a mass hypothesis of m a = 20 GeV by a factor of 2.5 when compared with the previous ATLAS result which uses the same integrated luminosity.

Figure 1 :
Figure 1: Representative tree-level Feynman diagram for the Z H production processes with the subsequent decays Z → ( = e, µ) and H → aa → (b b)(b b).

Figure 2 :
Figure 2: The identification efficiency of a → b b candidates with m a = 20 GeV as a function of the inverse of the t t background efficiency (rejection).For the signal sample, both b-quarks are required to lie within the reconstructed candidate, while for the background t t sample the reconstructed candidate contains a single b-quark.The left and right stars indicate the tight and the loose WPs, respectively, which are used to define, as described in the text, the HPC and LPC.

( 3 )
t track-jets are used exclusively for identification purposes.All band c-hadrons present in the event simulation with p T > 5 GeV are matched to the track-jets using the ghost-association method.If a track-jet has at least one simulated b-hadron matched to it, it is classified as B. If it does not contain a simulated b-hadron, but has a simulated c-hadron matched to it, it is classified as C. Otherwise it is classified as L. The flavor of an a → b b candidate is determined by the flavor of the two Exk (2) t jets.Most signal a → b b candidates are BB candidates, while most background candidates are BL candidates.A signal candidate can be classified as BL when the two b-quarks decay inside the same track-jet or when they have p T ≤ 5 GeV.The identification efficiencies for BB and BL a → b b candidates are measured separately in data for three transverse momentum ranges: b T < 90 GeV, 90 ≤ p a→b b T < 140 GeV, and p a→b b T ≥ 140 GeV.These three ranges were chosen based on the p a→b b T

Figure 3 :
Figure 3: Averaged signed impact parameter significance S jet d0 distributions of the track-jet (a) with a muon inside and (b) without a muon, after performing the fit of the flavor fractions to data.The total simulation yield is scaled to the same number of events observed in data.The fit is performed separately in p a→b b T ranges.The figure shows the 30 ≤ p a→b b T < 90 GeV range of a → b b candidates.

Figure 4 :
Figure 4: Efficiency of the a → b b identification criteria measured in data and simulated multijet events.The efficiency is measured in three transverse momentum ranges, separately for (a) BB and (b) BL candidates in both the HPC and LPC.The ratio of the measured values in data and simulation (bottom panels) are SFs used in the analysis when comparing simulation with data.The error bars in the top panels are statistical only, while the hashed bands on the ratios in the bottom panels include the full systematic uncertainties.

Figure 5 :
Figure 5: Averaged signed impact parameter significance S jet d0 distributions of the track-jet (a) with a muon inside and (b) without a muon, in the 30 ≤ p a→b b T < 90 GeV range of a → b b candidates in the HPC.The hashed area represents the total uncertainty in the predicted yields.

Figure 6 :
Figure 6: Efficiency measured using a simulated multijet pseudo-data where g → b b decays are replaced by a → b b decays with mass m a = 20 GeV.The efficiency is measured separately for BB and BL samples in the same p T ranges used in the data-to-simulation SF measurement.The values can be interpreted as the ratio between the SFs for a particle with mass m a = 20 GeV and the ones for a massless gluon.Only the statistical uncertainties are indicated.
< m < 100 GeV.Both a → b b candidates are required to satisfy p T > 30 GeV and |η| < 2.0.Two mass requirements are imposed to select events consistent with a cascade decay H → aa → (b b)(b b).Firstly, the mass difference ∆m a→b b = m a 1 − m a 2 between the two a → b b candidates is required to satisfy −25 < ∆m a→b b < 25 GeV.The ordering of a → b b candidates is based on their transverse momenta, with a 1 corresponding to the higher-p T a → b b candidate.Second, the mass of the pair of a → b b candidates is required to be consistent with the Higgs boson mass.The compatibility is assessed with the reduced mass:

Figure 7 :
Figure 7: Distribution of (a) expected signal events and (b) top-quark pair background in a plane defined by the two mass requirements described in the text, m red and ∆m a→b b .The mass requirements aim at selecting events where the two a → b b candidates have similar reconstructed masses and the mass of the pair of a → b b candidates is consistent with the Higgs boson mass.The signal events correspond to Z H, H → aa → (b b)(b b) with m a = 20 GeV and are normalized to the SM pp → Z H cross-section.

Figure 8 :
Figure 8: Expected composition of events in each of the signal and control regions defined for the search.Control regions (CR) have a negligible expected yield for the signal.Definitions of the regions are based on the dilepton mass and the purity of the two a → b b candidates.Regions labeled "on-Z" require the dilepton mass to be in the range 85 < m < 100 GeV, while regions "off-Z" require the dilepton mass to be outside this window.For a → b b candidates, the HPC and LPC are defined using ranges of the identification BDT score, as described in Section 5.

S
R o n -Z ,1 H P C 1 L P C S R o n -Z ,2 H P C C R o n -Z , 2 L P C C R o ff -Z , 2 H P C C R o ff -Z , 2 L P C C R o ff -Z , 1 H P C 1

Figure 9 :
Figure 9: Expected yields for the different background components in each signal and control region after the profile likelihood fit to data.The expected yield for signal with m a = 20 GeV is calculated before the profile likelihood fit and normalized to the observed limit in the cross-section times the branching ratio for Z H, H → aa → (b b)(b b).The data observed in each region is included for comparison.The hashed area represents the total uncertainty in the background.

Figure 10 :
Figure 10: Summary of the 95% C.L. upper limits on σ Z H B(H → aa → (b b)(b b)) from this result and a previous search targeting the higher masses (dash-dotted line)[17].For the results of this search, the observed limits are shown, together with the expected limits (dashed line).In the case of the expected limits, one-and two-standard-deviation uncertainty bands are also displayed.The SM NNLO cross-section for pp → Z H of 0.88 pb[66] is also shown.

A
search for Higgs bosons decaying into a pair of new spin-0 particles that subsequently decay into a final state with four b-quarks is presented.The search uses 36 fb −1 of 13 TeV proton-proton collision data collected by the ATLAS detector at the LHC.A dedicated strategy for reconstruction and identification of a → b b candidates in the mass range 15 ≤ m a ≤ 30 GeV is introduced.The measurement of the acceptance and efficiency of this strategy is described in detail and used to compare data with simulated events in regions with two a → b b candidates consistent with the cascade decay H → aa → (b b)(b b).The dominant background sources are measured in control regions defined by relaxing some of the identification criteria.No excess of data events consistent with H → aa → (b b)(b b) is observed, and upper limits at 95% C.L. on the production cross-section σ Z H B(H → aa → (b b)(b b)) are obtained as a function of the a-boson mass hypothesis.This novel search improves the expected limit on σ Z H B(H → aa → (b b)(b b))

Table 2 :
Impact of groups of systematic uncertainties on the expected upper limits for m a = 20 GeV.For comparison, the statistical uncertainty impact, defined as the 1σ uncertainty of the expected upper limit after removing all nuisance parameters from the profile likelihood, is also shown.
64] ATLAS Collaboration, Search for the Standard Model Higgs boson produced in association with top quarks and decaying into a b b pair in pp collisions at √ s = 13 TeV with the ATLAS detector, Phys.Rev. D 97 (2018) 072016, arXiv: 1712.08895[hep-ex].Also at Centro Studi e Ricerche Enrico Fermi; Italy.c Also at CERN, Geneva; Switzerland.d Also at CPPM, Aix-Marseille Université, CNRS/IN2P3, Marseille; France.e Also at Département de Physique Nucléaire et Corpusculaire, Université de Genève, Genève; Switzerland.f Also at Departament de Fisica de la Universitat Autonoma de Barcelona, Barcelona; Spain.g Also at Department of Financial and Management Engineering, University of the Aegean, Chios; Greece.h Also at Department of Physics and Astronomy, Michigan State University, East Lansing MI; United States of America.i Also at Department of Physics and Astronomy, University of Louisville, Louisville, KY; United States of America.j Also at Department of Physics, Ben Gurion University of the Negev, Beer Sheva; Israel.k Also at Department of Physics, California State University, East Bay; United States of America.l Also at Department of Physics, California State University, Fresno; United States of America.m Also at Department of Physics, California State University, Sacramento; United States of America.n Also at Department of Physics, King's College London, London; United Kingdom.o Also at Department of Physics, St. Petersburg State Polytechnical University, St. Petersburg; Russia.p Also at Department of Physics, University of Fribourg, Fribourg; Switzerland.q Also at Dipartimento di Matematica, Informatica e Fisica, Università di Udine, Udine; Italy.r Also at Faculty of Physics, M.V. Lomonosov Moscow State University, Moscow; Russia.s Also at Giresun University, Faculty of Engineering, Giresun; Turkey.t Also at Graduate School of Science, Osaka University, Osaka; Japan.u Also at Hellenic Open University, Patras; Greece.v Also at IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405, Orsay; France.w Also at Institucio Catalana de Recerca i Estudis Avancats, ICREA, Barcelona; Spain.x Also at Institut für Experimentalphysik, Universität Hamburg, Hamburg; Germany.y Also at Institute for Mathematics, Astrophysics and Particle Physics, Radboud University Nijmegen/Nikhef, Nijmegen; Netherlands.z Also at Institute for Nuclear Research and Nuclear Energy (INRNE) of the Bulgarian Academy of Sciences, Sofia; Bulgaria.aa Also at Institute for Particle and Nuclear Physics, Wigner Research Centre for Physics, Budapest; Hungary. b