Measurement of the CP-violating phase ϕs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _s$$\end{document} in Bs0→J/ψϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ B_{s}^{0} \rightarrow J/\psi \phi $$\end{document} decays in ATLAS at 13 TeV

A measurement of the Bs0→J/ψϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ B_{s}^{0} \rightarrow J/\psi \phi $$\end{document} decay parameters using 80.5fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 80.5\, \mathrm {fb^{-1}} $$\end{document} of integrated luminosity collected with the ATLAS detector from 13 Te\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Te}\text {V}$$\end{document} proton–proton collisions at the LHC is presented. The measured parameters include the CP-violating phase ϕs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{s} $$\end{document}, the width difference ΔΓs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta \Gamma _{s}$$\end{document} between the Bs0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{s}^{0}$$\end{document} meson mass eigenstates and the average decay width Γs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Gamma _{s}$$\end{document}. The values measured for the physical parameters are combined with those from 19.2fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 19.2\, \mathrm {fb^{-1}} $$\end{document} of 7 and 8 Te\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Te}\text {V}$$\end{document} data, leading to the following: ϕs=-0.087±0.036(stat.)±0.021(syst.)radΔΓs=0.0657±0.0043(stat.)±0.0037(syst.)ps-1Γs=0.6703±0.0014(stat.)±0.0018(syst.)ps-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \phi _{s}= & {} -0.087 \pm 0.036 ~\mathrm {(stat.)} \pm 0.021 ~\mathrm {(syst.)~rad} \\ \Delta \Gamma _{s}= & {} 0.0657 \pm 0.0043 ~\mathrm {(stat.)}\pm 0.0037 ~\mathrm {(syst.)~ps}^{-1} \\ \Gamma _{s}= & {} 0.6703 \pm 0.0014 ~\mathrm {(stat.)}\pm 0.0018 ~\mathrm {(syst.)~ps}^{-1} \end{aligned}$$\end{document}Results for ϕs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{s} $$\end{document} and ΔΓs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta \Gamma _{s}$$\end{document} are also presented as 68% confidence level contours in the ϕs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{s} $$\end{document}–ΔΓs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta \Gamma _{s}$$\end{document} plane. Furthermore the transversity amplitudes and corresponding strong phases are measured. ϕs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{s} $$\end{document} and ΔΓs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta \Gamma _{s}$$\end{document} measurements are in agreement with the Standard Model predictions.


Introduction
In the presence of new physics (NP) phenomena, sources of CP violation in b-hadron decays can arise in addition to those predicted by the Standard Model (SM) [1]. In the B 0 s → J/ψφ decay, CP violation occurs due to interference between a direct decay and a decay with B 0 s -B 0 s mixing. The oscillation frequency of B 0 s meson mixing is characterised by the mass difference, m s , of the heavy (B H ) and light (B L ) mass eigenstates. The CP-violating phase φ s is defined as the weak phase difference between the B 0 s -B 0 s mixing amplitude and the b → ccs decay amplitude. In the SM the phase φ s is small and is related to the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix elements via the relation φ s −2β s , with β s = arg[−(V ts V * tb )/(V cs V * cb )]. By combining beauty and kaon physics observables, and assuming no NP contributions to B 0 s mixing and decays, a value of −2β s = −0.03696 +0.00072 −0.00082 rad was predicted by the e-mail: atlas.publications@cern.ch CKMFitter group [2] and −2β s = −0.03700 ± 0.00104 rad according to the UTfit Collaboration [3]. While large NP enhancements of the mixing amplitude have been excluded by the precise measurement of the oscillation frequency [4], the NP couplings involved in the mixing may still increase the size of the observed CP violation by enhancing the mixing phase φ s with respect to the SM value.
Other physical quantities involved in B 0 s -B 0 s mixing are the decay width s = ( L + H )/2 and the width difference s = L − H , where L and H are the decay widths of the light and heavy mass eigenstates, respectively. The latest predictions for the width difference in the SM are s = 0.091±0.013 ps −1 [5] and s = 0.092±0.014 ps −1 [6]. A potential NP enhancement of φ s would also decrease the size of s , but it is not expected to be affected as significantly as φ s [7]. Nevertheless, extracting s from the data is an important test of theoretical predictions [7].
Theory predictions have been made for the lifetime ratios τ (B 0 s )/τ (B d ) and τ (B 0 s )/τ (B + ), with the latest update Ref. [8]. The lifetime τ (B 0 s ) has not been calculated in theory yet at a precision comparable with those obtained by experiments. The current world combined value of the decay width, s , obtained from experimental results is s = 0.6600 ± 0.0016 ps −1 [9].
The analysis presented here introduces a measurement of the B 0 s → J/ψφ decay parameters using 80.5 fb −1 of the LHC proton-proton ( pp) data collected by the ATLAS detector during 2015-2017, at a centre-of-mass energy, √ s, equal to 13 TeV. The analysis closely follows a previous ATLAS measurement [13] that was performed using 19.2 fb −1 of the data collected at 7 and 8 TeV, and introduces more precise signal and background models.

ATLAS detector and Monte Carlo simulation
The ATLAS detector 1 consists of three main components: an inner detector (ID) tracking system immersed in a 2 T axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer (MS). The inner tracking detector covers the pseudorapidity range |η| < 2.5, and consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. The ID is surrounded by a high-granularity liquid-argon (LAr) sampling electromagnetic calorimeter. A steel/scintillator tile calorimeter provides hadronic coverage in the central rapidity range. The endcap and forward regions are equipped with LAr calorimeters for electromagnetic and hadronic measurements. The MS surrounds the calorimeters and provides a system of tracking chambers and detectors for triggering. A full description can be found in Refs. [22][23][24].
The data were collected during periods with different instantaneous luminosity, so several triggers were used in the analysis. All triggers were based on the identification of a J/ψ → μ + μ − decay, with transverse momentum ( p T ) thresholds of either 4 GeV or 6 GeV for the muons. Data quality requirements are imposed on the data, notably on the performance of the MS, ID and calorimeter systems. The measurement uses 80.5 fb −1 of pp collision data. The uncertainty in the combined 2015-2017 integrated luminosity is 2.0% [25], obtained using the LUCID-2 detector [26] for the primary luminosity measurements.
To study the detector response, estimate backgrounds, and model systematic effects, 100M Monte Carlo (MC) simulated B 0 s → J/ψφ events were generated using Pythia 8.210 [27] tuned with ATLAS data, using the A14 set of parameter values [28] together with the CTEQ6L1 set of parton distribution functions [29]. The detector response was simulated using the ATLAS simulation framework based on Geant4 [30,31]. In order to account for the varying number of proton-proton interactions per bunch crossing (pile-up) and trigger configurations during data-taking, the MC events were weighted to reproduce the same pile-up and trigger conditions as in the data. Additionally, background samples of both the exclusive (B 0 d → J/ψ K 0 * and b → J/ψ pK − ) and inclusive (bb → J/ψ X and pp → J/ψ X ) decays were simulated, using the same simulation tools as in the case of 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point. The z-axis is along the beam pipe, the xaxis points to the centre of the LHC ring and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, r being the distance from the origin and φ being the azimuthal angle around the beam pipe. The pseudorapidity η is defined as η = − ln[tan(θ/2)] where θ is the polar angle. the signal events. For validation studies related to flavour tagging, detailed in Sect. 4, events with B ± → J/ψ K ± exclusive decays were also simulated.

