Measurement of CP violation in the $B_s^0\rightarrow\phi\phi$ decay and search for the $B^0\rightarrow \phi\phi$ decay

A measurement of the time-dependent CP-violating asymmetry in $B_s^0\rightarrow\phi\phi$ decays is presented. Using a sample of proton-proton collision data corresponding to an integrated luminosity of $5.0$fb$^{-1}$ collected by the $\mbox{LHCb}$ experiment at centre-of-mass energies $\sqrt{s} = 7$ TeV in 2011, 8 TeV in 2012 and 13 TeV in 2015 and 2016, a signal yield of around 9000 $B_s^0\rightarrow\phi\phi$ decays is obtained. The CP-violating phase $\phi_s^{s\bar{s}s}$ is measured to be $-0.073 \pm 0.115$(stat)$\pm 0.027$(syst) rad, under the assumption it is independent on the helicity of the $\phi\phi$ decay. In addition, the CP-violating phases of the transverse polarisations under the assumption of CP conservation of the longitudinal phase are measured. The helicity-independent direct CP-violation parameter is also measured, and is found to be $|\lambda|=0.99 \pm 0.05 $(stat)$ \pm 0.01 $(syst). In addition, $T$-odd triple-product asymmetries are measured. The results obtained are consistent with the hypothesis of CP conservation in $b\rightarrow\bar{s}s\bar{s}$ transitions. Finally, a limit on the branching fraction of the $B^0\rightarrow \phi\phi$ decay is determined to be $\mathcal{B}(B^0\rightarrow \phi\phi)<2.7\times 10^{-8}$ at 90% confidence level.


Introduction
In the Standard Model (SM) the B 0 s → φφ decay, where the φ(1020) is implied throughout this paper, is forbidden at tree level and proceeds predominantly via a gluonic b → sss loop (penguin) process. Hence, this channel provides an excellent probe of new heavy particles entering the penguin quantum loops [1][2][3]. In the SM, CP violation is governed by a single phase in the Cabibbo-Kobayashi-Maskawa quark mixing matrix [4]. Interference caused by the resulting weak phase difference between the B 0 s -B 0 s oscillation and decay amplitudes leads to a CP asymmetry in the decay-time distributions of B 0 s and B 0 s mesons. For B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decays, which proceed via b → scc transitions, the SM prediction of the weak phase is −2 arg (−V ts V * tb /V cs V * cb ) = −0.0369 +0.0010 −0.0007 rad according to the CKMfitter group [5], and −2 arg (−V ts V * tb /V cs V * cb ) = −0.0370 ± 0.0010 rad according to the UTfit collaboration [6]. The LHCb collaboration has measured the weak phase in several decay processes: B 0 s → J/ψ K + K − , B 0 s → J/ψ π + π − , B 0 s → J/ψ K + K − for the K + K − invariant mass region above 1.05 GeV/c, B 0 s → ψ(2S)φ and B 0 s → D + s D − s , corresponding to the combined result of −0.041 ± 0.025 rad [7]. These measurements are consistent with the SM prediction and place stringent constraints on CP violation in B 0 s -B 0 s oscillations [8]. The CP -violating phase, φ sss s , in the B 0 s → φφ decay is expected to be small in the SM. Calculations using quantum chromodynamics factorisation (QCDf) provide an upper limit of 0.02 rad for its absolute value [1][2][3]. The previous most accurate measurement is φ sss s = −0.17 ± 0.15 (stat) ± 0.03 (syst) rad [9]. CP violation can also be probed by time-integrated triple-product asymmetries. These are formed from T -odd combinations of the momenta of the final-state particles. These asymmetries complement the decay-time-dependent measurement [10] and are expected to be close to zero in the SM [11]. Previous measurements of the triple-product asymmetries in B 0 s decays from the LHCb and CDF collaborations [9,12] have shown no significant deviations from zero.
The B 0 s → φφ decay is a P → V V decay, where P denotes a pseudoscalar and V a vector meson. This gives rise to longitudinal and transverse polarisation of the final states with respect to their direction of flight in the B 0 s reference frame, the fractions of which are denoted by f L and f T , respectively. In the heavy quark limit, f L is expected to be close to unity at tree level due to the vector-axial structure of charged weak currents [2]. This is found to be the case for tree-level B decays measured at the B Factories [13][14][15][16][17][18]. However, the dynamics of penguin transitions are more complicated. Previously LHCb reported a value of f L ≡ |A 0 | 2 = 0.364 ± 0.012 in B 0 s → φφ decays [9]. The measurement is in agreement with predictions from QCD factorisation [2,3]. The observed value of f L is significantly larger than that seen in the B 0 s → K * 0 K * 0 decay [19,20]. In addition to the study of the B 0 s → φφ decay, a search for the as yet unobserved decay B 0 → φφ is made. In the SM this is an OZI suppressed decay [21], with an expected branching fraction in the range (0.1 − 3.0) × 10 −8 [1,2,22,23]. However, the branching fraction can be enhanced, up to the 10 −7 level, in extensions to the SM such as supersymmetry with R-parity violation [23]. The most recent experimental limit was determined to be 2.8 × 10 −8 at 90 % confidence level [24].
Measurements presented in this paper are based on pp collision data corresponding to an integrated luminosity of 5.0 fb −1 , collected with the LHCb experiment at centre-of-mass energies √ s = 7 TeV in 2011, 8 TeV in 2012, and 13 TeV from 2015 to 2016. This paper reports a time-dependent analysis of B 0 s → φφ decays, where the φ meson is reconstructed in the K + K − final state, that measures the CP -violating phase, φ sss s , and the parameter |λ|, that is related to the direct CP violation. Results on helicity-dependent weak phases are also presented, along with helicity amplitudes describing the P → V V transition and strong phases of the amplitudes. In addition, triple-product asymmetries for this decay are presented. The analysis also includes a search for the decay B 0 → φφ. Results presented here supersede the previous measurements based on data collected in 2011 and 2012 [9].

