Angular analysis and differential branching fraction of the decay $B^0_s\to\phi\mu^+\mu^-$

An angular analysis and a measurement of the differential branching fraction of the decay $B^0_s\to\phi\mu^+\mu^-$ are presented, using data corresponding to an integrated luminosity of $3.0\, {\rm fb^{-1}}$ of $pp$ collisions recorded by the LHCb experiment at $\sqrt{s} = 7$ and $8\, {\rm TeV}$. Measurements are reported as a function of $q^{2}$, the square of the dimuon invariant mass and results of the angular analysis are found to be consistent with the Standard Model. In the range $1<q^2<6\, {\rm GeV}^{2}/c^{4}$, where precise theoretical calculations are available, the differential branching fraction is found to be more than $3\,\sigma$ below the Standard Model predictions.


Introduction
The decay B 0 s → φµ + µ − is mediated by a b → s flavour changing neutral current (FCNC) transition. In the Standard Model (SM) it is forbidden at tree-level and proceeds via loop diagrams as shown in Fig. 1. In extensions of the SM, new heavy particles can appear in competing diagrams and affect both the branching fraction of the decay and the angular distributions of the final-state particles.
This decay channel was first observed and studied by the CDF collaboration [1,2] and subsequently studied by the LHCb collaboration using data collected during 2011, corresponding to an integrated luminosity of 1.0 fb −1 [3]. While the angular distributions were found to be in good agreement with SM expectations, the measured branching fraction differs from the recently updated SM prediction by 3.1 σ [4,5]. A similar trend is also seen for the branching fractions of other b → sµ + µ − processes, which tend to be lower than SM predictions [6][7][8].
This paper presents an updated analysis of the decay B 0 s → φ(→ K + K − )µ + µ − using data accumulated by LHCb in pp collisions, corresponding to an integrated luminosity of 1.0 fb −1 collected during 2011 at 7 TeV and 2.0 fb −1 collected during 2012 at 8 TeV centreof-mass energy. The differential branching fraction dB(B 0 s → φµ + µ − )/dq 2 is determined as a function of q 2 , the square of the dimuon invariant mass. In addition, a three-dimensional angular analysis in cos θ l , cos θ K and Φ is performed in bins of q 2 . Here, the angle θ K (θ l ) denotes the angle of the K − (µ − ) with respect to the direction of flight of the B 0 s meson in the K + K − (µ + µ − ) centre-of-mass frame, and Φ denotes the angle between the µ + µ − and the K + K − decay planes in the B 0 s meson centre-of-mass frame. Compared to the previously published fit of the one-dimensional projections of the decay angles [3], the full three-dimensional angular fit gives improved sensitivity and allows access to more angular observables.
The decay B 0 s → φµ + µ − is closely related to the decay B 0 → K * 0 µ + µ − , which has been studied extensively by LHCb [6,9,10]. Although B 0 s meson production is suppressed with respect to the B 0 meson by the fragmentation fraction ratio f s /f d ∼ 1/4, the narrow φ resonance allows a clean selection with low background levels. Furthermore, the contribution from the S wave, where the K + K − system is in a spin-0 configuration, is expected to be low [11]. Since the K + K − µ + µ − final state is not flavour-specific, the angular observables accessible in the decay B 0 s → φµ + µ − are the CP averages F L , S 3,4,7 and the CP asymmetries A 5,6,8,9 [12]. The flavour-averaged differential decay rate, as a function of the decay angles in bins of q 2 , is given by 1 4 (1 − F L ) sin 2 θ K cos 2θ l − F L cos 2 θ K cos 2θ l + S 3 sin 2 θ K sin 2 θ l cos 2Φ + S 4 sin 2θ K sin 2θ l cos Φ + A 5 sin 2θ K sin θ l cos Φ + A 6 sin 2 θ K cos θ l + S 7 sin 2θ K sin θ l sin Φ + A 8 sin 2θ K sin 2θ l sin Φ + A 9 sin 2 θ K sin 2 θ l sin 2Φ .
(1) The T-odd CP asymmetries A 8 and A 9 are predicted to be close to zero in the SM and are of particular interest, as they can be large in the presence of contributions beyond the SM [12].