Reconstruction and candidate selection
The reconstruction and candidate selection for the decay B 0 s → J/ψ(μ + μ − )φ(K + K − ) is described here. Events must pass the trigger selections described in Sect. 2. In addition, each event must contain at least one reconstructed primary vertex, formed from at least four ID tracks, and at least one pair of oppositely charged muon candidates that are reconstructed using information from the MS and the ID. The muons used in the analysis are required to meet the Tight 2 or Low-p T 3 working point identification criteria. The muon track parameters are determined from the ID measurement alone, since the precision of the measured track parameters is dominated by the ID track reconstruction in the p T range of interest for this analysis. Pairs of oppositely charged muon tracks are re-fitted to a common vertex and the pair is accepted if the quality of the fit meets the requirement χ 2 /ndof < 10. In order to account for varying mass resolution in different parts of the detector, the J/ψ candidates are divided into three subsets according to the pseudorapidity η of the muons. In the first subset, both muons have |η| < 1.05, where the values η = ±1.05 correspond to the edges of the barrel part of the MS. In the second subset, one muon has 1.05 < |η| < 2.5 and the other muon |η| < 1.05. The third subset contains candidates where both muons have 1.05 < |η| < 2.5. A maximum likelihood fit is used to extract the J/ψ mass and the corresponding mass resolution for these three subsets, and in each case the signal region is defined symmetrically around the fitted mass, so as to retain 99.7% of the J/ψ candidates identified in the fits.
The candidates for the decay φ → K + K − are reconstructed from all pairs of oppositely charged tracks, with p T > 1 GeV and |η| < 2.5, that are not identified as muons. Candidate events for B 0 s → J/ψ(μ + μ − )φ(K + K − ) decays are selected by fitting the tracks for each combination of J/ψ → μ + μ − and φ → K + K − to a common vertex. The fit is also constrained by fixing the invariant mass cal-culated from the two muon tracks to the J/ψ mass [33]. A quadruplet of tracks is accepted for further analysis if the vertex fit has χ 2 /ndof < 3. For the φ → K + K − candidate, the invariant mass of the track pairs (using a charged kaon mass hypothesis) must fall within the interval 1.0085 GeV < m(K + K − ) < 1.0305 GeV. The interval, chosen using MC simulation, is selected to retain 98% of true φ → K + K − decays. The B 0 s candidate with the lowest χ 2 /ndof is selected in events where more than one candidate passes all selections. In total, 2 977 526 B 0 s candidates are collected within the mass range of 5. 150-5.650 GeV. This range is chosen to give enough background events in the sidebands of the mass distributions to allow precise determination of the properties of the background events.
The mean number of interactions per bunch crossing is 30, necessitating a choice of the best candidate for the primary vertex at which the B 0 s meson is produced. Primary vertex positions are recalculated after removing any tracks used in the B 0 s meson reconstruction. The variable used to select the best candidate for the primary vertex is the three-dimensional impact parameter, a 0 , which is calculated as the minimum distance between each primary vertex candidate and the line extrapolated from the reconstructed B 0 s meson vertex in the direction of the B 0 s momentum. The chosen primary vertex is the one with the smallest a 0 . A simulated dataset is used to estimate the fraction of B 0 s candidates where the incorrect production vertex is selected (12%) and demonstrates that the mis-selection of reconstructed primary vertex does not bias the reconstructed proper decay time.
For each B 0 s meson candidate the proper decay time t is estimated using: where p T B is the reconstructed transverse momentum of the B 0 s meson candidate and m B denotes the mass of the B 0 s meson, taken from Ref. [33]. The transverse decay length, L xy , is the displacement in the transverse plane of the B 0 s meson decay vertex relative to the primary vertex, projected onto the direction of the B 0 s transverse momentum.

Flavour tagging
To identify, or tag, the flavour of a neutral B meson at the point of production, information is extracted using the decay of the other (or opposite) b-hadron that is produced from the pair production of b andb quarks. This method is called opposite-side tagging (OST). The OST algorithms each define a discriminating variable, based on charge information, which is sensitive to the flavour (i.e. b-orb-quark) of the opposite-side b-hadron.
The algorithms thus provide a probability that a signal B meson in a given event is produced in a given flavour. The calibration of the OST algorithms proceeds using data containing B ± → J/ψ K ± candidate decays, where the charge of the kaon determines the flavour of the B meson, providing a self-tagging sample of events. These OST algorithms are calibrated as a function of the discriminating variable, using yields of signal B ± mesons extracted from fits to the data. Once calibrated, the OST algorithms are applied to B 0 s → J/ψ(μ + μ − )φ(K + K − ) candidate events to provide a probability that each candidate was produced in a B 0 s or B 0 s meson state, which is used in the maximum likelihood fit (described in Sect. 5). This approach assumes invariance of the OST algorithm with respect to the specific signal bhadron type (i.e B ± meson or B 0 s meson), which is tested and the difference is considered as a systematic uncertainty.
Candidate B ± → J/ψ K ± decays are identified in a series of steps. First, J/ψ candidates are selected from oppositely charged muon pairs forming a good vertex, as described in Sect. 3. Each muon is required to have p T > 4 GeV and |η| < 2.5. Dimuon candidates with invariant mass 2.8 < m(μ + μ − ) < 3.4 GeV, as determined from the re-fitted track parameters of the vertex, are retained for further analysis. To form the B ± candidate, an additional track is required, which is not identified as an electron or muon. The track is assigned the charged-kaon mass hypothesis and combined with the dimuon candidate using a vertex fit, performed with the mass of the dimuon pair constrained to the J/ψ mass. Prompt background contributions are suppressed by a requirement on the proper decay time of the B ± candidate of t > 0.2 ps.
The tagging probabilities are determined from B + and B − signal events. These signal yields are derived from fits to the invariant mass distribution, m(J/ψ K ± ), and performed in intervals of the discriminating variables. To describe the B ± → J/ψ K ± signal, two Gaussian functions with a common mean are used. An exponential function is used to describe the combinatorial background and a hyperbolic tangent function to parameterise the low-mass contribution from incorrectly or partially reconstructed b-hadron decays. A Gaussian function is used to describe the B ± → J/ψπ ± contribution, with fixed parameters taken from simulation except for the normalisation, which is a free parameter. A fit to the overall mass distribution is used to define the shapes of signal and backgrounds. Subsequent fits are performed in the intervals of the tagging discriminating variables, separately for B + and B − candidate events, with the normalisations and also the slope of the exponential function left free. The B + and B − signal yields are extracted from these fits. Figure 1 shows the invariant mass distribution of B ± candidates over-  Fig. 1 The invariant mass distribution for selected B ± → J/ψ K ± candidates. Data are shown as points, and the overall result of the fit is given by the blue curve. The contributions from the combinatorial background component are indicated by the red dotted line, partially reconstructed b-hadron decays by the purple shaded area, and decays of B ± → J/ψπ ± , where the pion is misassigned as a kaon, by the green dashed line laid with a fit to all selected candidates, and including the individual fit components for the signal and backgrounds.

Flavour tagging methods
The flavour of the signal B meson at the point of production is inferred using several methods, which differ in their efficiency and discrimination power. The measured charge of a lepton (electron or muon) from the semileptonic decay of a B meson provides strong discrimination; however, the ATLAS sensitivity to b → transitions are diluted through processes that can change the charge of the observed lepton, such as through neutral B meson oscillations, or through the cascade decays b → c → . The separation power of lepton tagging is enhanced by considering a weighted sum of the charge of the tracks in a cone around the lepton. If no lepton is present, a weighted sum of the charge of the tracks in a jet associated with the opposite-side b-hadron decay is used to provide discrimination. This weighted sum, or cone charge, is defined as: where x = {μ, e, jet} refers to muon, electron, or jet charge, respectively, and the summation is made using the charge of the track, q i , and its p Ti , over a selected set of tracks, including the lepton, in a cone of size R = ( φ) 2 + ( η) 2 , around the lepton or jet direction. The value of the parameter κ is optimised on each OST method, by determining the value of κ that maximises the tagging power (defined in Sect. 4.3).
The requirements on the tracks and R are described below, dependent on the OST method. Two subcategories of Q x are considered: the first discrete category is used in the case where the cone charge is formed either from only one track or from more than one track of the same charge; this results in a cone charge of Q x = ±1. The second continuous category is used when more than one track is considered, and the sum contains tracks of both negative and positive charge. In the continuous case, Q x is divided into intervals within the range −1 < Q x < 1 for each OST algorithm.
A probability P(B|Q x ) is constructed, which is defined as the probability that a B meson is produced in a state containing ab-quark, given the value of the cone charge Q x . Since Q x is evaluated on the opposite side, a large, negative value of Q x tends to correspond to a higher value of P(B|Q x ). An equivalent probability for the b-quark case is defined as P(B|Q x ). Using the B ± calibration samples, P(Q x |B ± ) for each tagging method used can be defined. The probability to tag a B 0 s meson as containing ab-quark is therefore given as P(B|Q x ) = P(Q x |B + )/(P(Q x |B + ) + P(Q x |B − )), and correspondingly P(B|Q x ) = 1 − P(B|Q x ). If there is no OST information available for a given B 0 s meson, a probability of 0.5 is assigned to that candidate.

Muon tagging
For muon-based tagging, at least one additional muon is required in the event, with p T > 2.5 GeV, |η| < 2.5 and | z| < 5 mm, where | z| is the difference in z between the primary vertex and the longitudinal impact parameter of the ID track associated with the muon. Muons are classified and kept if their identification quality selection working point is either Tight or Lowp T ; these categories are subsequently treated as distinct flavour tagging methods. For muons with p T > 4 GeV, Tight muons are the dominant category, with the Low-p T requirement typically identifying muons of p T < 4 GeV. In the case of multiple muons satisfying selection criteria in one event, Tight muons are chosen over Lowp T muons. Within the same muon category, the muon with the highest p T that passes the selections is used.
A muon cone charge variable, Q μ , is constructed according to Eq. (1), with κ = 1.1 and the sum over the reconstructed ID tracks within a cone of size R = 0.5 around the muon direction. These tracks must have p T > 0.5 GeV, |η| < 2.5, and | z| < 5 mm. Tracks associated with the decay of a B meson signal candidate are excluded from the sum. In each interval of Q μ , a fit to the J/ψ K ± invariant mass spectrum is performed and the number of signal events extracted. The fit model used is described in Sect. 4.1. Figure 2 shows the distributions of the muon cone charge using B ± signal candidates for Tight muons, and includes the tag-ging probability as a function of the cone charge variable. The corresponding distributions for Lowp T muons are shown in Fig. 3.

Electron tagging
Electrons are identified using ID and calorimeter information, and must satisfy the Medium electron quality criteria [34]. The ID track associated with the electron is required to have p T > 0.5 GeV, |η| < 2.5, and | z| < 5 mm. To reject electrons from the signal-side of the decay, electrons with cos(ζ b ) > 0.93, where ζ b is the opening angle between the momentum of the signal B meson candidate and the electron momentum, are not considered. In the case of more than one electron passing the selection, the electron with the highest p T is chosen. Charged particle tracks within a cone of size R = 0.5 are used to form the electron cone charge Q e , constructed according to Eq. (1), with κ = 1.0. The resulting electron cone charge distributions are shown in Fig. 4, together with the corresponding tagging probability.

Jet tagging
In the absence of a muon or electron, a jet identified as containing a b-hadron is required. Jets are reconstructed from calorimetric information [35] using the anti-k t algorithm [36,37] with a radius parameter R = 0.4. The identification of a b-tagged jet uses a multivariate algorithm MV2c10 [38], utilising boosted decision trees (BDT), which output a classifier value. Jets are selected if this value exceeds 0.56. This value is chosen to maximise the tagging power of the calibration sample. In the case of multiple selected jets, the jet with the highest value of the BDT output classifier is used. Jets associated with the signal decay are not considered in this selection.
Tracks within a cone of size R = 0.5 around the jet axis are used to define a jet cone charge, Q jet , constructed according to Eq. (1), where κ = 1.1 and the sum is over the tracks associated with the jet, with | z| < 5 mm, and excluding tracks from the decay of the signal B meson candidate. Figure 5 shows the distribution of the opposite-side jet cone charge for B ± signal candidates.

