Search for Higgs boson decays into a pair of light bosons in the bb μμ ﬁnal state in pp collision at √ s = 13 TeV with the ATLAS detector

A search for decays of the Higgs boson into a pair of new spin-zero particles, H → aa , where the a -bosons decay into a b -quark pair and a muon pair, is presented. The search uses 36.1 fb − 1 of proton–proton collision data at √ s = 13 TeV recorded by the ATLAS experiment at the LHC in 2015 and 2016. No signiﬁcant deviation from the Standard Model prediction is observed. Upper limits at 95% conﬁdence level are placed on the branching ratio ( σ H / σ SM ) × B ( H → aa → bb μμ ) , ranging from 1 . 2 × 10 − 4 to 8 . 4 × 10 − 4 in the a -boson mass range of 20–60 GeV. Model-independent limits are set on the visible production cross-section times the branching ratio to the bb μμ ﬁnal state for new physics, σ vis ( X ) × B ( X → bb μμ ) , ranging from 0.1 fb to 0.73 fb for m μμ between 18 and 62 GeV. © 2018 The 3 .


Introduction
The discovery of the Standard Model (SM) Higgs boson [1,2] has opened up new avenues to search for physics beyond the SM (BSM) with perspectives to search for non-SM or "exotic" decays of the Higgs boson.Such searches could provide unique access to hidden-sector particles that are singlets under the SM gauge transformations [3].Exotic decays of the Higgs boson are predicted by many new-physics models [3,4], including those with an extended Higgs sector [5][6][7][8][9], dark matter (DM) models [10][11][12][13][14], models with a first-order electroweak phase transition [15,16] and theories with neutral naturalness [17,18].These models have also been used to explain the observations of a γ-ray excess from the galactic centre (GC) by the Fermi Large Area Telescope [19,20].For example, a model for the GC γ-ray excess was proposed in which 30 GeV DM particles pair-annihilate dominantly through a CP-odd scalar mediator that subsequently decays into SM fermions [13].If the mediator is sufficiently lighter than the SM Higgs boson (H) then H decay into the mediator pair can be observed at the LHC.
Existing measurements constrain the BSM or "exotic" branching ratio (B) of the 125 GeV Higgs boson decays to less than approximately 34% at 95% confidence level [21].Due to the narrow width (∼ 4 MeV) of the Higgs boson, even a small non-SM coupling of O(10 −2 ) can lead to O(10%) branching ratio into BSM states.This potentially large B(H → BSM states) motivates direct searches for exotic H decays. Searches for the Higgs boson with a mass of 125 GeV decaying into two spin-zero particles, H → aa, have been performed in various final states in ATLAS and CMS [22][23][24][25][26][27][28].The analysis presented in this Letter performs the search in the bbµµ final state.The a-boson can be either a scalar or a pseudoscalar under parity transformations, since the decay mode considered in this search is not sensitive to the difference in coupling.Assuming that the a-boson mixes with the SM Higgs boson and inherits its Yukawa couplings to fermions, the largest branching ratio is expected to be to the heaviest fermions accessible by kinematics (2m a < m H ), where m a and m H are the a-boson and Higgs boson masses.For m a 10 GeV this means the a-boson would decay preferentially into bb.However, in models with enhanced lepton couplings such as the Type-III 2HDM [29], the a → µµ branching ratio can also be relatively large.Additionally, the sensitivity of a given channel does not depend only on the expected signal rate in a particular model, but also on the efficiency for triggering and reconstructing events of interest.The presence of a clean dimuon resonance provides a distinctive signature that can be used for triggering and precision mass reconstruction, which helps to suppress background.

Data and simulation
The search presented in this Letter is based on the 36.1 fb −1 dataset of proton-proton collisions at a centreof-mass energy of √ s = 13 TeV recorded by the ATLAS experiment at the LHC during 2015 and 2016.
The ATLAS experiment [30] 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 hadronic calorimeters, and a muon spectrometer.Events are collected with single-muon triggers requiring the muon p T to be above 24 or 26 GeV, depending on the data-taking period.The trigger efficiency for the signal relative to the analysis selection is about 80%.
Simulated events are used to model the signal and SM backgrounds processes.Higgs boson production through the gluon-gluon fusion (ggF) and vector-boson fusion (VBF) processes was modelled at next-toleading order (NLO) using P -B v2 [31][32][33] interfaced with P 8.186 [34] using the AZNLO set of tuned parameters [35] for the simulation of the bbµµ decay of the Higgs boson, as well as for parton showering and hadronisation.The ggF Higgs boson production rate is normalised to the total cross-section predicted by a next-to-next-to-next-to-leading-order QCD calculation with NLO electroweak corrections applied [36][37][38][39][40].The VBF production rate is normalised to an approximate next-to-next-to-leading-order (NNLO) QCD cross-section with NLO electroweak corrections applied [41][42][43][44].Five mass points were simulated in the range m a = 20-60 GeV in steps of 10 GeV for both ggF and VBF production.S 2.2.1 [45] with the NNPDF3.0[46] set of parton distribution functions (PDF) was used for the generation of Drell-Yan, W+jets and diboson (WW, W Z, Z Z) backgrounds.Cross-sections were calculated at NNLO QCD accuracy for Z ( * ) /γ * +jets and W+jets production [47] and at NLO including LO contributions with two additional partons for the diboson processes [45,48,49].The t t and singletop-quark samples were generated with P -B v2 [32] using the CT10 PDF set [50] interfaced with P v6.428 [51] and the Perugia 2012 set of tuned parameters [52] for the parton shower.The mass of the top quark (m t ) was set to 172.5 GeV.The parameter h damp in P , used to regulate the high-p T radiation, was set to m t for improved agreement between data and simulation in the high p T region [53].The cross-section of t t was calculated at NNLO in QCD including resummation of next-to-next-to-leading logarithmic (NNLL) soft gluon terms [54,55].The cross-section for single-top-quark production was calculated with the prescriptions in Refs.[56,57].The production of t t pairs in association with W/Z bosons (denoted by t tV) was modelled with samples generated at LO using M G 5_ MC@NLO v2.2.2 [58] and showered with P v8.186.The samples are normalised to NLO cross-sections [59,60].
Additional pp collisions generated with P v8.186 were overlaid to model the effects of additional interactions in the same and neighbouring bunch crossings (pile-up) for all simulated events.The pile-up simulation used the A2 set of tuned parameters [61] and the MSTW2008LO PDF set [62].All the samples were processed through the full ATLAS detector simulation [63] based on GEANT4 [64] and processed with the same reconstruction algorithm as used for data.

Selection criteria
Interaction vertices from proton-proton collisions are reconstructed from at least two tracks with transverse momentum (p T ) larger than 0.4 GeV, and are required to be consistent with the beamspot envelope.The primary vertex (PV) is identified as the one with the largest p 2 T of associated tracks [65].Muon candidates are reconstructed using the information from the inner detector and the muon spectrometer [66].They are required to satisfy "medium" identification criteria [66], be matched to the PV and have p T > 7 GeV and |η| < 2.7.Additionally, the muons must satisfy the following criteria: the projected longitudinal impact parameter |z 0 sin θ| must be less than 0.5 mm and the ratio of the transverse impact parameter d 0 to its estimated uncertainty σ d 0 , |d 0 /σ d 0 |, must be less than 3. Finally, the selected muons must fulfil requirements on the scalar sum of p T of additional inner detector tracks and on the sum of the E T of calorimeter topological clusters [67] in a cone of size ∆R = 0.2 around the muon to ensure they satisfy "tight" isolation criteria [66].These requirements select signal muons with an identification efficiency of ∼94% and isolation efficiency ranging between ∼91% for m a = 20 GeV and ∼95% for m a = 60 GeV.Jets are reconstructed using the anti-k t algorithm [68] implemented in the F J package [69] with a radius parameter R = 0.4 applied to topological clusters of energy deposits in calorimeter cells.Jets from pile-up are suppressed with the use of tracking information as detailed in Ref. [70].All selected jets are required to have p T > 20 GeV, |η| < 2.5 and must pass quality requirements defined to minimise the impact of detector effects, beam backgrounds and cosmic rays.
Jets consistent with the hadronisation of a b-quark (b-jets) are identified using a multivariate discriminant [71,72].This analysis uses the 77% b-jet identification efficiency working point for which the purity of the b-tagged sample is approximately 95%, while the probability of misidentifying a jet initiated by a charm quark as a b-jet is approximately 16%, as determined from a sample of simulated t t events.
In order to reject non-prompt muons from the decay of hadrons within a jet, an overlap removal algorithm is applied.If a jet is found within ∆R = 0.4 of the muon candidate, the overlap is resolved in the following way: if there are more than two tracks with p T > 500 MeV associated with the jet then the muon is removed from the event, otherwise the muon is retained and the jet is removed.
The missing transverse momentum (E miss T ) used in the analysis is calculated as the magnitude of the negative vector sum ( − → p miss T ) of the transverse momenta of all selected and calibrated objects in the event and the additional "soft" term that takes into account tracks not associated with any of the these objects [73].The "soft" term is calculated from inner detector tracks matched to the PV and included to achieve a better E miss T resolution.
Events are required to have exactly two b-tagged jets with p T > 20 GeV and exactly two reconstructed muons of opposite charge, with the leading muon having p T > 27 GeV to be in the maximum-efficiency regime of the trigger and the subleading muon having p T > 7 GeV.The dimuon invariant mass (m µµ ) is required to be between 16 GeV and 64 GeV.The upper bound on m µµ is defined by the assumption that the 125 GeV Higgs boson decays into two on-shell particles of equal masses, while the lower bound is motivated by the kinematics of the a-boson decays.For lower values of m a , most of the signal jets fall below the reconstruction threshold and the jets tend to overlap geometrically in the detector so that the sensitivity of the analysis to the H → aa signal decreases.This set of selection criteria is referred to as the "preselection".
Signal events are characterised by the invariant mass of the two b-jets (m bb ) being equal, within the detector resolution, to the dimuon invariant mass and the four-object mass (m bbµµ ) being approximately 125 GeV.One side of the H → aa decay (a → µµ) is measured with approximately ten times better resolution than the other side of the decay (a → bb), as shown in Figures 1(a) and 1(b).
A kinematic-likelihood (KL) fit [74] exploiting the symmetry of H → aa decays is performed to test the compatibility of an event with the m bb m µµ hypothesis and improve the m bbµµ resolution in signal events.The KL fit finds the energies of the leading ( Êb 1 ) and subleading ( Êb 2 ) b-jets that maximise the likelihood for an event with measured leading and subleading b-jet energies E b 1 and E b 2 and with dimuon invariant mass m µµ .The likelihood is defined as follows,  bbµµ after the KL fit for events after the preselection stage, but removing the upper bound on m µµ .The t t contribution is modelled with the simulated sample normalised to the theoretical cross-section.The Drell-Yan contribution is taken from data templates (described in the text) and normalised to the total yield predicted by the Drell-Yan simulation.The signal distributions for all five simulated m a are also shown assuming the SM Higgs boson cross-section (including ggF, VBF and V H production) and B(H → aa → bbµµ) = 10%.The branching ratio in this and all subsequent figures is chosen so as to give good visibility on the plot.a width that is small compared to the m bb resolution.The transfer function W( Êb 1 , E b 1 ) is a double Gaussian probability density function derived from simulated events as a function of jet p T and η using the difference between true and reconstructed energies.The fit determines a maximum-likelihood value of L (denoted by ln(L max )), which quantifies how well an event fits to the constraints.The b-jet momenta determined by the fit are used to recompute the four-body mass denoted m KL bbµµ .As seen in Figure 1(d), the resolution of the m KL bbµµ distribution for the signal is improved by a factor of about two compared to the pre-fit m bbµµ shown in Figure 1(c), while the background shape within the m bbµµ signal peak remains almost unchanged.This allows the analysis to place tighter constraints on the difference between the reconstructed invariant mass of the bbµµ system and m H , rejecting more background events, while keeping the signal efficiency high.
Two criteria based on the kinematic likelihood fit are applied to select signal-like events and reject background events that do not fit the m bb = m µµ constraint well: |m KL bbµµ − m H | < 15 GeV and ln(L max ) > −8.Finally, the E miss T < 60 GeV requirement rejects a large portion of t t pairs where both top quarks decay semileptonically, while retaining most of the signal events.Adding these three requirements after the preselection stage defines the signal-enhanced region (SR).A search for a localised excess above the expected background is performed in multiple m µµ bins of the SR centred around the hypothesised m a .A bin width of 2 GeV is chosen for 16 < m µµ ≤ 40 GeV, 3 GeV for 40 < m µµ < 50 GeV and 4 GeV for 50 ≤ m µµ < 64 GeV respectively, in order to maximise the sensitivity.