Detector and simulation
The LHCb detector [13,14] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex, the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger [15], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. Simulated signal samples are used to determine the effect of the detector geometry, trigger, reconstruction and selection on the signal efficiency. In addition, simulated background samples are used to determine the pollution from specific background processes. In the simulation, pp collisions are generated using Pythia [16,17] with a specific LHCb configuration [18]. Decays of hadronic particles are described by EvtGen [19], in which final-state radiation is generated using Photos [20]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [21,22] as described in Ref. [23]. Data-driven corrections are applied to the simulated samples to account for imperfect modelling of particle identification performance, the B 0 s meson transverse momentum spectrum and B 0 s vertexing quality, as well as track multiplicity.

Selection of signal candidates
The B 0 s → φµ + µ − signal candidates are required to satisfy the hardware trigger requirement, which selects muons with p T > 1.48 GeV/c in the 7 TeV data and p T > 1.76 GeV/c in the 8 TeV data. In the subsequent software trigger, at least one of the final-state particles is required to have both p T > 0.8 GeV/c and IP larger than 100 µm with respect to all of the primary pp interaction vertices (PVs) in the event. The tracks of two or more of the final-state particles are also required to form a vertex that is significantly displaced from any PV.
Signal candidates are accepted if their reconstructed invariant mass is in the range 5267 < m(K + K − µ + µ − ) < 5800 MeV/c 2 and the invariant mass of the K + K − system is within 12 MeV/c 2 of the known φ mass [24]. The mass resolutions are 19 MeV/c 2 for the invariant K + K − µ + µ − mass and 4 MeV/c 2 for the K + K − invariant mass. The final-state particles are required to have significant χ 2 IP with respect to any PV in the event, where χ 2 IP denotes the change in the χ 2 of the PV when reconstructed with or without the considered track. The four final-state tracks are then fitted to a common vertex which is required to be of good quality and significantly displaced from any PV in the event. The signal candidate is required to have small χ 2 IP with respect to a PV in the event. Furthermore, the angle θ DIRA between the reconstructed B 0 s momentum and the vector connecting the PV with the decay vertex is required to be small.
To further reduce the combinatorial background, a boosted decision tree (BDT) [25] using the AdaBoost algorithm [26] is employed. The BDT is trained with a sample of B 0 s → J/ψ (→ µ + µ − )φ(→ K + K − ) decays as a signal proxy and events from the upper mass sideband, 5567 < m(K + K − µ + µ − ) < 5800 MeV/c 2 , as a proxy for the background. The discriminating variables of the BDT are the χ 2 IP of the B 0 s signal candidate and all final-state tracks, the B 0 s transverse momentum, the χ 2 of the vertex fit (χ 2 Vtx ), the flight distance significance of the signal candidate, and particle identification information for the final-state particles. The BDT selection has an efficiency of 96% for signal events with a background rejection of 95%. The total efficiency for signal events, including detector geometry, trigger and reconstruction effects is 1.1%.
The reconstructed B 0 s mass versus q 2 for signal candidates after the full selection is given in Fig. 2. The signal decay B 0 s → φµ + µ − is clearly visible as a vertical band. In the q 2 region around the J/ψ and ψ(2S) masses, the tree-level charmonium decays B 0 s → J/ψ φ and B 0 s → ψ(2S)φ dominate. The decay B 0 s → J/ψ φ is used as a control mode throughout the analysis.