Flavour tagging performance
In order to quantify and compare the performance of the various tagging methods, three figure-of-merit terms are constructed, which describe: the fraction of events used by a given tagging method, the purity of the method, and the overall power of the tagging method in the sample. The efficiency, x , of an individual tagging method is defined as the number of signal events tagged by that method divided by the total number of signal events in the sample. The purity of a particular flavour tagging method, called the dilution, is defined as D(Q x ) = 2P(B|Q x ) − 1. The tagging power of a particular tagging method is then defined as T x = i x i · D 2 (Q x i ), where the sum is over the probability distribution in intervals of the cone charge variable. An effective dilution, D x = √ T x / x , is calculated from the measured tagging power and efficiency.
By definition, there is no overlap between lepton-tagged and jet-charge-tagged events. The overlap between events with a muon (either Tight or Lowp T ) and events with an electron corresponds to around 0.6% of all tagged events. In the case of multiply tagged events, the OST method is selected in the following order: Tight muon, electron, Lowp T muon, jet. However, the ordering of muon-and electron-tagged events is shown to have negligible impact on the final results. A summary of the tagging performance for each method and the overall performance on the B ± sample is given in Table 1. 4.4 Using tag information in the B 0 s fit For the maximum likelihood fit performed on the B 0 s data, and described in detail in Sect. 5, the per-candidate probability, P(B|Q x ), that the B meson candidate was produced in a state B 0 s (versus aB 0 s ) is provided by the calibrations derived from the B ± → J/ψ K ± sample, described above, and shown in Figs. 2, 3, 4 and 5. Since the distributions of P(B|Q x ) from signal B 0 s mesons and backgrounds can be expected to be different, separate probability density functions (PDFs) are necessary to describe these distributions in the likelihood function. These PDFs are defined as P s (P(B|Q x )) and P b (P(B|Q x )), describing the probability distributions for signal and background, respectively, and are derived from the sample of B 0 s candidates. For the exclusive decays B d → J/ψ K 0 * and b → J/ψ pK − that are present in the sample of B 0 s candidates, P s (P(B|Q x )) is used to model the probability distributions for these contributions (described further in Sect. 5.2). The PDFs consist of the fraction of events that are tagged with a particular method (or are untagged), the fractions of those events categorised as discrete or continuous, and for those that are continuous, a PDF of the corresponding probability distribution.

Continuous PDF
The parameterisations of the continuous PDF components of P s,b (P(B|Q x )) for each OST method are defined as follows. In the sideband regions, 5.150 < m(J/ψ K K ) < 5.317 GeV and 5.417 < m(J/ψ K K ) < 5.650 GeV, unbinned maximum likelihood fits to the P(B|Q x ) distributions are performed to extract the background (continuous category) PDFs for P b (P(B|Q x )). For the Tight muon and electron methods, the parameterisation has the form of the sum of a second-order polynomial and two exponential functions. A   For each plot, in red (blue), the normalised B + (B − ) cone charge distribution is shown (corresponding to the right axis scale). A B + (B − ) candidate is more likely to have a large negative (positive) value of Q μ . Superimposed is the distribution of the tagging probability, P(B|Q μ ), as a function of the cone charge, derived from a data sample of B ± → J/ψ K ± decays, and defined as the probability to have a B + meson (on the signal-side) given a particular cone charge Q μ . The fitted parameterisation, shown in black, is used as the calibration curve to infer the probability to have a B 0 s orB 0 s meson at production in the decays to J/ψφ       Table 1 Summary of tagging performances for the different flavour tagging methods on the sample of B ± signal candidates, as described in the text. Uncertainties shown are statistical only. The efficiency ( x ) and tagging power (T x ) are each determined by summing over the individual bins of the cone charge distribution. The effective dilution (D x ) is obtained from the measured efficiency and tagging power. For the efficiency, effective dilution, and tagging power, the corresponding uncertainty is determined by combining the appropriate uncertainties in the individual bins of each charge distribution Tag method Gaussian function is used for the Lowp T muons. For the jet tagging algorithm an eighth-order polynomial is used. For the signal, fits are performed to the P(B|Q x ) distributions, using all events in the m(J/ψ K K ) distributions to extract the signal (continuous category) PDFs for P s (P(B|Q x )). In these fits, the parameters describing the background PDFs are fixed to their previously extracted values, as is the relative normalisation of signal and background, extracted from a fit to the m(J/ψ K K ) distribution. For the signal PDFs, the Tight muon tagging method uses the sum of two exponential functions and a constant function to describe the signal. For the electron tagging method, the signal function has the form of the sum of a second-order polynomial and two exponential functions, and for the Lowp T muon and jet tagging methods a Gaussian function is used.

Discrete PDF
In the case where the cone charge is discrete, the fractions of events f +1 ( f −1 ) with cone charges +1 (−1) are deter-mined separately for signal and background using events from the signal and sideband regions of the B 0 s mass distribution (as defined in Sect. 3). The remaining fraction of events, 1 − f +1 − f −1 , corresponds to the continuous parts of the distribution. Positive and negative charges are equally probable for background candidates formed from a random combination of a J/ψ and a pair of tracks, but this is not necessarily the case for background candidates formed from a partially reconstructed b-hadron. Table 2 summarises the fractions f +1 and f −1 obtained from each tagging method for signal and background events.
The fractions of signal and background events tagged using the different OST methods are found using a similar sideband-subtraction method, and are summarised in Table 3.
To account for possible deviations of the data from the selected fit models, variations of the procedure described here are used to determine systematic uncertainties, as described in Sect. 6.

Maximum likelihood fit
An unbinned maximum likelihood fit is performed on the selected events to extract the parameter values of the The fit uses information about the reconstructed mass, m, the measured proper decay time, t, the measured mass uncertainty, σ m , the measured proper decay time uncertainty, σ t , the measured transverse momentum, p T , the tagging probability, P(B|Q x ), and the transversity angles, , of each B 0 s → J/ψφ decay candidate. The measured value of the proper decay time uncertainty, σ t , is calculated from the covariance matrix associated with the vertex fit for each candidate event. The transversity angles = (θ T , ψ T , φ T ) are defined in Sect. 5.1. The like- Table 2 Fractions f +1 and f −1 of events with cone charges of +1 and −1, respectively, for signal and background events and for the different tagging methods. Only statistical uncertainties are given Tight muon 6.9 ± 0.3 7 .5 ± 0.3 4 .7 ± 0.1 4 .9 ± 0.1 lihood function is defined as a combination of the signal and background PDFs as follows: where N is the number of selected candidates, w i is a weighting factor to account for the trigger efficiency (described in Sect. 5.3). The terms F s , F B 0 , F b and F bkg are the PDFs modelling the signal, B 0 background, b background, and the other background distributions, respectively. The term f s is the fraction of signal candidates and f B 0 and f b are the background fractions of B 0 mesons and b baryons misidentified as B 0 s candidates, calculated relative to the number of signal events. These background fractions are fixed to their expectation from the MC simulation, and variations are applied as part of the evaluation of the effects of systematic uncertainties. The mass m i , the proper decay time t i and the decay angles i are the values measured from the data for each event i. A detailed description of the signal PDF terms in Eq.
(2) is given in Sect. 5.1. The three background functions are described in Sect. 5.2.

