Search for a high-mass dimuon resonance produced in association with b quark jets at $\sqrt{s}$ = 13 TeV

A search for high-mass dimuon resonance production in association with one or more b quark jets is presented. The study uses proton-proton collision data collected with the CMS detector at the LHC corresponding to an integrated luminosity of 138 fb$^{-1}$ at a center-of-mass energy of 13 TeV. Model-independent limits are derived on the number of signal events with exactly one or more than one b quark jet. Results are also interpreted in a lepton-flavor-universal model with Z$'$ boson couplings to a bb quark pair ($g_\mathrm{b}$), an sb quark pair ($g_\mathrm{b}\delta_\mathrm{bs}$), and any same-flavor charged lepton ($g_\ell$) or neutrino pair ($g_\nu$), with $\left|g_{\nu}\right| = \left|g_\ell\right|$. For a Z$'$ boson with a mass $m_{\mathrm{Z}'}$ = 350 GeV (2 TeV) and $\left|\delta_\mathrm{bs}\right|$ $\lt$ 0.25, the majority of the parameter space with 0.0057 $\lt$ $\left|g_\ell\right|$ $\lt$ 0.35 (0.25 $\lt$ $\left|g_\ell\right|$ $\lt$ 0.43) and 0.0079 $\lt$ $\left|g_\mathrm{b}\right|$ $\lt$ 0.46 (0.34 $\lt$ $\left|g_\mathrm{b}\right|$ $\lt$ 0.57) is excluded at 95% confidence level. Finally, constraints are set on a Z$'$ model with parameters consistent with low-energy b $\to$ s$\ell\ell$ measurements. In this scenario, most of the allowed parameter space is excluded for a Z$'$ boson with 350 $\lt m_{\mathrm{Z}'}$ $\lt$ 500 GeV, while the constraints are less stringent for higher $m_{\mathrm{Z}'}$ hypotheses. This is the first dedicated search at the LHC for a high-mass dimuon resonance produced in association with multiple b quark jets, and the constraints obtained on models with this signature are the most stringent to date.


Introduction
Over the last few years, a number of experimental results [1][2][3][4][5][6][7][8][9][10][11][12][13] have suggested that physics beyond the standard model (SM) could manifest itself in b → sℓ − ℓ + transitions. Experimental results for this process could be explained by the existence of a new neutral vector boson (Z ′ ) coupling to lepton (ℓ) pairs [14,15] with a mass near the TeV scale. In such a scenario, the new Z ′ boson would couple to b and s quarks, as represented by the Feynman diagrams in Fig. 1. While the most recent experimental results [16,17] on the ratio of the branching fractions for b → sµ − µ + to b → se − e + are in agreement with the SM, there remains some tension in other b → sℓ − ℓ + observables, such as the overall event rates [2, 3, 13] and angular distributions [1, 4, 5, 7-9], although theoretical predictions for these are less well established than those for the ratio of branching fractions. The current situation was recently discussed in Ref. [18].
bs bb b s g g g g Figure 1: Feynman diagrams of Z ′ → µ − µ + with a Z ′ boson produced via bb → Z ′ or sb → Z ′ , with at least one b quark in the final state. While a Z ′ bb coupling may be present in any generic model, a Z ′ sb coupling could arise through flavor mixing between the second-and third-generation quarks.
Inclusive searches for beyond-the-SM (BSM) Z ′ bosons have already been performed at the CERN LHC [19,20]. However, they are limited by a large Drell-Yan (DY) background and might not be sensitive to scenarios in which the Z ′ boson couples preferentially to secondor third-generation quarks. Therefore, the goal of this analysis is to search for a Z ′ → µ − µ + resonance with an explicit requirement on the presence of b quark jets, which strongly disfavors DY events.
The relevant interactions can be described through a lepton-flavor-universal (LFU) model with a Lagrangian simplified from the one in Ref. [21], ℓ=e,µ,τ ℓγ η P L ℓ + g ν ∑ ν=ν e ,ν µ ,ν τ ν γ η P L ν where P L denotes the left-handed projection operator. There are four coupling parameters in this model: a common g ℓ coupling for all charged leptons, a common g ν coupling for all neutrinos, the g b coupling that scales both Z ′ bb and Z ′ sb interactions, and the separate δ bs parameter that solely scales the Z ′ sb interaction. In this analysis, we probe values of δ bs between 0 and 0.25: The Z ′ sb interaction is fully suppressed at δ bs = 0, and δ bs values above 0.25 may yield scenarios incompatible with the measurement of the mass difference between the mass eigenstates of the neutral B s mesons [21].
The B 3 − L 2 model is a less generic Z ′ model introduced in 2017 [22, 23] to accommodate the experimental b → sℓ − ℓ + anomalies. It is based on a new B 3 −L 2 U(1) gauge symmetry in which B 3 and L 2 represent the third-generation baryon and muon lepton numbers, respectively. In this model, g Z ′ denotes the coupling of the Z ′ boson to SM fermions, and the angle θ 23 controls the mixing angle between the second-and third-generation quarks. This mixing angle is equivalent to arctan(−δ bs ) in the notation of Eq.
(1). The allowed parameter space for the B 3 −L 2 model has been recently reevaluated in Ref. [18].
Studies of dilepton invariant mass distributions in final states with b quark jets have previously been performed at the LHC [24][25][26]. However, these studies either focused on final states with exactly one b quark jet, or used inclusive samples with at least one such jet, or searched for nonresonant signatures in the high dimuon invariant mass tail. Some of the previous studies were based on data sets with only about one fourth of the integrated luminosity used in this analysis. Furthermore, they all suffered from large backgrounds arising from the pair production of t quarks (tt), each decaying to one b quark and one W boson, with W → ℓν. In this analysis, we reduce this background substantially, as noted below, and we further optimize the event categorization for the presence of one or more b quark jets. We focus on dimuon resonances because the experimental tensions that motivate this search are thought to be consistent with BSM contributions in the muon sector [18].
The analysis can be summarized as follows: • The search is performed for a narrow Z ′ → µ − µ + resonance with a mass m Z ′ > 350 GeV in the presence of at least one b quark jet. The width of the resonance is assumed to be narrow relative to the dimuon invariant mass (m µµ ) resolution. The validity of this assumption is ensured by restricting the search to the regions of parameter space where the Z ′ width does not exceed one half of the m µµ resolution. The resolution itself ranges from 6.6 GeV at m Z ′ = 350 GeV to 30 (106) GeV at m Z ′ = 1 (2.5) TeV. • The data sample is collected with the CMS detector during the LHC Run 2 datataking period of 2016-2018 at a center-of-mass energy of 13 TeV, and corresponds to an integrated luminosity of 138 fb −1 (2016: 36.3 fb −1 ; 2017: 41.5 fb −1 ; 2018: 59.8 fb −1 ). • Events in this data sample are categorized according to the multiplicity of b quark jets: N b = 1 and N b ≥ 2. Tight and relaxed requirements are used together in identifying these jets, as discussed in Section 2, in order to maximize the sensitivity to a possible signal. Events with N b = 0 are not considered because the background is about 150 times larger in this category due to the contribution from DY events. • The dominant backgrounds across the probed m µµ range arise from the DY process and tt production. The DY background is already reduced by requiring the detection of at least one b quark jet. The tt background is suppressed by requiring that the minimum invariant mass of any muon-b quark jet pairing (min(m µb )) be > 175 GeV, i.e., above the t quark mass. Such a requirement reduces the tt background by more than two orders of magnitude while retaining most of the predicted signal with large m Z ′ , thus offering an improved signal discrimination compared to other kinematic requirements used in previous similar studies. Other background sources, including tZ+X, tW+X, ttV (V = W, Z, γ * ), and ttH processes, as well as diboson (WW, WZ, and ZZ) production, are less important and further reduced by vetoing events with any additional lepton or isolated high transverse momentum (p T ) charged hadron. The veto for events with isolated charged hadrons is used to enhance the rejection of hadronically decaying τ leptons and unidentified electrons or muons. • We extract constraints on the total event yield as a function of the fraction of events in each of the N b = 1 and N b ≥ 2 categories. The constraints are determined using unbinned maximum likelihood fits of the m µµ distributions with analytic functions within mass ranges sliding coherently with the probed value of m Z ′ , as further discussed in Section 4. The functional forms used in the fit and the general statistical procedure are also described in Section 4 and are similar to those in Ref. [27]. The constraints on the total event yield are then reinterpreted in terms of the parameters of specific models with a Z ′ boson coupled to b quarks.
All of the results for this analysis can also be found on HEPData [28].

Experimental setup
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker (up to a pseudorapidity coverage of |η| < 2.5), a lead tungstate crystal electromagnetic calorimeter (|η| < 3), and a brass and scintillator hadron calorimeter (HCAL, with |η| < 3), each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors up to |η| = 5. Muons are measured in gas-ionization detectors (|η| < 2.4) embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [29].
Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [30]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [31].
A particle-flow (PF) algorithm [32] aims to reconstruct and identify each individual particle in an event, using an optimized combination of all subdetector information. The particles reconstructed with this algorithm are hereafter referred to as PF candidates, and their isolation is measured from the flux of other charged PF candidates within a cone of ∆R = √ (∆ϕ) 2 + (∆η) 2 = 0.3 around their direction, where ϕ is the azimuthal angle.
The primary proton-proton (pp) interaction vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone (as described in Section 9.4.1 of Ref. [33]). Muons originating from the PV are identified using a set of requirements [34] that are optimized to maximize their selection efficiency at high p T , i.e., roughly above 200 GeV. The muon p T assignment is performed with the 'TuneP' algorithm [35], which mainly relies on the tracker information instead of data from the PF algorithm. This avoids the potential bias arising from electromagnetic showers in the muon detectors, caused by extra particles radiated by the high-p T muon. The difference between the p T value assigned using only tracker information and that obtained from the PF algorithm is taken into account when applying requirements on other PF-based quantities in the events. As the analysis targets highp T muons, which may radiate photons when traversing detector material, muons are required to be isolated with respect to nearby charged-particle tracks, i.e., the presence of neutral PF candidates is not considered in assessing the level of isolation. The tracker-only isolation of muons, which is defined as the scalar p T sum of all tracks, excluding that of the muon, within a longitudinal distance of 0.2 cm from the PV and a cone of ∆R = 0.3 relative to the muon track, is required to be less than 5 GeV and less than 5% of the muon p T as measured in the tracker.
Events are collected with triggers that require the presence of at least one muon at HLT with p T > 50 GeV and |η| < 2.4. The trigger selection efficiency is measured using the "tag-andprobe" method [36] and independent data sets. The efficiency varies from about 70% in the detector endcaps, i.e., at |η| ≳ 2.1, to about 95% in the detector barrel, yielding an overall efficiency of about 90% with only a slight dependence on the data-taking period [37].
Events are further required to contain two oppositely charged muons with p T > 53 GeV and |η| < 2.4 satisfying the identification and isolation criteria, with an overall muon reconstruction and selection efficiency larger than 90%. The two muons are then used to form a µ − µ + resonance candidate. Events with additional muons of p T > 10 GeV and |η| < 2.4 satisfying the same identification and isolation criteria are vetoed. We also veto events that contain isolated PF muon candidates originating from the PV with p T > 5 GeV. These vetoes, together with those described below, reduce the SM background arising from tZ+X, tW+X, ttV, ttH, and diboson production.
In the analysis, electrons are identified using loose requirements [38], for the purpose of vetoing events with electrons originating from the PV. The electron isolation is measured from the flux of PF photons and hadron candidates within a cone of ∆R = 0.3 around its direction, including corrections for contributions from additional pp interactions within the same or nearby bunch crossings (pileup). This PF-based isolation is required to be below a threshold that varies as a function of the electron p T . Events with any such electron of p T > 10 GeV and |η| < 2.5 are vetoed, together with events that contain any isolated PF electron candidate originating from the PV with p T > 5 GeV.
The largest fraction of all visible τ lepton decay products consists of a single charged hadron, often accompanied by multiple neutral pions. For this reason, the veto for events with extra leptons is expanded to include isolated PF charged-hadron candidates originating from the PV with p T > 10 GeV.
Jets are reconstructed using the anti-k T algorithm [39, 40] with a distance parameter of 0.4. Jet energies are corrected for instrumental effects and contributions from pileup [41]. We select jets with p T > 20 GeV and |η| < 2.5, that must also be separated by ∆R > 0.4 from each muon candidate. The identification of b quark jets uses the DEEPJET algorithm [42,43]. At least one b quark jet must satisfy tight identification criteria (∼58% tagging efficiency for ∼0.1% misidentification probability of jets from light quarks and gluons), and all other b quark jets are identified using relaxed requirements (∼76% tagging efficiency for ∼1% misidentification probability). In order to reject µb pairs that originate from t quark decays, we calculate all possible values of m µb among the selected muons and b quark jets, and require min(m µb ) > 175 GeV, i.e., above the value of the t quark mass. This requirement is used to suppress the tt background, which is reduced by more than two orders of magnitude, while retaining most of the possible signals with large m Z ′ , as shown in Fig. 2.
The missing transverse momentum vector (⃗ p miss T ) is estimated from the negative of the vector p T sum of all PF candidates, where the pileup-per-particle-identification algorithm [41, 44] weights each PF candidate for its probability to originate from the PV in order to reduce the pileup dependence of this observable. The ⃗ p miss T vector is corrected further for residual incon-sistencies that may be introduced from the non-PF high-p T muon identification criteria and the muon p T assignment through the 'TuneP' algorithm. While no source of significant genuine p miss T is expected except for neutrinos in the case of semileptonic decays of the heavy-flavor quarks, large p miss T may nevertheless be observed if the momentum of a muon or b quark jet is mismeasured. Thus, events are vetoed if p miss T > 250 GeV and ⃗ p miss T is aligned (|∆ϕ| < 0.3) or antialigned (|∆ϕ| > π − 0.3) with any selected muon or b quark jet. This requirement retains almost the totality of the possible signals. Other anomalous high-p miss T events, due to reconstruction issues, detector malfunctions, or noncollision backgrounds, are rejected by dedicated filters. These filters identify more than 85-90% of these spurious high-p miss T events with a mistagging rate less than 0.1% [45]. In addition, the three-dimensional angle between the two muons is required to be smaller than π − 0.02 to further suppress cosmic ray muon contributions.
In 2018, a sector of the HCAL (−3.2 < η < −1.3 and −1.57 < ϕ < −0.87) was not operational for a data-taking period corresponding to an integrated luminosity of 39 fb −1 , resulting in the misreconstruction of hadronic jets in that η-ϕ range. Therefore, events collected during this period are rejected if jets or electrons are found in the affected region, with upper and lower boundaries for jets enlarged by 0.2 in both η and ϕ in order to account for the size of the jet cone. The simulation of events, discussed in Section 3, is also adjusted by applying the same requirement in a subset of events corresponding to the fraction of data affected by this issue. The efficiency loss over the full data set is about 5% regardless of the value of m µµ .

Event simulation
We use simulated samples of various BSM Z ′ models to motivate the event selection and provide model-dependent limits. On the other hand, the SM background is estimated directly from data as a continuum background in the m µµ spectrum, parametrized by analytic functions. Thus, the background estimation does not rely on simulation. The SM simulated samples are used only to validate and optimize the event selection as well as to visually compare the observed dimuon mass distribution to the expected SM background from simulation.
The simulated samples for the dominant DY and tt backgrounds, and the tW+X process, are generated using the POWHEG 2 [46-49] program at next-to-leading order (NLO) in perturbative quantum chromodynamics (QCD). Samples for the tZ+X, ttW, ttZ, and tt γ * processes, and for the ttH process are generated at NLO in QCD using MADGRAPH5 aMC@NLO [50] v2.6.5 and v2.6.1, respectively, with the FxFx [51] scheme to match jets from matrix element calculations and parton showers. Samples for diboson production are generated using either MAD-GRAPH5 aMC@NLO v2.6.5 or POWHEG 2 at NLO in QCD.
As discussed in Section 1, we consider a signal Z ′ model involving Z ′ bb and Z ′ sb couplings and Z ′ → µµ decays. We generate such interactions using MADGRAPH5 aMC@NLO v2.9.9 at LO in QCD with m Z ′ ranging from 350 to 2500 GeV. For all the generated samples, the total width Γ(Z ′ ) of the Z ′ is always smaller than one half of the m µµ resolution to ensure the validity of the narrow width approximation intrinsic to the search strategy. Predictions for any other model with a likewise narrow width value can be obtained through reweighting by the double ratio of the quantity Γ(Z ′ → q i q j )Γ(Z ′ → µ − µ + )/Γ(Z ′ ) between the target model and the generated sample. Here, Γ(Z ′ → q i q j ) denotes the partial width for Z ′ → bb, Z ′ → sb and Z ′ → bs, and Γ(Z ′ → µ − µ + ) is the partial width for Z ′ → µ − µ + . The reweighting is performed by separating the simulated event sample into contributions from matrix elements with different q i q j → Z ′ interactions, as shown in Fig. 1, and the use of partial widths ensures generality (1), for several Z ′ mass hypotheses. For illustrative purposes, we choose couplings |g ℓ | = |g ν | = |g b | = 0.03 and δ bs = 0. The efficiency of the requirement on the Z ′ signal process is about 30%, 70%, and 85% for m Z ′ = 400, 1000, and 2000 GeV, respectively. The contribution of background processes other than DY and tt is so small that it is only barely visible at the bottom of the stacked histogram. The hatched region indicates the statistical uncertainty arising from the limited size of the SM simulated samples (Section 3). Histograms are normalized to unit area.
across different ways to construct target models.
In all samples, the parton shower and hadronization are modeled with PYTHIA 8.230 [52] using the CP5 [53] tune. The parton distribution functions are taken from NNPDF 3.1 [54]. For the SM simulated samples, they are taken at next-to-NLO in QCD. For the signal simulation, they are taken at LO in QCD to match the signal cross section calculations. All the simulated events are generated with a distribution of additional pp interactions per bunch crossing that is adjusted to match the corresponding pileup distribution measured in data. Finally, the detector response is simulated with the GEANT4 [55] package.

Analysis strategy
The events selected as described in Section 2 are categorized according to the multiplicity of b quark jets: N b = 1 and N b ≥ 2. The search is performed by simultaneously fitting the unbinned m µµ distributions in these two categories with analytic functions, within mass ranges that vary coherently with the value of m Z ′ to be probed.
We parametrize the signal m µµ distribution as the sum of a Gaussian distribution with a doublesided Crystal Ball function [56,57] with varying proportions, and a common resolution parameter σ mass . The choice of the functional form and the dependence of the corresponding pa-rameters on m Z ′ are determined from simulation. The fits are performed using a finite grid of simulated signal samples, and the fit parameters are subsequently interpolated across m Z ′ . As mentioned in Section 1, the values of σ mass range from 6.6 GeV at m Z ′ = 350 GeV to 30 (106) GeV at m Z ′ = 1 (2.5) TeV. Because of the narrow-width approximation intrinsic to the search strategy, the values of σ mass determined from simulation as a function of m Z ′ are independent of the signal coupling values. Fits to the data are performed within an m µµ mass window of ±10 σ mass around the probed m Z ′ value with the restriction m µµ > 275 GeV. The fit window is chosen to be relatively wide to ensure that the background is reliably estimated, without any significant bias that may otherwise arise in the presence of a potential signal at its center. The window slides coherently with the probed value of m Z ′ , which ranges between 350 and 2500 GeV, in steps reflecting the size of σ mass . This approach is used, instead of fitting the full m µµ distribution at once, so as not to rely on prior knowledge of the SM background shape across the full range of interest.
With a relatively wide mass window, the background m µµ distribution can be modeled analytically using Bernstein polynomials, exponential, or power-law functions. The choice of the functional form for the background model is treated as a discrete nuisance parameter in the fit. An envelope of functional forms is provided to the fit. Then, in the minimization process of the negative logarithm of the likelihood, the discrete profiling method [58] is used to select the best fit background model as a function of the signal strength. The statistical procedures to determine the goodness of fit for the backgrounds and extract the final fit results are similar to those described in Ref. [27]. For Bernstein polynomials, the best order in each fit mass window and event category is selected by means of a Fisher test [59]. Typically, first order Bernstein polynomials are selected, since higher order functions do not significantly improve the description of the data because of the small number of observed events. For the same reason, only Bernstein polynomials of the first order are used to model the backgrounds in windows with less than ten observed events. To avoid artifacts from fluctuations due to limited event samples, models that would yield an increasing background as a function of m µµ within the fit mass window are rejected. Only background models that fit the data appropriately are considered: a goodnessof-fit test is performed based on a χ 2 test statistic. This is converted into a p-value through a set of pseudo-experiments, and models with p < 0.01 are rejected.
The potential presence of a bias in the measurement of a signal due to the choice of the background functional forms is assessed by means of pseudo-experiments. First, a varying amount of signal is injected on top of the background generated according to the selected functional form. Then, a background+signal fit is performed allowing the signal yield to float freely. The bias is quantified as the difference between the measured and injected signal yields relative to the statistical uncertainty in the measured signal yield. No evidence of a statistically significant bias is found.
A systematic uncertainty of 1.6% in the expected signal yields arises from the integrated luminosity measurement [60][61][62]. In addition, uncertainties in signal event yields from jet energy scale, trigger, muon reconstruction and b tagging efficiencies, and finite simulated sample size are included in the fits. We include a 5% uncertainty in the signal acceptance to account for possible mismodeling of the muon identification efficiency in the simulation. We also assess an uncertainty ≲5% arising from the estimation of the size of the fit mass window. Finally, we account for possible effects on the signal shape from the uncertainties in the muon momentum scale (≲0.1% at m Z ′ = 1 TeV) and resolution (≲10% σ mass ). A summary of the sources of signal uncertainty and their sizes is given in Table 1. Uncertainties from other sources are found to be negligible. Table 1: Summary of signal uncertainties together with their sizes or size ranges. The uncertainties are grouped based on whether they affect the normalization or the shape of the signal, and any variations for the two categories of N b are shown. The uncertainties in the signal normalization are expressed as relative uncertainties (%) with respect to the nominal expected yields. The fit parameter m µµ defines the position of the maximum of the reconstructed m µµ distribution for a given m Z ′ hypothesis. The fit resolution parameter σ mass is distinguished from the values of σ mass extracted from simulation.

Results
The m µµ distributions in the two event categories are shown in Fig. 3 for data, SM background simulation, and several representative signal hypotheses. As already mentioned, the SM background simulation is only used for illustrative purposes since the background is estimated directly from data across the full m µµ range of interest.
The results of a single set of background-only fit to data in the event categories with N b = 1 and N b ≥ 2 are shown in Fig. 4 for a sample signal mass hypothesis, m Z ′ = 500 GeV, within the corresponding m µµ fit window. No significant excess is observed over the background-only expectation in any of the mass windows explored.
The results are used to set model-independent limits at 95% confidence level (CL) on the number of signal events with N b ≥ 1. The relative fraction of events in the N b ≥ 2 category, f 2b , is varied to probe a range of hypotheses of signal production in association with b quarks. These limits are shown in Fig. 5. They are extracted by performing unbinned fits to the m µµ distribution in data, using windows that slide in steps reflecting the size of σ mass . As different data events are included in each of the fit windows, fluctuations are present in both the expected and observed exclusion limits. The frequency of these fluctuations follows the size of the steps used to probe different m Z ′ hypotheses, and their amplitude is typically less than the 68% expected band. The statistical procedure to set the limits follows a modified frequentist approach employing the CL s criterion [63][64][65][66]. (1) with the assumed values |g ν | = |g ℓ |, and either |δ bs | = 0 or 0.25. As |g ℓ | or |g b | increases, the expected value of Γ(Z ′ ) may approach and exceed the value of σ mass . Since Γ(Z ′ ) is assumed to be smaller than σ mass in obtaining these results, we restrict our exclusion ranges to the regions where Γ(Z ′ ) < σ mass /2. For a Z ′ boson with a mass of 350 GeV (2 TeV) and |δ bs | < 0.25, the majority of the parameter space with 0.0057 < |g ℓ | < 0.35 (0.25 < |g ℓ | < 0.43) and 0.0079 < |g b | < 0.46 (0.34 < |g b | < 0.57) is excluded at 95% CL.   shown together with the corresponding selected background functional forms used as input to the discrete profiling method [58] when probing the m Z ′ = 500 GeV hypothesis. The expected signal distribution for the LFU model described in Eq. (1), with couplings |g ℓ | = |g ν | = |g b | = 0.03 and δ bs = 0, is overlaid. The displayed mass range corresponds to the fit window used for this m Z ′ hypothesis, which is ±10 σ mass around the probed m Z ′ value. While the likelihood fits are performed on unbinned data, here we present the data in binned histograms with binning chosen to reflect the size of σ mass .
The results are also used to set constraints on the B 3 −L 2 model from Ref. [18]. As in the case of the LFU model, we restrict our exclusion ranges to regions with Γ(Z ′ ) < σ mass /2. As shown in Fig. 7, most of the allowed parameter space is excluded for a Z ′ boson with 350 < m Z ′ < 500 GeV, while the constraints are less stringent for higher m Z ′ . Since exclusion limits are nearly independent of the θ 23 parameter in the B 3 −L 2 model or, equivalently, δ bs in the LFU model, we also show in Fig. 8 exclusion regions in the |g Z ′ |-m Z ′ plane for a fixed value of θ 23 = 0. For larger values of the couplings, the narrow width approximation intrinsic to the search strategy is not considered valid. For a given mass, the region enclosed between the solid black (dashed red) and dotted curves is (expected to be) excluded. The dotted curve for m Z ′ = 500 GeV lies beyond the displayed |g Z ′ | range and is, therefore, not shown. The shaded blue area represents the region preferred from the global fit in Ref.
[18] at 95% CL. The region above the green dash-dotted curve is incompatible at 95% CL with the measurement of the mass difference between the mass eigenstates of the neutral B s mesons [21]. For larger values of the couplings, the narrow width approximation intrinsic to the search strategy is not considered valid. The region enclosed between the solid black (dashed red) and the dotted curves is (expected to be) excluded. The shaded blue area represents the |g Z ′ | range preferred from the global fit in Ref.

Summary
A search for high-mass dimuon resonance production in association with one or more b quark jets has been presented, using data collected with the CMS experiment at the LHC that correspond to an integrated luminosity of 138 fb −1 at a center-of-mass energy of 13 TeV.
The limits are derived on the total number of signal events with N b = 1 and ≥2, where N b denotes the multiplicity of b quark jets, and these limits are model independent to the extent that the signal arises from the direct production and subsequent decay of a narrow dimuon resonance. The relative fraction of events with N b ≥ 2 is varied to probe a range of hypotheses of signal production in association with b quarks. The limits are presented as a function of the analyzed dimuon resonance mass values.
Results are also interpreted in terms of a lepton-flavor-universal model that involves Z ′ boson couplings to b quarks (g b ) and muons, where the Z ′ boson couplings to all neutrinos (g ν ) and to all charged leptons (g ℓ ) are assumed to be equal, the g b coupling scales both Z ′ bb and Z ′ sb interactions, and the separate δ bs coupling solely scales the Z ′ sb interaction. The exclusions in this model are presented in terms of the coupling strengths g ℓ and g b , and the mass of the Z ′ boson (m Z ′ ). For a Z ′ boson with m Z ′ = 350 GeV (2 TeV) and |δ bs | < 0.25, the majority of the parameter space with 0.0057 < |g ℓ | < 0.35 (0.25 < |g ℓ | < 0.43) and 0.0079 < |g b | < 0.46 (0.34 < |g b | < 0.57) is excluded at 95% confidence level.
Constraints are also set on a specific Z ′ model (B 3 −L 2 ), constructed to accommodate possible contributions to b → sℓ − ℓ + transitions beyond the standard model. In this scenario, most of the allowed parameter space is excluded for a Z ′ boson with 350 < m Z ′ < 500 GeV, while the constraints are less stringent for higher m Z ′ hypotheses.
This is the first dedicated search at the LHC for a high-mass, narrow dimuon resonance produced in association with multiple b quark jets, and the constraints obtained on models with this signature are the most stringent to date.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: