Search for scalar resonances decaying into $\mu^{+}\mu^{-}$ in events with and without $b$-tagged jets produced in proton-proton collisions at $\sqrt{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$^{-1}$ of $\sqrt{s}= 13$ TeV proton-proton collision data recorded by the ATLAS experiment at the LHC. No significant 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$\%$ 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 1.6 and 44 fb respectively, which is consistent with expectations.


Introduction
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.  Figure 1: Feynman diagrams for the b-quark associated production mode 1(a) and the gluon-gluon fusion production mode 1(b). 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.0 TeV 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-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. For the generation of tt 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 tt 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 tt system recoils. The tt 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] [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 tt 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 tt) 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 tt, 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.0 TeV 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.0 TeV, 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 statisticq µ 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.0 TeV 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 tt. 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.
Background Asimov data Observed data Z+ HF 1.00 ± 0.23 1.47 ± 0.26 Z+ LF 1.00 ± 0.02 1.02 ± 0.02 tt 1.00 ± 0.04 1.02 ± 0.04 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.

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  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.

Theoretical uncertainties 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 tt 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 tt 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 tt 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 tt 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 tt 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  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.0 TeV. 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.  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(a) and 5(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.  Figure 5: The 5(a) observed and 5(b) expected 95% CL upper limit on the production cross section times branching ratio for Φ → µµ as a function of the fractional contribution from b-quark associated production and the scalar boson mass.

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.