Signal PDF
The PDF used to describe the signal events, F s , has the following composition: The mass term P s (m i |σ m i ) is modelled in the following way: The term P s (m i |σ m i ) uses per-candidate mass errors, σ m i , calculated for each J/ψφ candidate from the covariance matrix associated with the four-track vertex fit. Each measured candidate mass is convolved with a Gaussian function with a width equal to σ m i multiplied by a scale factor S m , introduced to account for any mismeasurements. Both S m and the mean value m B s , which is the B 0 s meson mass, are free parameters determined in the fit.
The PDF term P s (t i , i |σ t i , P i (B|Q x )) takes into account the lifetime resolution, so each time element in Table 4 is convolved with a Gaussian function defined as: S t is a scale factor (a parameter of the fit) and σ t i is the per-candidate uncertainty on proper decay time t i . This convolution is performed numerically on an event-by-event basis and the value σ t i is measured for each B 0 s candidate, based on the tracking error matrix of the four final state particles. The probability term P s (σ t i |p Ti ) is introduced to account for differences between signal and background events for the values of the per-candidate time errors. Distributions of this variable for signal and background described by gamma functions are shown in Fig. 6. The average value of the time error for signal events is 69 fs.  (t) and the functions of the transversity angles g (k) (θ T , ψ T , φ T ). The amplitudes |A 0 (0)| 2 and |A (0)| 2 are for the CP-even components of the B 0 s → J/ψφ decay, |A ⊥ (0)| 2 is the CP-odd amplitude; they have corresponding strong phases δ 0 , δ and δ ⊥ . By convention, δ 0 is set to be zero. The S-wave amplitude |A S (0)| 2 gives the fraction of B 0 s → J/ψ K + K − ( f 0 ) and has a related strong phase δ S . The factor α is described in the text of Sect. 5.1. The ± and ∓ terms denote two cases: the upper sign describes the decay of a meson that was initially a B 0 s meson, while the lower sign describes the decays of a meson that was initiallyB 0  Fig. 6 The proper decay time uncertainty distribution for data (black), and the fits to the background (blue) and the signal (purple) contributions. The total fit is shown as a red curve The same approach was applied for the probability terms P s (σ m i |p T i ) and P s (p Ti ) accounting for differences between signal and background events for the values of the percandidate mass error and p T i values, respectively. The tagging probability term for signal P s (P i (B|Q x )) is described in Sect. 4.4.
The term P s (t i , i |σ t i , P i (B|Q x )) is a joint PDF for the decay time t and the transversity angles for the B 0 s → J/ψ(μ + μ − )φ(K + K − ) decay. Ignoring detector effects, the distribution for the time t and the angles is given by the differential decay rate [39]: where O (k) (t) are the time-dependent functions corresponding to the contributions of the four different amplitudes (A 0 , A || , A ⊥ , and A S ) and their interference terms, and g (k) (θ T , ψ T , φ T ) are the angular functions. Table 4 shows the time-dependent and the angular functions of the transversity angles. The formulae for the time-dependent functions have the same structure for B 0 s andB 0 s but with a sign reversal in the terms containing m s , which is a fixed parameter of the fit (using Ref. [33]). The formalism used throughout this analysis assumes no direct CP violation, i.e. λ = (q/ p)(Ā/A) = 1. In Table 4, the parameter A ⊥ (t) is the time-dependent amplitude for the CP-odd final-state configuration while A 0 (t) and A (t) correspond to CP-even final-state configurations. The amplitude A S (t) gives the contribution from the CP-odd non-resonant B 0 s → J/ψ K + K − S-wave state [40] (which includes the f 0 meson). The corresponding functions are given in the last four lines of Table 4 (k = 7-10). The amplitudes are parameterised by |A i |e iδ i , where i = {0, ||, ⊥, S}, δ 0 = 0 and are normalised such that |A 0 (0)| 2 + |A ⊥ (0)| 2 + |A (0)| 2 = 1. The amplitude |A ⊥ (0)| is determined according to this condition, while the remaining three amplitudes are parameters of the fit. The value |A S | 2 is the ratio of the S-wave yield to the φ → K + K − yield in the interval of m(K + K − ) used in the analysis. In the sum over the mass interval, the interference terms (lines 8-10 in Table 4) are corrected by a factor α that takes into account the mass-dependent differences in absolute amplitude and phase between the φ → K + K − and the S-wave amplitudes. The correction is based on the Breit-Wigner description of the φ and on model assumptions for the shape and the phase variations of the S-wave amplitude. The phase δ S is the phase difference between A 0 (0) and the S-wave amplitudes at the φ → K + K − peak. The values of α and the related systematic uncertainty are discussed in Sect. 6.
The angles (θ T , ψ T , φ T ), are defined in the rest frames of the final-state particles. The x-axis is determined by the direction of the φ meson in the J/ψ rest frame, and the K + K − system defines the x-y plane, where p y (K + ) > 0. The three angles are defined as: • θ T , the angle between p(μ + ) and the normal to the x-y plane, in the J/ψ meson rest frame, • φ T , the angle between the x-axis and p xy (μ + ), the projection of the μ + momentum in the x-y plane, in the J/ψ meson rest frame, • ψ T , the angle between p(K + ) and − p(J/ψ) in the φ meson rest frame.
The angular acceptance of the detector and the kinematic cuts on the angular distributions are included in the likelihood function through A( i , p Ti ). This is calculated using a four-dimensional binned acceptance method, applying an event-by-event efficiency according to the transversity angles (θ T , ψ T , φ T ) and the p T of the candidate. The p T binning is necessary, because the angular acceptance is influenced by the p T of the B 0 s candidate. The acceptance is calculated from the B 0 s → J/ψφ MC events with additional weighting for p T and η distributions. In the likelihood function, the acceptance is treated as an angular acceptance PDF, which is multiplied with the time-and angle-dependent PDF describing the B 0 s → J/ψ(μ + μ − )φ(K + K − ) decays. As both the acceptance and time-and angle-dependent decay PDFs depend on the transversity angles they must be normalised together. This normalisation is done numerically during the likelihood fit. The PDF is normalised over the entire B 0 s mass range, 5.150-5.650 GeV.

