Branching fraction measurements of the rare $B^0_s\rightarrow\phi\mu^+\mu^-$ and $B^0_s\rightarrow f_2^\prime(1525)\mu^+\mu^-$ decays

The branching fraction of the rare $B^0_s\rightarrow\phi\mu^+\mu^-$ decay is measured using data collected by the LHCb experiment at center-of-mass energies of $7$, $8$ and $13\,\rm{TeV}$, corresponding to integrated luminosities of $1$, $2$ and $6\,{\rm fb}^{-1}$, respectively. The branching fraction is reported in intervals of $q^2$, the square of the dimuon invariant mass. In the $q^2$ region between $1.1$ and $6.0\,{\rm Ge\kern -0.1em V}^2\!/c^4$, the measurement is found to lie $3.6$ standard deviations below a Standard Model prediction based on a combination of Light Cone Sum Rule and Lattice QCD calculations. In addition, the first observation of the rare $B^0_s\rightarrow f_2^\prime(1525)\mu^+\mu^-$ decay is reported with a statistical significance of nine standard deviations and its branching fraction is determined.

Recent studies of rare semileptonic b → s + − decays exhibit tensions between experimental results and Standard Model (SM) predictions of branching fractions [1][2][3][4][5], angular distributions [6][7][8][9][10][11], and lepton universality [11][12][13][14][15][16][17][18][19].Since these decays are only allowed via higher-order electroweak (loop) diagrams in the SM, they constitute powerful probes for non-SM contributions.One of the most significant discrepancies appears in the branching fraction of the B 0 s → φµ + µ − decay [1,2].Using 3 fb −1 of data collected with the LHCb experiment at center-of-mass energies of 7 and 8 TeV, the branching fraction was measured below the SM prediction at the level of three standard deviations (σ) [1].This Letter presents an updated measurement using data taken at center-of-mass energies of 7, 8 and 13 TeV during the 2011, 2012 and 2015-2018 data-taking periods, with integrated luminosities corresponding to 1, 2 and 6 fb −1 , respectively.Compared to the 3 fb −1 sample alone, this represents an increase of about a factor of four in the number of produced B 0 s mesons.The branching fraction is determined in intervals of q 2 , the squared invariant mass of the dimuon system.In addition, the observation of the B 0 s → f 2 (1525)µ + µ − decay and a determination of its branching fraction are reported.This constitutes the first observation of a rare semileptonic decay involving a spin-2 meson in the final state, and provides complementary information to transitions involving pseudoscalar or vector mesons.In the following, the shorthand notation f 2 is used to refer to the f 2 (1525) meson.The inclusion of charge-conjugate processes is implied throughout.
The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, detailed in Refs.[20,21].The online event selection is performed by a trigger [22] that consists of hardware and software stages.The former selects signal candidates containing a muon with significant transverse momentum with respect to the beam axis.At the software stage, a full event reconstruction is applied.Simulated events are used in this analysis to determine the reconstruction and selection efficiency of signal candidates, and to estimate contamination from residual background.The simulated samples are produced using the software described in Refs.[23][24][25].Residual mismodeling in simulation is corrected for using control samples from data.The B 0 s → φµ + µ − and B 0 s → f 2 µ + µ − decays are reconstructed in the K + K − µ + µ − final state.Particle identification criteria are applied to the kaon and muon candidates.The muons (kaons) are further required to have χ 2 IP > 9(6) with respect to any primary pp interaction vertex (PV) in the event, where χ 2 IP denotes the difference in the vertex-fit χ 2 of the PV when reconstructed with or without the considered track.The four finalstate tracks are fit to a common vertex that is required to have good quality and to be significantly displaced from any PV in the event.Signal candidates are retained if the lies between 5270 and 5700 MeV/c 2 .The invariant mass of the dikaon system, m(K + K − ), is required to be within 12 MeV/c 2 of the known φ mass for the B 0 s → φµ + µ − decay, or within 225 MeV/c 2 of the known mass of the wider f 2 resonance for the B 0 s → f 2 µ + µ − decay [26].The q 2 regions between 8.0 and 11.0 GeV 2 /c 4 , and between 12.5 and 15.0 GeV 2 /c 4 , are dominated by tree-level B 0 s decays into final states with a J/ψ or ψ(2S) meson.While these regions are vetoed in the selection of the signal modes, the decays to charmonium are used as high-yield control modes.The B 0 s → J/ψφ decay is used for normalization.The q 2 region from 0.98 to 1.1 GeV 2 /c 4 is also vetoed to remove contributions from B 0 s → φ(→ µ + µ − )φ decays.
To reduce combinatorial background, formed from random track combinations, a boosted decision tree (BDT) algorithm [27,28] is applied.The BDT classifier is trained on data using cross-validation techniques [29], with B 0 s → J/ψφ events as signal proxy and candidates from the upper mass sideband m(K + K − µ + µ − ) > 5567 MeV/c 2 as background proxy.The classifier combines the B 0 s transverse momentum and χ 2 IP , the angle between the B 0 s momentum and the vector connecting the PV and the decay vertex of the B 0 s candidate, the fit quality of the B 0 s vertex and its displacement from the associated PV, particle identification information and χ 2 IP of the final-state particles.The criterion on the BDT output is optimized by maximizing the expected significance of the B 0 s → φµ + µ − and B 0 s → f 2 µ + µ − signals separately, due to different levels of background contamination.The requirement on the BDT classifier yields a signal efficiency of 96% (85%) and a background rejection of 96% (95%) for the B 0 s → φµ + µ − (B 0 s → f 2 µ + µ − ) decay mode.Finally, information from particle identification is combined with invariant mass variables, constructed under the relevant particle-hypotheses, to reject background from Λ 0 b → pK − µ + µ − decays, where the proton is misidentified as a kaon, and from B 0 s → J/ψφ, B 0 s → ψ(2S)φ and B 0 → J/ψK * 0 decays, where a final state hadron is misreconstructed as a muon and vice-versa.
As the relative efficiencies vary according to the data-taking conditions, the data are split into the 2011-2012, 2015-2016 and 2017-2018 periods.The yields of the normalization mode for the different data-taking periods are determined using extended unbinned maximum-likelihood fits to the m(K + K − µ + µ − ) distribution.The B 0 s → J/ψφ decay is modeled using the sum of two Gaussian functions with a common mean and a power-law tail towards upper and lower mass.The combinatorial background is modeled using an exponential function.The m(K + K − µ + µ − ) distribution of the normalization mode for the full data sample, overlaid with the fit projections, is shown in Fig. 1 (left).The yields of the normalization mode, N J/ψφ , are determined to be 62 980 ± 270, 70 970 ± 290, and 148 490 ± 410 for the three different data-taking periods, where the uncertainties are statistical only.
For the rare B 0 s → φµ + µ − decay, a simultaneous extended maximum-likelihood fit of the data samples for the different periods is performed in intervals of q 2 , where the signal yields are parameterized using Eq. ( 1) and the differential branching fraction is shared between the samples.The model used to describe the m(K + K − µ + µ − ) distribution is the same as for the B 0 s → J/ψφ normalization mode.The model parameters for the signal component are fixed to those from the fit of the normalization mode, where the q 2 dependence of the mass resolution is accounted for with scaling factors determined from simulation.
Negligible contributions from physical background, including B 0 s → K + K − µ + µ − decays with the K + K − system in an S-wave configuration, are not considered in the fit and a systematic uncertainty is assigned.Integrated over the full q 2 range, signal yields, s → φµ + µ − signal candidates, integrated over q 2 and overlaid with the fit projections.
N φµ + µ − , of 458 ± 12, 484 ± 13, and 1064 ± 28 are found from the simultaneous fit to the different data sets.Figure 1 (right) shows the m(K + K − µ + µ − ) distribution of the full data sample, integrated over q 2 and overlaid with the fit projections.Figures for the different data-taking periods are available as Supplemental Material.
The relative branching fraction measurement is affected by systematic uncertainties on the fit model and the efficiency ratio, where the latter is determined using SM simulation.A summary of the systematic uncertainties is provided in the Supplemental Material.The dominant systematic uncertainty on the absolute branching fraction (Eq. 1) originates from the model used to simulate B 0 s → φµ + µ − events (0.04-0.10 × 10 −8 GeV −2 c 4 ).The model depends on ∆Γ s , the decay width difference in the B 0 s system [31], and the specific form factors used.The effect of the model-choice on the relative efficiency is assessed by varying ∆Γ s by 20%, corresponding to the difference in ∆Γ s between the default value [32] and that of Ref. [26], and by comparing the form factors in Ref. [33] with the older calculations in Ref. [34].The observed differences are taken as a systematic uncertainty.Other leading sources of systematic uncertainty arise from the limited size of the simulation sample (0.02-0.07 × 10 −8 GeV −2 c 4 ) and the omission of small background contributions from the fit model (0.01-0.04 × 10 −8 GeV −2 c 4 ).
The resulting relative and total branching fractions are given in Table 1.In addition, the differential branching fraction is shown in Fig. 2, overlaid with SM predictions.These predictions are based on form factor calculations using Light Cone Sum Rules (LCSRs) [33,35] at low q 2 and Lattice QCD (LQCD) [36,37] at high q 2 , which are implemented in the flavio software package [38].In the q 2 region between 1.1 and 6.0 GeV 2 /c 4 , the measured branching fraction of (2.88 ± 0.22) × 10 −8 GeV −2 c 4 , lies 3.6 σ below a precise SM prediction of (5.37 ± 0.66) × 10 −8 GeV −2 c 4 which uses both LCSR and LQCD calculations.A less precise SM prediction of (4.77 ± 1.01) × 10 −8 GeV −2 c 4 based on LCSRs alone lies 1.8 σ above the measurement.To determine the total branching fraction, the branching fractions of the individual q 2 intervals are summed and corrected for the vetoed q 2 regions using q 2 veto = (65.47± 0.27) %.This efficiency is determined using SM simulation, and its uncertainty originates from the comparison of form factors s → φµ + µ − )/dq 2 branching fraction, both relative to the normalization mode and absolute, in intervals of q 2 .The uncertainties are, in order, statistical, systematic, and due to the uncertainty on the branching fraction of the normalization mode.
The B 0 s → f 2 µ + µ − decay is searched for using the combined q 2 region [0.1, 0.98] ∪ [1.1, 8.0] ∪ [11.0, 12.5] GeV 2 /c 4 .The branching fraction of the signal decay is determined relative to the B 0 s → J/ψφ normalization mode, according to where the ratio of branching fractions B(φ [26] is used.To separate the f 2 signal from S-and P-wave contributions to the wide m(K + K − ) mass window, a two-dimensional fit to the m(K ) using the sum of two Gaussian functions with a power-law tail towards upper and lower mass, and in m(K + K − ) using a relativistic spin-2 Breit-Wigner function.The model parameters are determined from data using fits to the B 0 s → J/ψf 2 control mode and are fixed for the signal mode.Contributions from the S-wave and P-wave resonances, e.g. the φ and the φ(1680) mesons, are combined and described with a linear function in m(K + K − ) and use the same model as the signal in m(K + K − µ + µ − ).Interference effects are neglected as these were found to be small in the study of B 0 s → J/ψK + K − decays in Ref. [39].The combinatorial background is modeled using an exponential function in both the reconstructed B 0 s mass and the mass of the dikaon system.Background from B 0 → K + π − µ + µ − and Λ 0 b → pK − µ + µ − decays is found to be non-negligible in the wide m(K + K − ) window.These background components are included in the fit model, with their yields constrained to the expected values and line shapes determined on simulated events.
The branching fraction of the B 0 s → f 2 µ + µ − decay is determined using a simultaneous fit to the three data samples.The branching fraction of the signal and the S-and P-wave contributions are shared between the data samples.From this fit, the signal yields, N f 2 µ + µ − , are found to be 62 ± 8, 67 ± 8, and 161 ± 20 for the different data-taking periods.Figure 3 shows the m(K + K − µ + µ − ) and m(K + K − ) mass distributions, where the latter is shown within 50 MeV/c 2 of the known B 0 s mass [26], overlaid with the fit projections.The significance of the signal is determined using Wilks' theorem [40], comparing the log-likelihood with and without the signal component.The B 0 s → f 2 µ + µ − decay is observed with a statistical significance of 9 σ.Systematic effects on the significance due to the choice of fit model are negligible.
The dominant systematic uncertainties on the relative branching fraction of the B 0 s → f 2 µ + µ − decay originate from the uncertainty of the branching fraction ratio ), the modeling of the parameters of the Breit-Wigner function describing the f 2 resonance, and the simplified fit model for the m(K + K − ) distribution (0.03 × 10 −7 ).The effect of the simplified fit model is evaluated using pseudoexperiments, in which events are generated using the amplitude model in Ref. [39] and fit with the default model.The observed difference in the determined yield is taken as a systematic uncertainty.Further details on the systematic uncertainties associated with B(B 0 s → f 2 µ + µ − ) are given in the Supplemental Material.The fraction of signal events within the considered q 2 region is calculated using the q 2 -differential distribution in Ref. [41] and found to be q 2 veto = (73.8± 2.8) %.Accounting for this factor, the relative and total branching fractions are determined to be where the given uncertainties are, in order, statistical, systematic, from the extrapolation to the full q 2 range and, for the absolute branching fraction, from the uncertainty on the branching fraction of the normalization mode.The total B 0 s → f 2 µ + µ − branching fraction is found to be in agreement with SM predictions [41][42][43].
In summary, the most precise measurement of the branching fraction of the rare B 0 s → φµ + µ − decay is presented, using LHCb data corresponding to an integrated luminosity of 9 fb −1 .Consistent with earlier measurements [1,2], the data are found to lie below SM expectations.In the q 2 region between 1.1 and 6.0 GeV 2 /c 4 the measurement deviates by 3.6 σ with respect to a precise SM prediction [33,[35][36][37][38].These results supersede, and are consistent with, those of Refs.[1,2].In addition, the first observation of the rare B 0 s → f 2 µ + µ − decay is reported with a statistical significance of nine standard deviations and the resulting branching fraction is found to be in agreement with SM predictions [41][42][43].
IFIN-HH (Romania), CBPF (Brazil), PL-GRID (Poland) and NERSC (USA).We are indebted to the communities behind the multiple open-source software packages on which we depend.Individual groups or members have received support from ARC and ARDC (Australia); AvH Foundation (Germany); EPLANET, Marie Sk lodowska-Curie Actions and ERC (European Union); A*MIDEX, ANR, IPhU and Labex P2IO, and Région Auvergne-Rhône-Alpes (France); Key Research Program of Frontier Sciences of CAS, CAS PIFI, CAS CCEPP, Fundamental Research Funds for the Central Universities, and Sci.& Tech.Program of Guangzhou (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal and GENCAT (Spain); the Leverhulme Trust, the Royal Society and UKRI (United Kingdom).

Supplemental material Supplemental figures
Figure 4 shows the K + K − µ + µ − invariant mass versus q 2 for selected (top) B 0 s → φµ + µ − and (bottom) B 0 s → f 2 µ + µ − candidates.The signal modes are clearly visible as a vertical band around the known mass of the B 0 s meson.
across all data-taking periods. i The K + K − µ + µ − invariant mass of the selected (left) B 0 s → J/ψφ and (right) B 0 s → φµ + µ − candidates, integrated over q 2 , is shown in Fig. 5 for the different data-taking periods.The total fit projection (black line) is overlaid on the data, along with the signal component (blue line) and background component describing combinatorial background (red dotted line).

Systematic uncertainties
The systematic uncertainties associated with the measurement of the branching fractions of the B 0 s → φµ + µ − and B 0 s → f 2 µ + µ − decays are summarized in Table 2.The Physics model in Table 2 refers to the model used for the generation of B 0 s → φµ + µ − and B 0 s → f 2 µ + µ − decays in simulation.Studied variations to the B 0 s → φµ + µ − physics model are detailed in the Letter.The model used for B 0 s → f 2 µ + µ − decays accounts only for phase-space effects.The q 2 distribution is therefore not an exact description of data and is corrected to the predictions in Ref. [41].The difference in the relative efficiency with and without this correction is assigned as a systematic uncertainty.
The Residual background in Table 2 refers to contamination of the signal modes from residual physical background.Background contributions to B 0 s → φµ + µ − decays are neglected in the fit and a systematic uncertainty is assigned.For B 0 s → f 2 µ + µ − decays, a systematic uncertainty is associated with the choice of lineshape used to describe background from B 0 → K + π − µ + µ − and Λ 0 b → pK − µ + µ − decays in the fit.The systematic uncertainty associated with the B 0 s → φµ + µ − signal fit model in Table 2 is obtained using an alternative description for the radiative tail of the B 0 s meson.For the B 0 s → f 2 µ + µ − decay, the lineshape in m(K + K − ) of both the non-signal B 0 s → K + K − µ + µ − contributions and the signal decay are varied.For the B 0 s → K + K − µ + µ − contributions, an alternative description is taken from Ref. [39], as detailed in the Letter.For the f 2 lineshape, the input values for the Blatt-Weisskopf barrier functions [44] are varied, namely the barrier radius of the f 2 and B 0 s mesons, along with the orbital angular momentum of the B 0 s meson.The simulation corrections in Table 2 refer to the uncertainties associated with applied corrections to simulated events.Corrections are recalculated using alternative binning schemes or accounting for the finite statistics of the control modes used to derive the corrections.The uncertainty associated with small levels of mismodelling in distributions which are not directly corrected for in the default approach (e.g.tracking efficiencies) are indicated under residual mismodelling.

Figure 1 :
Figure1: Reconstructed invariant mass of the K + K − µ + µ − system for (left) the B 0 s → J/ψφ normalization mode and (right) the B 0 s → φµ + µ − signal candidates, integrated over q 2 and overlaid with the fit projections.

Figure 3 :
Figure3: Reconstructed invariant mass of (left) the K + K − µ + µ − system and (right) the K + K − system for B 0 s → f 2 µ + µ − candidates, overlaid with the fit projections.The m(K + K − ) distribution is shown in the B 0 s signal region ±50 MeV/c 2 around the known B 0 s mass.

Figure 6
Figure6shows the (left) K + K − µ + µ − and (right) K + K − invariant mass distributions of selected B 0 s → f 2 µ + µ − candidates for the different data-taking periods.The total fit projection (black line) is overlaid on the data along with projections of individual fit components describing: the signal (blue line), other B 0 s → K + K − µ + µ − decays (green dash-dotted line), combinatorial background (red dotted line) and Λ 0 b (magenta long dashed line) and B 0 (cyan medium size dashed line) decays.