Backgrounds
The dominant backgrounds in the signal region are Drell-Yan (DY) dimuon events in association with b-quarks and pair production of top quarks where both W bosons from top quarks decay into muons.Each of the dominant backgrounds amounts to approximately 50% of the total background in the SR.Two control regions (CR) are defined to constrain the contributions of the dominant backgrounds in the signal region.They are chosen such that they have negligible signal contamination, but are kinematically close to the SR to reduce model dependence.The top control region (TCR) is defined by applying the same selection criteria as for the signal region, but inverting the requirement on the missing transverse momentum to E miss T > 60 GeV.According to the simulation, approximately 95% of the events in TCR originate from t t production.The Higgs boson mass sidebands of the signal region are used as the Drell-Yan control region (DYCR): the constraint on the bbµµ invariant mass after the KL fit is inverted to 80 < m KL bbµµ < 110 GeV or 140 < m KL bbµµ < 170 GeV.The DYCR consists of about 50% DY events and about 50% t t events.
The shapes of the t t kinematic variables are modelled using simulated events, while the distributions for the Drell-Yan process are taken from data templates as described below.The t t simulated sample and the DY templates are normalised in profile likelihood fits to the data.In one fit variant, the two background normalisations are simultaneously determined from the event yields in the TCR and DYCR assuming no presence of signal.In a second variant, the two background normalisations and the signal strength are determined using the event yields measured in the TCR, DYCR, and a given signal window.Two validation regions are defined to compare the number of observed events with the number of SM events predicted by the fit.One validation region (VR1) is defined in the high tail of the bbµµ invariant mass distribution, 170 < m KL bbµµ < 300 GeV, while for the second validation region (VR2) only the requirement on the ln(L max ) is changed relative to the SR, −11 < ln(L) < −8.All the analysis regions are illustrated in Figure 2.
The DY templates for each of the kinematic variables considered in a particular region of the analysis (SR, CR or VR) are taken from the data in a corresponding template region (DYTR).For each analysis region the associated DYTR is defined by changing the two-b-tag requirement (present in every SR, CR and VR) to a zero-b-tag requirement, while keeping all other selection requirements the same.All the DYTR are > 90% pure in DY events.The small contribution from non-DY backgrounds, namely t t, dibosons, W+jets, single-top and t tV, is subtracted from the data in a DYTR using the simulated samples and the remaining data events are assigned to the DY template.To construct b-jet-based variables, such as m bb and m bbµµ , in a DYTR the two leading non-tagged jets are taken and used in the computation instead of the b-jets.It is verified in both the simulation and the data that the shapes of all the muon-based variables (most importantly m µµ ) are consistent between the sample with no b-tagged jets and the sample with two btagged jets.To account for differences in jet kinematics between the DYTR dominated by light-flavour jets and the corresponding analysis region dominated by heavy-flavour jets, an event-reweighting based on the leading jet p T is applied to the events in the DYTR.The event weights are derived in the data after the preselection as the ratio of the leading b-tagged jet p T in the two-b-tag sample to the leading jet p T in the sample with zero b-tags.An improvement in the modelling of jet-based kinematic variables after the reweighting is verified both in simulation and in data in the DYCR, while the shape of the m µµ distribution remains unchanged.
Minor backgrounds include diboson production, W boson production in association with b-jets (with one non-prompt muon satisfying the isolation criteria) and production of a single top quark or t t pair in association with a vector boson.The contribution of the minor backgrounds in the signal region is at the percent level.They are estimated using simulation normalised to the best available theory prediction.

