Search for scalar resonances decaying into (cid:22) + (cid:22) (cid:0) in events with and without b -tagged jets produced in proton-proton collisions at p s = 13 TeV with the ATLAS detector

: A search for a narrow scalar resonance decaying into an opposite-sign muon pair produced in events with and without b -tagged jets is presented in this paper. The search uses 36.1 fb (cid:0) 1 of p s = 13 TeV proton-proton collision data recorded by the ATLAS experiment at the LHC. No signi(cid:12)cant excess of events above the expected Standard Model background is observed in the investigated mass range of 0.2 to 1.0 TeV. The observed upper limits at 95% con(cid:12)dence level on the cross section times branching ratio for b -quark associated production and gluon-gluon fusion are between 1.9 and 41 fb and 1.6 and 44 fb respectively, which is consistent with expectations.


Introduction
This paper presents a search, in proton-proton (pp) collision data at a centre-of-mass energy of 13 TeV, targeting a scalar particle Φ with a mass in the range 0.2-1.0TeV decaying into two opposite-sign muons.The natural width of the particle is assumed to be much narrower than the experimental resolution.To maximize the sensitivity to the presence of b-quarks, the analysis is performed in two data categories.The first search category requires at least one jet tagged as containing b-hadrons (b-tagged) in the final state.The second category requires exactly zero b-tagged jets in the event.Simultaneous fits to these two categories are used to search for both the gluon-gluon fusion and b-associated production mechanisms of the scalar particle.This work is primarily motivated as a signature-based search; the event selection is not designed to target any specific model and instead facilitates the comparison of the data with predictions from various theoretical models.
The discovery of the 125 GeV Higgs boson at the Large Hadron Collider (LHC) [1,2] was a significant step in understanding the mechanism of electroweak symmetry breaking (EWSB).Measurements of the properties of this particle [3][4][5][6][7] are so far consistent with expectations for the Standard Model (SM) Higgs boson [8][9][10][11][12][13].However, the observed Higgs boson could be only a part of an extended scalar sector, as predicted by several theoretical models beyond the SM.For example, such an extended scalar sector is predicted by a class of extensions of the SM known as the two-Higgs-doublet models (2HDM) [14,15].The 2HDM posit the existence of three neutral bosons with properties and coupling strengths differing from those of the SM Higgs particle.One of the higher-mass neutral Higgs bosons may couple to muons at a higher rate than to τ-leptons, in contradiction to the SM Yukawa ordering.This is the case in the Flavourful Higgs model [16], for example.The search described in this paper sets limits on the cross-section times branching ratio of Φ decaying into muon pairs which are more stringent than the equivalent limits for Φ decaying into τ-lepton pairs [17].This is due to higher identification efficiency and lower background rate for muons than for τ-leptons, and the invariant-mass resolution being better for muon pairs than for τ-lepton pairs.Similarly, the coupling of these higher-mass neutral Higgs bosons to b-quarks could be enhanced relative to the SM Higgs boson coupling.Hence, the production of Φ in association with b-quarks (bbΦ), shown in Figure 1(a), could be more abundant than the production of Φ by gluon-gluon fusion (ggF), shown in Figure 1(b).The Φ → b b decay mode, where Φ is produced in associations with b-quarks, is investigated by CMS in Ref. [18].
Additional interest for the search for an excess of events with respect to the SM predictions in the dimuon plus b-tagged jets final state arises from Z and leptoquark models [19], which could result in either resonant or non-resonant contributions to the pp → tt µ + µ − → W bW bµ + µ − final state.This work complements the search for Z candidates decaying into muon pairs presented in Ref. [20]: by classifying events in terms of the presence or absence of a b-tagged jet, it may be sensitive to signatures not observable with an inclusive selection.
Previous similar searches focused explicitly on beyond-the-SM Higgs bosons decaying into muon pairs, in the mass range 90-500 GeV [21,22] using √ s = 7 and 8 TeV data.This work extends the exclusion limits for new scalar resonances with masses up to 1.0 TeV.
The analysis described in this paper makes use of a fit to the observed dimuon invariant mass (m µµ ) spectra.
The presence of a signal produced either in association with b-quarks or via gluon-gluon fusion is tested with simultaneous fits to the two data regions (with and without a b-tagged jet) separately for bbΦ and ggF production modes.No assumptions are made about the relative contributions of the bbΦ and ggF cross-sections.
This paper is structured as follows.Section 2 describes the ATLAS detector.Section 3 discusses the data and the simulated event samples used to model the signal and the background processes.The event reconstruction is discussed in Section 4, while Section 5 describes the background estimate and introduces the signal interpolation procedure used to model the m µµ distribution for resonance mass hypotheses for which no simulated sample was generated.The event yields and the description of the statistical analysis are discussed in Section 6, followed by Section 7, which provides an overview of the systematic uncertainties.Section 8 summarizes the results.