Detector description
The LHCb detector [25,26] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region [27], a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [28] placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [29]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [25].
The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware trigger stage, events are required to contain a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. In the software trigger, B 0 s → φφ candidates are selected either by identifying events containing a pair of oppositely charged kaons with an invariant mass within 30 MeV/c 2 of the known φ meson mass, m φ = 1019.5 MeV/c 2 [30], or by using a topological b-hadron trigger. This topological trigger requires a three-track secondary vertex with a large sum of the p T of the charged particles and significant displacement from the PV. At least one charged particle should have p T > 1.7 GeV/c and χ 2 IP with respect to any primary vertex greater than 16, where χ 2 IP is defined as the difference in χ 2 of a given PV fitted with and without the considered track. A multivariate algorithm [31] is used for the identification of secondary vertices consistent with the decay of a b hadron.
Simulation samples are used to optimise the signal candidate selection, to derive the angular acceptance and the correction to the decay-time acceptance. In the simulation, pp collisions are generated using Pythia [32] with a specific LHCb configuration [33]. Decays of hadronic particles are described by EvtGen [34], in which final-state radiation is generated using Photos [35]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [36], as described in Ref. [33].

Selection and mass model
For decay-time-dependent measurements and the T -odd asymmetries presented in this paper, the previously analysed data collected in 2011 and 2012 [9] is supplemented with the additional data taken in 2015 and 2016, to which the selection described below is applied. For the case of the B 0 → φφ search, a wider invariant-mass window is required, along with more stringent background rejection requirements.
Events passing the trigger are required to satisfy loose criteria on the fit quality of the four-kaon vertex, the χ 2 IP of each track, the transverse momentum of each particle, and the product of the transverse momenta of the two φ candidates. In addition, the reconstructed mass of the φ candidates is required to be within 25 MeV/c 2 of the known φ mass [30].
In order to separate further the B 0 s → φφ signal candidates from the background, a multilayer perceptron (MLP) [37] is used. To train the MLP, simulated B 0 s → φφ candidates satisfying the same requirements as the data candidates are used as a proxy for signal, whereas the four-kaon invariant-mass sidebands from data are used as a proxy for background. The invariant-mass sidebands are defined to be inside the region is the four-kaon invariant mass. Separate MLP classifiers are trained for each data taking period. The variables used in the MLP comprise the minimum and the maximum p T and η of the kaon and φ candidates, the p T and η of the B 0 s candidate, the quality of the four-kaon vertex fit, and the cosine of the angle between the momentum of the B 0 s and the direction of flight from the PV to the B 0 s decay vertex, where the PV is chosen as that with the smallest impact parameter χ 2 with respect to the B 0 s candidate. For measurements of CP violation, the requirement on each MLP is chosen to maximise N S / √ N S + N B , where N S (N B ) represents the expected signal and background yields in the signal region, defined as m B 0 s ± 3σ, where m B 0 s is the known B 0 s mass [30]. The signal yield is estimated using simulation, whereas the number of background candidates is estimated from the data sidebands. For the search of the B 0 → φφ decay, the figure of merit is chosen to maximise ε/(a/2 + √ N B ) [38], where a = 3 corresponds to the desired significance, and ε is the signal efficiency, determined from simulation. This figure of merit does not depend on the unknown B 0 → φφ decay rate.
The presence of peaking backgrounds is studied using simulation. The decay modes considered include B 0 → φK * 0 , Λ 0 b → φpK − , B 0 → φπ + π − and B + → φK + , where the last decay mode could contribute if an extra kaon track is added. The B 0 → φπ + π − and B + → φK + decays do not contribute significantly. The B 0 → φK * 0 decay, resulting from a misidentification of a pion as a kaon, is vetoed by rejecting candidates which simultaneously have K + π − (K + K − K + π − ) invariant masses within 50 (30) MeV/c 2 of the known K * 0 (B 0 ) masses. The K + π − and K + K − K + π − invariant masses are computed by taking the kaon with the highest probability of being misidentified as a pion and assigning it the pion mass. These vetoes reduce the number of B 0 → φK * 0 candidates to a negligible level. Similarly, the number of Λ 0 b → φpK − decays, resulting from a misidentification of a proton as a kaon, is estimated from data by assigning the proton mass to the final-state particle that has the largest probability to be a misidentified proton based on the particle-identification information. This method yields 241 ± 30 Λ 0 b → φpK − decays in the total data set.
In order to determine the B 0 s → φφ yield in the final data sample, the four-kaon    invariant-mass distributions are fitted with the sum of the following components: a B 0 s → φφ signal model, which comprises the sum of a Crystal Ball [39] and a Student's t-function; the peaking background contribution modelled by a Crystal Ball function, with the shape parameters fixed to the values obtained from a fit to simulated events, and the combinatorial background component, described using an exponential function. The yield of the Λ 0 b → φpK − peaking background contribution is fixed to the number previously stated. Once the MLP requirements are imposed, an unbinned extended maximum-likelihood fit to the four-kaon invariant mass gives a total yield of 8843 ± 102 B 0 s → φφ decays and 2813 ± 67 combinatorial background candidates in the total data set. The fits to the four-kaon invariant-mass distributions, after the selection optimised for the CP -violation measurement, separately for each data taking year, are shown in Fig. 1.