Backgrounds
The decays B 0 s → J/ψ φ and B 0 s → ψ(2S)φ, primarily originating from b → ccs tree-level processes, are vetoed by rejecting candidates in the q 2 regions 8.0 < q 2 < 11.0 GeV 2 /c 4 and 12.5 < q 2 < 15.0 GeV 2 /c 4 . The B 0 s → J/ψ φ decay can also constitute a peaking background if one of the final-state muons is misidentified as a kaon and vice-versa. This background is vetoed by rejecting candidates for which the invariant mass of the K ± µ ∓ system, with the kaon reconstructed under the muon mass hypothesis, is within 45 MeV/c 2 of the known J/ψ meson mass [24], unless the final-state particles fulfil stringent particle identification requirements. After the veto is applied, this background contribution is found to be negligible.
The rare baryonic decay Λ 0 b → Λ(1520)(→ pK − )µ + µ − can mimic the signal decay if the proton in the final state is misidentified as a kaon. This potential background is vetoed by rejecting events with invariant mass close to the known Λ 0 b baryon mass [24] where one kaon has the proton mass hypothesis assigned, unless the kaon passes stringent particle identification requirements. Assuming a q 2 dependence following Ref. [27] and using B(Λ 0 b → Λµ + µ − ) = (0.96 ± 0.29) × 10 −6 [28] as an estimate for the unknown Λ 0 b → Λ(1520)µ + µ − branching fraction, a yield of 2.0 ± 0.8 Λ 0 b → Λ(1520)µ + µ − background events is expected in the signal region, within 50 MeV/c 2 of the known B 0 s mass [24], after the veto. The rare decay B 0 → K * 0 µ + µ − can be a peaking background if the pion in the final state is reconstructed as a kaon. After suppressing this background using particle identification information, a yield of 1.7 ± 0.4 events is expected in the signal region. The background pollution from Λ 0 b → Λ(1520)µ + µ − and B 0 → K * 0 µ + µ − decays is neglected in the fit and treated as a systematic uncertainty. Backgrounds from semileptonic b → c(→ sµ −ν µ )µ + ν µ cascade decays and fully hadronic decays such as where hadrons are misidentified as muons, are estimated to be small and can therefore be neglected. Figure 3 shows the K + K − µ + µ − invariant mass distribution for the B 0 s → φµ + µ − signal decay integrated over q 2 , as well as for the control mode B 0 s → J/ψ φ. To determine the B 0 s → φµ + µ − signal yields in bins of q 2 , extended maximum likelihood fits are performed. The combinatorial background is described by an exponential function, whilst the signal component is modelled with the sum of two Gaussian functions with a common mean and a radiative power-law tail toward smaller invariant mass values. The parameters describing the signal mass shape are determined from a fit to the B 0 s → J/ψ φ control mode. The q 2 dependence of the signal mass resolution is accounted for by using scale factors, which are determined from simulation. The K + K − µ + µ − invariant mass distributions for the signal decay in bins of q 2 are given in Appendix A; the yields and corresponding uncertainties are listed in Table 1. Integrating over the q 2 bins, the signal yield is found to be 432 ± 24. A fit to the control mode B 0 s → J/ψ φ, which is used for normalisation, gives N J/ψ φ = 62 033 ± 260 decays.