Background PDF
The background PDF has the following composition: The proper decay time function P b (t i |σ t i ) is parameterised as a peak modelled by a Gaussian distribution, two positivetime exponential functions and two negative-time exponen-tial functions. These functions are convolved with the same resolution function, defined in Eq. (4) as the signal decay time-dependence. The prompt peak models the combinatorial background events, which are expected to have reconstructed lifetimes distributed around zero. The two positivetime exponential functions represent a fraction of longerlived backgrounds with non-prompt J/ψ, combined with hadrons from the primary vertex or from a B/D meson in the same event. The two negative-time exponential functions take into account events with poor vertex resolution. The probability terms P b (σ m i |p T i ), P b (σ t i |p T i ) and P b (p T i ) are described by gamma functions. They are unchanged from the analysis described in Ref.
[41] and explained in detail there. The tagging probability term for background events The shape of the background angular distribution, P b ( i ) arises primarily from detector and kinematic acceptance effects. The best description is achieved by Legendre polynomial functions: where l max = 14, k max = 14 and the coefficients a k,l,m are adjusted to give the best fit to the angular distributions for events in the sidebands of the B 0 s mass distribution. These parameters are then fixed in the main fit, defined by Eq. (2). The B 0 s mass interval used for the background fit is between 5.150 and 5.650 GeV excluding the signal mass region |(m(B 0 s )−5.366| < 0.110 GeV. Higher-order Legendre polynomial functions were tested as a systematic check, described in Sect. 6. The background mass model, P b (m i ) is the sum of an exponential and a constant, with the exponential slope and the relative normalisation left as free parameters of the fit.
Contamination from B d → J/ψ K 0 * and b → J/ψ pK − events misreconstructed as B 0 s → J/ψφ is accounted for in the fit through the F B 0 and F b terms in the PDF described in Eq. (2). The PDFs are determined using MC simulation of these decays. The fractions of these contributions, f B 0 = (4.3 ± 0.5)% and f b = (2.1 ± 0.6)%, are defined relative to the number of the B 0 s → J/ψφ signal events and are evaluated from MC simulation using production cross sections and branching fractions from Refs. [33,[42][43][44][45][46]. MC simulated events are also used to determine the shape of the mass and transversity angle distributions. The 3D angular distributions of B 0 d → J/ψ K 0 * and of the conjugate decay are modelled using input from Ref.
[47], while angular distributions for b → J/ψ pK − and the conjugate decay are considered flat. These distributions are sculpted for detector acceptance effects and then described by Legendre polynomial functions with l max = 10 and k max = 10, Eq. (5). These shapes are used as templates in the fit. The B d and b lifetimes are accounted for in the fit by adding additional exponential terms, scaled by the ratio of B d /B 0 s or b /B 0 s masses as appropriate, where the lifetimes and masses are taken from Ref. [33]. The PDF terms that describe each of the tagging, mass, decay time and p T probability distributions are taken from the same PDFs used to describe the B 0 s → J/ψφ signal (described in Sects. 4.4 and 5.1) to account for the fact that these dedicated background events are fully reconstructed b-hadron decays. Systematic uncertainties due to the background from B d → J/ψ K 0 * and b → J/ψ pK − decays are described in Sect. 6. The contribution of the S-wave B d → J/ψ K π decays as well as their interference with the P-wave B d → J/ψ K 0 * decays are included in the PDF of the fit, using the parameters measured in Ref. [47].

Proper decay time dependence of muon trigger efficiency
In the triggers used in this analysis, there is no minimum cut applied on the transverse impact parameter d 0 of muons. On the other side, trigger muons with values of d 0 > 10 mm are not accepted. This results in inefficiency at large values of the proper decay time. This inefficiency is estimated using MC simulated events, by comparing the B 0 s proper decay time distribution obtained before and after applying the trigger selection. To account for this inefficiency in the fit, the events are reweighted by a factor w, inversely proportional to the trigger efficiency: where Erf denotes the error function and p 0 , p 1 , p 2 and p 3 are parameters determined in the fit to MC events. For more than 99% of the B 0 s candidates the inefficiency is below 2%. For the most affected candidates the inefficiency reaches up to 50% at high decay time. No significant bias or inefficiency due to offline track reconstruction, vertex reconstruction, or track quality selection criteria is observed.

Summary of the fit parameters
The joint PDF of proper decay time and decay angles includes the main physics parameters of interest: • the CP-violating phase φ s , • the average decay width s and the decay width difference s , • the size of the CP-state amplitudes at t = 0: |A (0)| 2 , |A 0 (0)| 2 and their corresponding strong phases δ ⊥ and δ • and the size of the S-wave amplitude at t = 0: |A S (0)| 2 and corresponding strong phase δ S .
The size of the remaining amplitude |A ⊥ (0)| 2 is constrained by the normalisation condition, phase δ 0 is set to zero and m s is fixed as mentioned above.
The likelihood function also includes other parameters referred to as "nuisance parameters" such as: the B 0 s signal fraction f s , parameters describing the invariant mass and decay time-angular distributions of combinatorial background events and scale factors of mass and decay time uncertainties. In addition, there are also other nuisance parameters describing: acceptance functions, parametrisations of the angles of dedicated backgrounds B d → J/ψ K 0 * and b → J/ψ pK − and their fractions f B 0 and f b , the probability density functions of time error distributions P(σ t i | p T i ), mass error distributions P(σ m i | p T i ), p T distributions P( p T i ) and tagging parameters and calibrations. These parameter values are mainly fixed in the fit to the values extracted from the B 0 s mass signal and sideband regions or from MC simulations.

Systematic uncertainties
Systematic uncertainties are evaluated for effects that are described below.
• Flavour tagging: The effects on the main physics parameters from the fit, due to uncertainties introduced by the flavour tagging procedure, are assessed as follows: The statistical uncertainty due to the size of the sample of B ± → J/ψ K ± decays is included in the overall statistical uncertainty. The systematic uncertainty arising from the precision of the OST calibration, described in Sect. 4.2, is estimated by changing the models used to parameterise the probability distribution, P(B|Q x ), as a function of the cone charge from the function used by default (a third-order polynomial for muons and a sinusoid for electrons) to one of several alternative functions: a linear function; a fifth-order polynomial; or two thirdorder polynomials that describe the positive and negative regions and have common constant and linear terms, but independent quadratic and cubic terms. The B 0 s fit is repeated using the alternative models and the largest deviation from the nominal fit is assigned as the systematic uncertainty. To validate the calibration procedure, calibration curves are derived from simulated samples of B ± and B 0 s signals. The variations between the curves from these two samples are propagated to the calibration curves derived from data. The differences in the parameter values between the nominal fit and that with the varied calibration curves are included in the systematic uncertainty. An additional systematic uncertainty is assigned to account for potential dependencies on the pile-up distribution. The calibration data are split into subsets of approximately equal size, separated according to the estimated pile-up of the event, and separate calibrations are made for each subset. For the B 0 s fit, the fit is repeated using the calibrations corresponding to the estimated pile-up of that event. Differences between the nominal and the modified fit for the parameters of interest are taken as the systematic uncertainty. For the terms P b (P(B|Q x )) and P s (P(B|Q x )), variations of the parameterisation are considered (including using histograms in place of a parameterisation). The resulting changes in the parameter values of the B 0 s fit are similarly included in the systematic uncertainties.
• ID alignment: The changes of the fit parameters due to residual misalignments of the ID were studied and the observed deviations are included in the systematic uncertainties. • Angular acceptance method: The angular acceptance of the detector and the kinematic cuts, A( i , p T i ), described in Sect. 5.1, is calculated from a binned fit to MC simulated data. In order to estimate the systematic uncertainty introduced by the choice of binning, different acceptance functions are calculated using different numbers of p T bins as well as different widths and central values of the bins. • Time efficiency: To correct for the proper decay time dependence of trigger inefficiencies, the events are reweighted according to Eq. (5). To estimate systematic uncertainties connected with this procedure, several alternative fits are performed. Firstly, using different sets of p T binning in the MC sample used to determine the efficiency. Secondly, to assess the effects of mis-modelling of the B 0 s vertex χ 2 /ndof in simulated data, an alternative fit is done with MC χ 2 /ndof reweighted according to data. Finally, the groups of similar triggers used to determine the efficiencies were subdivided further by their common features. These smaller groups allowed to address even more details of trigger varieties in time-efficiency determination. Deviations from the default fit result are included in the systematic uncertainties of the measurement.
• Best candidate selection: After applying all selection cuts for B 0 s signal events, approximately 5% of the events are found to contain multiple candidates. In the default fit, the B 0 s candidate with the lowest χ 2 /ndof is selected. To assess the systematic uncertainty due to this selection, an equivalent sample is created where all candidates in the event are retained, each weighted by a factor of 1/N cand , where N cand is number of B 0 s candidates in the event. Deviations from the default fit are included in the systematic uncertainties of the measurement. • Background angles model: The shape of the background angular distribution, P b (θ T , ϕ T , ψ T ), is described by the fourteenth-order Legendre polynomial functions, given in Eq. (5). Alternatively, higher-order Legendre polynomial functions with l max = 16 and k max = 16 were tested, and the changes in the fit parameter values relative to the default fit are taken as systematic uncertainties.
The shapes are primarily determined by detector and kinematic acceptance effects and are sensitive to the p T of the B 0 s meson candidate. For this reason, the parameterisation using the Legendre polynomial functions is performed in six p T intervals: [10][11][12][13][14][15][15][16][17][18][19][20][20][21][22][23][24][25][25][26][27][28][29][30][30][31][32][33][34][35] GeV and >35 GeV. The systematic uncertainties due to the choice of p T intervals are estimated by repeating the fit, with these intervals enlarged and reduced by 1 GeV and by 2 GeV. The largest changes in the fit results are taken to represent the systematic uncertainties. The sensitivity of the fit to the choice of the invariant mass window is tested by varying its size. The difference to the default fit results is assigned as a systematic uncertainty. The parameters of the Legendre polynomial functions given in Eq. (5)  → J/ψφ is accounted for in the final fit. Studies are performed to evaluate the systematic uncertainties due to trigger effects, the B d → J/ψ K 0 * fraction, and the distributions of the mass, transversity angles, and lifetime PDFs. In the MC events the angular distribution of the B d → J/ψ K 0 * decay is modelled using parameters taken from Ref. [47]. The contribution of the S-wave B d → J/ψ K π decays as well as its interference with the P-wave B d → J/ψ K 0 * decays are also included in the PDF of the fit, following the parameters measured in Ref. [47]. The uncertainties of these parameters are taken into account in the estima-tion of the systematic uncertainty. After applying the B 0 s signal selection cuts, the angular distributions are fitted using Legendre polynomial functions. The uncertainties of this fit are included in the systematic uncertainty. • b contribution: The contamination from b → J/ψ pK − events misreconstructed as B 0 s → J/ψφ is accounted for in the final fit. Studies are performed to evaluate the effect of the uncertainties in the b → J/ψ pK − fraction f b , and the shapes of the distributions of the mass, transversity angles, and lifetime. Additional studies are performed to determine the effect of the uncertainties in the b → J/ψ * branching ratios used to reweight the generated MC sample.
• Alternate m s : The systematics due to fixing the parameter m s to the PDG value were estimated by running an alternative fit where the default model was altered by releasing m s within a Gaussian constraint of which the width was taken from uncertainties assigned to m s in [33]. The resulting changes to all the fit parameter values are found to be negligible except for δ ⊥ . The effects on the other variables are small for the parameters φ s and δ . • Fit model mass and lifetime: To estimate the systematic uncertainties due to the signal B 0 s mass model, the default model was altered by adding a second Gaussian function in Eq. (3), which has the same structure as the first Gaussian function but a different scale factor, S 1 m , which is an additional free parameter of the fit. The resulting changes to other fit parameter values are found to be negligible. To test the sensitivity of the part of the fit model describing the lifetime, two systematic tests are performed. The determination of signal and background lifetime errors is sensitive to the choice of p T bins, in which the relative contributions of these two components are evaluated. To estimate the systematic uncertainty, the fit is repeated varying the intervals of the default p T binning. The determination of signal and background lifetime errors is also sensitive to the determination of the signal fraction. The fit is repeated by varying this fraction within one standard deviation of its uncertainty and differences are included in the systematic uncertainty.
• Fit model S-wave phase: As explained in Sect. 5.1, the model for the interference between B 0 s → J/ψφ(K + K − ) and S-wave B 0 s → J/ψ K + K − is corrected by a factor α to account for the mass-dependent variations in absolute amplitude and phase of the two amplitudes within the interval 1.0085-1.0305 GeV in m(K + K − ). The value of α is 0.51± 0.08. The central value is obtained under the hypothesis of uniformity of the S-wave amplitude. The uncertainty is systematic and it is due to: detector mass resolution and mass scale uncertainties, uncertainties in the description of the φ resonance (mass, width, and shape descriptions as relativistic Breit-Wigner or Flatté parameterisation [48]), and to uncertainties in the description of the shape and phase variation of the Swave amplitude. For the last effect, which is dominant, the S-wave amplitude is assumed to be due to the f 0 (980) resonance and it is described as in Ref. [49], similarly to what was done in Ref. [16]. The variation from the value of α obtained with the default assumption is taken as a systematic uncertainty. This procedure uses descriptions of the f 0 (980) based on other measurements [33]. To account for the uncertainty in α, the fit was repeated with α = 0.51 + 0.08 and α = 0.51 − 0.08 values. The variations of the parameter values relative to those from the default fit using the central value of α are included in the systematic uncertainties.
• Fit bias: Due to its complexity, the fit model can be sensitive to some nuisance parameters. This limited sensitivity could potentially lead to a bias in the measured physics parameters, even when the model describes the fitted data well. To test the stability of the results obtained from the chosen default fit model, a set of pseudo-experiments is conducted using the default model in both the generation and fit. The systematic uncertainties are determined from the mean of the pull distributions of the pseudoexperiments scaled according to the statistical uncertainty of that parameter in the fit to data. The observed deviations are included in the systematic uncertainties.
The systematic uncertainties are listed in Table 5. For each parameter, the total systematic uncertainty is obtained by adding all of the contributions in quadrature.

Fit results
The results of the likelihood fit are shown in Table 6. The total number of B 0 s meson candidates is 453 570 ± 740. The fitted value of the B 0 s mass agrees well with the world average value [33]. Fit projections, including ratio plots, are shown in Fig. 7 for the mass and proper decay time and in Fig. 8 for the angles. The ratio plots show the difference between each data point and the total fit line divided by the statistical and systematic uncertainties summed in quadrature (σ ) for that point. The deviations of ratio plots are within 2σ , which shows that the total uncertainties cover any discrepancy between data and fit model.
While for most of the physics parameters, including φ s , s and s , the fit determines a single solution with Gaussian behaviour of the projection of the log-likelihood (see Fig. 10 in Sect. 7.2), for the strong-phases δ and δ ⊥ two well separated local maxima of the likelihood are found, and shown as solution (a) and (b) in Table 6. The difference in 10 −3 rad) (10 −3 ps −1 ) (10 −3 ps −1 ) (10 −3 ) (10 −3 ) (10 −3 ) (10 −3 rad) (10 −3 rad) (10 −3 rad) −2 ln(L) between the two solutions is 0.03. As discussed in Sect. 7.2, the two-fold behaviour of the likelihood in the strong phases is the result of an approximate symmetry of the signal PDF. The effect is completely negligible for all other variables, for which the fit values and uncertainty ranges overlap accurately. The correlation between statistical uncertainties of the parameters are also shown in Table 7 for solution (a) and in Table 8 for solution (b). The correlations do not change significantly between the physics variables which remain stable between the solutions (a) and (b). The correlations with the two strong phases: δ and δ ⊥ flip the sign, as expected, except for the correlations that are smaller than 0.01. The correlations between δ and δ ⊥ themselves keep the same sign in both solutions (a) and (b).
Releasing the direct CP related λ parameter (see Sect. 5.1) in the fit results in a λ value compatible within the statistical precision with unity as well as with results of other experiments. The effect on the values of the other parameters of interest is negligible, due to small correlations with λ.

Fit to strong phases
As shown in Table 6, the likelihood fit has determined two solutions with well separated values for the strong phases δ || and δ ⊥ . Figure 9 shows results of the 2D log-likelihood scan in the δ || , δ ⊥ plane, revealing two minima, identified at (δ || = 3.35, δ ⊥ = 3.12) and (δ || = 2.94, δ ⊥ = 2.91). These minima are represented by two-dimensional contours at the level of − 2 ln(L) = 2.30, 6.18, and 11.83, where − 2 ln(L) = 2(ln(L i ) − ln(L a )) is the difference between the likelihood values (L i ) of the fit in which the two strong phases are fixed to the values shown on the horizontal and vertical axis, and L a which is the likelihood value for solution (a) of the fit.
An approximate symmetry in the signal PDF is at the origin of this duality. The strong phases are determined by the six interference terms 4) − 6) and 8) − 10) in Table 4. Four of them, namely terms 4), 5), 8) and 9), are invariant under the transformation This transformation is proportional to π −δ , which from our data is equal to 0.21. The two local maxima of the likelihood determined in this analysis satisfy the relation of Eq. (6) very accurately for δ || and (δ ⊥ -δ S ), and within the statistical uncertainty in the transformation for δ ⊥ . The difference in the log-likelihoods, − 2 ln(L), between the two solutions is equal to 0.03, favouring (a) but without ruling out (b).
As discussed in Sect. 7.1, the two-fold nature of the likelihood maxima for the strong phases has a negligible effect on Table 6 Fitted values for the physical parameters of interest with their statistical and systematic uncertainties. For variables δ ⊥ and δ the values are given for the two solutions (a) and (b). The difference in -2 ln(L) between solutions (b) and (a) is 0.03. For the rest of the variables the values for the two minima are consistent. The same is true for the statistical and systematic uncertainties of all the variables   Figure 10 shows the one-dimensional likelihood scans for all other physics parameters, separately for the two solutions (a) and (b). For each scan, the other physics parameters and all nuisance parameters are optimised in a profile-likelihood fit. The 1D scans are almost identical for the two solutions.

