Main

Neutral mesons with strange, charm or beauty quantum numbers can mix with their antiparticles, as these quantum numbers are not conserved by the weak interaction. The neutral meson comprising an antibeauty quark and a strange quark, namely, the \({B}_{\mathrm{s}}^{0}\) meson, and its antiparticle, namely, the \({\overline{B}}_{\mathrm{s}}^{0}\) meson, are one such example. In the \({B}_{\mathrm{s}}^{0}\)\({\overline{B}}_{\mathrm{s}}^{0}\) system, the observed particle and antiparticle states are linear combinations of the heavy (H) and light (L) mass eigenstates. These mass eigenstates have corresponding masses mH and mL and decay widths ΓH and ΓL (ref. 1). As a consequence, the \({B}_{\mathrm{s}}^{0}\)\({\bar{B}}_{\mathrm{s}}^{0}\) system oscillates with a frequency given by the mass difference, Δms = mH − mL. This oscillation frequency is an important parameter of the standard model of particle physics. In combination with the B0 –\({\overline{B}}^{0}\) oscillation frequency, Δmd, it provides a powerful constraint on the Cabibbo–Kobayashi–Maskawa quark-mixing matrix2,3,4,5,6. A precise measurement of Δms is also required to reduce the systematic uncertainty associated with measurements of matter–antimatter differences in the \({B}_{\mathrm{s}}^{0}\)\({\overline{B}}_{\mathrm{s}}^{0}\) system7.

In this Article, we present a measurement of Δms using \({B}_{\mathrm{s}}^{0}\) mesons that decay to a charmed–strange \({D}_{\mathrm{s}}^{-}\) meson and a pion, \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\), and decays with the opposite charge, \({\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\). We refer to both charge combinations as \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) throughout the paper, and similarly for decays of the \({D}_{\mathrm{s}}^{-}\) meson. The measurement is performed using data collected between 2015 and 2018, denoted as Run 2 of the Large Hadron Collider (LHC), corresponding to an integrated luminosity of 6 fb−1 of proton–proton (pp) collisions at a centre-of-mass energy of 13 TeV.

The first measurement in which the significance of the observed \({B}_{\mathrm{s}}^{0}\)\({\overline{B}}_{\mathrm{s}}^{0}\) oscillation signal exceed five standard deviations was obtained by the CDF collaboration8. More recently, the LHCb collaboration has performed several measurements of Δms using data collected at the LHC: a measurement using \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decays9; two measurements using \({B}_{\mathrm{s}}^{0} \to J/\psi {K}^{+}{K}^{-}\) decays10,11; and a measurement using \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{\mp }{\pi }^{\pm }{\pi }^{\pm }{\pi }^{\mp }\) decays12. Theoretical predictions for Δms are available6,1317, with the most precise prediction provided in ref. 18. The prediction is consistent with—but considerably less precise than—existing experimental results.

The \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decay time distribution, in the absence of detector effects, can be written as

$$P(t) \approx {\mathrm{e}}^{-{{{\varGamma }}}_{\mathrm{s}}t}\left[\cosh \left(\frac{{{\Delta }}{\varGamma }_{\mathrm{s}}t}{2}\right)+C \cos ({{\Delta }}{m}_{\mathrm{s}}t)\right]\ ,$$
(1)

