Measurement of CP violation in $B_s^0 \to \phi \phi$ decays

A measurement of the decay time dependent CP-violating asymmetry in $B_s^0 \to \phi\phi$ decays is presented, along with measurements of the $T$-odd triple-product asymmetries. In this decay channel, the CP-violating weak phase arises from the interference between $B_s^0$-$\bar{B}_s^0$ mixing and the loop-induced decay amplitude. Using a sample of proton-proton collision data corresponding to an integrated luminosity of $3.0\, fb^{-1}$ collected with the LHCb detector, a signal yield of approximately 4000 $B_s^0 \to \phi\phi$ decays is obtained. The CP-violating phase is measured to be ${\phi_s =-0.17\pm0.15\mathrm{\,(stat)}\pm0.03\mathrm{\,(syst)}}$ rad. The triple-product asymmetries are measured to be ${A_U=-0.003\pm0.017\mathrm{\,(stat)}\pm0.006\mathrm{\,(syst)}}$ and ${A_V=-0.017\pm0.017\mathrm{\,(stat)}\pm0.006\mathrm{\,(syst)}}$. Results are consistent with the hypothesis of CP conservation.


Introduction
The B 0 s → φφ decay is forbidden at tree level in the Standard Model (SM) 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 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, which is characterised by a CP -violating weak phase. Due to different decay amplitudes the actual value of the weak phase is dependent on the B 0 s decay channel. 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 given by −2 arg (−V ts V * tb /V cs V * cb ) = −0.0364 ± 0.0016 rad [5]. The LHCb collaboration has measured the weak phase in the combination of B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decays to be 0.07 ± 0.09 (stat) ± 0.01 (syst) rad [6]. A recent analysis of B 0 s → J/ψ π + π − decays using the full LHCb Run I dataset of 3.0 fb −1 has measured the CP -violating phase to be 0.070 ± 0.068 (stat) ± 0.008 (syst) rad [7]. These measurements are consistent with the SM and place stringent constraints on CP violation in B 0 s -B 0 s oscillations [8]. The CP -violating phase, φ 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 |φ s | [1][2][3].
Triple-product asymmetries are formed from T -odd combinations of the momenta of the final state particles. Such asymmetries provide a method of measuring CP violation in a decay time integrated method that complements the decay time dependent measurement [9]. These asymmetries are calculated from functions of the angular observables and are expected to be close to zero in the SM [10]. Particle-antiparticle oscillations reduce non-zero triple-product asymmetries due to CP -conserving strong phases, known as "fake" triple-product asymmetries by a factor Γ/(∆m), where Γ and ∆m are the decay rates and oscillation frequencies of the neutral meson system in question. Since one has Γ s /(∆m s ) ≈ 0.04 for the B 0 s system, "fake" triple-product asymmetries are strongly suppressed, allowing for "true" CP -violating triple-product asymmetries to be calculated without the need to measure the initial flavour of the B 0 s meson [9]. Theoretical calculations can be tested further with measurements of the polarisation fractions, where the longitudinal and transverse polarisation fractions are denoted by f L and f T , respectively. In the heavy quark limit, f L is expected to be the dominant polarisation 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 [11][12][13][14][15][16]. However, the dynamics of penguin transitions are more complicated. In the context of QCDf, f L is predicted to be 0.36 +0. 23 −0.18 for the B 0 s → φφ decay [3]. In this paper, a measurement of the CP -violating phase in B 0 s → φ(→ K + K − )φ(→ K + K − ) decays, along with a measurement of the T -odd triple-product asymmetries is presented. The results are based on pp collision data corresponding to an integrated luminosity of 1.0 fb −1 and 2.0 fb −1 collected by the LHCb experiment at centre-of-mass energies √ s = 7 TeV in 2011 and 8 TeV in 2012, respectively. Previous measurements of the triple-product asymmetries from the LHCb [17] and CDF [18] collaborations, together with the first measurement of the CP -violating phase in B 0 s → φφ decays [19], have shown no evidence of deviations from the SM. The decay time dependent measurement improves on the previous analysis [19] through the use of a more efficient candidate selection and improved knowledge of the B 0 s flavour at production, in addition to a data-driven determination of the efficiency as a function of decay time.
The results presented in this paper supersede previous measurements of the CP -violating phase [19] and T -odd triple-product asymmetries [17], made using 1.0 fb −1 of data collected at a √ s = 7 TeV.

Detector description
The LHCb detector [20] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [21] placed downstream. The combined tracking system provides a momentum measurement with relative uncertainty that varies from 0.4% at low momentum to 0.6% at 100 GeV/c, and impact parameter resolution of 20 µm for tracks with large transverse momentum, p T . Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov (RICH) detectors [22]. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. The trigger [23] 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. The hardware trigger selects B 0 s → φφ candidates by requiring large transverse energy deposits in the calorimeters from at least one of the final state particles. 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 close to that of the φ meson or by using a topological b-hadron trigger. The topological software trigger requires a two-, three-or four-track secondary vertex with a large sum of the p T of the charged particles and a significant displacement from the primary pp interaction vertices (PVs). At least one charged particle should have p T > 1.7 GeV/c and χ 2 IP with respect to any primary interaction 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 [24] is used for the identification of secondary vertices consistent with the decay of a b-hadron.
In the simulation, pp collisions are generated using Pythia [25] with a specific LHCb configuration [26]. Decays of hadronic particles are described by EvtGen [27], in which final state radiation is generated using Photos [28]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [29] as described in Ref. [30].

Selection and mass model
Events passing the trigger are initially required to pass loose requirements on the fit quality of the four-kaon vertex fit, 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 φ meson candidates is required to be within 25 MeV/c 2 of the known φ mass [31].
In order to further separate the B 0 s → φφ signal from the background, a boosted decision tree (BDT) is implemented [32,33]. To train the BDT, simulated B 0 s → φφ events passing the same loose requirements as the data events are used as signal, whereas events in the fourkaon invariant mass sidebands from data are used as background. The signal mass region is defined to be less than 120 MeV/c 2 from the known B 0 s mass, m B 0 s [31]. The invariant mass sidebands are defined to be inside the region 120 < |m is the four-kaon invariant mass. Separate BDTs are trained for data samples collected in 2011 and 2012, due to different data taking conditions in the different years. Variables used in the BDT consist of the minimum and maximum kaon p T and η, the minimum and the maximum p T and η of the φ candidates, the p T and η of the B 0 s candidate, the minimum probability of the kaon mass hypothesis using information from the RICH detectors, the quality of the four-kaon vertex fit, and the χ 2 IP of the B 0 s candidate. The BDT also includes kaon isolation asymmetries. The isolation variable is calculated as the scalar sum of the p T of charged particles inside a region defined as ∆ϕ 2 + ∆η 2 < 1, where ∆ϕ (∆η) is the difference in azimuthal angle (pseudorapidity), not including the signal kaon from the B 0 s decay. The asymmetry is then calculated as the difference between the isolation variable and the p T of the signal kaon, divided by the sum. After the BDT is trained, the optimum requirement on each BDT is chosen to maximise N S / √ N S + N B , where N S (N B ) represent the expected number of signal (background) events in the signal region of the data sample.
The presence of peaking backgrounds is extensively studied. The decay modes considered include B + → φK + , B 0 → φπ + π − , B 0 → φK * 0 , and Λ 0 b → φpK − , of which only the last two are found to contribute, and are the result of a mis-identification of a pion or proton as a kaon, respectively. The number of B 0 → φK * 0 events present in the data sample is determined from scaling the number of B 0 → φK * 0 events seen in data through a different dedicated selection with the relative efficiencies between the two selections found from simulated events. This method yields values of 7.3 ± 0.4 and 17.8 ± 0.9 events in the 2011 and 2012 datasets, respectively. The amount of Λ 0 b → φpK − decays is estimated directly from data by changing the mass hypothesis of the final-state particle most likely to have the mass of the proton from RICH detector information. This method yields 52 ± 19 and 51 ± 29 Λ 0 b → φpK − events in the 2011 and 2012 datasets, respectively. In order to correctly determine the number of B 0 s → φφ events in the final data sample, the four-kaon invariant mass distributions are fitted with the B 0 s → φφ signal described by a double Gaussian model, and the combinatorial background component described using an exponential function. The peaking background contributions are fixed to the shapes found in simulated events. The yields of the peaking background contributions The use of the four-kaon invariant mass to assign signal weights allows for a decay time dependent fit to be performed with only the signal distribution explicitly described. The method for assigning the signal weights is described in greater detail in Sec. 8.1.

Phenomenology
The B 0 s → φφ decay is composed of a mixture of CP eigenstates, that are disentangled by means of an angular analysis in the helicity basis, defined in Fig. 2.

Decay time dependent model
The B 0 s → φφ decay is a P → V V decay, where P denotes a pseudoscalar and V a vector meson. However, due to the proximity of the φ resonance to that of the f 0 (980), there will also be contributions from S-wave (P → V S) and double S-wave (P → SS) processes, Figure 2: Decay angles for the B 0 s ! decay, where the K + momentum in the 1,2 rest frame and the parent 1,2 momentum in the rest frame of the B 0 s meson span the two meson decay planes, ✓ 1,2 is the angle between the K + track momentum in the 1,2 meson rest frame and the parent 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.
where S denotes a spin-0 meson or a pair of non-resonant kaons. Thus the total amplitude is a coherent sum of P -, S-, and double S-wave processes, and is accounted for during fitting by making use of the di↵erent functions of the helicity angles associated with these terms. The choice of which meson is used to determine ✓ 1 and which is used to determine ✓ 2 is randomised. The total amplitude (A) containing the P -, S-, and double S-wave components as a function of decay time, t, can be written as [36] where A 0 , A k , 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. The di↵erential decay rate may be found through the square of the total amplitude leading to the 15 terms [36] The K i (t) term can be written as 5 Figure 2: Decay angles for the B 0 s → φφ decay, where the K + momentum in the φ 1,2 rest frame and the parent φ 1,2 momentum in the rest frame of the B 0 s meson span the two φ meson decay planes, θ 1,2 is the angle between the K + track momentum in the φ 1,2 meson rest frame and the parent φ 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.
where S denotes a spin-0 meson or a pair of non-resonant kaons. Thus the total amplitude is a coherent sum of P -, S-, and double S-wave processes, and is accounted for during fitting by making use of the different functions of the helicity angles associated with these terms. The choice of which φ meson is used to determine θ 1 and which is used to determine θ 2 is randomised. The total amplitude (A) containing the P -, S-, and double S-wave components as a function of decay time, t, can be written as [35] 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. The differential decay rate may be found through the square of the total amplitude leading to the 15 terms [35] The K i (t) term can be written as where the coefficients are shown in Table 1, ∆Γ 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 three CP -violating terms introduced in Table 1 are defined as where φ s measures CP violation in the interference between the direct decay amplitude and that via mixing, λ ≡ (q/p)(A/A), q and p are the complex parameters relating the B 0 s flavour and mass eigenstates, and A (A) is the decay amplitude (CP conjugate decay amplitude). Under the assumption that |q/p| = 1, |λ| measures direct CP violation. The CP violation parameters are assumed to be helicity independent. The association of φ s and |λ| with S-wave and double S-wave terms implies that these consist solely of contributions with the same flavour content as the φ meson, i.e. an ss resonance.

Triple-product asymmetries
Scalar triple products of three momentum or spin vectors are odd under time reversal, T . Non-zero asymmetries for these observables can either be due to a CP -violating phase or a CP -conserving phase and 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. [9].
The two independent terms in the time dependent decay rate that contribute to a T -odd asymmetry are the K 4 (t) and K 6 (t) terms, defined in Eq. 3. The triple products that allow access to these terms are 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 non-zero triple-product asymmetries implies CP violation or final state interactions (in the case of B 0 s meson decays), the measurements of triple-product asymmetries consistent with zero do not rule out the presence of CP -violating effects, as strong phase differences can cause suppression [9].
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 negative sign otherwise.
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 Extraction of the triple-product asymmetries is then reduced to a simple counting exercise.

Decay time resolution
The sensitivity to φ s is affected by the accuracy of the measured decay time. In order to resolve the fast B 0 s -B 0 s oscillation period of approximately 355 fs, it is necessary to have a decay time resolution that is much smaller than this. To account for decay time resolution, all decay time dependent terms are convolved with a Gaussian function, with width σ t i that is estimated for each event, i, based upon the uncertainty obtained from the vertex and kinematic fit. In order to apply an event-dependent resolution model during fitting, the estimated per-event decay time uncertainty must be calibrated. This is done using simulated events that are divided into bins of σ t i . For each bin, a Gaussian function is fitted to the difference between reconstructed decay time and the true decay time to determine the resolution σ t true . A first-order polynomial is then fitted to the distribution of σ t i versus σ t true , with parameters denoted by q 0 and q 1 . The calibrated 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 time dependent fit. Cross-checks, consisting of the variation of an effective single Gaussian resolution far beyond the observed differences in data and simulated events yield negligible modifications to results, hence no systematic uncertainty is assigned. The results are verified to be largely insensitive to the details of the resolution model, as supported by tests on data and observed in similar measurements [6].
The effective single Gaussian resolution is found from simulated datasets to be 41.4 ± 0.5 fs and 43.9 ± 0.5 fs for the 2011 and 2012 datasets, respectively. Differences in the resolutions from 2011 and 2012 datasets are expected due to the independent selection requirements.

Acceptances
The four observables used to analyse B 0 s → φφ events consist of the decay time and the three helicity angles, which require a good understanding of efficiencies in these variables. It is assumed that the decay time and angular acceptances factorise.

Angular acceptance
The geometry of the LHCb detector and the momentum requirements imposed on the final-state particles introduce efficiencies that vary as functions of the helicity angles. Simulated events with the same selection criteria as those applied to B 0 s → φφ data events are used to determine this efficiency correction. Efficiencies as a function of the three helicity angles are shown in Fig. 3.
Acceptance functions are included in the decay time dependent fit through the 15 integrals (Ω)f k (Ω)dΩ, where f k are the angular functions given in Table 1 and (Ω) is the efficiency as a function of the set of helicity angles, Ω. The inclusion of the integrals in the normalisation of the probability density function (PDF) is sufficient to describe the  angular acceptance as the acceptance factors for each event appear as a constant in the log-likelihood, the construction of which is described in detail in Sec. 8.1, and therefore do not affect the fitted parameters. The method for the calculation of the integrals is described in detail in Ref. [36]. The integrals are calculated correcting for the differences between data and simulated events. This includes differences in the BDT 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. Simulated events are then 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 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 B 0 s → D − s (→ K + K − π − )π + data events, with an upper limit of 1 ps applied to the D − s decay time to ensure topological similarity to the B 0 s → φφ decay. After the same decay time-biasing selections are applied to the B 0 s → D − s π + decay as used in the B 0 s → φφ decay, B 0 s → D − s π + events are reweighted according to the minimum track transverse momentum to ensure the closest For the case of the decay time dependent fit, the efficiency as a function of the decay time is modelled as a histogram, with systematic uncertainties arising from the differences in B 0 s → φφ and B 0 s → D − s π + simulated events. Figure 4 shows the comparison of the efficiency as a function of decay time calculated using B 0 s → D − s π + data in 2011 and 2012. Also shown is the comparison between B 0 s → φφ and B 0 s → D − s π + simulated events. In the fit to determine the triple-product asymmetries, the decay time acceptance is treated only as a systematic uncertainty, which is based on the acceptance found from B 0 s → D − s π + data events.

Flavour tagging
To maximise the sensitivity on φ s , the determination of the initial flavour of the B 0 s meson is necessary. This results from the terms in the differential decay rate with the largest sensitivity to φ s requiring the identification (tagging) of the flavour at production. At LHCb, tagging is achieved through the use of different algorithms described in Refs. [6,38]. This analysis uses both the opposite side (OS) and same side kaon (SSK) flavour taggers.
The OS flavour tagging algorithm [39] makes use of the b(b)-quark produced in association with the signal b(b)-quark. In this analysis, the predicted probability of an incorrect flavour assignment, ω, is determined for each event 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 π +  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 [38]. The SSK tagger 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 [19,40].
Parameter estimation is achieved from a minimisation of the negative log likelihood. The likelihood, L, is weighted using the sPlot method [41, 42], with the signal weight of an event e calculated from the equation where j sums over the number of fit components to the four-kaon invariant mass, with PDFs F , associated yields N , and V sj is the covariance between the signal yield and the yield associated with the j th fit component. The log-likelihood then takes the form where α = e W e / e W 2 e is used to account for the weights in the determination of the statistical uncertainties, and S TD is the signal model of Eq. 2, accounting also for 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 where ω e is the calibrated probability of an incorrect flavour assignment, and R denotes the Gaussian resolution function. In Eq. 14, q e = 1(−1) for a B 0 s (B . The Gaussian constraints applied to the Γ s and ∆Γ s parameters use the combination of the measured values from B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decays. Constraints are therefore applied taking into account a correlation of 0.1 for the statistical uncertainties [6]. The systematic uncertainties are taken to be uncorrelated between the B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decay modes. The events selected in this analysis are within the two-kaon invariant mass range 994.5 < m K + K − < 1044.5 MeV/c 2 , and are divided into three regions. These correspond to both φ candidates with invariant masses smaller than the known φ mass, one φ candidate

Parameter
Best fit value φ s ( rad) −0.17 ± 0.15 |λ| 1.04 ± 0.07 0.102 ± 0.012 ∆m s ( ps −1 ) 17.774 ± 0.024 with an invariant mass smaller than the known φ mass and one larger, and a third region in which both φ candidates have invariant masses larger than the known φ mass. Binning the data in this way allows the analysis to become insensitive to correction factors that must be applied to each of the S-wave and double S-wave interference terms in the differential cross section. These factors modulate 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. Their parameterisations are denoted by g(m K + K − ) and h(m K + K − ), respectively. The spin-1 configuration is described by a Breit-Wigner function and the spin-0 configuration is assumed to be approximately uniform. The correction factors, denoted by C SP , are defined from the relation [6] where m h and m l are the upper and lower edges of a given m K + K − bin, respectively. Alternative assumptions on the P -wave and S-wave lineshapes are found to have a negligible effect on the parameter estimation. A simultaneous fit is then performed in the three m K + K − invariant mass regions, with all parameters shared except for the fractions and strong phases associated with the S-wave and double S-wave, which are allowed to vary independently in each region. The correction factors are calculated as described in Ref. [6]. The correction factor used for each region is calculated to be 0.69.

Results
The results of the fit to the parameters of interest are given in Table 3. The S-wave and double S-wave parameter estimations for the three regions defined in Sec. 8.1 are given in Table 4. The fraction of S-wave is found to be consistent with zero in all three mass regions. Table 4: S-wave and double S-wave results of the decay time dependent fit for the three regions identified in Sec. 8.1, where M −− indicates the region with both two-kaon invariant masses smaller than the known φ mass, M −+ the region with one smaller and one larger, and M ++ indicates the region with both two-kaon invariant masses larger than the known φ mass.

Region
|A 0.001 ± 0.003 −2.58 ± 2.08 0.020 ± 0.022 0.53 ± 0.55 Table 5: Correlation matrix associated with the result of the decay time dependent fit. Correlations with a magnitude greater than 0.5 are shown in bold. The correlation matrix is shown in Table 5. The largest correlations are found to be between the amplitudes themselves and the CP -conserving strong phases themselves. The observed correlations have been verified with simulated datasets. Cross-checks are performed on simulated datasets generated with the same number of events as observed in data, and with the same physics parameters, to ensure that generation values are recovered with negligible 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 are event weighted to yield the signal distribution and include acceptance effects.
The scan of the natural logarithm of the likelihood for the φ s parameter is shown in Fig. 6. At each point in the scan, all other parameters are re-minimised. A parabolic minimum is observed, and a point estimate provided. The shape of the profile log-likelihood is replicated in simplified simulations as a cross-check.

Systematic uncertainties
The most significant systematic effects arise from the angular and decay time acceptances. Minor contributions are also found from the mass model used to construct the event weights, the uncertainty on the peaking background contributions, and the fit bias. An uncertainty due to the angular acceptance arises from the limited number of simulated events used to determine the acceptance correction. This is accounted for by varying the normalisation weights within their statistical uncertainties accounting for correlations. The varied weights are then used to fit simulated datasets. This process is repeated and the width of the Gaussian distribution is used as the uncertainty. A further uncertainty arises from the assumption that the angular acceptance does not depend on the algorithm used for the initial flavour assignment. Such a dependence can be expected due to the kinematic correlations of the tagging particles with the signal particles. This introduces a tagging efficiency based on the kinematics of the signal particles. The difference between the nominal data result and the result with angular acceptances calculated independently for the different flavour tagging algorithms leads to a non-negligible uncertainty on the polarisation amplitudes. Further checks are performed to verify that the angular acceptance does not depend on the way in which the event was triggered.
[rad] The systematic uncertainty on the decay time acceptance is evaluated from the difference in the decay time acceptance evaluated from B 0 s → φφ and B 0 s → D − s π + simulated events. The simulated datasets are generated with the decay time acceptance of B 0 s → φφ simulation and then fitted with the B 0 s → D − s π + decay time 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 signal weights derived from a single Gaussian B 0 s → φφ model, rather than the nominal double Gaussian. The uncertainty due to peaking background contributions is found through the recalculation of the signal weights with peaking background contributions varied according to the statistical uncertainties on the yields of the Λ 0 b → φpK − and B 0 → φK * 0 contributions. Fit bias arises in likelihood fits when the number of events used to determine the free parameters is not sufficient to achieve the Gaussian limit. This uncertainty is evaluated by generating and fitting simulated datasets and taking the resulting bias as the uncertainty.
Uncertainties due to flavour tagging 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 phase.
A summary of the systematic uncertainties is given in Table 6.

Likelihood
In order to determine the triple-product asymmetries, a separate likelihood fit is performed. This is based around the simultaneous fitting of separate datasets to the four-kaon invariant mass, which are split according to the sign of U and V observables. Simultaneous mass Parameter fits are performed for the U and V observables separately. The set of free parameters in fits to determine the U and V observables consist of the asymmetries of the B 0 s → φφ signal and combinatoric background (A U (V ) and A B U (V ) ), along with their associated total yields (N S and N B ). The mass model is the same as that described in Sec. 3. The total PDF, S TP is then of the form where j indicates the sum over the background components with corresponding PDFs, P j , and G S is the double Gaussian signal PDF as described in Sec. 3. The parameters f k i found in Eq. 16 are related to the asymmetry, A k U (V ) , through where k denotes a four-kaon mass fit component, as described in Sec. 3. Peaking backgrounds are assumed to be symmetric in U and V .

Results
The background-subtracted distributions of the U and V observables are shown in Fig. 7 for the mass range 5246.8 < m K + K − K + K − < 5486.8 MeV/c 2 . Distributions are found to agree between 2011 and 2012 datasets and show qualitatively symmetric distributions. The triple-product asymmetries found from the simultaneous fit described in Sec. 9.1 are measured to be Statistical uncertainties are therefore to have approximately halved with respect to the previous LHCb measurements [17], due to more efficient selection requirements and a larger data sample, and are verified through fits to simulated datasets. No evidence for CP violation is found.

Systematic uncertainties
As for the case of the decay time dependent fit, the largest contributions to the systematic uncertainty arise from the decay time and angular acceptances. Minor uncertainties also result from the mass model and peaking background knowledge.
The effect of the decay time acceptance is determined through the generation of simulated samples including the decay time acceptance obtained from B 0 s → D − s π + data, and fitted with the method described in Sec. 9.1. The resulting bias is used to assign a systematic uncertainty.
The effect of the angular acceptance is evaluated by generating simulated datasets with and without the inclusion of the angular acceptance. The resulting bias found on the fit results of the triple-product asymmetries is then used as a systematic uncertainty.
Uncertainties related to the mass model are evaluated by taking the difference between the nominal fit results and those using a single Gaussian function to model the B 0 s → φφ decay. The effect of the peaking background is evaluated by taking the largest difference between the nominal fit results and the fit results with the peaking background yields varied according to their uncertainties, as given in Sec. 3.
The total systematic uncertainty is estimated by choosing the larger of the two individual
Results are found to agree with the theoretical predictions [1][2][3]. When compared with the CP -violating phase measured in B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decays [6], these results show that no large CP violation is present either in B 0 s -B 0 s mixing or in the b → sss decay amplitude.
Values of the polarisation amplitudes are found to agree well with the previous measurements [17][18][19]. Measurements in other B → V V penguin transitions at the B factories generally give higher values of f L ≡ |A 0 | 2 [11][12][13][14][15][16]. The value of f L found in the B 0 s → φφ channel is almost equal to that in the B 0 s → K * 0 K * 0 decay [44]. As reported in Ref. [17], the results are in agreement with QCD factorisation predictions [2,3], but disfavour the perturbative QCD estimate given in Ref. [45]. The fractions of S-wave and double S-wave are found to be consistent with zero in all three regions of m K + K − mass.
[40] G. A. Krocker, Development and calibration of a same side kaon tagging algorithm and measurement of the B 0 s -B 0 s oscillation frequency ∆m s at the LHCb experiment, PhD thesis, CERN-THESIS-2013-213.