ATLAS detector
ATLAS [23], [24] is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and near 4π coverage in solid angle.1It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer.
The inner tracking detector (ID) covers the pseudorapidity range |η| < 2.5, and is surrounded by a superconducting solenoid providing a 2 T magnetic field.At small radii, a high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track.It is followed by the silicon microstrip tracker, which provides eight measurement points per track.The silicon detectors are complemented by a gas-filled straw-tube transition radiation tracker, which extends the tracking capability radially with typically 35 measurements per track for particles at |η| = 2.0.
Electromagnetic (EM) calorimetry, within the region |η| < 3.2, is provided by barrel and endcap highgranularity lead/liquid-argon (LAr) sampling calorimeters, with an additional thin LAr presampler covering |η| < 1.8 to correct for energy loss in the upstream material.For |η| < 2.5 the EM calorimeter is divided into three layers in depth.A steel/scintillator-tile calorimeter provides hadronic calorimetry for |η| < 1.7.LAr ionization, with copper as absorber, is used for the hadronic calorimeters in the endcap 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 z-axis along the beam pipe.The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upward.Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis.The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2).The angular distance is measured in units of ∆R = (∆η) 2 + (∆φ) 2 region 1.5 < |η| < 3.2.The solid-angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules in 3.1 < |η| < 4.9, optimized for EM and hadronic measurements, respectively.
The muon spectrometer surrounds the calorimeters and comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field provided by one barrel and two end-cap toroid magnets.The precision chamber system covers the region |η| < 2.7 with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region.The muon trigger system covers the range |η| < 2.4 with resistive-plate chambers in the barrel and thin-gap chambers in the endcap regions.
A two-level trigger and data acquisition system is used to record events for offline analysis [25].The level-1 trigger is implemented in hardware and uses a subset of the detector information to reduce the event rate to a value of at most 100 kHz.It is followed by a software-based high-level trigger which filters events using the full detector information and outputs events for permanent storage at an average rate of 1 kHz.

Data and simulated event samples
This search is performed with a sample of pp collision data recorded at a centre-of-mass energy √ s =13 TeV corresponding to an integrated luminosity of 36.1 fb −1 collected during 2015 and 2016.Events are considered for further analysis only if they were collected under stable LHC beam conditions and the relevant detector components were fully operational.
The data used in this analysis were collected using a combination of muon triggers.For the 2015 dataset an isolated muon with transverse momentum p µ T greater than 20 GeV was required, while for 2016 this requirement was raised to 26 GeV.To avoid inefficiencies due to the isolation requirement, these triggers were complemented by a trigger requiring a reconstructed muon with p µ T greater than 50 GeV, without any isolation requirements.Simulated signal events were generated similarly to the SM Higgs boson H, where H was replaced by a higher-mass scalar boson Φ with a constant width of 4 MeV.Nine MC samples were generated in the mass range 0.2-1.0TeV at intervals of 100 GeV.Both the ggF and bbΦ production modes were considered.Events were generated at next-to-leading order (NLO) with P -B 2 [26][27][28][29] and M G 5_ MC@NLO [30,31], for ggF and bbΦ respectively.The CT10 [32] set of parton distribution functions (PDFs) was used in the generation of ggF events while CT10 _ 4 [33] PDFs were used to produce the b-quark associated signal samples in the four-flavour scheme.In the gluon-gluon fusion production mode, P 8.210 [34] with the AZNLO [35] set of tuned parameters was used together with the CTEQ6L1 [36] PDF set for the parton shower, underlying event and hadronization simulation.For the b-quark associated production, P 8.210 with the A14 [37] set of tuned parameters was used together with the NNPDF2.3LO[38] PDF set for the parton shower, underlying event, and hadronization simulation.The EvtGen v1.2.0 program [39] was used to model the properties of the bottom and charm hadron decays in all signal samples.Since the generated signal samples are spaced equally in m Φ while the experimental resolution increases with m Φ , signal distributions for intermediate mass hypotheses were obtained using an interpolation procedure, described in Section 5.
Events containing Z/γ * + jets were simulated using P -B 2 [28,40] interfaced to the P 8.186 parton shower model and the CT10 PDF set.The AZNLO tune was used, with PDF set CTEQ6L1 for the modelling of non-perturbative effects.The EvtGen v1.2.0 program was used to model the properties of the bottom and charm hadron decays.Photos++ 3.52 [41] was used for photons radiated from electroweak vertices and charged leptons.Generator-level filters are employed to enhance the fraction of simulated events in the phase-space region that is most relevant for the analysis.Event yields were corrected with a mass-dependent rescaling at next-to-next-leading order (NNLO) in the QCD coupling constant, computed with VRAP 0.9 [42] and the CT14NNLO PDF set [33].Mass-dependent EW corrections were computed at NLO with mcsanc-1.20[43], along with photon-induced contributions (γγ → via t-and u-channel processes) computed with the MRST2004QED PDF set [44].
Additional Z/γ * + jets events, used to evaluate the systematic uncertainty from the modelling of the fitted m µµ distribution, were simulated with S 2.2.1 [45] and with M G 5_ MC@NLO interfaced with P 8.186.In S , Z/γ * + jets matrix elements were evaluated with up to two extra partons at NLO, and three or four extra partons were included at LO in QCD.The merging of different parton multiplicities was achieved through a matching scheme based on the CKKW-L [46] [47] merging technique using a merging scale of Q cut = 20 GeV.The models used for the parton shower and underlying event are the ones provided internally by S .The S 2.2.1 generator adopts a full five-flavour scheme, with massless b-quarks and c-quarks in the matrix elements.These quarks are given mass in the final state and massive quarks can be produced directly via the parton shower.The M G 5_ MC@NLO v2 generator merges the LO (QCD) matrix-element calculations for Z/γ * + jets with up to four additional partons (higher jet multiplicities are modelled by the parton shower algorithm).The merging scheme applied to combine different parton multiplicities is the CKKW-L scheme with a merging scale of Q cut = 30 GeV.For the LO matrix-element calculation the NNPDF2.3LO PDFs were used (with α S = 0.13).The A14 tune, together with the NNPDF2.3LO PDFs, was used for the parton shower.Similarly to S 2.2.1, M G 5_ MC@NLO adopts a full five-flavour scheme.
For the generation of t t and single top-quark production in the Wt-and t-channel, the P -B 2 generator with the CT10 PDF set was used for the matrix element calculations.The parton shower, fragmentation, and the underlying event were simulated using P 6.428 with the CTEQ6L1 PDF set and the corresponding Perugia 2012 tune [48].The top-quark mass was set to 172.5 GeV.For the generation of t t events, the h damp parameter of P -B , which controls the transverse momentum of the first additional emission beyond the Born configuration, was set to the mass of the top quark.The main effect of this parameter is to regulate the high transverse momentum emission against which the t t system recoils.The t t production sample was normalized to the predicted production cross-section as calculated with the T ++2.0 program to NNLO in perturbative QCD, including soft-gluon resummation to next-to-next-to-leading-log (NNLL) order [49].The normalization of the single top-quark event samples used an approximate calculation at NLO in QCD for the t-channels [50,51] while NLO+NNLL predictions were used for the Wt-channel [51].The EvtGen v1.2.0 program was used to describe the properties of the bottom and charm hadron decays.
Diboson processes were modelled using the S 2.1.1 [45] generator and they were calculated for up to one (Z Z) or zero (WW, W Z) additional partons at NLO and up to three additional partons at leading order (LO) using the Comix [52] and OpenLoops [53] matrix-element generators.The matrix element calculation was merged with the S internal parton shower [54] using the ME+PS@NLO prescription [55].The CT10 PDF set was used in conjunction with dedicated parton shower tuning developed by the S authors.The generator cross-sections, calculated at NLO, were used in this case.
The generated signal samples were simulated with the fast simulation [56,57] framework of ATLAS, which replaces the full simulation of the electromagnetic and hadronic calorimeters by a parameterized model.The background samples were simulated using the full Geant4-based simulation of the detector [58].
Finally, the simulated events were processed through the same reconstruction software as the data.The effects of pile-up from multiple proton-proton interactions in the same and neighbouring bunch crossings were accounted for by overlaying minimum-bias events simulated using P 8 with the A2 tune [59] and interfaced with the MSTW2008LO PDFs.Simulated events were then reweighted to match the pile-up conditions observed in the data.

Event reconstruction
Interaction vertices are reconstructed [60] from tracks measured by the inner detector.The primary vertex is defined as the reconstructed vertex with at least two associated tracks, each with transverse momentum p T > 400 MeV, that has the largest p 2 T for the associated tracks.Muon candidates are reconstructed from an ID track combined with a track or track segment detected in the muon spectrometer [61].These two tracks are used as inputs to a combined fit (for p T less than 300 GeV) or to a statistical combination (for p T greater than 300 GeV) [61].The combined fit takes into account the energy loss in the calorimeter and multiple-scattering effects.The statistical combination for high transverse momenta is performed to mitigate the effects of relative ID and MS misalignments.High-p T muon candidates are required to have at least three hits, one in each of the three layers of precision chambers in the MS.For these muons, specific regions of the MS where the alignment is suboptimal are vetoed as a precaution, as well as the transition region between the MS barrel and endcap (1.01 < |η| < 1.10).Further quality requirements are applied to the consistency of the ID and MS momentum measurements, and to the measured momentum uncertainty for the MS track.
The reconstructed muon candidates are required to have transverse momentum, p µ T , greater than 30 GeV and pseudorapidity |η µ | < 2.5.They are further required to be consistent with the hypothesis that they originate from the primary vertex by applying selections on the transverse (d 0 ) and longitudinal (z 0 ) impact parameters, defined relative to the primary vertex position: |d 0 /σ d 0 | < 3 and |z 0 sin θ| < 0.5 mm where σ d 0 denotes the uncertainty in the transverse impact parameter.A loose isolation requirement is applied, based on the sum of the momenta of inner detector tracks which lie within a variable-sized cone, with ∆R max = 0.3, around the muon track.This isolation requirement is tuned to yield a 99% efficiency over the full range of muon p T .
Jets are reconstructed from noise-suppressed energy clusters in the calorimeter [62] with the anti-k t algorithm [63, 64] with radius parameter R = 0.4.The energies of the jets are calibrated using a jet energy scale (JES) correction derived from both simulation and in situ studies using data [65].All jets are required to have p T > 25 GeV and |η| < 4.5.A multivariate selection that reduces contamination from pile-up [66] is applied to jets with p T < 60 GeVand with |η| < 2.4 utilizing calorimeter and tracking information to separate hard-scatter jets from pile-up jets.
Selected jets in the central region can be tagged as containing b-hadrons (b-tagged) by using a multivariate discriminant (MV2c10) [67, 68] that combines information from an impact-parameter-based algorithm, from the explicit reconstruction of a secondary vertex and from a multi-vertex fitter that attempts to reconstruct the full b-to c-hadron decay chain.At the chosen working point, the algorithm provides nominal light-flavour (u,d,s-quark and gluon) and c-jet misidentification rates of 3% and 32%, respectively, for an average 85% b-jet tagging efficiency, as estimated from simulated t t events for jets with p T > 20 GeV and |η| < 2.5.The flavour-tagging efficiencies from simulation are corrected separately for b-, c-and light-flavour jets, based on the respective data-based calibration analyses [69].Simulated events were corrected to reflect the muon and jet momentum scales and resolutions, b-jet identification algorithm calibration, as well as the muon trigger, identification, and isolation efficiencies measured in data.
An overlap removal procedure is applied to avoid a single particle being reconstructed as two different objects as follows.Jets not tagged as b-jets, but which are reconstructed within ∆R = 0.2 of a muon, are removed if they have fewer than three associated tracks or if the muon energy constitutes most of the jet energy.Muons reconstructed within a cone of size ∆R = min 0.4, 0.04 + 10 GeV/p µ T around the jet axis of any surviving jet are removed.Jets are also discarded if they are within a cone of size ∆R = 0.2 around an electron candidate.Electrons are reconstructed from clusters of energy deposits in the electromagnetic calorimeter that match ID tracks, and are identified as described in Ref. [70].For electrons with transverse energy E T > 10 GeV and pseudorapidity of |η| < 2.47, a likelihood-based selection [70] at the "loose" operating point is used.
Following the overlap removal, jet cleaning criteria are applied to identify jets arising from non-collision sources or noise in the calorimeters, and any event containing such a jet is removed [71].
The missing transverse momentum, E miss T , is defined as the negative vector sum of the transverse momenta of muons, electrons and jets associated with the primary vertex.A soft term [72] is added to include well-reconstructed tracks matched to the primary vertex that are not associated with any of the objects.
At least two reconstructed muons are required in the event, and the two highest-p T muons of the event are used to form a dimuon candidate.If these muons do not have opposite charges, the event is rejected.

