First measurement of the $CP$-violating phase $\phi_s^{d\overline{d}}$ in $B_s^0\to(K^+\pi^-)(K^-\pi^+)$ decays

A flavour-tagged decay-time-dependent amplitude analysis of $B_s^0\to(K^+\pi^-)(K^-\pi^+)$ decays is presented in the $K^{\pm}\pi^{\mp}$ mass range from 750 to 1600 MeV$/c^2$. The analysis uses $pp$ collision data collected with the LHCb detector at centre-of-mass energies of $7$ and $8$ TeV, corresponding to an integrated luminosity of $3.0$ fb$^{-1}$. Several quasi-two-body decay modes are considered, corresponding to $K^{\pm}\pi^{\mp}$ combinations with spin 0, 1 and 2, which are dominated by the $K_0^*(800)^0$ and $K_0^*(1430)^0$, the $K^*(892)^0$ and the $K_2^*(1430)^0$ resonances, respectively. The longitudinal polarisation fraction for the $B_s^0\to K^*(892)^0\overline{K}^*(892)^0$ decay is measured as $f_L=0.208 \pm 0.032 \pm 0.046$, where the first uncertainty is statistical and the second is systematic. The first measurement of the mixing-induced $CP$-violating phase, $\phi_s^{d\overline{d}}$, in $b\to d\overline{d}s$ transitions is performed, yielding a value of $\phi_s^{d\overline{d}}=-0.10$ $\pm$ $0.13$ (stat) $\pm$ $0.14$ (syst) rad.


Introduction
The CP -violating weak phases φ s arise in the interference between the amplitudes of B 0 s mesons directly decaying to CP eigenstates and those decaying to the same final state after B 0 s -B 0 s oscillation. The B 0 s → K * 0 K * 0 decay, 1 which in the Standard Model (SM) is 1 Throughout this article, charge conjugation is implied and K * 0 refers to the K * (892) 0 resonance, unless otherwise stated. dominated by the gluonic loop diagram shown in figure 1, has been discussed extensively in the literature as a benchmark test for the SM and as an excellent probe for physics beyond the SM [1][2][3][4][5][6][7]. New heavy particles entering the loop would introduce additional amplitudes and modify properties of the decay from their SM values. In general, the weak phase φ s depends on the B 0 s decay channel under consideration, and can be different between channels as it depends on the contributions from tree-and loop-level processes. The notation φ dd s is used when referring to the weak phase measured in b → dds transitions. For b → ccs transitions, e.g. B 0 s → J/ψ K + K − and B 0 s → J/ψ π + π − decays, the weak phase φ cc s has been measured by several experiments [8][9][10][11]. The world average reported by HFLAV, φ cc s = −0.021 ± 0.031 rad [8], is dominated by the LHCb measurement φ cc s = −0.010 ± 0.039 rad [9]. The LHCb collaboration has also measured the φ ss s phase in B 0 s → φφ transitions [12], reporting a value of φ ss s = −0.17 ± 0.15 rad. The decay B 0 s → K * 0 K * 0 , with K * 0 → K + π − and K * 0 → K − π + , was first observed by the LHCb collaboration, based on pp collision data corresponding to an integrated luminosity of 35 pb −1 at a centre-of-mass energy √ s = 7 TeV [13]. A branching fraction and a final-state polarisation analysis were reported. An updated analysis of the B 0 s → (K + π − )(K − π + ) decay was performed by LHCb using 1.0 fb −1 of data at √ s = 7 TeV [14]. In both analyses, the invariant mass of the two Kπ pairs 2 was restricted to a window of ±150 MeV/c 2 around the known K * 0 mass. This publication reports the first decay-time-dependent amplitude analysis of B 0 s → (K + π − )(K − π + ) decays using a Kπ mass window that extends from 750 to 1600 MeV/c 2 , approximately corresponding to the region between the Kπ production threshold and the D 0 → K − π + resonance. At the current level of sensitivity, the assumption of common CP -violating parameters for the contributing amplitudes is appropriate. Consequently, such a wide window provides a four-fold increase of the signal sample size with respect to the narrow window of 150 MeV/c 2 around the K * 0 mass. The analysis uses pp collision data collected by LHCb in 2011 and 2012 at √ s = 7 and 8 TeV, corresponding to an integrated luminosity of 3.0 fb −1 . In this study, nine different quasi-two-body decay 2 Hereafter the notation Kπ will stand for both K + π − and K − π + pairs.  Table 1. Quasi-two-body decay channels and corresponding polarisation amplitudes contributing to the B 0 s → (K + π − )(K − π + ) final state in the Kπ mass window from 750 to 1600 MeV/c 2 . The different contributions are identified by the spin j 1 (j 2 ) of the K + π − (K − π + ) pair and the helicity h. In cases where more than one amplitude contributes, the polarisations are defined as being longitudinal, parallel, or perpendicular, which are then denoted by 0, and ⊥ respectively, following the definitions given in ref. [15]. The subscripts 1 and 2 in the parallel and perpendicular helicities of the tensor-tensor component denote different spin states leading to a parallel or a perpendicular configuration, as discussed in appendix A. channels are considered, corresponding to the different possible combinations of Kπ pairs with spin 0, 1 or 2. Additional contributions were studied and found to be negligible in the phase-space region considered in this analysis. The Kπ spectrum is dominated by the K * 0 (800) 0 , K * 0 (1430) 0 , K * (892) 0 and K * 2 (1430) 0 resonances. Angular momentum conservation in the decay allows for one single amplitude in modes involving at least one scalar Kπ pair, three amplitudes for vector-vector or vector-tensor decays and five amplitudes for a tensor-tensor decay. These possibilities are listed in table 1. There is a physical difference between decay pairs of the form scalar-vector and vector-scalar. Namely, in the used convention, the spectator quark from the B 0 s decay (see figure 1) always ends up in the second Kπ pair. The CP -averaged fractions of the contributing amplitudes, f i , as well as their strong-phase differences, δ i , are determined together with the CP -violating weak phase φ dd s and a parameter that accounts for the amount of CP violation in decay, |λ|. This is the first time that the weak phase in b → dds transitions has been measured. It is also the first time that the tensor components in the (K + π − )(K − π + ) system have been studied.