Systematic uncertainties
Dominant sources of experimental systematic uncertainty are the calibration and resolution of jet energies and muon momenta, the measurement of the b-tagging efficiency and the measurement of the scale and resolution of the soft term of the missing transverse momentum.Each of these uncertainties affects the t t yields by up to 14% in any of the m µµ bins of the signal region.Other experimental uncertainties have a sub-percent effect on the expected yields.These include the uncertainties in the measurement of muon identification and isolation efficiencies and the uncertainties associated with the integrated luminosity and the simulation of pile-up interactions.The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%.It is derived, following a methodology similar to that detailed in Ref. [75], from a calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016.
Four sources of theoretical uncertainty in the modelling of the t t process are considered in the analysis.As the t t simulation is normalised to the data in TCR, all of these uncertainties are applied to the acceptance ratio between TCR and SR.Hadronisation and parton-showering model uncertainties are estimated using a sample generated with P and showered by H ++ v2.7.1 and comparing it with the nominal P sample showered with P v6.428.The uncertainty due to the choice of the event generator is estimated by comparing the expected yields obtained using a t t sample generated with MC@NLO and one that is generated with P .Both samples are showered with H ++ v2.7.1.The event generator and hadronisation/parton-showering uncertainties are found to have the largest effect among all the uncertainties affecting the total t t expectation in the signal region: ∼18% and ∼16%, respectively.Systematic uncertainties in the modelling of initial-and final-state radiation (ISR and FSR) are assessed with P samples showered with two alternative settings of P 6.428.The first of these uses the PERUGIA2012radHi tune and has the renormalisation and factorisation scales set to twice the nominal value, resulting in more radiation in the final state.In addition, it has h damp set to 2 × m t .The second sample, using the PERUGIA2012radLo tune, has h damp = m t and the renormalisation and factorisation scales are set to half of their nominal values, resulting in less radiation in the event.This uncertainty has about a 5% effect on the final t t yields.Finally, the uncertainties due to the choice of PDF are evaluated by taking the maximum difference in the acceptance ratio between TCR and SR obtained with the nominal CT10 set and the alternative PDF4LHC15 set [76].The PDF uncertainty has up to a 2% effect on the final t t yields.
The uncertainties in the theoretical cross-sections (described earlier in this Letter) are assigned to the minor backgrounds whose yields are taken directly from the simulation: dibosons (10%), single top (5%) and t tV (13%).A 100% uncertainty is applied to the W+jets process to account for the limited precision of the simulation when modelling the non-prompt muons satisfying the isolation criteria.Due to the minor contribution of the W+jets background to the analysis, this uncertainty has negligible effect.As these backgrounds have very small contributions to the SR, no theoretical uncertainties affecting the acceptance have been applied.
The systematic uncertainties applied to the data-driven DY template include the uncertainties in the shape of the template due to the background subtraction and different jet-flavour composition between the DYTR and SR.The uncertainty in the background subtraction is estimated from a comparison of the nominal template after the non-DY backgrounds are subtracted and the template where no subtraction is performed.The effect of this systematic uncertainty on the DY yields in the signal region is up to 4%.The uncertainty in the template shape due to the jet-flavour composition is assessed by comparing the nominal template extracted from the DYTR with zero b-tagged jets to the template extracted from the corresponding region, but with exactly one b-tagged jet.The average per-bin difference between the two templates in the m µµ distribution is taken as an overall uncertainty in the shape, amounting to 14%.
The systematic uncertainties affecting the acceptance of the H → aa signal that correspond to the QCD scale uncertainties, the process of parton showering and hadronization and the choice of PDF set are evaluated.The renormalisation and factorisation scales are independently varied up and down from their nominal value by a factor of two and the largest resulting change is taken as the overall uncertainty due to the QCD scale.The parton-shower uncertainties are derived by independently shifting up and down the P internal parameters that control the amount of ISR and FSR.Uncertainties due to the PDF are evaluated by taking the maximum difference between the yields obtained with the nominal PDF set and the alternative PDF4LHC15 and NNPDF3.0PDF sets.The uncertainties due to the missing higher-order QCD corrections are applied to the ggF and VBF Higgs boson production cross-sections, amounting to 3.9% and 2.1%, respectively [36,77].The uncertainties due to the choice of PDF and α S are also applied to the Higgs boson cross-section, amounting to 3.2% for ggF and 0.4% for VBF production [36,77].
Additionally, the ggF signal sample is compared with the alternative sample generated using the NNLOPS approach [78].The Higgs boson rapidity distribution in the original P signal sample is found to be consistent with the one predicted by the NNLOPS calculations, while the Higgs boson transverse momentum (p T (H)) distribution is found to be harder than the one obtained using the NNLOPS approach.A reweighting is derived as a function of p T (H) by fitting the ratio of the two generated p T distributions with a continuous function.The ggF signal sample is then reweighted with this function to obtain the nominal signal prediction.A 2.5% difference in the SR event yields observed between the weighted and unweighted sample is applied as a systematic uncertainty in the modelling of p T (H).
The signal contribution of the Higgs boson produced in the association with a vector boson (V H) is taken into account by increasing the total cross-section of the ggF and VBF processes by an estimated 3.5% V H contribution.A 100% uncertainty is applied to this procedure to account for kinematic differences between the estimated V H contribution and the generated ggF and VBF processes.The contribution from other Higgs boson production processes is minor and therefore not included.
Table 1 shows a summary of the dominant post-fit systematic uncertainties in the total background and signal yields across multiple m µµ bins of the signal region.All of the uncertainties shown in Table 1, except the normalisation and cross-section uncertainties, affect the shapes of the signal and background distributions and therefore the extrapolation of the predicted yields from the CRs to the SR.