Signal and background estimate
Events satisfying the preselection criteria of Section 4 are considered for further analysis.Depending on the value of the reconstructed dimuon invariant mass and the number of b-tagged jets, they are classified as being in either a control region, used to measure the rate of the dominant backgrounds (Z/γ * + jets and t t) or a signal region, used to search for the Φ → µµ signal.All signal and control regions are designed to be orthogonal.
Selected events are retained in the signal region if the dimuon invariant mass, m µµ , exceeds 160 GeV, while they are retained in the control region if 100 GeV < m µµ < 160 GeV.The control regions do not include the Z-boson peak; this avoids constraints on systematic uncertainties in a region of phase space which may behave differently with respect to the high mass tails where the signal search is focused.
Events in the signal region are further classified as "SRbTag", if at least one b-tagged jet is identified in the event, or "SRbVeto", otherwise.No further optimization of the signal region definition is considered, in order to minimize the assumptions made about the signal kinematics.
Events in the control region are further classified as "CRbVeto" if there is no b-tagged jet identified in the event, while remaining control region events are classified as "CRbTag", if the missing transverse momentum is less than 100 GeV, or "CRttbar" otherwise.
The dominant background sources for this signature are muon pairs through Z/γ * + jets, top-quark pair, and diboson production.Contributions from W+jets and QCD multi-jet events were estimated and found to be negligible and are therefore not included.All relevant background contributions were modelled using simulated samples as described in Section 3. Simulated Z/γ * + jets events were categorized depending on the generator-level "truth" labels of the jets in the event.Simulated jets were truth-labelled according to which hadrons with p T > 5 GeV were found within a cone of size ∆R = 0.3 around the jet axis.If a b-hadron was found, the jet was truth-labelled as a b-jet.If no b-hadron was found, but a c-hadron was present, then the jet was truth-labelled as a c-jet.Otherwise the jet was truth-labelled as a light-flavour (i.e., u,d,s-quark, or gluon) jet.The Z/γ * + jets events were classified as Z+ heavy flavour (Z+ HF) if a b-jet or c-jet was found at generator level and Z+ light flavour (Z+ LF) otherwise.These two categories of simulated events were treated as different background components.
The predictions for the expected numbers of t t, Z+ HF and Z+ LF events are adjusted to match the number of events in data in the control regions via a global fit (see Section 6), while the other backgrounds are normalized to their expected cross-section, discussed in Section 3.
The expected invariant mass distribution of the narrow resonance signal m Φ is modelled with a double-sided Crystal Ball (DSCB) [73] function in the mass range 0.2 TeV < m Φ < 1.0 TeV for all nine simulated MC samples.The width of the m µµ distribution is dominated by the experimental resolution, which can be described by a Gaussian distribution.The power-law asymmetric terms of the DSCB have enough degrees of freedom to model the m µµ tails for signals in the 0.2-1.0TeV mass range.In order to interpolate the signal parameterization to any mass value in this fit range, second-order polynomial parameterizations of all six signal-shape parameters as a function of m Φ are obtained from a simultaneous fit to all the generated mass points m Φ , separately for events with and without at least one b-tagged jet.Since all signal samples are scaled to the same arbitrary cross-section, any differences in the number of events are due to differences in the acceptances of the selection.The fitted normalization is parameterized with a second-order polynomial in a similar way to the other DSCB parameters.As a result of this procedure, the m µµ distribution for an arbitrary hypothesis mass can be found by constructing a DSCB from parameters calculated using the fitted coefficients and the hypothesized mass m Φ .The interpolation is performed separately for each systematic variation, and for the gluon-gluon fusion and b-quark associated production signal samples.Binned templates are then generated from the DSCB.To validate the interpolation procedure, one simulated sample at a time is removed from the simultaneous fit, and the template generated from the DCSB is compared with the m µµ distribution from the corresponding simulated sample and they agree well with each other.Nine MC samples were generated in the mass range 0.2-1.0TeV, every 100 GeV.The interpolation procedure was used to generate signal templates every 10 GeV between 0.2 and 0.3 TeV, every 20 GeV between 0.3 and 0.6 TeV, and every 50 GeV between 0.6 and 1.0 TeV.
The acceptance of b-quark associated production in the SRbTag region varies in the range 11-19% depending on m Φ .The acceptance in the SRbVeto region varies in the range 20-15%.For this reason, both SRbTag and SRbVeto are included in the global fit (see Section 6).The acceptance of gluon-gluon fusion in the SRbVeto region varies in the range 31-35%, but it is less than 2% for the SRbTag.