Phenomenology
The phenomenon of quark mixing means that a B 0 s meson can oscillate into its antiparticle equivalent, B 0 s . Consequently, the physical states, B where p and q are complex coefficients that satisfy |p| 2 + |q| 2 = 1. The time evolution of the initially pure flavour eigenstates at t = 0, |B 0 s (0) and |B 0 s (0) , is described by where the decay-time-dependent functions g ± (t) are given by with m s and Γ s being the average mass and width of the B 0 s,H and B 0 s,L states. Negligible CP violation in mixing is assumed in this analysis, leading to the parameterisation q/p = e −iφ M , where φ M is the B 0 s -B 0 s mixing phase. The total decay amplitude of the flavour eigenstates at t = 0 into the final state f = (K + π − )(K − π + ), denoted by f |B 0 s (0) and f |B 0 s (0) , is a coherent sum of scalar-scalar (SS), scalar-vector (SV), vector-scalar (VS), scalar-tensor (ST), tensor-scalar (TS), vector-vector (VV), vector-tensor (VT), tensor-vector (TV) and tensor-tensor (TT) contributions. The quantum numbers used to label the (K + π − )(K − π + ) final states are the spin j 1 (j 2 ) of the K + π − (K − π + ) pair and the helicity h. The vector component is represented in this analysis by the K * 0 meson, since this resonance is found to be largely dominant in this spin configuration. Potential contributions from the K * 1 (1410) 0 and K * 1 (1680) 0 resonances are considered as sources of systematic uncertainty. For the tensor case, only the K * 2 (1430) 0 resonance contributes in the considered Kπ mass window. The scalar component, denoted in this paper by (Kπ) * 0 requires a more careful treatment. It can have contributions from the K * 0 (800) 0 and K * 0 (1430) 0 resonances and from a nonresonant Kπ component. The parameterisation of the Kπ invariant mass spectrum for the scalar contribution is explained later in this section. All of the considered decay modes, together with the quantum numbers for the corresponding amplitudes, are shown in table 1. In order to separate components with different CP eigenvalues, η j 1 j 2 h = ±1, the differential decay rate is expressed as a function of three angles and the two Kπ invariant masses. The angles θ 1 , θ 2 and ϕ, are written in the helicity basis and defined according to the diagram shown in figure 2. The invariant mass of the K + π − pair is denoted as m 1 , while that of the K − π + pair as m 2 . The symbol Ω is used to represent all three angles and the two invariant masses, Ω = (m 1 , m 2 , cos θ 1 , cos θ 2 , ϕ). Summing over the possible states and using the partial wave formalism, the decay amplitudes at t = 0 can be written as Figure 2. Graphical definition of the angles in the helicity basis. Taking the example of a B 0 s → Q 1 Q 2 decay (this analysis uses B 0 , with each final-state quasi-two-body meson decaying to pseudoscalars (Q 1 → K + π − and Q 2 → K − π + ), θ 1 (θ 2 ) is defined as the angle between the directions of motion of K + (K − ) in the Q 1 (Q 2 ) rest frame and Q 1 (Q 2 ) in the B 0 s rest frame, and ϕ as the angle between the plane defined by K + π − and the plane defined by K − π + in the B 0 s rest frame.
The complex parameters A j 1 j 2 h and A j 1 j 2 h contain the physics of the decays to the final states with j 1 , j 2 and h as defined in table 1. The angular terms, Θ j 1 j 2 h , are built from combinations of spherical harmonics as shown in appendix A. The η j 1 j 2 h factor is equal to The mass-dependent terms are parameterised as where F j 1 j 2 h (m 1 , m 2 ) is the Blatt-Weisskopf angular-momentum centrifugal-barrier factor [16] and M j describes the shape of the Kπ invariant mass of a Kπ pair with spin j. Relativistic Breit-Wigner functions of spin 1 and 2, parameterising the K * 0 and the K * 2 (1430) 0 resonances, are used for M 1 and M 2 , respectively. The parameterisation of M 0 is based on the phenomenological S-wave scattering amplitude of isospin 1/2 presented in ref. [17]. Since only the phase evolution of M 0 is linked to that of the scattering amplitude (by virtue of Watson's theorem [18]), its modulus is parameterised with a fourth-order polynomial whose coefficients are determined in the final fit to data. Details of this parameterisation can be found in appendix B. The normalisation condition for the mass-dependent terms is where Φ 4 is the four-body phase-space factor. The phase of H j 1 j 2 h (m 1 , m 2 ) is set to 0 at m 1 = m 2 = M (K * 0 ), where M (K * 0 ) is the mass of the K * 0 state [15], in order to normalise the relative global phases of the Kπ mass-dependent amplitudes. The CP -violating effects are assumed to be the same for all of the modes under study. Consequently, the value of φ dd s and |λ| determined in this article is effectively an average over the various channels considered in table 1. Within this approach, the physical amplitudes A j 1 j 2 h and A j 1 j 2 h in eq. (2.4) can be separated into a CP -averaged complex amplitude, A j 1 j 2 h , a direct CP 3 and a CP -violating weak phase in the decay, φ D , as (2.7) In the expressions above the CP transformation also changes j 1 j 2 to j 2 j 1 . The total CP -violating phase associated to the interference between mixing and decay is given by φ dd s = φ M − 2φ D and its determination is the main goal of this analysis. In the SM the size of φ dd s is expected to be small due to an almost exact cancellation in the values of φ M and 2φ D [5]. The parameter |λ| is defined in terms of the direct CP asymmetry by (2.8)

Detector and simulation
The LHCb detector [19,20] is a single-arm forward spectrometer covering the pseudorapidity range between 2 and 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with 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 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 (RICH) detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger, 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. The software trigger requires a two-, three-or four-track secondary vertex with significant displacement from the primary pp interaction vertices. At least one charged particle must have transverse momentum p T > 1.7 GeV/c and be inconsistent with originating from a PV. A multivariate algorithm [21] is used for the identification of 3 The direct CP asymmetry is often notated elsewhere as ACP . secondary vertices consistent with the decay of a b hadron. Simulated samples of resonant B 0 s → K * 0 K * 0 , B 0 s → K * 0 K * 0 (1430) 0 and B 0 s → K * 0 (1430) 0 K * 0 (1430) 0 decays, as well as phase-space B 0 s → K + π − K − π + decays, are used to study the signal. Simulated samples are created to study peaking backgrounds. In the simulation, pp collisions are generated using Pythia [22] with a specific LHCb configuration [23]. Decays of particles are described by EvtGen [24], in which final-state radiation is generated using Photos [25]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [26,27] as described in ref. [28].

Signal candidate selection
Events passing the trigger are required to satisfy requirements on the fit quality of the B 0 s decay vertex as well as the p T and χ 2 IP of each track, where χ 2 IP is defined as the difference between the χ 2 of the secondary vertex reconstructed with and without the track under consideration. The tracks are assigned as kaon or pion candidates using particle identification information from the RICH detectors by requiring that the likelihood for the kaon hypothesis is larger than that for the pion hypothesis and vice versa. In addition, the p T of each Kπ pair is required to be larger than 500 MeV/c, the reconstructed mass of each Kπ pair is required to be within the range 750 ≤ m(Kπ) ≤ 1600 MeV/c 2 and the reconstructed mass of the B 0 s candidate is required to be within the range 5000 ≤ m(K + π − K − π + ) ≤ 5800 MeV/c 2 . A boosted decision tree (BDT) algorithm [29,30] is trained to reject combinatorial background, where at least one of the final-state tracks originates from a different decay or directly from the PV. The signal is represented in the BDT training with simulated B 0 s → K * 0 K * 0 candidates, satisfying the same requirements as the data, while selected data candidates in the four-body invariant mass sideband, 5600 ≤ m(K + π − K − π + ) ≤ 5800 MeV/c 2 , are used to represent the background. The input variables employed in the training are kinematic and geometric quantities associated with the four final-state tracks, the two Kπ candidates and the B 0 s candidate. The features used to train the BDT response are chosen to minimise any correlation with the B 0 s and two Kπ pair invariant masses. Separate trainings are performed for the data samples collected in 2011 and 2012, due to the different data-taking conditions. The k-fold cross-validation method [31], with k = 4, is used to increase the training statistics while reducing the risk of overtraining. The requirement on the BDT response is optimized by maximising the metric where N S is the estimated number of signal candidates after selection and N B is the estimated number of combinatorial background candidates within ±60 MeV/c 2 of the known B 0 s mass [15]. The BDT requirement is 95% efficient for simulated signal candidates and rejects 70% of the combinatorial background. After applying the BDT requirement, specific background contributions containing two real oppositely charged kaons and two real oppositely charged pions are removed by mass vetoes on the two-and three-body invariant masses. Candidates are removed if they fulfill either m(K + K − π ± ) < 2100 MeV/c 2 or m(K + K − ) within 30 MeV/c 2 of the known D 0 mass [15]. Sources of peaking background in which one of the final-state tracks is misidentified are suppressed by introducing further -7 -JHEP03(2018)140 particle identification requirements. The particle identification quantities make use of information from the RICH detectors and are calibrated using D * + → D 0 π + and Λ 0 b → Λ + c µ − ν µ decays in data. These requirements significantly reduce contributions from B 0 → ρ 0 K * 0 , B 0 → φK * 0 and Λ 0 b → pπ − K − π + in which a pion or proton is misidentified as a kaon, or a kaon is misidentified as a pion. In addition, there are specific extra particle identification requirements for candidates whose reconstructed mass falls within ±30 MeV/c 2 of the known B 0 or Λ 0 b mass under the relevant mass-hypothesis change (K → p, K → π or π → K). These requirements remove 40% of the simulated signal but almost all of the simulated background: 80% of B 0 → φK * 0 , 96% of B 0 → ρK * 0 and 88% of Λ 0 b → pπ − K − π + events. Subsequently, each of these background components is found to have a small effect on the signal determination. After all of the selection criteria have been imposed, 1.4% of selected events contain multiple candidates, from which one is randomly selected. A fit to the four-body invariant mass distribution is performed in order to determine a set of signal weights, obtained using the sPlot procedure [32], which allows the decay-timedependent CP fit to be performed on a sample that represents only the signal. For the invariant mass fit the B 0 s → (K + π − )(K − π + ) signal and the peaking background com- are modelled as Ipatia functions [33] in which the tail parameters are fixed to values obtained from fits to the simulated samples. The mass difference between the B 0 s and B 0 mesons is fixed to its known value whilst the mean of the B 0 s component, as well as the width of both the B 0 and B 0 s components, are allowed to vary freely. The yields of the B 0 components are allowed to freely vary, whilst the yields of the other components are Gaussian constrained to values relative to the known B 0 → (K + K − )(K − π + ) branching fraction taking into account the relevant production fractions [34] and reconstruction efficiencies. There is an additional background contribution in the low-mass region from partially reconstructed b-hadron decays in which a pion is missed in the final state. This component is modelled as an ARGUS function [35] convolved with a Gaussian mass resolution function. The ARGUS cutoff parameter is fixed to the fitted B 0 s mass minus the neutral pion mass, with the other parameters and yield allowed to vary. The combinatorial background is modelled as an exponential function whose shape parameter and yield are allowed to vary. The result of the four-body invariant mass fit, which is used to obtain the sPlot signal weights, is shown in figure 3. The two Kπ pair invariant masses, with the signal weights applied, are shown in figure 4. The resulting yields of the various fit components are shown in table 2.

Flavour tagging
At the LHC, b quarks are predominantly produced in bb pairs. This analysis focuses on events where one of the quarks hadronises to produce the B 0 s meson while the other quark hadronises and decays independently. Taking advantage of this effect, two types of tagging algorithms aimed at identifying the b-quark flavour at production time are used in this analysis: same-side (SS) taggers, based on information from accompanying particles associated with the signal B 0 s hadronisation process; and opposite-side (OS) taggers, based Partially reconstructed 2580 ± 151 0 Combinatorial 2810 ± 214 372 Table 2. Yields of the signal decay and the various background components considered in the four-body invariant mass fit. The uncertainties are statistical only. The signal region is defined as ±60 MeV/c 2 from the known B 0 s meson mass [15].  on particles produced in the decay of the other b quark. This analysis uses the neuralnetwork-based SS-kaon tagging algorithm presented in ref. [36]; and the combination of OS tagging algorithms explained in ref. [37], based on information from b-hadron decays to electrons, muons or kaons and the total charge of tracks that form a vertex. Both the SS and OS tagging algorithms provide for each event a tagging decision, q, and an estimated mistag probability, η tag . The tagging decision takes the value 1 for B 0 s , −1 for B 0 s and 0 for untagged. To obtain the calibrated mistag probability for a B 0 s (B 0 s ) meson, ω (ω), the estimated probability is calibrated on several flavour-specific control channels. The following linear functions are used in the calibration , where X ∈ {OS, SS}, η X tag is the mean η X tag of the sample, p X 0,1 correspond to calibration parameters averaged over B 0 s and B 0 s , and ∆p X 0,1 account for B 0 s and B 0 s asymmetries in the calibration. Among other modes, the portability of the SS tagger calibration was checked on B 0 s → φφ decays [36], which are kinematically similar to the considered signal mode. The tagging efficiency, tag , denotes the fraction of candidates with a nonzero tagging decision. The tagging power of the sample, eff = tag (1 − 2 ω ) 2 , characterises the tagging performance. Information from the SS and OS algorithms is combined on a per-event basis (see eq. (7.3)) for the decay-time-dependent amplitude fit discussed in section 7. The overall effective tagging power is found to be (5.15 ± 0.14)%. The flavourtagging performance is shown in table 3. When separating the B 0 s and B 0 s components at t = 0, the value of the production asymmetry )) is the production cross-section for the B 0 s (B 0 s ) meson, also has to be incorporated in the model. This asymmetry was measured by LHCb in pp collisions at √ s = 7 TeV by means of a decay-time-dependent analysis of B 0 s → D − s π + decays [38]. To correct for the different kinematics of B 0 s → D − s π + and B 0 s → (K + π − )(K − π + ) decays, a weighting in bins of B 0 s transverse momentum and pseudorapidity is performed, yielding a value of A p = −0.005 ± 0.019. No detection asymmetry need be considered in this analysis since the final state under consideration is charge symmetric.