Formalism
The final state of the B 0 s → φφ decay comprises a mixture of CP eigenstates, which are disentangled by means of an angular analysis in the helicity basis. In this basis, the decay is described by three angles, θ 1 , θ 2 and φ, defined in Fig. 2

Decay-time-dependent model
As discussed in Sec. 1, the B 0 s → φφ decay is a P → V V decay. However, due to the proximity of the φ resonance to the scalar f 0 (980) resonance, there are irreducible Figure 2: Decay angles for the B 0 s → φφ decay, where θ 1,2 is the angle between the K + momentum in the φ 1,2 meson rest frame and the φ 1,2 momentum in the B 0 s rest frame, Φ is the angle between the two φ meson decay planes andn V 1,2 is the unit vector normal to the decay plane of the φ 1,2 meson.
contributions to the four-kaon mass spectrum from P → V S (S-wave) and P → SS (double S-wave) processes, where S denotes a scalar meson, or a nonresonant pair of kaons. Thus, the total amplitude is a coherent sum of P -, S-, and double S-wave processes, and is modelled by making use of the different dependence on the helicity angles associated with these terms, where the helicity angles are defined in Fig. 2. A randomised choice is made for which φ meson is used to determine θ 1 and which is used to determine θ 2 . The total amplitude (A) containing the P -, S-, and double S-wave components as a function of time, t, can be written as [40] where A 0 , A , and A ⊥ are the CP -even longitudinal, CP -even parallel, and CP -odd perpendicular polarisations of the B 0 s → φφ decay. The P → V S and P → SS processes are described by the A S and A SS amplitudes, respectively, where P → V S is CP -odd and P → SS is CP -even. The resulting differential decay rate is proportional to the square of the total amplitude and consists of 15 terms [40] dΓ dt d cos where the f i terms are functions of the angular variables and the time-dependence is contained in The coefficients N i , a i , b i , c i and d i , which are functions of the CP observables, are defined in Appendix A. ∆Γ s ≡ Γ L − Γ H is the decay-width difference between the light and heavy B 0 s mass eigenstates, Γ s ≡ (Γ L + Γ H )/2 is the average decay width, and ∆m s is the B 0 s -B 0 s oscillation frequency. The differential decay rate for a B 0 s meson produced at t = 0 is obtained by changing the sign of the c i and d i coefficients. The amplitudes of helicity state k are expressed as where g + (t) and g − (t) describe the time evolution of B 0 s and B 0 s mesons, respectively. CP violation is parameterised through where, q and p relate the light and heavy mass eigenstates to the flavour eigenstates and η k is the CP eigenvalue of the polarisation being considered. Defining the amplitude in this way leads to the forms of N i , a i , b i , c i and d i , listed in Table 7 (Appendix A). The CP -violating asymmetry in B 0 s mixing, which can be characterised by the semileptonic asymmetry, a s sl is small [41]. Thus, to good approximation |q/p| = 1, and |λ k | quantifies the level of CP violation in the decay. Two different fit configurations are performed, one in which the CP -violation parameters are assumed to be helicity independent and the other in which CP -violation parameters are allowed to differ as a function of helicity. The helicity independent fit assumes one CP -violating phase, φ sss s , which takes the place of all φ s,k contained in the coefficients of Appendix A, and likewise one parameter that describes direct CP violation, |λ|, which takes the place of all λ k coefficients. Due to the small sample size, the number of degrees of freedom is reduced for the case of the helicity-dependent CP -violation fit. This involves assuming CP conservation for the case of the direct CP -violation parameters, λ = 1, and also for the phase of the longitudinal polarisation, φ sss s,0 = 0. The longitudinal polarisation has been theoretically calculated as close to zero in the B 0 s → φφ decay [1]. The φ sss s and |λ| parameters are measured with respect to contributions with the same flavour content as the φ meson, i.e. ss. Regarding the S-wave and double S-wave terms, the impact of the non-ss component of the φ wavefunction is negligible in this analysis.