Statistical analysis
The test statistic qµ as defined in Ref.
[74] is used to determine the probability that the background-only model is compatible with the observed data, to extract the local p-value, and, if no hint of a signal is found in this procedure, to derive exclusion intervals using the CL s method.The binned likelihood function is built as the product of Poisson probability terms associated with the bins in the m µµ distribution in the 0.16-1.5 TeV range.It depends on the parameter of interest, on the normalization factors of the dominant backgrounds and on additional nuisance parameters (representing the estimates of the systematic uncertainties) that are each constrained by a Gaussian prior.The bin-by-bin statistical uncertainties of the simulated backgrounds are also considered.Separate fits are performed to test the different hypotheses: only bbΦ signal, only ggF signal, and signal composed of different mixtures of the two production modes.Data are fitted in the 0.16-1.5 TeV range to check that the background description is correct and to constrain systematic uncertainties, while the signal presence is tested in a narrower 0.2-1.0TeV range.
Logarithmic binning in m µµ is chosen to scale with experimental dimuon mass resolution (defined as full-width at half-maximum), which varies from 5% at m µµ = 200 GeV to 14% at m µµ = 1 TeV, while ensuring that the number of simulated background events in each bin is sufficient to reliably predict the background.The numbers of events in the CRbVeto, CRbTag, and CRttbar control regions are also included in the maximum-likelihood fit.Data in these regions are used to constrain the normalizations of the dominant background processes: Z+ LF, Z+ HF, and t t.The normalization factors used to scale the expected number of events for each process to match the data are free parameters in the fit.Their expected uncertainty is estimated using an Asimov dataset [74], a pseudo-data distribution equal to the signal-plus-background expectation for a given value of the signal cross-section times branching ratio.The expected and measured normalization factors are shown in Table 1.They are measured in a fit to data under the background plus signal hypothesis, where the signal corresponds to bbΦ production with mass m Φ = 480 GeV.
Table 1: Normalization factors of the dominant backgrounds as measured in a fit to data under the background plus signal hypothesis (480 GeV bbΦ signal).The uncertainty includes both the statistical and systematic sources.The reduction of the observed relative uncertainty of the Z+ HF normalization in the data fit is due to the large increase in the number of Z+ HF events when the data in signal regions are included.As discussed in Section 7, the experimental uncertainties are treated as fully correlated among different processes and regions, while the theoretical uncertainties are mostly uncorrelated among processes.
The numbers of observed events in the signal regions together with the predicted event yields from signal and background processes are shown in Table 2.The numbers are the results of the maximum-likelihood fit to data under the background plus signal hypothesis, where the signal corresponds to bbΦ production with m Φ = 480 GeV.
The observed numbers of events are compatible with those expected from SM processes, within uncertainties.The m µµ distribution in the SRbTag region is shown in Figure 2(a) for a bbΦ-only fit.The m µµ distribution in the SRbVeto region is shown in Figure 2(b), for a ggF-only fit.

Systematic uncertainties
The sources of systematic uncertainty can be divided into three groups: those of experimental origin, those related to the modelling of the simulated backgrounds, and those associated with signal simulation.The finite size of the simulated background samples is also an important source of uncertainty.
Table 2: Post-fit numbers of events from the combined fit under the background plus signal hypothesis (480 GeV bbΦ signal).Here "Total Bkg" represents the sum of all backgrounds.The quoted uncertainties are the combination of statistical and systematic uncertainties.The uncertainty in the total background determined by the fit is smaller than the sum in quadrature of the individual components due to the normalization factors and systematic uncertainties that introduce correlation between them.The number of signal events corresponds to the expected upper limit (UL) times the cross-section.The uncertainties in the expected signal yields are from fits using the Asimov dataset.

Experimental uncertainties
The dominant experimental uncertainties originate from residual mismodelling of the muon reconstruction and selection after the simulation-to-data efficiency correction factors have been applied.They include the uncertainty obtained from Z → µµ data studies and a high-p T extrapolation uncertainty corresponding to the decrease in the muon reconstruction and selection efficiency with increasing p T which is predicted by the Monte Carlo (MC) simulation.The degradation of the muon reconstruction efficiency was found to be approximately −3% per TeV as a function of muon p T .Uncertainties in the isolation and trigger efficiencies of muons [61], along with the uncertainty in their energy scale and resolution, are estimated from Z → µµ 13 TeV data.These are found to have only a small impact on the result.Other sources of experimental uncertainties are the residual differences in the modelling of the flavour-tagging algorithm, the jet energy scale and the jet energy resolution after the simulation-to-data correction factors are applied.The uncertainties in the energy scale and resolution of the jets and leptons are propagated to the calculation of E miss T , which also has additional uncertainties from the scale, resolution, and efficiency of the tracks used to define the soft term, along with the modelling of the underlying event.
The pile-up modelling uncertainty is assessed by varying the number of pile-up interactions in simulated events.The variations are designed to cover the uncertainty in the ratio of the predicted and measured cross-section of non-diffractive inelastic events producing a hadronic system of mass m X > 13 GeV [75].
The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%.It is derived, following a methodology similar to that detailed in Ref. [76], and using the LUCID-2 detector for the baseline luminosity measurements [77], from calibration of the luminosity scale using x-y beam-separation scans.