where Γs = (ΓH + ΓL)/2 is the \({B}_{\mathrm{s}}^{0}\) meson decay width and ΔΓs = ΓH − ΓL is the decay-width difference between the heavy and light mass eigenstates. The parameter C takes the value C = 1 for unmixed decays, that is, \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\), and C = − 1 for decays in which the initially produced meson mixed into its antiparticle before decaying, that is, \({B}_{\mathrm{s}}^{0}\to {\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\). The mixed decay is referred to as \({\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) throughout the paper. The mass difference Δms corresponds to a frequency in natural units, and is measured in inverse picoseconds (ps–1).

The LHCb detector19,20 is designed to study the decays of beauty and charm hadrons produced in pp collisions at the LHC. It instruments a region around the proton beam axis, covering the polar angles between 10 and 250 mrad, in which approximately a quarter of the b-hadron decay products are fully contained. The detector includes a high-precision tracking system with a dipole magnet, providing measurements of the momentum and decay-vertex position of particles. Different types of charged particle are distinguished using information from two ring-imaging Cherenkov detectors, a calorimeter and a muon system.

Simulated samples of \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decays and data control samples are used to verify the analysis procedure and to study systematic effects. The simulation provides a detailed model of the experimental conditions, including the pp collision, decays of the produced particles, their final-state radiation and response of the detector. Simulated samples are corrected for residual differences in relevant kinematic distributions to improve the agreement with data. The software used is described in refs. 21,22,23,24,25,26,27.

The \({B}_{\mathrm{s}}^{0}\) mesons travel a macroscopic distance at LHC energies (1 cm on average) before decaying and are considerably heavier than most other particles produced directly in pp collisions. Thus, their decay products have large displacement relative to the pp collision point and a larger momentum transverse to the beam axis compared with other particles. The candidate selection exploits these fundamental properties. Two fast real-time selections use partial detector information to reject LHC bunch crossings likely to be incompatible with the presence of the signal, before a third selection uses fully aligned and calibrated data in real time to reconstruct and select topologies consistent with the signal28. Selected collisions are recorded to permanent storage. All selections, except the first real-time selection, are based on multivariate classifiers. Two subsequent selections fully reconstruct the decays with the \({D}_{\mathrm{s}}^{-}\) meson reconstructed in both KK+π and ππ+π final states. After the real-time stages, the initial ‘offline’ selection is based on track kinematic quantities and displacement relative to the pp collision point that favours the signal, followed by a multivariate classifier trained on properties of the full signal decay. These selections sequentially improve the signal purity of the sample to the final value of 84%, which is optimized using a simulation to maximize the product of signal significance and signal efficiency. This criterion gives the optimal sensitivity to the oscillation frequency.

The remaining sources of background after selection consist of random track combinations (combinatorial background); \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{* -}{\pi }^{+}\) decays, where the photon from the \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}\gamma\) decay is not reconstructed; and contributions from b-hadron decays with similar topologies to the signal, namely, B0Dπ+, \({\overline{\varLambda }}_{\mathrm{b}}^{0}\to {\overline{\varLambda }}_{{{{\rm{c}}}}}^{-}{\pi }^{+}\) and \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{\mp }{K}^{\pm }\) decays. Decays with similar topology are suppressed by applying kinematic vetoes and additional particle identification requirements.

To measure Δms, a \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decay time distribution is first constructed in the absence of a background. This is achieved by performing an unbinned two-dimensional likelihood fit to the observed \({D}_{\mathrm{s}}^{-}{\pi }^{+}\) and KK+π or ππ+π invariant mass distributions. This fit is used to determine the signal yield and a set of weights29 used to statistically subtract the background in the subsequent fit to the decay time distribution. The invariant mass distributions of the selected decays are shown in Fig. 1. The non-peaking contribution in the combinatorial background distribution, as shown in Fig. 1 (right), is due to events in which a fake \({D}_{\mathrm{s}}^{-}\) candidate is produced from a combination of random tracks. The peaking contribution is due to genuine \({D}_{\mathrm{s}}^{-}\) candidates combined with a random track resulting in a fake \({B}_{\mathrm{s}}^{0}\) candidate.

Fig. 1: Invariant mass distributions.
figure 1

Distributions of the \({D}_{s}^{-}{\pi }^{+}\) (left) and K+Kπ± or π+ππ± (right) invariant mass for the selected candidates, namely, \(m({D}_{\mathrm{s}}^{-}{\pi }^{+})\) and m(K+Kπ±, π+ππ±), respectively. The mass fit described in the text is overlaid. The different contributions are shown as coloured areas (for the background) or by dashed lines (for the signal). The vertical bars, typically visible only in regions with low numbers of candidates, correspond to the statistical uncertainty in the number of observed candidates in each bin. The horizontal bin width is indicated on the vertical axis legend.

The probability density functions describing the invariant mass distributions of the signal and background are obtained using a mixture of control samples in the data and simulation. The \({D}_{s}^{-}{\pi }^{+}\) and KK+π or ππ+π invariant-mass signal shapes are described by the sum of Hypatia30 and Johnson SU (ref. 31) functions. The combinatorial background contribution for both invariant mass distributions is described by an exponential function in each of them with the parameters determined in the fit. The B0Dπ+, \({\overline{\varLambda }}_{\rm{b}}^{0}\to {\overline{\varLambda }}_{{{{\rm{c}}}}}^{-}{\pi }^{+}\) or \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{\mp }{K}^{\pm }\) background components constitute less than 2% of the signal yield and are accounted for in the fit to the invariant mass distributions using yields obtained from known branching fractions and relative efficiencies, as determined from the simulated samples, which are weighted to account for differences between the data and simulation. The \({B}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) and \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{* -}{\pi }^{+}\) background components are also obtained from the simulated samples and included in the mass fit. The combined \({B}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) and \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{* -}{\pi }^{+}\) yield is a free parameter of the fit. The signal yield obtained from the invariant mass fit is 378,700 ± 700.

Decay time parameterization (equation (1)) is modified to account for the following detector effects: decay-time-dependent reconstruction efficiency; time-dependent decay time resolution; imperfect knowledge of the initial flavour of the reconstructed \({B}_{\mathrm{s}}^{0}\) or \({\overline{B}}_{\mathrm{s}}^{0}\) meson; asymmetry in \({B}_{\mathrm{s}}^{0}\) or \({\overline{B}}_{\mathrm{s}}^{0}\) production in pp collisions; and asymmetry in the reconstruction of final-state particles due to interactions in the detector material32.

Due to the lifetime biasing effect of the selections, the reconstruction efficiency is low at small decay times and increases to a plateau after 2 ps. The decay-time-dependent reconstruction efficiency is modelled with cubic b-spline curves as described elsewhere33. The spline coefficients are allowed to vary in the fit to the observed decay time distribution.

The decay time resolution is measured using a data sample of \({D}_{\mathrm{s}}^{-}\) mesons originating from pp interactions without being required to come from an intermediate \({B}_{\mathrm{s}}^{0}\) meson decay. These ‘prompt’ candidates pass the same real-time selection procedure as the signal sample. After real-time selection, additional requirements are applied to ensure the selection of a \({D}_{\mathrm{s}}^{-}\) signal peak with high background rejection but without any requirement on displacement from the pp collision point. The multivariate classifier trained using the full signal decay is not applied. The reconstructed decay time in this control sample is proportional to the distance between the \({D}_{\mathrm{s}}^{-}\) production vertex and an artificial \({B}_{\mathrm{s}}^{0}\) decay vertex, formed by combining the prompt \({D}_{\mathrm{s}}^{-}\) meson with a π+ track from the same pp collision. It is, therefore, compatible with a zero decay time up to bias and resolution effects. A linear relationship is observed between the decay time resolution measured at zero decay time and the decay time uncertainty estimated in the vertex fit. This relationship is used to calibrate the \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decay time uncertainty. Simulated prompt \({D}_{\mathrm{s}}^{-}\) and \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decays, for which the generated decay time is known, are used to check the suitability of this method, which determines a 0.005 ps bias in the reconstructed decay time due to residual detector misalignments. This bias is corrected in the analysis. The uncertainty in Δms due to these residual detector misalignments is evaluated using simulated samples that were intentionally misaligned. This uncertainty is reported in Table 1.

Table 1 Systematic uncertainties affecting the measurement of Δms

To determine if a neutral meson oscillated into its antiparticle, knowledge of the \({B}_{s}^{0}\) or \({\overline{B}}_{s}^{0}\) flavour at production and decay is required. In \({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\) decays, the \({B}_{s}^{0}\) flavour at decay is identified by the charge of the pion as the \({D}_{s}^{+}{\pi }^{-}\) decay cannot be directly produced. To determine whether the \({B}_{s}^{0}\) oscillated before decay, the flavour at production is inferred from the hadronization of the \({B}_{s}^{0}\) meson or the decay of other beauty hadrons produced in the collision using a combination of several flavour-tagging algorithms34,35,36,37. Each algorithm estimates the probability that a candidate has been assigned the wrong flavour tag. The algorithms that use information independent of signal fragmentation are calibrated using B+ meson decays, and a combined wrong-tag estimate is used in the fit. The tagging efficiency is measured to be ε = (80.30 ± 0.07)% with a probability to tag a candidate as the wrong flavour of ω = (36.21 ± 0.17)%, where the uncertainties account for the calibration.

In the unbinned maximum likelihood fit to the decay time distribution used to extract Δms, the calibration parameters for the combined wrong-tag estimate are allowed to vary. Additional free parameters are the values of the spline coefficients used to describe the decay-time-dependent reconstruction efficiency and the \({B}_{\mathrm{s}}^{0}\)\({\overline{B}}_{\mathrm{s}}^{0}\) production and detection asymmetries.

The parameters Γs and ΔΓs are fixed in the fit to their known values38. The other fixed parameters are the estimate of the wrong-tag fraction and efficiency of the tagging algorithms, decay-time bias correction and decay-time resolution calibration parameters. The decay time distribution of the tagged–mixed (\({\overline{B}}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\)), tagged–unmixed (\({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\)) and untagged (where the initial flavour is unknown) samples are shown in Fig. 2 (left). The corresponding fit projection is overlaid. To highlight the oscillation phenomenon, the asymmetry distribution between the tagged–unmixed and tagged–mixed samples is defined as

$$A(t)=\frac{N({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)-N({\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)}{N({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)+N({\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)},$$
(2)

with t modulo 2π/Δms, as shown in Fig. 2 (right). Here \(N({\overline{B}}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+},t)\) and \(N({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+},t)\) indicate the tagged–mixed and tagged–unmixed decays observed at time t, respectively. For this distribution, each event—in addition to the weight used to statistically subtract the background—is also weighted by the product of two factors. The first is a flavour-tagging dilution factor related to the probability that the flavour tag is indeed correct. The second is an effective decay-time uncertainty dilution factor depending on the reconstructed decay time per-event resolution and on Δms, for which the central value of the decay time fit is being used. The overlaid continuous line corresponds to the fit result. The result of the fit for Δms is 17.7683 ± 0.0051 ps−1, where this uncertainty is related to the sample size.

Fig. 2: Decay time distribution of the signal decays.
figure 2

Distribution of the decay time of the \({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\) signal decays (left) and decay time asymmetry between the mixed and unmixed signal decays (right). The vertical bars correspond to the statistical uncertainty in the number of observed candidates in each bin. The horizontal bars represent the bin width. In the left plot, the horizontal bin width is indicated on the vertical axis legend. The three components—unmixed, mixed and untagged—are shown in blue, red and grey, respectively. The inset shows a zoomed-in view of the region delineated in grey. The fit described in the text is overlaid.

Several sources of systematic uncertainty have been investigated and those with a non-negligible contribution are listed in Table 1. These include uncertainty in the momentum scale of the detector, obtained by comparing the reconstructed masses of known particles with the most accurate available values38; residual detector misalignment and length-scale uncertainties; and uncertainties due to choice of mass and decay-time fit models, determined using alternate parameterizations and pseudo-experiments. To verify the robustness of the measurement to variations in Δms as a function of decay kinematics, the data sample is split into mutually disjoint subsamples, each having about the same number of signal events, in relevant kinematic quantities, such as the \({B}_{\mathrm{s}}^{0}\) momentum, and the Δms values obtained from each subsample are compared. The largest observed variation is included as a systematic uncertainty. The total systematic uncertainty is 0.0032 ps−1, with the leading contribution due to residual detector misalignment and detector length-scale uncertainties.

The value of the \({B}_{\mathrm{s}}^{0}\)\({\overline{B}}_{\mathrm{s}}^{0}\) oscillation frequency determined in this article is as follows.

$${{\Delta }}{m}_{\mathrm{s}}\ \ =17.7683\pm 0.0051 \, ({{{\rm{stat}}}})\pm 0.0032 \, ({{{\rm{syst}}}})\,{{\mathrm{ps}}}^{-1}$$

This is the most precise measurement to date. The precision is further enhanced by combining this result with the values determined in refs. 9,12. One study9 used \({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\) decays collected in 2011. The other study12 used a sample of \({B}_{s}^{0} \to {D}_{s}^{-} {\pi }^{+} {\pi }^{+} {\pi }^{-}\) decays selected from the combined 2011–2018 dataset, corresponding to 9 fb−1. The measurements are statistically independent. The systematic uncertainties related to the momentum scale, length scale and residual detector misalignment are assumed to be fully correlated. Due to aging of the detector and different alignment procedures used in Run 1 and Run 2, the effect of residual detector misalignment is larger in measurements using the Run 2 data. Given the precision of the measurement described in this paper, a detailed study of the detector misalignment effects is performed, and the related uncertainty due to the decay time bias has been substantially reduced compared with previous measurements using the Run 2 data. The values of the fixed parameters ΔΓs and Γs used as inputs to the previous analyses have evolved over time as additional measurements have been made. However, as the correlation between Δms and ΔΓs and Γs is negligible, these small differences have been ignored in the combination procedure. A covariance matrix is constructed by adding statistical and systematic uncertainties in quadrature for each input, including correlations between systematic uncertainties. The results are averaged by minimizing χ2 from the full covariance matrix. The value of Δms obtained is 17.7666 ± 0.0057 ps−1. Additionally, these results are combined with those from refs. 10,11 where Δms is determined using \({B}_{s}^{0} \to {{{\rm{J}}}} /\psi {K}^{+}{K}^{-}\) decays in the 2011–2012 (3 fb−1) and 2015–2016 (2 fb−1) datasets. The decay time resolution for the measurements used in the combination912, including the analysis presented here, varies from 35 to 45 fs, depending on the decay mode. The result for Δms is 17.7656 ± 0.0057 ps−1. The different measurements, as well as the resulting combination, are shown in Fig. 3.

Fig. 3: Summary of LHCb measurements.
figure 3

Comparison of LHCb Δms measurements obtained from refs. 912, the result presented in this article, and their average. The measurement described in this paper is labelled as \({D}_{\mathrm{s}}^{-}{\pi }^{+}\) 6 fb−1. The horizontal bars correspond to the total uncertainty reported for each measurement. The band indicates the size of the uncertainty on average for comparison purposes. The combination procedure and inputs are described in the text.

In summary, this paper presents the most precise measurement of the Δms oscillation frequency, namely, 17.7683 ± 0.0051 (stat) ± 0.0032 (syst) ps−1, where the first uncertainty is statistical and the second is systematic. The result is obtained using a sample of \({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+}\) decays collected with the LHCb detector during Run 2 of the LHC. Combining the result of this paper with previous measurements by the LHCb collaboration yields Δms = 17.7656 ± 0.0057 ps−1. This value is compatible with, and considerably more precise than, the predicted value from lattice quantum chromodynamics (refs. 13,14,15) and sum-rule calculations16,17 of \(18.{4}_{-1.2}^{+0.7}\,{{\mathrm{ps}}}^{-1}\) (ref. 18). The combined result represents a considerable improvement over previous measurements, and is a legacy measurement of the original LHCb detector. The experiment is currently undergoing a major upgrade to operate at five times the instantaneous luminosity from 2022 onwards39. The largest sources of systematic uncertainty for this measurement, that is, those related to the detector length scale and misalignment, will be a focal point to further improve on this result for future data-taking periods.

Methods

The LHCb detector

The LHCb detector19,20 is a single-arm forward spectrometer covering the pseudo-rapidity range of 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 region40, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 T m, and three stations of silicon-strip detectors and straw drift tubes41 placed downstream of the magnet. The tracking system provides a measurement of momentum p of the charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV c–1. The minimum distance of a track to a primary vertex (PV), that is, the impact parameter, is measured with a resolution of (15 + 29/pT) μm, where pT is the component of the momentum transverse to the beam (unit, GeV c–1). Different types of charged hadron are distinguished using information from two ring-imaging Cherenkov detectors42. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors and an electromagnetic and hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers43.

Simulation of the LHCb detector response is required to model the effects of detector acceptance and imposed selection requirements. In the simulation, pp collisions are generated using PYTHIA21 with a specific LHCb configuration22. Decays of unstable particles are described by EvtGen23, in which the final-state radiation is generated using PHOTOS27. The interaction of the generated particles with the detector and its response are implemented using the GEANT4 toolkit24,25 described in ref. 26.

Selection

A fast decision about which pp collisions are of interest is made by a trigger system44. It consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage that reconstructs the pp collision based on all the available detector information. The software trigger selects candidates consistent with a b-hadron decay topology, with tracks originating from a vertex detached from the primary pp collision point, known as the PV. The mean \({B}_{\mathrm{s}}^{0}\) lifetime is 1.5 ps (ref. 38), which corresponds to an average flight distance of 1 cm in the LHCb detector.

After being accepted by the trigger, a further selection is applied that forms \({D}_{s}^{-}\to {{{{\rm{K}}}}}^{-}{{{{\rm{K}}}}}^{+}{\pi }^{-}\) and \({D}_{s}^{-}\to {\pi }^{-}{\pi }^{+}{\pi }^{-}\) candidates from the reconstructed charged tracks and subsequently combines them with a fourth track to form \({B}_{s}^{0}\to {D}_{s}^{-}{\pi }^{+}\) candidates. Particle identification information is used to assign mass hypotheses to each of the final-state tracks.

The \({B}_{\mathrm{s}}^{0}\) invariant mass resolution is improved by constraining the \({D}_{\mathrm{s}}^{-}\) invariant mass to its known value38. The K+Kπ± or π+ππ± and \({D}_{s}^{-}{\pi }^{+}\) invariant mass ranges considered in this analysis are [1,930, 2,015] and [5,300, 5,800] MeV c–2, respectively.

To suppress the \({B}_{s}^{0}\) candidates formed from random track combinations, a gradient boosted decision tree (BDT) is used, implemented using the XGBoost library45. Training uses data from both signal and background samples to avoid mismatches between the data and simulation. This classifier uses information on the fit quality of the \({D}_{s}^{-}\) and \({B}_{s}^{0}\) decay vertices; \({D}_{s}^{-}\) and \({B}_{s}^{0}\,\)\({\chi }_{{{{\rm{IP}}}}}^{2}\) defined as the difference in χ2 of the vertex fit for a given PV reconstructed with and without the considered particle; angles between their momentum vector and vector connecting their production and decay vertices; and pT and impact parameter \({\chi }_{{{{\rm{IP}}}}}^{2}\) of the final-state tracks. The BDT classifier threshold is chosen to maximize the product of the signal significance and signal efficiency. This choice optimizes sensitivity to the oscillation frequency.

Flavour tagging

The initial flavour of the \({B}_{\mathrm{s}}^{0}\) meson must be known to determine if it has oscillated before decay. Flavour-tagging algorithms are used to determine the initial flavour from properties of the b-hadron production in the pp collision.

Beauty quarks are predominantly produced in pairs. Opposite-side (OS) tagging algorithms35 determine the initial flavour of the \({B}_{s}^{0}\) meson based on information from the other beauty quark decay. These include the OS muon and OS electron taggers, which identify the flavour from the charge of leptons produced in the other b-hadron decay. The OS kaon tagger identifies bcs transitions, the OS charm quark tagger identifies bc transitions, and the OS vertex charge tagger calculates the effective charge of an OS displaced vertex36. In addition, a same-side kaon tagger exploits the charge information of the kaon originating from the \(\bar{s}\) or s quark leftover from the \({B}_{s}^{0}\) or \({\overline{B}}_{s}^{0}\) meson fragmentation37. Each algorithm determines the initial flavour of the \({B}_{s}^{0}\) meson from the charge of the reconstructed tagging particle or the reconstructed vertex in case of the OS vertex tagger.

The tagging information is incorporated in the decay time description. The amplitude of oscillation is reduced by a dilution factor D = (1 − 2ω), where ω is the average fraction of incorrect tags (known as mistag rate in the literature). Different machine-learning algorithms provide an estimate of the mistag rate that is calibrated with the data to match the true mistag distribution. A linear calibration of the average mistag as a function of the predicted mistag for the combined OS tag and same-side kaon tag information is then implemented in the decay time fit with freely varying calibration parameters. The combined tagging efficiency of the sample is ε = (80.30 ± 0.07)% with a mistag fraction of ω = (36.21 ± 0.02 ± 0.17)%, where the first uncertainty is due to the finite size of the calibration sample and the second one is due to the calibration procedure. This results in a combined effective performance of (6.10 ± 0.02 ± 0.15)% with respect to a perfect tagging algorithm that would have 100% tagging efficiency and zero mistag rate.

Decay time fit

The observed decay time distribution is fitted using an unbinned maximum likelihood fit in which all the combinations of initial-state flavours (\({B}_{\mathrm{s}}^{0}\), \({\overline{B}}_{\mathrm{s}}^{0}\) or untagged) and final-state charges (\({D}_{\mathrm{s}}^{-}\)π+ or \({D}_{\mathrm{s}}^{+}{\pi }^{-}\)) are simultaneously fit. The decay time distribution of each measured final state is described by the sum of all the processes contributing to that state. Experimental effects are taken into account with several adjustments to the theoretical prediction in equation (1), namely,

$$\begin{array}{l}{P}_{\exp}({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t) \approx (1+{a}_{\det })\\ \left[(1-\omega )(1-{a}_{{{{\rm{prod}}}}})P({B}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)+\omega (1+{a}_{{{{\rm{prod}}}}})P({\overline{B}}_{\mathrm{s}}^{0}\to {D}_{\mathrm{s}}^{-}{\pi }^{+},t)\right]\ .\end{array}$$
(3)

Production and detection asymmetries are parameterized by factors aprod and adet, respectively, which are allowed to deviate from unity. The decay time distribution of both flavours contain a fraction (1−ω) of the correctly tagged decay time parameterization plus a fraction ω of the incorrectly tagged decay time parameterization. The mistag rate ω is obtained from a per-event estimation, after linear calibration. Different calibration parameters are used for the \({B}_{\mathrm{s}}^{0}\) and \({\overline{B}}_{\mathrm{s}}^{0}\) initial flavours.

The experimental decay time distributions of both flavours are convolved with a Gaussian function to account for the finite detector resolution. The mean of this function is shifted by the decay-time bias correction factor, and the width is obtained from a per-event estimate of the decay time uncertainty after linear calibration.

A decay-time-dependent efficiency is finally modelled by a time-dependent cubic spline function, which multiplies the decay time distribution obtained from the previous step.

Systematic uncertainties

The following sources of systematic uncertainty have been found to give a non-negligible contribution to the Δms measurement. These are summarized in Table 1.

The measured decay time is inversely proportional to \({B}_{s}^{0}\) momentum and therefore depends on an accurate determination of the momentum-scale uncertainty of the tracking system. The uncertainty is determined by varying the momentum of the \({B}_{s}^{0}\) meson by ±0.03% (obtained from a comparison of masses of different particles with their known values) in the simulated signal samples. The corresponding uncertainty in Δms is 0.0007 ps−1.

The measured decay time is also proportional to the distance travelled by the \({B}_{\mathrm{s}}^{0}\) meson between production and decay, which is affected by the precise knowledge of the position of the vertex detector elements along the proton beam axis. The measured uncertainty is 100 μm over a length of 1 m (ref. 40). The corresponding uncertainty in Δms is 0.0018 ps−1.

The relative alignment of the tracking detector elements is a source of bias in decay time and contributes to resolution effects. The uncertainty in Δms due to imprecise knowledge of this alignment has been obtained from the analysis of the simulated signal samples in which the detector elements have been deliberately misaligned. Different misalignments, namely, translations, rotations and their combinations—have been investigated. The leading effect is due to translation along the x axis—the axis perpendicular to the beam direction pointing towards the centre of the LHC ring. As a consequence, the simulated signal samples have been misaligned with x-axis translations in the range between 0 and 9 μm, as determined from the survey results. Each misaligned simulated sample is then corrected for decay time bias in the same manner as the data, and the extracted Δms value is compared with the value obtained in simulation without any misalignment. This comparison produces a corresponding uncertainty in the bias correction procedure of 0.0020 ps−1.

Alternative parameterizations of the background contributions to the invariant mass fit have been obtained by using different weighting methods; the difference between these parameterizations corresponds to an uncertainty of 0.0002 ps−1.

For the specific \({B}_{s}^{0}\to {D}_{s}^{* -}{\pi }^{+}\) and \({B}^{0}\to {D}_{s}^{-}{\pi }^{+}\) background contributions, the relative fraction of these components cannot be reliably determined from the data. Their relative contributions are nominally set to an equal mixture and varied between 0 (pure \({B}^{0}\to {D}_{s}^{-}{\pi }^{+}\)) and 1 (pure \({B}_{s}^{0}\to {D}_{s}^{* -}{\pi }^{+}\)) to determine the maximum deviation in Δms corresponding to an uncertainty of 0.0005 ps−1.

The decay time resolution is obtained from data using a sample of \({D}_{s}^{-}\) mesons that are directly produced in a pp collision. These are combined with a π+ meson coming from the same collision to produce a fake \({B}_{s}^{0}\) candidate with a decay time equal to zero, ignoring the resolution effects. Different parameterizations of the measured decay time distribution are applied to a simulated signal sample. The largest deviation of the extracted Δms value with respect to nominal parameterization is found to be 0.0011 ps−1.

The procedure used to subtract background contributions in the fit to the decay time distribution assumes no large correlations between the decay time and the reconstructed \({B}_{s}^{0}\) and \({D}_{s}^{-}\) invariant masses. This is studied by analysing the simulated signal and background samples where any residual correlations between these observables are removed. The difference in measured value of Δms between the decorrelated and nominal samples is found to be 0.0011 ps−1.

The data sample was split into mutually disjoint subsamples to study the effect of potential correlations between the kinematic ranges, data-taking periods, flavour-tagging categories, BDT-based selection and measured value of Δms. The measured values obtained from each subsample are compared and the largest observed variation is found to be 0.0003 ps−1.

Several additional effects have been considered consisting of possible biases introduced by the fit procedure, changes to the signal and background parameterizations, and changes in the reweighting procedure used when obtaining the invariant mass shapes of partially reconstructed backgrounds constituting less than 2% of the signal yield. Their impact has been found to be negligible with respect to the sources listed in Table 1.

The largest sources of systematic uncertainty are found to be due to the imprecise knowledge of position and alignment of the tracking detector closest to the nominal pp collision region.