Triple-product asymmetries
Scalar triple products of three-momentum or spin vectors are odd under time reversal, T . Nonzero asymmetries for these observables can either be due to a CP-violating phase or from CP-conserving strong final-state interactions. Four-body final states give rise to three independent momentum vectors in the rest frame of the decaying B 0 s meson. For a detailed review of the phenomenology the reader is referred to Ref. [10].
Two triple products can be defined: wheren V i (i = 1, 2) is a unit vector perpendicular to the vector meson (V i ) decay plane andp V 1 is a unit vector in the direction of V 1 in the B 0 s rest frame, defined in Fig. 2. This then provides a method of probing CP violation without the need to measure the decay time or the initial flavour of the B 0 s meson. It should be noted, that while the observation of nonzero triple-product asymmetries implies CP violation or final-state interactions (in the case of B 0 s meson decays), measurements of triple-product asymmetries consistent with zero do not rule out the presence of CP -violating effects, as the size of the asymmetry also depends on the differences between the strong phases [10].
In the B 0 s → φφ decay, two triple products are defined as U ≡ sin Φ cos Φ and V ≡ sin(±Φ) where the positive sign is taken if cos θ 1 cos θ 2 ≥ 0 and the negative sign otherwise [10]. The T -odd asymmetry corresponding to the U observable, A U , is defined as the normalised difference between the number of decays with positive and negative values of sin Φ cos Φ, Similarly, A V is defined as Here, A ⊥ , A and A 0 correspond to the three transversity amplitudes. The determination of the triple-product asymmetries is then reduced to a simple counting experiment. Comparing these formulae with Eq. 3 and Appendix A it can be seen that the triple products are related to the K 4 (t) and K 6 (t) terms in the decay amplitude.

Decay-time resolution
The sensitivity to φ sss s is affected by the accuracy of the measured decay time. In order to resolve the fast B 0 s -B 0 s oscillations, it is necessary to have a decay-time resolution that is much smaller than the oscillation period. To account for the resolution of the measured decay-time distribution, all decay-time-dependent terms are convolved with a Gaussian function, with width σ t i that is estimated for each candidate, i, based upon the uncertainty obtained from the vertex and kinematic fit [42].
In order to apply a candidate-dependent resolution model during fitting, the estimated per-event decay time uncertainty is needed. This is calibrated using the fact the decay time resolution for the B 0 s → φφ mode is dominated by the secondary vertex resolution. A sample of good-quality tracks, which originate from the primary interaction vertex is selected. Due to the small opening angle of the kaons in the decay of a φ meson, it is sufficient to use a single prompt track and assign it the mass of a φ meson. When combining this with another pair of tracks, the invariant mass of the three-body combination is required to be within 250 MeV/c 2 of the known B 0 s mass. That the decay-time resolution of the signal B 0 s decays can be described well using three tracks has been validated using simulation.
A linear function is then fitted to the distribution of σ t i versus σ t true , with parameters q 0 and q 1 . Here, σ t true denotes the difference between reconstructed decay time and the exact decay time of simulated signal. The per-event decay-time uncertainty used in the decay-time-dependent fit is then calculated as σ cal i = q 0 + q 1 σ t i . Gaussian constraints are used to account for the uncertainties on the calibration parameters in the decay-timedependent fit. The effective single-Gaussian decay-time resolution is found to be between 41 and 44 fs, depending on the data-taking year, in agreement with the expectation from the simulation.

Acceptances
The B 0 s → φφ differential decay rate depends on the decay time and three helicity angles as shown in Eq. 2. Good understanding of the efficiencies in these variables is required. The decay-time and angular acceptances are assumed to factorise. Control channels show this assumption has a negligible systematic uncertainty on the physics parameters.

Angular acceptance
The geometry of the LHCb detector and the momentum requirements imposed on the final-state particles introduce distortions of the helicity angles, giving rise to acceptance effects. Simulated signal events, selected with the same criteria as those applied to data are used to determine these efficiency corrections. The angular acceptances as a function of the three helicity angles are shown in Fig. 3.
The efficiency is parameterised in terms of the decay angles as where Ω depends on the decay angles, cos θ 1 , cos θ 2 and φ, the c ijk are coefficients, P i (cos θ 1 ) are Legendre polynomials, and Y jk (cos θ 2 , Φ) are spherical harmonics. The procedure followed to calculate the coefficients is described in detail in Ref. [43] and exploits the orthogonality of Legendre polynomials. The coefficients are given by This integral is calculated by means of a Monte Carlo technique, which reduces the integral to a sum over the number of accepted simulated events (N obs ) where P gen is the probability density function (PDF) without acceptance where the parameters are set to values used in the Monte Carlo generation. In order to easily incorporate the angular acceptance, it is convenient to write angular functions of Eq. 2 in the same basis as the efficiency parameterisation, i.e.
where P ij (cos θ 1 ) are the associated Legendre polynomials, κ ijkl,a are coefficients and a numerates the 15 terms outlined earlier. The parameterisation for each angular function is given in Table 1.
The normalisation of the angular component in the decay-time dependent fit occurs through the 15 integrals ζ k = (Ω)f k (Ω)dΩ, where (Ω) is the efficiency as a function of the helicity angles as shown in Eq. 10 and f k (Ω) are the angular functions as defined in Eq. 13.
The angular acceptance is calculated correcting for the differences in kinematic variables between data and simulation. This includes differences in the MLP training variables that can affect acceptance corrections through correlations with the helicity angles.
The fit to determine the triple-product asymmetries assumes that the U and V observables are symmetric in the acceptance corrections. Simulation is used to assign a systematic uncertainty related to this assumption.

Decay-time acceptance
The impact-parameter requirements on the final-state particles efficiently suppress the background from the numerous pions and kaons originating from the PV, but introduce a decay-time dependence in the selection efficiency.
The efficiency as a function of the decay time is taken from the B 0 s → D − s (→ K + K − π − )π + decay, in the case of data taken between 2011 and 2012, and from the B 0 → J/ψ (→ µ + µ − )K * 0 (→ K + π − ) decay in the case of data taken between 2015 and 2016. The reason for the change in control channel is related to changes  to the software-trigger selection between the two data-taking periods. The decay-time acceptances of the control modes are weighted by a multivariate algorithm based on simulated kinematic and topological information, in order to match more closely those of the signal B 0 s → φφ decay. Cubic splines are used to model the acceptance as a function of decay time in the PDF. The PDF can then be computed analytically with the inclusion of the decay-time acceptance following Ref. [44]. Example decay-time acceptances are shown for the case of the B 0 s → D − s π + and B 0 → J/ψ K * 0 decays in Fig. 4. To simplify the measurement of the triple-product asymmetries, the decay-time acceptance is not applied in the fit to determine the triple-product asymmetries. The time acceptance correction has an impact on the asymmetry of 0.3% and is treated as a source of systematic uncertainty, as further described in Sec. 9.3.

Flavour tagging
To obtain sensitivity to φ sss s , the flavour of the B 0 s meson at production must be determined. At LHCb, tagging is achieved through the use of various algorithms described in Refs. [45,46]. With these algorithms, the flavour-tagging power, defined as tag D 2 can be evaluated. Here, tag is the flavour-tagging efficiency defined as the fraction of candidates with a flavour tag with respect to the total, and D ≡ (1 − 2ω) is the dilution, where ω is the average fraction of candidates with an incorrect flavour assignment. This analysis uses opposite-side (OS) and same-side kaon (SSK) flavour taggers.
The OS flavour-tagging algorithm [45] makes use of the b (b) hadron produced in association with the signal b (b) hadron. In this analysis, the predicted probability of an incorrect flavour assignment, ω, is determined for each candidate by a neural network that is calibrated using B + → J/ψ K + , B + → D 0 π + , B 0 → J/ψ K * 0 , B 0 → D * − µ + ν µ , and B 0 s → D − s π + data as control modes. Details of the calibration procedure can be found in Ref. [47].
When a signal B 0 s meson is formed, there is an associated s quark formed in the first branches of the fragmentation that about 50 % of the time forms a charged kaon, which is likely to originate close to the B 0 s meson production point. The kaon charge therefore allows for the identification of the flavour of the signal B 0 s meson. This principle is exploited by the SSK flavour-tagging algorithm [46], which is calibrated with the B 0 s → D − s π + decay mode. A neural network is used to select fragmentation particles, improving the flavour-tagging power quoted in the previous decay-time-dependent measurement [9]. Table 2 shows the tagging power for the candidates tagged by only one of the algorithms and those tagged by both. Uncertainties due to the calibration of the flavour tagging algorithms are applied as Gaussian constraints in the decay-time-dependent fit. The initial flavour of the B 0 s meson established from flavour tagging is accounted for during fitting.