Results
The expected SM background in each of the analysis regions is determined by a profile likelihood fit to the data.The numbers of observed and predicted events in each of the bins included in the likelihood are described by Poisson probability density functions.The systematic uncertainties are implemented as nuisance parameters constrained by Gaussian distributions with widths corresponding to the sizes of the uncertainties.
The background-only version of the fit is performed to verify that the post-fit background yields agree with the data in the VRs and SR.In this version of the fit, only the data in TCR and DYCR are used to constrain the t t and DY backgrounds and determine their normalisation factors.Both TCR and DYCR are considered as one bin each.The free fit parameters are the overall normalisation factors for the t t and Drell-Yan backgrounds.The derived t t (DY) normalisation factors are then applied to the number of t t (DY) events as predicted by the simulation (template) in any of the VR or SR bins.The post-fit distributions are shown in Figures 3-5.Both the normalisation and the shapes of the predicted background distributions describe the data well in all of the analysis control and validation regions, as well as in the SR.The post-fit yields in five m µµ bins of the SR, for which the signal sample was simulated, are shown in Table 2.
Events / 10 GeV Since no significant deviation from the predicted background is observed in the signal region, upper limits on signal yields at 95% confidence level (CL) are set as a function of m µµ using the CL s prescription [79,80].A series of profile likelihood fits is applied to the data in order to test 36 hypotheses for the m a value in steps half the size of the mass-bin width optimised in each m µµ region.In each fit the likelihood function is based on the observed and predicted yields in a SR m µµ bin corresponding to the m a hypothesis under test and on the expected and measured yields in the TCR and DYCR.The profile likelihood is maximised to extract the best-fit values for the signal strength and the t t and DY normalisation factors.Model-dependent limits are set on (σ H /σ SM ) × B(H → aa → bbµµ) assuming the signal acceptance × efficiency as given by the simulation.The signal acceptance × efficiency varies between 1.3% and 2.5% for ggF production and between 0.94 % and 3.2% for VBF Higgs boson production.To obtain the signal yield for masses for which no events were simulated, the acceptance × efficiency is interpolated with spline functions between the five simulated points.All signal-related uncertainties are taken into account in the likelihood, with an additional 3% interpolation uncertainty applied to the intermediate masses.The limits are set in the 20 ≤ m a ≤ 60 GeV range for which the signal samples were simulated and range between 2 × 10 −4 and 10 −3 (see Figure 6(a)).
A model-independent fit that does not include any prediction for the signal yields in SRs and CRs is also performed.The upper limit on the number of BSM events for each mass bin of the SR is translated to a 95% CL upper bound on the visible cross-section for new physics times branching ratio into bbµµ final state (including the KL fit constraint on m bb ∼ m µµ ), σ vis (X)×B(X → bbµµ).The visible cross-section is defined as the product of the production cross-section and acceptance × efficiency (σ vis (X) = σ prod (X)× X ) of a potential signal after all the analysis selection criteria have been applied.The limits range from 0.1 fb to 0.73 fb, depending on the dimuon mass, and are shown in Figure 6(b).The most significant excess of data over the SM prediction is found at m µµ = 38 GeV, with a local significance of 1.6 standard deviations.
[GeV] Figure 6: The (a) observed and expected upper limits at the 95% confidence level on B(H → aa → bbµµ) given the SM Higgs boson production cross-section in the ggF, VBF and V H modes and (b) model-independent upper limits on the visible cross-section for new physics times branching ratio to the bbµµ final state σ vis (X) × B(X → bbµµ).

Conclusions
In summary, a search for exotic decays of the Higgs boson into two spin-zero particles in the bbµµ final state is presented.The analysis uses 36.1 fb −1 of pp collision data collected by ATLAS during the 2015 and 2016 runs of the LHC at √ s = 13 TeV.The search for a narrow dimuon resonance is performed over the range 18 GeV ≤ m µµ ≤ 62 GeV using mass bins that are 2, 3 or 4 GeV wide depending on m µµ .No significant excess of the data above the SM prediction is observed.Upper limits are set on (σ H /σ SM ) × B(H → aa → bbµµ) and range between 2 × 10 −4 and 10 −3 , depending on m a .Modelindependent limits are set on the visible cross-section for new physics times branching ratio to the bbµµ final state (σ vis (X) × B(X → bbµµ)), and range from 0.1 fb to 0.73 fb, depending on the dimuon mass.
where m KL bb is the dijet invariant mass computed from the b-jet four-momenta corresponding to Êb 1 and Êb 2 , W is the transfer function of the b-jets, and F BW is a Breit-Wigner function centred on m µµ with

Figure 1 :
Figure 1: The (a) m µµ , (b) m bb before the KL fit, (c) m bbµµ before and (d) m KLbbµµ after the KL fit for events after the preselection stage, but removing the upper bound on m µµ .The t t contribution is modelled with the simulated sample normalised to the theoretical cross-section.The Drell-Yan contribution is taken from data templates (described in the text) and normalised to the total yield predicted by the Drell-Yan simulation.The signal distributions for all five simulated m a are also shown assuming the SM Higgs boson cross-section (including ggF, VBF and V H production) and B(H → aa → bbµµ) = 10%.The branching ratio in this and all subsequent figures is chosen so as to give good visibility on the plot.

37 Figure 2 :
Figure 2: Illustration of the signal, control and validation regions used in the analysis.The E miss T requirement in the VR2 region is the same as in the SR.The corresponding DY-template regions are defined in the same way except that the two-b-tag requirement is changed to a zero-b-tag requirement.

Figure 3 :
Figure 3: The predicted and observed m KL bbµµ distributions (a) after the preselection and the KL-fit ln(L max ) > −8 constraint and (b) across DYCR, SR and VR1 (shown separated by vertical dashed lines).Both are shown after the background-only fit and differ only in the E miss T < 60 GeV criterion being applied in (b).The signal distribution for m a = 40 GeV is also shown assuming the SM Higgs boson cross-section (including ggF, VBF and V H production) and (a) B(H → aa → bbµµ) = 0.5% and (b) B(H → aa → bbµµ) = 0.15%.The hashed bands show the total statistical and systematic uncertainties of the backgrounds.

Figure 4 : 1 Figure 5 :
Figure 4: The predicted and observed KL-fit ln(L max ) distribution across VR2 and SR (shown separated by a vertical dashed line) after the background-only fit.The signal distribution for m a = 40 GeV is also shown assuming the SM Higgs boson cross-section (including ggF, VBF and V H production) and B(H → aa → bbµµ) = 0.1%.The hashed bands show the total statistical and systematic uncertainties of the backgrounds.

Table 1 :
Summary of the dominant post-fit systematic uncertainties on the background and signal yields.The uncertainties are expressed as a percentage of the total background (middle column) and signal (rightmost column) yields per m µµ bin of the signal region.Shown are the uncertainties that exceed 2% in at least one m µµ bin.

Table 2 :
Total and individual background yields in five representative m µµ bins of the signal region.The yields are the post-fit values as determined by the background-only fit.The uncertainties shown include all systematic uncertainties and the statistical MC uncertainty.W+jets contribution in the SR is found to be negligible and is therefore not shown in the table.