Combination with 7 TeV and TeV results
The measured values of solution (a) and solution (b) are consistent with those obtained in the previous analysis [13] using 19.2 fb −1 of data collected at √ s = 7 TeV and 8 TeV. A best linear unbiased estimator (BLUE) [50,51] is used to combine the current results with those from the previous analysis. The measured values, uncertainties, and correlations are all background components. Below each figure is a ratio plot that shows the difference between each data point and the total fit line divided by the statistical and systematic uncertainties summed in quadrature (σ ) of that point taken from the measurements performed at each centre-ofmass energy. The statistical correlation between these three measurements is zero as the events are different. The correlations of the systematic uncertainties between the three measurements are estimated and tested in several categories depending on whether the given systematic effect changed significantly between the measurements. Solution (a) and solution (b) are treated separately, leading to the two sets of combined results as shown in Table 9. The correlation matrices of these two combinations are shown in Tables 10 and 11.
The two-dimensional likelihood contours in the φ ss plane for the ATLAS result based on 7 and 8 TeV data, the solution (a) of the 13 TeV measurement, and the combined result for the solution (a) are shown in Fig. 11. The statistical and systematic uncertainties are combined in quadrature and correlations are taken into account in the construction of Gaussian contours. Because there is no significant difference in the φ s and s values between solution (a) and solution (b), only contours for solution (a) are shown.
Two-dimensional likelihood contours in the φ ss plane are shown in Fig. 12