Differential branching fraction
The differential branching fraction for a given q 2 bin [q 2 min , q 2 max ] is calculated according (2) where N φµµ and N J/ψ φ denote the yield of the signal and normalisation mode, and φµµ and J/ψ φ their respective efficiencies. The branching fractions are given by B(J/ψ → µ + µ − ) = (5.961 ± 0.033) × 10 −2 [24] and B(B 0 s → J/ψ φ) = (10.76 ± 0.81) × 10 −4 . For the branching fraction of the normalisation channel B 0 s → J/ψ φ the LHCb measurement [11] is recalculated using an updated measurement of f s /f d = 0.259 ± 0.015 [29]. A weighted average is calculated by combining this updated measurement with the measurements by Belle [30] and CDF [31]. The resulting relative and absolute differential branching fractions are given in Table 1.
The differential branching fraction is also shown in Fig. 4, overlaid with SM predictions from Refs. [4,5]. In the q 2 region 1.0 < q 2 < 6.0 GeV 2 /c 4 the measured differential branching fraction lies 3.3 σ below the SM expectation of (4.81 ± 0.56) × 10 −8 GeV −2 c 4 [4,32]. For the SM predictions, the form factors are determined in a combined fit to the results of light-cone sum rule calculations at low q 2 [5] and lattice QCD calculations at high q 2 [33,34]. Standard Model predictions for the branching fraction at high q 2 that exclusively use the results from lattice calculations [35] are found to be larger than the results from the combined fit. Owing to their proximity to the charmonium resonances, no predictions are available corresponding to the q 2 bins 5.0 < q 2 < 8.0 GeV 2 /c 4 and 11.0 < q 2 < 12.5 GeV 2 /c 4 .
The total branching fraction of the signal decay is given by the integral over the six q 2 bins. To account for the fraction of signal events in the vetoed q 2 regions, a correction factor f veto = 1.520 ± 0.003 ± 0.043 is applied, which is determined using the calculation in Ref. [36] with updated form factors from Ref. [37]. The first given uncertainty is statistical, the second is systematic.
The resulting relative and total branching fractions are where the uncertainties are (from left to right) statistical, systematic, and from the extrapolation to the full q 2 region. For the total branching fraction, a further uncertainty originates from the uncertainty on the branching fraction of the normalisation mode.

Systematic uncertainties
For the branching fraction ratio B(B 0 s → φµ + µ − )/B(B 0 s → J/ψ φ), systematic uncertainties are mostly due to uncertainties on the efficiency ratio J/ψ φ / φµµ , which is taken from simulation. To evaluate the size of uncertainties affecting the efficiency ratio, it is recalculated after applying the corresponding systematic variation to the simulated samples. The observed deviation is taken as systematic uncertainty. The procedure to correct the tracking efficiency in simulation introduces a systematic uncertainty on the efficiency ratio of less than 0.6%. The correction to particle identification performance in simulation has a systematic uncertainty of 0.5%. The relative efficiency is further affected by the data-driven corrections to the simulation in the distribution of the variables p T (B 0 s ) and χ 2 Vtx (B 0 s ), as well as the track multiplicity, which have a combined systematic effect of 1.0%. The non-uniform angular acceptance detailed in Sec. 5 introduces a dependence of the signal efficiency on the underlying physics model. Its effect on the branching fraction Table 1: The signal yields for B 0 s → φµ + µ − decays, as well as the differential branching fraction relative to the normalisation mode and the absolute differential branching fraction, in bins of q 2 . The given uncertainties are (from left to right) statistical, systematic, and the uncertainty on the branching fraction of the normalisation mode.  measurement is evaluated by varying the Wilson coefficient C 9 used in the generation of simulated signal events. By allowing a New Physics contribution of −1.5, which is motivated by the global fit results in Ref. [38], the resulting systematic uncertainty is found to be less than 1.6%. The selection requirements introduce a decay-time dependence of the efficiencies which can, due to the sizeable lifetime difference in the B 0 s system [39], affect the measured branching fraction [40]. The systematic uncertainty is determined with simulated B 0 s → φµ + µ − signal events, generated using time-dependent decay amplitudes as described in Ref. [12]. When varying the Wilson coefficients, the size of the effect is found to be at most 1.6%, which is taken as the systematic uncertainty. The statistical uncertainty due to the limited size of the simulated signal samples leads to a systematic uncertainty of 1.9%.
The systematic uncertainties due to the parametrisation of the mass shapes are evaluated using pseudoexperiments. For the signal mass model, events are generated using a double Gaussian mass shape, and then fitted using both the double Gaussian as well as the nominal signal mass shape, taking the observed deviation as the systematic uncertainty. For the parametrisation of the combinatorial background, the nominal exponential function is compared with a linear mass model. The systematic uncertainties due to the modelling of the signal and background mass shape are 2.1% and 1.6%, respectively. Peaking backgrounds are neglected in the fit for determination of the signal yields. The main sources of systematic uncertainty are caused by contributions from the decays Λ 0 b → pK − µ + µ − and B 0 → K * 0 µ + µ − , resulting in systematic uncertainties of 0.2 − 2.2%, depending on the q 2 bin. Finally, the uncertainty on the branching fraction of the decay J/ψ → µ + µ − amounts to a systematic uncertainty of 0.6%. The complete list of systematic uncertainties is given in Table 2.
For the total branching fraction of the signal decay, the uncertainty on the branching fraction of the normalisation channel is the dominant systematic uncertainty, at the level of 7.5%. The uncertainty on the correction factor f veto to account for signal events that are rejected by the charmonium vetoes is estimated by varying the Wilson coefficients and form-factor parameters, leading to a systematic uncertainty of 2.9%.

Angular analysis
For the determination of the four CP averages F L , S 3,4,7 and the four CP asymmetries A 5,6,8,9 an unbinned maximum likelihood fit to the three-dimensional angular distribution and the K + K − µ + µ − invariant mass distribution is performed in each q 2 bin. The models described in Sec. 4 are used to parametrise the mass line shapes for signal and background. The angular distribution of the signal component is given by Eq. 1. The angular background distribution is described by the product of second-order Chebyshev polynomials in the three decay angles.
The non-uniform efficiency due to the reconstruction, triggering and selection of signal candidates distorts the angular distributions of the final-state particles, as well as the q 2 distribution. This acceptance effect is parametrised using Legendre polynomials, according to (cos θ l , cos θ K , Φ, q 2 ) = klmn c klmn P k (cos θ l )P l (cos θ K )P m (Φ)P n (q 2 ), where P i (x) denote Legendre polynomials of order i and c klmn the coefficients that are determined by performing a moments analysis using a large sample of simulated B 0 s → φµ + µ − signal events generated according to a phase-space model. The maximum order of the polynomials that is included is four for cos θ l , two for cos θ K , six for the angle Φ and five for q 2 . In addition, the acceptance is assumed to be symmetric in the decay angles. This choice corresponds to the lowest orders of polynomials that describe the acceptance effect. The acceptance description is cross-checked using the control mode B 0 s → J/ψ φ. An angular analysis of the control mode is performed and the angular observables are found to be in good agreement with the previous measurement [39].
Appendix B gives the one-dimensional angular distributions of the signal decay in each q 2 bin, overlaid with the projections of the likelihood fit. For the q 2 bins with the lowest number of signal candidates, pseudoexperiments show the likelihood estimator to be biased for certain observables due to boundary effects which arise from the requirement of Eq. 1 being positive for all values of the decay angles. Therefore, the Feldman-Cousins method [41] is used to determine confidence regions for each observable, which guarantees correct coverage for low signal yields. The remaining signal observables are treated as nuisance parameters following the plugin method [42]. The Feldman-Cousins scans for the angular observables in bins of q 2 are given in Appendix C. In some q 2 bins, the shape of the obtained confidence level is dominated by effects arising from the boundary condition. Table 3 gives the minima of the Feldman-Cousins scans and the 68% confidence intervals. The linear correlations between the angular observables in the different q 2 bins are given in Appendix D. The angular observables are shown in Fig. 5, overlaid with SM predictions from Refs. [4,5]. No predictions are given for S 7 and A 5,6,8,9 ; they are expected to be close to zero in the SM.