Acceptance and resolution effects
The LHCb geometrical coverage and selection procedure induce acceptance effects that depend on the three decay angles, the Kπ two-body invariant masses and the decay time. In addition, imperfect reconstruction gives rise to resolution effects. Any deviations caused by imperfect angular and mass resolution are small and are accounted for within the eval-   Table 3. The flavour-tagging performance of the SS and OS tagging algorithms, as well as the combination of both, for the signal data sample used in the analysis. The quoted uncertainty includes both statistical and systematic contributions.
uation of systematic uncertainties (see section 8). However, knowledge of the decay-time resolution is of key importance in the determination of φ dd s and is consequently included in the decay-time-dependent fit. In this analysis, both acceptance and resolution effects are studied using samples of simulated events which have been weighted to match the data distributions in several important kinematic variables. In the description of the acceptance, the decay-time-dependent part is factorised with respect to the part that depends on the kinematic quantities, since they are found to be only 5% correlated. The acceptance and the decay-time resolutions are determined from simulated events that contain an appropriate combination of the vector-vector B 0 s → K * 0 K * 0 component with a sample of B 0 s → K + π − K − π + decays generated according to a phase-space distribution. This combination sufficiently populates the phase-space regions to represent the signal decay. To obtain the acceptance function, the simulated events are weighted by the inverse of the probability density function (PDF) used for generation (defined in terms of angles, masses and decay time). The decay-time acceptance is treated analytically and parameterised using cubic spline functions, following the procedure outlined in ref. [39], with the number of knots chosen to be six. The effect of this choice is addressed as a systematic uncertainty in section 8. The decay-time acceptance is shown in figure 5 (bottom right). The fivedimensional kinematic acceptance in angles and masses is included by using normalisation weights in the denominator of the PDF used in the fit to the data, following the procedure described in ref. [40]. When visualising the fit results (see figure 7), the simulated events are weighted using the matrix element of the amplitude fit model. For illustrative purposes, some projections of the kinematic acceptance are shown in figure 5. In order to obtain the best possible sensitivity for the measurement of the φ dd s phase, the time resolution is evaluated event by event, using the estimated decay-time uncertainty, δ t , obtained in the track reconstruction process. This variable is calibrated using the simulation sample described above to provide the per-event decay-time resolution, σ t , using a linear relationship where δ t is the mean δ t of the sample and p σt (0,1) are the calibration parameters. During fitting, σ t is taken to be the width of a Gaussian resolution function which convolves the decay-time-dependent part of the total amplitude model. Figure 6 shows the relationship between the estimated decay-time uncertainty, δ t , and the calibrated per-event decay-time resolution, σ t .   s → K * 0 K * 0 and pure phase-space B 0 s → (K + π − )(K − π + ) candidates scaled by the mean acceptance. In the bottom right plot the decay-time acceptance obtained from the simulated sample is shown as the black points and the parametric form of the acceptance obtained with cubic splines is shown as the red curve. In the other three plots the black points show the acceptance distribution for the masses and angles. The two cos θ variables and the two m(Kπ) masses have been averaged for the purpose of illustration. In the fit, the kinematic acceptance enters via the normalisation weights.  The model used to fit the data is built by taking the squared moduli of the amplitudes f |B 0 s (t) and f |B 0 s (t) introduced in section 2, multiplying them by the four-body phasespace factor, incorporating the relevant flavour-tagging and production-asymmetry parameters, and including the acceptance and resolution factors obtained in section 6. The observables η SS tag , η OS tag and δ t (introduced in section 5 and section 6) are treated as conditional variables. The effective 4 normalised PDF can be written as where the subscript α (β) represents the state labels {j 1 , j 2 , h} ({j 1 , j 2 , h }), K αβ (t) parameterises the decay-time dependence and is defined in eq. (7.2), and F αβ (Ω) are terms that parameterise the angular and mass dependence. Both the numerator and the denominator of eq. (7.1) are constructed as a sum over 190 real terms, which arise when squaring the amplitudes decomposed in the combination of the nineteen contributing polarisation states. The decay-time-dependent factors are constructed as where R(t, δ t ) is the decay-time resolution function and the factors ζ ± contain the flavourtagging and production-asymmetry information. These factors are where with X ∈ {OS, SS}. The complex quantities a αβ , b αβ , c αβ and d αβ are defined in terms of the CP -averaged amplitudes, the CP -violating parameters and the η j 1 j 2 h factors, as where the bars on the amplitude indices α and β denote the CP transformation of the considered final state, i.e. the change of quantum numbers j 1 j 2 → j 2 j 1 . The functions K untag αβ are obtained by summing K αβ over the tagging decisions. The angular-and massdependent terms are constructed as where δ αβ is the Kronecker delta and the other terms have been introduced in section 2. The decay-time acceptance function, t (t), and the normalisation weights, ξ αβ , are included in the denominator of eq. (7.1). The normalisation weights correspond to angular and mass integrals that involve the five-dimensional kinematic acceptance, Ω (Ω), and are obtained by summing over the events in the simulated sample where G(Ω) is the model used for generation. The CP -conserving amplitudes, A j 1 ,j 2 h , the direct CP -asymmetry parameter, |λ|, and the mixing induced CP -violating phase, φ dd s , are allowed to vary during the fit. Gaussian constraints are applied to ∆m s , Γ s and ∆Γ s from their known values [8], and to the flavour-tagging and decay-time resolution calibration parameters, introduced in section 5 and section 6. The CP -averaged amplitudes are characterised in the fit by wave fractions, f w , polarisation fractions, f w h , and strong phases, δ w h , given by