.1 Likelihood fit
The fit parameters in the polarisation-independent fit are the CP violation parameters, φ sss s and |λ|, the squared amplitudes, |A 0 | 2 , |A ⊥ | 2 , |A S | 2 , and |A SS | 2 , and the strong phases, δ ⊥ , δ , δ 0 , δ S , and δ SS , as defined in Sec. 4.1. The P -wave amplitudes are defined such that |A 0 | 2 + |A ⊥ | 2 + |A | 2 = 1, hence only two of the three amplitudes are free parameters. This normalisation effectively means the S and SS components are measured relative to the P -wave. The polarisation-dependent fit allows for a perpendicular, parallel and longitudinal component of φ sss s and |λ|. The measurement of the parameters of interest is performed through an unbinned negative log likelihood minimisation. The log-likelihood, L, of each candidate is weighted using the sPlot method [48,49], to remove partly reconstructed and combinatorial background. The negative log-likelihood then takes the form where W e are the signal sPlot weights calculated using the four-kaon invariant mass as the discriminating variable. The correlations between the angular and decay-time variables used in the fit with the four-kaon mass are small enough for this technique to be appropriate. The factor α = e W e / e W 2 e accounts for the sPlot weights in the determination of the statistical uncertainties. The parameter S e TD is the differential decay rate of Eq. 2, modified to the effects of decay-time and angular acceptance, in addition to the probability of an incorrect flavour tag. Explicitly, this can be written as where ζ k are the normalisation integrals used to describe the angular acceptance described in Sec. 6.1 and The calibrated probability of an incorrect flavour assignment is given by ω e , R denotes the Gaussian time-resolution function, and the ⊗ denotes a convolution operation. In Eq. 16, q e = 1 (−1) for a B 0 s (B 0 s ) meson at t = 0 or q e = 0 where no flavour-tagging information is assigned. The data samples corresponding to the different years of data taking are assigned independent signal weights, decay-time and angular acceptances, and separate Gaussian constraints are applied to the decay-time resolution parameters, as defined in Sec. 5. The B 0 s -B 0 s oscillation frequency is constrained to the value measured by LHCb of ∆m s = 17.768 ± 0.023 (stat) ± 0.006 (syst) ps −1 [50], with the assumption that the systematic uncertainties are uncorrelated with those of the current measurement. The values of the decay width and decay-width difference are constrained to the current best known values of Γ s = 0.6646 ± 0.0020 ps −1 and ∆Γ s = 0.086 ± 0.006 ps −1 [51].
Correction factors must be applied to each of the S-wave and double S-wave interference terms in the differential decay width. These factors modulate the sizes of the contributions of the interference terms in the angular PDF due to the different line-shapes of kaon pairs originating from spin-1 and spin-0 configurations. This takes the form of a multiplicative factor for each time a S-wave pair of kaons interferes with a P -wave pair. Their K + K − invariant-mass parameterisations are denoted by g(m K + K − ) and h(m K + K − ), respectively. The P -wave configuration is described by a Breit-Wigner function and the S-wave configuration is assumed to be uniform. The correction factors, denoted by C SP , are defined in Ref. [47] where m h and m l are the upper and lower edges of the m K + K − window and the phase of C SP is absorbed in the measurements of δ S − δ ⊥ . The factor |C SP |, is calculated to be 0.36. In order to determine systematic uncertainties due to the model dependence of the S-wave, C SP factors are recalculated based on the S-wave originating from an f 0 (980) resonance and incorporating the effects of the m K + K − resolution. These alternative assumptions on the P -wave and S-wave lineshapes yield a |C SP | value of 0.34, which is found to have a negligible effect on the parameter estimation.