Background processes
As discussed earlier, the rates of major backgrounds are measured using data control regions, thus theoretical uncertainties in the predicted cross-section for Z/γ * + jets and t t processes are not considered.Instead, the effects of theoretical uncertainties in the modelling of the m µµ distribution and in the ratio of the event yields in signal regions to those in the control regions are evaluated.The dominant theoretical uncertainties are the shape uncertainty in the dimuon mass spectra of the t t and Z+ HF contributions.
For both the Z+ LF and Z+ HF processes, the following sources of theoretical and modelling uncertainties are considered: PDF uncertainties (estimated via eigenvector variations and by comparing different PDF sets), limited accuracy of the fixed-order calculation (estimated by QCD scale variations), variations in the choice of strong coupling constant value (α S (M Z )), EW corrections, and photon-induced corrections.These variations are treated as fully correlated between Z+ LF and Z+ HF processes.
The PDF variation uncertainty is obtained using the 90% confidence level (CL) CT14NNLO PDF error set and by following the procedure described in Refs.[78,79].Rather than a single nuisance parameter to describe the 28 eigenvectors of this PDF error set, which could lead to an underestimation of its effect, a re-diagonalized set of 7 PDF eigenvectors was used [33].This represents the minimal set of PDF eigenvectors that maintains the necessary correlations, and the sum in quadrature of these eigenvectors matches the original CT14NNLO error envelope well.They are treated as separate nuisance parameters.The uncertainties due to the variation of PDF scales and α S are derived using VRAP.The former is obtained by varying the renormalization and factorization scales of the nominal CT14NNLO PDF up and down simultaneously by a factor of two.The value of α S used (0.118) is varied by ±0.002.The EW correction uncertainty was assessed by comparing the nominal additive (1 + δ EW + δ QCD ) treatment with the multiplicative approximation ((1 + δ EW )(1 + δ QCD )) treatment of the EW correction in the combination of the higher-order EW and QCD effects.The uncertainty in the photon-induced correction is calculated from the uncertainties in the quark masses and the photon PDF.Following the recommendations of the PDF4LHC forum [79], an additional uncertainty due to the choice of nominal PDF set is derived by comparing the central values of CT14NNLO with those from other PDF sets, namely MMHT14 [80] and NNPDF3.0[81].The maximum width of the envelope of these comparisons is used as the PDF choice uncertainty, but only if it is larger than the width of the CT14NNLO PDF eigenvector variation envelope.
An additional modelling uncertainty is considered for the Z+ HF process.The m µµ spectrum simulated by P was compared with those simulated with S 2.2.1 [45] and with M G 5_ MC@NLO interfaced with P 8.186.A functional form was chosen to describe the envelope of the differences in the m µµ distribution shape between events with at least one b-quark simulated by S 2.2.1 and M G 5_ MC@NLO v2.
These systematic uncertainties not only imply uncertainties in the shape of the m µµ distribution in the signal region, but they also affect the ratio of the expected number of events in the signal region to the expected number in the control region.The effect on the Z+ LF normalization in the bVeto signal region is 2%, while the size of the uncertainty for Z+ HF in the bTag signal region is 7%.
Moreover, the ratio of the expected number of Z+ LF events in the bVeto control region to the events in bTag control and signal regions was estimated with the nominal P +P Z+jet sample, and with the alternative M G 5_ MC@NLO v2 sample.To cover the observed difference between the two, an additional uncertainty of 27% in the normalization of the Z+ LF process in the bTag signal and control region was applied.
For the t t process, theoretical uncertainties in the modelling of the m µµ spectrum and in the extrapolation from control region to signal region are also considered.These are estimated by comparing the nominal prediction with alternative ones.To estimate the impact of initial and final state radiation modelling, two alternative P +P samples were generated with the following parameters.In the first one, the renormalization (µ ) and factorization (µ ) scale were varied by a factor of 0.5, the value of h damp was doubled (2m t ) and the corresponding Perugia 2012 radiation tune variation was used.In the second sample, the renormalization and factorization scales were varied by a factor of 2, while the h damp parameter was not changed.The associated variation from the P2012 tune was used.These choices of parameters have been shown to encompass the cases where µ and µ are varied independently, and covered the measured uncertainties of the data for unfolded t t distributions [82,83].Differences in parton shower and hadronization models were investigated using a sample where P -B was interfaced with H ++ 2.7.1 [84] with the UE-EE-5 tune [37] and the corresponding CTEQ6L1 PDFs.The nominal sample was also compared with a sample generated with M G 5_ MC@NLO 2.2.1 [30] interfaced with H ++. A NLO matrix element and CT10 PDF were used for the t t hard-scattering process.The parton shower, hadronization and the underlying events were modelled using the H ++ 2.7.1 generator.The UE-EE-5 tune and the corresponding CTEQ6L1 PDF were used.A functional form was chosen to encompass the differences in the m µµ distribution shape between the nominal sample and all these alternative samples.These samples were also used to estimate the 3.5% extrapolation uncertainty from CRttbar to the SRbTag.Theoretical uncertainties that arise from higher-order contributions to the cross-sections and PDFs and affect the values of the predicted cross-sections for the diboson and single-top backgrounds are considered.

Signal processes
Uncertainties related to signal modelling include the uncertainties associated with the initial-and final-state radiation, the modelling of underlying events, the choice of the renormalization and factorization scales, and the parton distribution functions.None of them has a sizeable effect on the shape of the m µµ spectra, but some of them do affect the acceptance in the signal regions.Uncertainties for the different mass hypotheses were evaluated, but for simplicity, only the largest of the values is used.In the calculation of the uncertainties, appropriate requirements on the muon and jets kinematics are applied at hadron level for the SRbTag and SRbVeto regions.
The factorization and renormalization scales were varied by a factor of two up and down, including correlated and anti-correlated variations, both in P -B and in MG5_ MC@NLO.For the bbΦ process, the largest deviation from the value of the nominal acceptance is considered as the final scale uncertainty (2% in SRbTag, and 1% in SRbVeto).A 25% uncertainty for the acceptance of ggF events in the SRbTag is considered, as well as its anti-correlated effect in the SRbVeto region, following the procedure adopted in Ref. [85].
The PDF uncertainties are estimated by taking the envelope of the changes in the acceptance associated with the PDF4LHC15_nlo_100 (PDF4LC15_nlo_nf4_30) [79] error set for the ggF (bbΦ) signal process.They correspond to changes in acceptances not larger than 1%.Systematic variations of the parameters of the A14 (for bbΦ) and AZNLO (for ggF) tunes are used to account for the uncertainties associated with the initial and final state radiation and the modelling of underlying events.For bbΦ, these uncertainties correspond to changes of 3.8% and 3.2% of the acceptance in the SRbTag and SRbVeto regions.They are associated with a migration of events from one signal region to the other, hence they are treated as anti-correlated between the two signal regions.For ggF, the uncertainty in the SRbVeto is negligible, while it is 3.8% in the SRbTag.
The current result is dominated by the statistical uncertainty of the data (which accounts for 66% of the total uncertainty on the fitted value of the signal cross-section), and the finite size of the simulated background samples (which accounts for 25% of the total uncertainty).Amongst the systematic uncertainties, the modelling of the shape of the m µµ distribution for the Z+jets and top-antitop backgrounds dominates (3% of the total uncertainty) followed by the muon identification efficiency (1.5%).All other considered sources of experimental and theoretical uncertainties have a combined impact on the upper limit of less than 5%.