JHEP03(2018)140
so not all the fractions are independent of each other, for example The phase of the longitudinal polarisation amplitude of the vector-vector component is set to zero to serve as a reference.

Systematic uncertainties
The decay-time-dependent amplitude model and the fit procedure are cross-checked in several independent ways: using purely simulated decays, fitting in a narrow window around the dominant K * 0 resonance, fitting only in the high-mass region above the K * 0 resonance, considering higher-spin contributions (whose effect is found to be negligible), ensuring that there is no bias when repeating the fit procedure on ensembles of pseudo-experiments and by repeating the fit on subsamples of the data set split by the year of data taking, the magnet polarity and using a different mass range. These checks give compatible results. Several sources of systematic uncertainty are considered for each of the physical observables extracted in the decay-time-dependent fit. These are described in this section. A summary of the systematic uncertainties is given in table 4.

Fit to the four-body invariant mass distribution
The uncertainty on the yield of each of the partially reconstructed components used in the four-body invariant mass fit is propagated to the decay-time-dependent amplitude fit by recalculating the sPlot signal weights after varying each of the yields by one standard deviation. Sources of systematic uncertainty which arise from mismodelling the shapes of both the background and signal components are calculated by performing the full fit procedure using alternative parameterisations. The signal is replaced with a double-sided Crystal Ball function [41] instead of the nominal Ipatia shape described in section 4 and the combinatorial-background shape is replaced with a first-order polynomial instead of the nominal exponential function.

Weights derived from the sPlot procedure
The sPlot procedure assumes that there is no correlation between the fit variable used to determine the weights, in this case the four-body invariant mass, m(K + π − K − π + ), and the projected variables in which the signal distribution is unfolded, in this case the three angles and two masses, Ω. This is checked to be valid to a close approximation for signal decays. In order to assess the impact of any residual correlations in the signal weights, the fourbody mass fit is performed by splitting the data into different bins of cos θ for each (Kπ) pair. For each subcategory the four-body fit is repeated and the resulting model is used to compute a new set of signal weights for the full sample. The largest difference between each subcategory value and the nominal fit value is taken as the systematic uncertainty.

Decay-time-dependent fit procedure
An ensemble of pseudoexperiments is generated to estimate the bias on the parameters of the decay-time-dependent fit. For each experiment, a sample with a similar size to the selected signal is generated using the matrix element of the nominal model (employing -15 -JHEP03(2018)140 the measured amplitudes) and then refitted to determine the deviation induced in the fit parameters. The systematic uncertainty is calculated as the mean of the deviation over the ensemble.

Decay-time-dependent fit parameterisation
Several sources of systematic uncertainty originating from the decay-time-dependent fit model have been studied. These include the parameterisations of the angular momentum centrifugal-barrier factors, the mean and width of the Breit-Wigner functions and the model for the S-wave propagator. An alternative model-independent approach is used, as described in appendix B. The systematic uncertainties are obtained for each of these cases by comparing the fitted parameter values of the alternative model with the fitted values from the nominal model. Additional contributions from higher mass (Kπ) vector resonances, namely the K * 1 (1410) 0 and the K * 1 (1680) 0 states, are also considered. In this case, the size of these components is first estimated on data through a simplified fit. Afterwards, an ensemble of pseudoexperiments is generated including these resonances in the model and then refitting with the nominal PDF. The total systematic uncertainty for the decay-time-dependent fit model is taken as the sum in quadrature of these alternatives.

Acceptance normalisation weights
The kinematic acceptance weights, explained in section 7, are computed from simulated samples of limited size, which induces an uncertainty. This systematic uncertainty is calculated using an ensemble of pseudoexperiments in which the acceptance weights are randomly varied according to their covariance matrix (evaluated on the simulated sample). The root-mean-square of the distribution of the differences between the nominal fitted value and the value obtained in each pseudoexperiment is taken as the size of the systematic uncertainty. This effect is found to be the largest systematic uncertainty impacting the measurement of the φ dd s phase.

Other acceptance and resolution effects
Various other acceptance and resolution effects for the decay angles, the two Kπ pair masses and the decay-time are accounted for. Most of these quantities are nominally computed in the decay-time-dependent fit using simulation samples. Any differences between data and simulation are accounted for by the systematic uncertainties described in this section. Furthermore, various other effects originating from mismodelling of the decay-time acceptance and decay-time resolution functions are considered. Each of these effects are summed in quadrature to provide the value listed in table 4. The kinematic and decay-time acceptances, shown in figure 5, are computed from samples of simulated signal events. Small systematic effects can arise due to differences between the data and the simulated samples.
In particular, mismodelling of the B 0 s and the four-track momentum distributions can impact the acceptance in cos θ. This effect is checked by producing a data-driven correction for the simulation in several relevant physical quantities. 5 This correction is produced using an 5 The variables used to correct the distributions of the simulation are the momentum and pseudorapidity of the kaons and pions, the transverse momentum of the B 0 s and the number of tracks in the event.

JHEP03(2018)140
iterative procedure that removes any effects arising from differences between the model used in the event generation and the actual decay kinematics of B 0 s → (K + π − )(K − π + ) decays. The systematic uncertainty is computed as the difference in the fit parameters before and after the iterative correction has been applied. Systematic effects due to the possible mismodelling of the decay-time-dependent acceptance are studied by generating ensembles of pseudoexperiments in two different configurations: one in which the decay-time acceptance spline coefficients are randomised and one in which the configuration of the decay-time acceptance knots is varied. The nominal decay-time-dependent fit procedure is repeated for each pseudoexperiment and the systematic uncertainty for each of these two effects is computed as the average deviation of the fit parameters from their generated values over each ensemble. Sources of systematic uncertainty which affect the decay-time resolution are studied by modifying the calibration function in eq. (6.1) that is used to obtain the perevent decay-time resolution. First, the nominal function is substituted by an alternative quadratic form, to asses the effect of nonlinearity in the calibration. Second, the nominal function is multiplied by a scale factor that accounts for possible remaining differences between data and simulation. This scale factor is taken from the analysis of B 0 s → J/ψφ decays performed by LHCb in ref. [42]. In the both cases, the systematic uncertainties are obtained by comparing the values resulting from the alternative configurations with the nominal values. The effect of the resolution on the masses and angles is studied by generating ensembles of pseudoexperiments for which the masses and angles are smeared using a multi-dimensional Gaussian resolution function, obtained from simulation. The systematic uncertainty is computed as the mean deviation between the fitted and generated values.

Production asymmetry
The uncertainty of the production asymmetry for the B 0 s meson is studied by computing the maximum difference between the nominal conditions and when the production asymmetry is shifted to ±1σ of its nominal value.

Fit results
An unbinned maximum likelihood fit is applied to the background-subtracted data using the PDF defined in eq. (7.1). The large computational load due to the complexity of the fit motivates the parallelisation of the process on a Graphics Processing Unit (GPU), for which the Ipanema software package [43,44]  |λ| Parameter Parameter  Table 4. Summary of the systematic uncertainties on the two CP parameters, the CP -averaged fractions and the strong phase differences (in radians) for each of the components listed in table 1.

Summary
A flavour-tagged decay-time-dependent amplitude analysis of the B 0 s → (K + π − )(K − π + ) decay, for (K ± π ∓ ) invariant masses in the range from 750 to 1600 MeV/c 2 , is performed on a data set corresponding to an integrated luminosity of 3.0 fb −1 obtained by the LHCb experiment with pp collisions at √ s = 7 TeV and √ s = 8 TeV. Several quasi-two-body decay components are considered, corresponding to (K ± π ∓ ) combinations with spins of 0, 1 and 2. The longitudinal polarisation fraction for the B 0 s → K * 0 K * 0 vector-vector decay is determined to be f V V L = 0.208 ± 0.032 ± 0.046, where the first uncertainty is statistical and the second one systematic. This confirms, with improved precision, the relatively low value reported previously by LHCb [14]. The first determination of the CP asymmetry of the (K + π − )(K − π + ) final state and the best, sometimes the first, measurements of 19 CP -averaged amplitude parameters corresponding to scalar, vector and tensor final states, are also reported. This analysis determines for the first time the mixing-induced CPviolating phase φ s using a b → dds transition. The value of this phase is measured to be φ dd s = −0.10 ± 0.13 ± 0.14 rad, which is consistent with both the SM expectation [7] and the corresponding LHCb result of φ ss s = −0.17 ± 0.15 ± 0.03 rad measured using B 0 s → φφ decays [12]. The statistical uncertainty of the two measurements is at a similar level although the systematic uncertainty of this measurement is larger, which is mainly due to the treatment of the multi-dimensional acceptance. It is expected that this can be reduced by increasing the size of the simulation sample used to determine the acceptance effects. Most other sources of systematic uncertainty are expected to scale with larger data samples.

A Angular distributions
The angular dependence of the decay amplitudes introduced in eq. (2.4) is shown in table 6.

B Scalar Kπ mass-dependent amplitude
The variation of the phase with m(Kπ) in the nominal model used for the scalar Kπ massdependent amplitude is taken from ref. [17]. The modulus line-shape is parameterised with  Table 7. Parameters used in the nominal model for the scalar Kπ mass-dependent amplitude. The correlations among them are found to be small, the largest ones been of the order of 50%.  Table 8. Coefficients used in the model-independent parameterisation of the scalar Kπ massdependent amplitude.
The coefficients measured in the decay-time-dependent fit for this case are given in table 8.
The line-shapes of the two scalar mass amplitude models are shown in figure 8. Both approaches are found to be qualitatively compatible with each other.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.