Results
The resulting parameters are given in Table 3. A polarisation-independent fit is performed to calculate values for φ sss s and |λ|. A negligible fraction of S-wave and double S-wave is observed.
In addition, the CP -violating phases are also determined in a polarisation-dependent manner. Due to limited size of data samples, the phases φ s, and φ s,⊥ are measured under the assumptions that the longitudinal weak phase is CP -conserving and that there is no direct CP violation. In addition, all S-wave and double S-wave components of the fit are set to zero. The results of the polarisation dependent fit are shown in Table 4. The results for |A 0 | 2 , |A ⊥ | 2 , δ ⊥ and δ are not shown but are in agreement with the results reported in Table 3.
The correlation matrices for the two fit configurations are provided in Appendix B. Correlations with such decay-time-dependent measurements depend on the central values of the parameters. No large correlation is expected between the CP -violating parameters when the central values are consistent with CP conservation. The largest correlations are found to be between the different decay amplitudes. Cross-checks are performed on   simulated data sets generated with the same yield as observed in data, and with the same physics parameters, to establish that the generated values are recovered without biases. Figure 5 shows the distributions of the B 0 s decay time and the three helicity angles. Superimposed are the projections of the fit result. The projections include corrections for acceptance effects. Pseudoexperiments were used to confirm that the deviation of the data around cos θ 2 = ±0.5 from the resulting distribution of the fit is compatible with a statistical fluctuation.

Systematic uncertainties
Various sources of systematic uncertainty are considered in addition to those applied as Gaussian constraints in the fit. These arise from the angular and decay-time acceptances, the mass model used to describe the mass distribution, the determination of the time resolution calibration, and the fit bias. A summary of the systematic uncertainties is given in Table 5.
An uncertainty due to the angular acceptance arises from the choice of weighting scheme described in Sec. 6. This is accounted for by performing multiple alternative weighting schemes for the weighting procedure, based on different kinematic variables in the decay. The largest variation is then assigned as the uncertainty. Further checks are performed to verify that the angular acceptance does not depend on the way in which the event was triggered.
Two sources of systematic uncertainty are considered concerning the decay-time acceptance. These are the statistical uncertainty from the spline coefficients, and also the residual disagreement between the weighted control mode and the signal decay acceptances (see Sec. 6.2). The former is evaluated through fitting the signal data set with 30 different spline functions, whose coefficients are varied according to the corresponding uncertainties. This study is performed with the three different choices of the knot points. The RMS of the fitted parameters is then assigned as the uncertainty. The residual disagreement between the control mode and the signal mode is accounted for with a simulation-based correction. Simplified simulation is used with the corrected acceptance and then fitted with the nominal acceptance. This process is repeated and the resulting bias on the fitted parameters is used as an estimate of the systematic uncertainty.
The uncertainty on the mass model is found by refitting the data with various alternative signal models, consisting of the sum of two Crystal Ball models, the sum of a double-sided Crystal Ball and a Gaussian model. In addition, a Chebyshev polynomial is used to describe the combinatorial background. The signal weights are recalculated and the largest deviation from the nominal fit results is used as the corresponding uncertainty.
Fit biases can arise in maximum-likelihood fits where the number of candidates is small compared to the number of free parameters. The effect of such a bias is taken as a systematic uncertainty which is evaluated by generating and fitting simulated data sets and taking the resulting bias as the uncertainty.
The uncertainties of the effective flavour-tagging power are included in the statistical uncertainty through Gaussian constraints on the calibration parameters, and amount to 10 % of the statistical uncertainty on the CP -violating phases.

.1 Likelihood
To determine the triple-product asymmetries, the data sets are divided according to the sign of the observables U and V . Simultaneous likelihood fits to the four-kaon mass distributions are preformed for the U and V variables separately. The set of free parameters in the fits to determine the U and V observables consists of their total yields and the asymmetries A U (V ) . The mass model is the same as that described in Sec. 3. The total PDF, D TP , is then of the form where k indicates the sum over the background components with corresponding PDFs, P j , and G S is the signal PDF, as described in Sec. 3. The parameters f S i found in Eq. 18 are related to the asymmetry, where S denotes the signal component of the four-kaon mass fit, as described in Sec. 3. Peaking backgrounds are assumed to be symmetric in U and V .

Results
The triple-product asymmetries found from the simultaneous fit described in Sec. 9.1 are measured separately for the 2015 and 2016 data. The results are combined by performing likelihood scans of the asymmetry parameters and summing the two years. This gives where the uncertainties are statistical only.