Results
The observed p-values as a function of m Φ , obtained applying the statistical analysis described in Section 6, are shown in Figure 3(a) for the bbΦ-only fit and in Figure 3(b) for the ggF-only fit.The lines in each figure correspond to separate fits.The dotted line represents the p-value for a fit including only SRbTag and all control regions (CRs), the dashed line represents the p-value for a fit including only SRbVeto+CRs, while the solid line corresponds to the combined SRbTag+SRbVeto+CRs fit.The dotted lines in the two figures show similar behaviour and so do the dashed lines, while the combined p-values are calculated for separate signal production modes implying very different contributions from the SRbVeto and SRbTag regions, which lead to substantially different limits.The largest excess of events above the expected background is observed for the b-quark associated production at about m Φ = 480 GeV and amounts to a local significance of 2.3σ, which becomes 0.6σ after considering the look-elsewhere effect over the mass range 0.2-1.0TeV.The look-elsewhere effect is estimated using the method described in Ref. [86].Since the data are in agreement with the predicted backgrounds the results are given in terms of exclusion limits.These are set using the modified frequentist CL s method [87].Upper limits on the cross section times branching ratio for a massive scalar particle are set at the 95% confidence level (CL) as a function of the particle mass.
The upper limits at 95% CL on σ Φ × B(Φ → µµ) (where σ Φ is the cross section and B is the branching ratio) assume the natural width of the resonance to be negligible compared to the experimental resolution, and they cover the mass range 0.2-1.0TeV.In Figure 4(a) and 4(b) the 95% CL upper limits are shown for b-quark associated production and gluon-gluon fusion, respectively.The expected limits are in the ranges 1.3-25 fb and 1.8-25 fb respectively.The observed limits are in the ranges 1.9-41 fb and 1.6-44 fb respectively.Both the expected and observed limits are calculated in the asymptotic approximation [74], which was verified with MC pseudo-experiments for the m Φ = 1000 GeV hypothesis.The upper limit on the ggF production mode is dominated by the SRbVeto, while for the bbΦ production mode both the SRbTag and SRbVeto regions contribute to the result because of the relative acceptance of the bbΦ production mode in the two regions.Limits shown in Figure 4(a) and 4(b) are obtained in simultaneous fits of the SRbVeto and SRbTag data regions and correspond to the different signal production modes.Usage of the same data events in both the bbΦ-only and the ggF-only fits explains correlations in the behaviour of the observed limits.

Conclusion
The ATLAS detector at the LHC has been used to search for a massive scalar resonance decaying into two opposite-sign muons, produced via b-quark associated production or via gluon-gluon fusion, assuming that the natural width of the resonance is negligible compared to the experimental resolution.The search is conducted with 36.1 fb −1 of pp collision data at √ s = 13 TeV, recorded during 2015 and 2016.
The observed dimuon invariant mass spectrum is consistent with the Standard Model prediction within uncertainties for events both with and without a b-tagged jet over the 0.2-1.0TeV range.The observed upper limits at 95% confidence level on the cross-section times branching ratio for b-quark associated production and gluon-gluon fusion are between 1.9 and 41 fb and

Figure 1 :
Figure 1: Feynman diagrams for the b-quark associated production mode 1(a) and the gluon-gluon fusion production mode 1(b).

Figure 2 :
Figure 2: Distributions of the dimuon invariant mass, m µµ , after the combined fit to data are shown normalized: 2(a) in the SRbTag, under the background plus signal hypothesis (480 GeV bbΦ signal) and 2(b) in the SRbVeto, under the background plus signal hypothesis (480 GeV ggF signal).The fit for a m Φ =480 GeV signal corresponds to the largest excess observed above the background expectation.Each background process is normalized according to its post-fit cross-section.The templates for the m Φ =200 GeV, m Φ =480 GeV and m Φ =1 TeV mass hypotheses are normalized to the expected upper limit.The data are shown by the points, while the size of the statistical uncertainty is shown by the error bars.The blue arrows represent the data points outside of the frame.The hatched band shows the total systematic uncertainty of the post-fit yield.

Figure 3 :
Figure 3: Observed p-values as a function of m Φ for 3(a) a bbΦ-only fit and 3(b) a ggF-only fit.The dotted line corresponds to the data in the SRbTag only, the dashed line to the data in SRbVeto only, and the solid curve to the combination of the two data categories.All three fits include data from all control regions.

Figure 4 :
Figure 4: The observed and expected 95% CL upper limits on the production cross section times branching ratio for a massive scalar resonance produced via 4(a) b-quark associated production and 4(b) gluon-gluon fusion.