Systematic uncertainties
The systematic uncertainties on the angular observables are evaluated using large numbers of pseudoexperiments where simulated events are generated to reflect the measured angular distributions and event yield. To reduce statistical effects, each sample is fitted twice, once with and once without systematic variations. The angular acceptance correction, which is determined from simulation, is a significant source of systematic uncertainties for the angular observables. The data-driven corrections of the distributions of p T (B 0 s ), χ 2 Vtx (B 0 s ) and track multiplicity, as well as particle identification performance and tracking efficiency, amount to a systematic uncertainty of less than 0.01 in total. Furthermore, the kinematic distributions of the final-state particles are cross-checked using the control mode B 0 s → J/ψ φ. Correcting for kinematic differences between data and simulation amounts to a systematic deviation of less than 0.01. The effect of the limited size of the simulated signal samples on the acceptance  description is evaluated by varying the Legendre coefficients in Eq. 3 according to their corresponding covariance matrix. The resulting systematic uncertainty is smaller than 0.02 for all observables and q 2 bins. The four-dimensional acceptance correction is evaluated at the centre of each q 2 bin. To estimate the systematic effect due to this, an alternative acceptance description is used, where a separate three-dimensional acceptance is used for each q 2 bin. The resulting systematic deviation is negligible. The systematic effect of neglecting peaking backgrounds is evaluated by performing toy studies, where simulated Λ 0 b → pK − µ + µ − and B 0 → K * 0 µ + µ − background events are added, according to their expected yields for the specific q 2 bin. The resulting systematic deviations of the angular observables are smaller than 0.01 for all observables and q 2 bins. The S-wave pollution for the decay B 0 s → φµ + µ − is expected to be similar to that of the B 0 s → J/ψ φ decay, at the level of 1.1% [11]. The effect of the S wave in the K + K − system on the angular observables is determined to be smaller than 0.01 using toy studies. The combinatorial background is described using second-order Chebyshev polynomials determined from the upper mass sideband. The systematic uncertainty associated with this model choice is estimated by using first-order polynomials as an alternative. With a systematic effect of up to 0.04 on the angular observables, depending on q 2 bin, this constitutes the dominant systematic uncertainty for the angular analysis. In addition, the effect of fixing the angular background parameters in the nominal fit is evaluated using toy studies. The systematic deviation is found to be smaller than 0.02 for all observables and q 2 bins.
The total systematic uncertainty, given by the quadratic sum over all systematic effects, is found to be small compared to the statistical uncertainties for all angular observables in all q 2 bins.

Conclusions
Measurements of the differential branching fraction and the first full three-dimensional angular analysis of the decay B 0 s → φµ + µ − are presented, using data collected by the LHCb experiment in pp collisions, corresponding to an integrated luminosity of 3.0 fb −1 . The results are given in Tables 1 and 3 and are the most precise measurements of these quantities to date. The CP -averaged angular observables S 4 and S 7 are determined for the first time for this decay. The determination of the CP asymmetries A 5 and A 8 constitutes the first measurement of these quantities for any rare b → s decay, providing additional constraints in global fits. All angular observables are found to be compatible with SM predictions.
The B 0 s → φµ + µ − branching fraction relative to the normalisation mode B 0 s → J/ψ φ is measured to be and the resulting total absolute branching fraction is measured to be where the uncertainties are (from left to right) statistical, systematic, and from the extrapolation to the full q 2 region. For the total branching fraction, a further uncertainty originates from the uncertainty on the branching fraction of the normalisation mode. The measured branching fraction is compatible with the previous measurement [3] and lies below SM expectations. For the q 2 region 1.0 < q 2 < 6.0 GeV 2 /c 4 the differential branching fraction of (2.58 +0.33 −0.31 ± 0.08 ± 0.19) × 10 −8 GeV −2 c 4 is more than 3 σ below the SM prediction of (4.81 ± 0.56) × 10 −8 GeV −2 c 4 [4,5,32].    Figure 16: Confidence level obtained from a likelihood scan (shaded blue) and from a Feldman-Cousins method (solid black). The shaded blue and solid red vertical lines indicate the corresponding 68% CL intervals obtained from the likelihood scan and the Feldman-Cousins method, respectively. Table 4: Correlation matrices for the q 2 bins 0.1 < q 2 < 2.0 GeV 2 /c 4 , 2.0 < q 2 < 5.0 GeV 2 /c 4 and 5.0 < q 2 < 8.0 GeV 2 /c 4 .