Systematic uncertainties
As for the case of the decay-time-dependent fit, significant contributions to the systematic uncertainty arise from the decay-time and angular acceptances. Minor uncertainties also result from the knowledge of the mass model of the signal and the composition of peaking backgrounds. The effect of the decay-time acceptance is determined through the generation of simulated samples including the decay-time acceptance and fitted with the method described in Sec. 9.1. The resulting deviation from the nominal fit results is used to assign a systematic uncertainty. The effect of the angular acceptance is evaluated by generating simulated data sets with and without the inclusion of the angular acceptance. The difference between the nominal fit results and the results obtained using the simulated samples including the decay-time acceptance is then used as a systematic uncertainty.
Uncertainties related to the mass model are evaluated using a similar approach to that described in Sec. 8.3. Additional uncertainties arise from the assumption that the peaking background is symmetric in U and V . The deviation observed without this assumption is then added as a systematic uncertainty. Similarly, the assumption that the combinatorial background has no asymmetry yields an identical uncertainty. The systematic uncertainties are summarised in Table 6.  10 Search for the B 0 → φφ decay The selection criteria for the B 0 → φφ mode are based on the B 0 s → φφ selection, with some modifications. The Punzi figure of merit [38] is used for the B 0 → φφ search, resulting in a more stringent MLP requirement. Furthermore, the uncertainty on the four-kaon mass is required to be less than 25 MeV/c 2 , corresponding to roughly 3σ separation between the B 0 s and B 0 mass peaks. The B 0 s → φφ decay is used as normalisation decay mode. The signal PDF for the mass of the B 0 meson is assumed to be the same as that of the B 0 s decay, with the modification of the resolution according to a scaling factor, which is defined as

Combination of Run 1 and Run 2 results
where m K is the known K + mass. Figure 6 shows the fit to the full data set. The Λ 0 b → φpK contribution is fixed to 109 candidates, following the same method described in Sec. 3. The fit returns a yield of 4.9 ± 9.2 B 0 → φφ decays.
The Confidence Levels (CL s ) method [52] is used to set a limit on the B 0 → φφ branching fraction. A total of 10,000 pseudoexperiments are used to calculate each point of the scan. Figure 7 shows the results of the CL s scan. At 90 % CL, N B 0 < 23.7. These limits are converted to a branching fraction using where N B 0 is the limit on the B 0 → φφ yield, and N B 0 s →φφ is the B 0 s yield from the fit displayed in Fig. 6. The relative reconstruction and selection efficiency of the B 0 s → φφ and B 0 → φφ decays, B 0 →φφ / B 0 s →φφ , is determined to be 0.986 using simulation. The ratio of the fragmentation functions has been measured at 7 and 8 TeV to be f s /f d = 0.259 ± 0.015 within the LHCb acceptance [53]. The production fraction at 13 TeV has been shown to be consistent with that of the 7 and 8 TeV data [54]. The B(B 0 s → φφ) = (1.84±0.05 (stat)±0.07 (syst)±0.11(f s /f d )±0.12 (norm))×10 −5 branching fraction is an external input taken from Ref. [24]. To set the limit, the uncertainties on the B 0 s → φφ branching fraction are propagated to the limit, where the uncertainty on the B 0 s → φφ branching fraction arising from f s /f d is already included in the uncertainty on the normalisation mode, B 0 → φK * . The maximum value of B(B 0 s → φφ) including the systematic contribution is found to be 1.99 × 10 −5 and is used in Eq. 22. This therefore translates to a limit of B(B 0 → φφ) < 2.7 (3.0) × 10 −8 at 90 % (95 %) CL, which supersedes the previous best limit.

Summary and conclusions
Measurements of CP violation in the B 0 s → φφ decay are presented, based on a sample of proton-proton collision data corresponding to an integrated luminosity of 5.0 fb −1 collected with the LHCb detector. The CP -violating phase, φ sss s , and CP violation parameter, |λ|, are determined in a helicity-independent manner to be φ sss s = −0.073 ± 0.115 (stat) ± 0.027 (syst) rad, |λ| = 0.99 ± 0.05 (stat) ± 0.01 (syst).
The polarisation amplitudes and strong phases measured in the polarisation-dependent fit are in agreement with the results listed here. In addition, values of the polarisation amplitudes are found to agree well with previous measurements [9,12,55,56] and with predictions from QCD factorisation [2,3]. The most precise measurements to date of the triple-product asymmetries are determined from a separate time-integrated fit to be A U = −0.003 ± 0.011 (stat) ± 0.004 (syst), A V = −0.014 ± 0.011 (stat) ± 0.004 (syst), in agreement with previous measurements [9,12,55]. The measured values of the CPviolating phase and triple-product asymmetries are consistent with the hypothesis of CP conservation in b → sss transitions.
In addition, the most stringent limit on the branching fraction of the B 0 → φφ decay is presented and it is found to be B(B 0 → φφ) < 2.7 × 10 −8 at 90 % CL.

A Time-dependent terms
In Table 7, δ S and δ SS are the strong phases of the P → V S and P → SS processes, respectively. The P -wave strong phases are dependent on δ and δ 0 .