Figures 5 (Figure 5 :
Figures 5(a) and5(b)  show the observed and expected 95% CL upper limits on the production cross section times branching ratio for Φ → µµ as a function of the fractional contribution from b-quark associated production (σ bbΦ /[σ bbΦ + σ ggF ]) and the scalar resonance mass.These fractions are derived assuming that other production mechanisms do not affect the result.
1.6 and 44 fb respectively, which is consistent with expectations.Sweden; SERI, SNSF and Cantons of Bern and Geneva, Switzerland; MOST, Taiwan; TAEK, Turkey; STFC, United Kingdom; DOE and NSF, United States of America.In addition, individual groups and members have received support from BCKDF, CANARIE, CRC and Compute Canada, Canada; COST, ERC, ERDF, Horizon 2020, and Marie Skłodowska-Curie Actions, European Union; Investissements d' Avenir Labex and Idex, ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; CERCA Programme Generalitat de Catalunya, Spain; The Royal Society and Leverhulme Trust, United Kingdom.Topological cell clustering in the ATLAS calorimeters and its performance in LHC Run 1, Eur.Phys.J. C 77 (2017) 490, arXiv: 1603.02934[hep-ex].JHEP 09 (2016) 001, arXiv: 1606.03833[hep-ex].[74] G. Cowan, K. Cranmer, E. Gross and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics, Eur.Phys.J. C 71 (2011) 1554, arXiv: 1007.1727[physics.data-an].[75] ATLAS Collaboration, Measurement of the Inelastic Proton-Proton Cross Section at √ s = 13 TeV with the ATLAS Detector at the LHC, Phys.Rev. Lett.117 (2016) 182002, arXiv: 1606.02625[hep-ex].Luminosity determination in pp collisions at √ s = 8 TeV using the ATLAS detector at the LHC, Eur.Phys.J. C 76 (2016) 653, arXiv: 1608.03953[hep-ex]. 176Waseda University, Tokyo; Japan. 177Department of Particle Physics, Weizmann Institute of Science, Rehovot; Israel. 178Department of Physics, University of Wisconsin, Madison WI; United States of America. 179Fakultät für Mathematik und Naturwissenschaften, Fachgruppe Physik, Bergische Universität Wuppertal, Wuppertal; Germany. 180Department of Physics, Yale University, New Haven CT; United States of America. 181Yerevan Physics Institute, Yerevan; Armenia.Also at Borough of Manhattan Community College, City University of New York, NY; United States of America.b Also at California State University, East Bay; United States of America.c Also at Centre for High Performance Computing, CSIR Campus, Rosebank, Cape Town; South Africa.d Also at CERN, Geneva; Switzerland.e Also at CPPM, Aix-Marseille Université, CNRS/IN2P3, Marseille; France.f Also at Département de Physique Nucléaire et Corpusculaire, Université de Genève, Genève; Switzerland.g Also at Departament de Fisica de la Universitat Autonoma de Barcelona, Barcelona; Spain.h Also at Departamento de Física Teorica y del Cosmos, Universidad de Granada, Granada (Spain); Spain.i Also at Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa, Lisboa; Portugal.j Also at Department of Applied Physics and Astronomy, University of Sharjah, Sharjah; United Arab Emirates.Also at II.Physikalisches Institut, Georg-August-Universität Göttingen, Göttingen; Germany.z Also at Institucio Catalana de Recerca i Estudis Avancats, ICREA, Barcelona; Spain.aa Also at Institut für Experimentalphysik, Universität Hamburg, Hamburg; Germany.ab Also at Institute for Mathematics, Astrophysics and Particle Physics, Radboud University Nijmegen/Nikhef, Nijmegen; Netherlands.ac Also at Institute for Particle and Nuclear Physics, Wigner Research Centre for Physics, Budapest; Hungary.ad Also at Institute of Particle Physics (IPP); Canada.ae Also at Institute of Physics, Academia Sinica, Taipei; Taiwan.a f Also at Institute of Physics, Azerbaijan Academy of Sciences, Baku; Azerbaijan.
[69] ATLAS Collaboration, Measurements of b-jet tagging efficiency with the ATLAS detector using t t events at √ s = 13 TeV, JHEP 08 (2018) 089, arXiv: 1805.01845[hep-ex].[70]ATLASCollaboration,Electronefficiencymeasurements with the ATLAS detector using the 2015LHC proton-proton collisiondata, ATLAS-CONF-2016-024, 2016, : https://cds.cern.ch/record/2157687.[71]ATLASCollaboration, Selection of jets produced in 13 TeV proton-proton collisions with the ATLAS detector, ATLAS-CONF-2015-029, 2015, : https://cds.cern.ch/record/2037702.[72]ATLASCollaboration, Performance of missing transverse momentum reconstruction with the ATLAS detector using proton-proton collisions at √ s = 13 TeV, Eur.Phys.J. C 78 (2018) 903, arXiv: 1802.08168[hep-ex].[73]ATLASCollaboration, Search for resonances in diphoton events at √ s = 13 TeV with the ATLAS detector, a k Also at Department of Financial and Management Engineering, University of the Aegean, Chios; Greece.lAlso at Department of Physics and Astronomy, University of Louisville, Louisville, KY; United States of America.mAlso at Department of Physics and Astronomy, University of Sheffield, Sheffield; United Kingdom.nAlso at Department of Physics, California State University, Fresno CA; United States of America.oAlso at Department of Physics, California State University, Sacramento CA; United States of America.pAlso at Department of Physics, King's College London, London; United Kingdom.qAlso at Department of Physics, St. Petersburg State Polytechnical University, St. Petersburg; Russia.rAlso at Department of Physics, Stanford University; United States of America.sAlso at Department of Physics, University of Fribourg, Fribourg; Switzerland.tAlsoat Department of Physics, University of Michigan, Ann Arbor MI; United States of America.uAlso at Giresun University, Faculty of Engineering, Giresun; Turkey.vAlso at Graduate School of Science, Osaka University, Osaka; Japan.wAlso at Hellenic Open University, Patras; Greece.xAlso at Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest; Romania.y ag Also at Institute of Theoretical Physics, Ilia State University, Tbilisi; Georgia.ah Also at Instituto de Física Teórica de la Universidad Autónoma de Madrid; Spain.