Summary
This paper presents a measurement of the time-dependent CP asymmetry parameters in B 0 s → J/ψφ decays from an 80.5 fb −1 data sample of pp collisions collected with the ATLAS detector during the 13 TeV LHC run. The values from the 13 TeV analysis are consistent with those obtained in the previous ATLAS analysis using 7 TeV and 8 TeV data. The two measurements are statistically combined.
The CP-violating phase φ s is measured to be −0.087 ± 0.036 (stat.) ± 0.021 (syst.) rad, the decay width difference between heavy and light B 0 s mass eigenstates, s = 0.0657 ± 0.0043 (stat.) ± 0.0037 (syst.) ps −1 and the measurement for their average decay width, s = 0.6703 ± 0.0014 (stat.) ± 0.0018 (syst.) ps −1 . The measurement of the CP-violating phase φ s is consistent with the Standard Model prediction, and it improves on the precision of previous ATLAS measurements, while for the average decay width, s , the comparison of the 13 TeV analysis with the current world combined value reveals a tension at the level of 3σ .

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes. All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN. Data associated with journal publications are also made available: tables and data from plots (e.g. cross section values, likelihood profiles, selection efficiencies, cross section limits, ...) are stored in appropriate repositories such as HEPDATA (http:// hepdata.cedar.ac.uk/). ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models. For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.org/). This information is taken from the ATLAS Data Access Policy, which is a public document that can be downloaded from http://opendata.cern.ch/record/413[opendata. cern.ch].] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .