Amplitude analysis of the $B_{(s)} \to K^{*0} \overline{K}^{*0}$ decays and measurement of the branching fraction of the $B \to K^{*0} \overline{K}^{*0}$ decay

The $B^0 \to K^{*0} \overline{K}^{*0}$ and $B^0_s \to K^{*0} \overline{K}^{*0}$ decays are studied using proton-proton collision data corresponding to an integrated luminosity of 3fb$^{-1}$. An untagged and time-integrated amplitude analysis of $B^0_{(s)} \to (K^+\pi^-)(K^-\pi^+) $ decays in two-body invariant mass regions of 150 MeV$/c^2$ around the $K^{*0}$ mass is performed. A stronger longitudinal polarisation fraction in the ${B^0 \to K^{*0} \overline{K}^{*0}}$ decay, ${f_L = 0.724 \pm 0.051 \,({\rm stat}) \pm 0.016 \,({\rm syst})}$, is observed as compared to ${f_L = 0.240 \pm 0.031 \,({\rm stat}) \pm 0.025 \,({\rm syst})}$ in the ${B^0_s\to K^{*0} \overline{K}^{*0}}$ decay. The ratio of branching fractions of the two decays is measured and used to determine $\mathcal{B}(B^0 \to K^{*0} \overline{K}^{*0}) = (8.0 \pm 0.9 \,({\rm stat}) \pm 0.4 \,({\rm syst})) \times 10^{-7}$.


Introduction
The B 0 → K * 0 K * 0 decay is a Flavour-Changing Neutral Current (FCNC) process. 1 In the Standard Model (SM) this type of processes is forbidden at tree level and occurs at first order through loop penguin diagrams. Hence, FCNC processes are considered to be excellent probes for physics beyond the SM, since contributions mediated by heavy particles, contemplated in these theories, may produce effects measurable with the current sensitivity.
Evidence of the B 0 → K * 0 K * 0 decay has been found by the BaBar collaboration [1] with a measured yield of 33.5 +9. 1 −8.1 decays. An untagged time-integrated analysis was presented finding a branching fraction of B = (1.28 +0. 35 −0.30 ±0.11)×10 −6 and a longitudinal polarisation fraction of f L = 0.80 +0.11 −0.12 ± 0.06. In untagged time-integrated analyses the distributions for B 0 and B 0 decays are assumed to be identical and summed, so that they can be fitted with a single amplitude. However, if CP -violation effects are present, the distribution is given by the incoherent sum of the two contributions. The Belle collaboration also searched for this decay [2] and a branching fraction of B = (0.26 +0.33+0.10 −0.29−0.07 ) × 10 −6 was measured, disregarding S-wave contributions. There is a 2.  [3]. Perturbative QCD predicts B = (0.64 +0. 24 −0.23 ) × 10 −6 [4]. 2 These theoretical predictions agree with the experimental results within the large uncertainties. The measurement of f L agrees with the naïve hypothesis, based on the quark helicity conservation and the V −A nature of the weak interaction, that charmless decays into pairs of vector mesons (V V ) should be strongly longitudinally polarised. See, for example, the Polarization in B Decays review in Ref. [5].
The B 0 s → K * 0 K * 0 decay was first observed by the LHCb experiment with early LHC data [4]. A later untagged time-integrated study, with data corresponding to 1 fb −1 of integrated luminosity, measured B = (10.8 ± 2.1 ± 1.5) × 10 −6 and f L = 0.201 ± 0.057 ± 0.040 [6]. More recently, a complete CP -sensitive time-dependent analysis of B 0 s → (K + π − )(K − π + ) decays in the (Kπ) mass range from 750 to 1600 MeV/c 2 has been published by LHCb [7], with data corresponding to 3 fb −1 of integrated luminosity. A determination of f L = 0.208 ± 0.032 ± 0.046 was performed as well as the first measurements of the mixing-induced CP -violating phase φ dd s and of the direct CP asymmetry parameter |λ|. These LHCb analyses of B 0 s → (K + π − )(K − π + ) decays lead to three conclusions: firstly, within their uncertainties, the measured observables are compatible with the absence of CP violation; secondly, a low polarisation fraction is found; finally, a large S-wave contribution, as much as 60%, is measured in the 150 MeV/c 2 window around the K * 0 mass. The low longitudinal polarisation fraction shows a tension with the prediction of QCDF (f L = 0.63 +0.42 −0.29 [3]) and disfavours the hypothesis of strongly longitudinally polarised V V decays. Theoretical studies try to explain the small longitudinal polarisation with mechanisms such as contributions from annihilation processes [3,8]. It is intriguing that the two channels B 0 → K * 0 K * 0 and B 0 s → K * 0 K * 0 , which are related Figure 2: Definition of the helicity angles, employed in the angular analysis of the B 0 (s) → K * 0 K * 0 decays. Each angle is defined in the rest frame of the decaying particle.
of the K +(−) meson and the direction opposite to the B-meson momentum in the rest frame of the K * 0 (K * 0 ) resonance, and φ, the angle between the decay planes of the two vector mesons in the B-meson rest frame. From angular momentum conservation, three relative polarisations of the final state are possible for V V final states that correspond to longitudinal (0 or L), or transverse to the direction of motion and parallel ( ) or perpendicular (⊥) to each other. For the two-body invariant mass of the (K + π − ) and (K − π + ) pairs, noted as m 1 ≡ M (K + π − ) and m 2 ≡ M (K − π + ), a range of 150 MeV/c 2 around the known K * 0 mass [5] is considered. Therefore, (Kπ) pairs may not only originate from the spin-1 K * 0 meson, but also from other spin states. This justifies that, besides the helicity angles, a phenomenological description of the two-body invariant mass spectra, employing the isobar model, is adopted in the analytic model. In the isobar approach, the decay amplitude is modelled as a linear superposition of quasi-two-body amplitudes [14].
For the S-wave (J = 0), the K * 0 (1430) 0 resonance, the possible K * 0 (700) 0 (or κ) and a non-resonant component, (Kπ) 0 , need to be accounted for. This is done using the LASS parameterisation [15], which is an effective-range elastic scattering amplitude, interfering with the K * 0 (1430) 0 meson, where represents the K * 0 (1430) 0 width. In Eq. (1) and Eq. (2) q is the (Kπ) centre-of-mass decay momentum, and M 0 , Γ 0 and q 0 are the K * 0 (1430) 0 mass, width and centre-of-mass decay momentum at the pole, respectively. The effective-range elastic scattering amplitude component depends on where a is the scattering length and b the effective range. For the P-wave (J = 1), only the K * (892) 0 resonance is considered. Other P-wave resonances, such as K * (1410) 0 or K * (1680) 0 , with pole masses much above the fit region, are neglected. Resonances with higher spin, for instance the D-wave K * 2 (1430) 0 meson, are negligible in the considered two-body mass range [7] and are also disregarded. The K * 0 amplitude is parameterised with a spin-1 relativistic Breit-Wigner amplitude, The mass-dependent width is given by where M 1 and Γ 1 are the K * 0 mass and width, r is the interaction radius parameterising the centrifugal barrier penetration factor, and q 1 corresponds to the centre-of-mass decay momentum at the resonance pole. The values of the mass propagator parameters are summarised in Table 1.
The differential decay rate for B 0 (s) mesons 3 at production is given by [6,17], where Φ 4 is the four-body phase space factor. The index i runs over the first column of Table 2 where the different decay amplitudes, A i ≡ |A i |e iδ i , and the angular-mass functions, g i , are listed. The angular dependence of these functions is obtained from spherical harmonics as explained in Ref. [17]. For CP -studies, the CP -odd, A + S , and CP -even, A − S , eigenstates of the S-wave polarisation amplitudes are preferred to the vector-scalar (V S) and scalar-vector (SV ) helicity amplitudes, to which they are related by The remaining amplitudes, except for A ⊥ , correspond to CP -even eigenstates. The contributions can be quantified by the terms F ij , defined as which are normalised according to This condition ensures that 6 i=1 |A i | 2 = 1. The polarisation fractions of the V V amplitudes are defined as where A 0 , A and A ⊥ are the longitudinal, parallel and transverse amplitudes of the P-wave. Therefore, f L is the fraction of B 0 (s) → K * 0 K * 0 longitudinally polarised decays. The polarisation fractions are preferred to the amplitude moduli since they are independent of the considered (Kπ) mass range. The P-wave amplitudes moduli can always be recovered as The phase of all propagators is set to be zero at the K * 0 mass. In addition, a global phase can be factorised without affecting the decay rate setting δ 0 ≡ 0. The last two requirements establish the definition of the amplitude phases (δ , δ ⊥ , δ − S , δ + S and δ SS ) as the phase relative to that of the longitudinal P-wave amplitude at the K * 0 mass.
Since B 0 (s) mesons oscillate, the decay rate evolves with time. The time-dependent amplitudes are obtained replacing A i → A i (t) andĀ i →Ā i (t) in Eq. (5) being where In this analysis, no attempt is made to identify the flavour of the initial B 0 (s) meson and time-integrated spectra are considered. Consequently, the selected candidates correspond to untagged and time-integrated decay rates and there is no sensitivity to direct and mixing-induced CP violation. Moreover, since the origin of phases is set in a CP -even eigenstate (δ 0 = 0), for the CP -odd eigenstates, the untagged time-integrated decay is only sensitive to the phase difference δ ⊥ − δ + S . The present experimental knowledge is compatible with small CP violation in mixing [18] and with the absence of direct CP violation in the B 0 s → (K + π − )(K − π + ) system [7]. The dependence of the decay rate in an untagged and time-integrated analysis of a B 0 (s) meson can be expressed as where the A i amplitudes account for the the average of B 0 (s) and B 0 (s) decays and N is a normalisation constant. For the B 0 meson, a further simplification of the decay rate is considered, since ∆Γ/Γ = −0.002 ± 0.010 [18] the light and heavy mass eigenstate widths can be assumed to be equal, and this factor can be extracted as part of the normalisation constant in Eq. (8). For the B 0 s meson the central values Γ H = 0.618 ps −1 and Γ L = 0.708 ps −1 [18] are considered.

Detector and simulation
The LHCb detector [19,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 siliconstrip 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 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. 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.
The magnetic field deflects oppositely charged particles in opposite directions and this can lead to detection asymmetries. Periodically reversing the magnetic field polarity throughout the data-taking almost cancels the effect. The configuration with the magnetic field pointing upwards (downwards), MagUp (MagDown), bends positively (negatively) charged particles in the horizontal plane towards the centre of the LHC ring.
The online event selection is performed by a trigger [21], 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. In the offline selection, trigger signatures are associated with reconstructed particles. Since the trigger system uses the p T of the charged particles, the phase-space and time acceptance is different for events where signal tracks were involved in the trigger decision (called trigger-on-signal or TOS throughout) and those where the trigger decision was made using information from the rest of the event only (noTOS).
Simulated samples of the B 0 → K * 0 K * 0 and B 0 s → K * 0 K * 0 decays with longitudinal polarisation fractions of 0.81 and 0.64, respectively, are primarily employed in these analyses, particularly for the acceptance description as explained in Sect. 6. Simulated samples of the main peaking background contributions, In the simulation, pp collisions are generated using Pythia [22] with a specific LHCb configuration [23]. Decays of hadronic 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 selection
Both data and simulation are filtered with a preliminary selection. Events containing four good quality tracks with p T > 500 MeV/c are retained. In events that contain more than one PV, the B 0 (s) candidate constructed with these four tracks is associated with the PV that has the smallest χ 2 IP , where χ 2 IP is defined as the difference in the vertex-fit χ 2 of the PV reconstructed with and without the track or tracks in question. Each of the four tracks must fulfil χ 2 IP > 9 with respect to the PV and originate from a common vertex of good quality (χ 2 /ndf < 15, where ndf is the number of degrees of freedom of the vertex).
To identify kaons and pions, a selection in the difference of the log-likelihoods of the kaon and pion hypothesis (DLL Kπ ) is applied. This selection is complemented with fiducial constraints that optimise the particle identification determination: the pion and kaon candidates are required to have 3 < p < 100 GeV/c and 1.5 < η < 4.5 and be inconsistent with muon hypothesis. The final state opposite charge (Kπ) pairs are combined into K * 0 and K * 0 candidates with a mass within 150 MeV/c 2 of the K * 0 mass. The K * 0 and K * 0 candidates must have p T > 900 MeV/c and vertex χ 2 /ndf < 9. The intermediate resonances must combine into B 0 (s) candidates within 500 MeV/c 2 of the B 0 s mass, with a distance of closest approach between their trajectories of less than 0.3 mm. To guarantee that the B 0 (s) candidate originates in the interaction point, the cosine of the angle between the B 0 (s) momentum and the direction of flight from the PV to the decay vertex is required to be larger than 0.99 and the χ 2 IP with respect to the PV has to be smaller than 25. A multivariate selection based on a Boosted Decision Tree with Gradient Boost [29] (BDTG) is employed. It relies on the aforementioned variables and on the B 0 (s) candidate flight distance with respect to the PV and its p T . Simulated B 0 → K * 0 K * 0 decays with tracks matched to the generator particles and filtered with the preliminary selection are used as signal sample, whereas the four-body invariant-mass sideband 5600 < M (K + π − K − π + ) < 5800 MeV/c 2 , composed of purely combinatorial (K + π − )(K − π + ) combinations, is used as background sample for the BDTG training. The number of events in the signal training sample of the BDTG is determined using the ratio between the B 0 s and the B 0 yields from Ref. [6] and the B 0 s yield obtained with a four-body mass fit to the data sample after the preliminary selection.
decays, are strongly suppressed by the requirement in the (Kπ) mass. Resonances in three-body combinations (K + K − π + ) and (K + π + π − ) are also explored. In the case of the former, the three-body invariant mass in the data sample is above all known charm resonances. For the latter, no evidence of candidates originated in is found. Three-body combinations with a pion misidentified as a kaon are reconstructed, mainly searching for All of them are suppressed to a negligible level by the applied selection. A search of three-body combinations with a proton misidentified as a kaon is performed, finding no relevant contribution from decays involving a Λ + c baryon. Decays into five final-state particles are also investigated. Contributions of the B 0 → η (γπ + π − )K * 0 decay can be neglected due to the small misidentification probability and the four-body mass distribution whereas the B 0 s → φ(π 0 π + π − )φ(K + K − ) decay is negligible due to the requirement on the (Kπ) mass.

Four-body mass spectrum
The signal and background yields are determined by means of a simultaneous extended maximum-likelihood fit to the invariant-mass spectra of the four final-state particles in the 2011 and 2012 data samples. The B 0 (s) → (K + π − )(K − π + ) signal decays are parameterised with double-sided Hypatia distributions [30] with the same parameters except for their means that are shifted by the difference between the B 0 and B 0 s masses, 87.13 MeV/c 2 [5].
contributions are described with the sum of a Crystal Ball [31] function and a Gaussian distribution which shares mean with the Crystal Ball core. The parameters of these distributions are obtained from simulation, apart from the mean and resolution values which are free to vary in the fit. Whereas the distribution mean values are constrained to be the same in the 2011 and 2012 data, the resolution is allowed to have different values for the two samples. The small contributions from decays have a broad distribution in the four-body mass and are the object of specific treatment. The contribution from B 0 → ρ 0 K * 0 decays has an expected yield of 3.5 ± 1.3 (6.6 ± 2.3) in the 2011 (2012) sample. It is estimated from the detection and selection efficiency measured with simulation, the collected luminosities, the cross section for bb production, the hadronisation fractions of B 0 and B 0 s mesons and the known branching fraction of the mode. Simulated events containing this decay mode are added with negative weights to the final data sample to subtract its contribution. The contribution of Λ 0 b → (pπ − )(K − π + ) decays in the 2011 (2012) sample is determined to be 36 ± 16 (120 ± 28) from a fit to the (pπ − K − π + ) four-body mass spectrum of the selected data. In this study the four-body invariant mass is recomputed assigning the proton mass to the kaon with the largest DLL pK value. In these fits the Λ 0 b component is described with a Gaussian distribution and the dominant B 0 s → (K + π − )(K − π + ) background is described with a Crystal Ball function. The parameters of both lineshapes are obtained from simulation. The remaining contributions, mainly B 0 → (K + π − )(K − K + ) and partially reconstructed events, are parameterised with a decreasing exponential with a free decay constant. The Λ 0 b → (pπ − )(K − π + ) decay angular distribution is currently unknown and its contribution can not be subtracted with negatively weighted simulated events. Its subtraction is commented further below.
Finally, contributions from partially reconstructed b-hadron decays and combinatorial background are also considered. The former is composed of Band B 0 s -meson decays containing neutral particles that are not reconstructed. Because of the missing particle, the measured four-body invariant mass of these candidates lies in the lower sideband of the spectrum. All contributions to this background are jointly parameterised with an ARGUS function [32] convolved with a Gaussian resolution function, with the same width as the signal. The endpoint of the distribution is also fixed to the B 0 s mass minus the π 0 mass. The combinatorial background is composed of charged tracks that are not originating from the signal decay chain. It is modelled with a linear distribution, with a free slope parameter, separate for 2011 and 2012 data samples.
The results of the fit to the four-body mass spectrum are shown in Fig. 3 and the yields are reported in Table 3. In total, about three hundred B 0 → (K + π − )(K − π + ) signal candidates are found, a factor seven larger than previous analyses [1,2]. To perform a background-subtracted amplitude analysis, the sPlot technique [33,34] is applied to isolate , the dotted yellow line to B 0 → (K + π − )(K − K + ) and the dotted cyan line represents the partially reconstructed background. The tiny combinatorial background contribution is not represented. The black points with error bars correspond to data to which the B 0 → ρ 0 K * 0 contribution has been subtracted with negatively weighted simulation, and the overall fit is represented by the thick blue line.
, for which the yield is fixed, is treated using extended weights according to Appendix B.2 of Ref. [33]. The sPlot method suppresses the background contributions using their relative abundance in the four-body invariant mass spectrum and, therefore, no assumption is required for their phase-space distribution.

Amplitude analysis
Each of the background-subtracted samples of B 0 → (K + π − )(K − π + ) and B 0 s → (K + π − )(K − π + ) decays is the object of a separate amplitude analysis based on the model described in Sect. 2. As a first step, the effect of a non-uniform efficiency, depending on the helicity angles and the two-body invariant masses, is examined. For this purpose, four categories are defined according to the hardware trigger decisions (TOS or noTOS) and data-taking period (2011 and 2012). The efficiency is accounted for through the complex integrals [35] where ε is the total phase-space dependent efficiency, k is the sample category and F ij are defined in Eq. (6). The integrals of Eq. (9) are determined using simulated signal samples of each of the four categories, selected with the same criteria applied to data. A single set of integrals is used for both the B 0 s and the B 0 amplitude analyses. A probability density function (PDF) for each category is built where A i and η i are given in Table 2.
Candidates from all categories are processed in a simultaneous unbinned maximumlikelihood fit, separately for each signal decay mode, using the PDFs in Eq. (10). To avoid nonphysical values of the parameters during the minimisation, some of them are redefined as , where x f , x |A + S | 2 and x |A SS | 2 are used in the fit, together with f L , |A − S | 2 ,δ , δ ⊥ − δ + S , δ − S and δ SS . The former three variables are free to vary within the range [0,1], ensuring that the sum of all the squared amplitudes is never greater than 1. The fit results are corrected for a small reducible bias, originated in discrepancies between data and simulation, as explained in Sect. 7. The final results are shown in Table 4. Figs. 4 and 5 show the one-dimensional projections of the amplitude fit to the B 0 → (K + π − )(K − π + ) and B 0 s → (K + π − )(K − π + ) signal samples in which the background is statistically subtracted by means of the sPlot technique. Three contributions are shown: V V , produced by (K + π − ) (K − π + ) pairs originating in a K * 0 K * 0 decay; V S, accounting for amplitudes in which only one of the (Kπ) pairs originates in a K * 0 decay; and SS, where none of the two (Kπ) pairs originate in a K * 0 decay.
The fraction of V V decays, or purity at production, of the B 0 → K * 0 K * 0 signal, f P B 0 , is estimated from the amplitude analysis and found to be .050 (stat) ± 0.017 (syst). Table 4: Results of the amplitude analysis of B 0 → (K + π − )(K − π + ) and B 0 s → (K + π − )(K − π + ) decays. The observables above the line are directly obtained from the maximum-likelihood fit whereas those below are obtained from the former, as explained in the text, with correlations accounted for in their estimated uncertainties. For each result, the first quoted uncertainty is statistical and the second systematic. The estimation of the latter is described in Sect. 7. 0.023 ± 0.014 ± 0.004 0.087 ± 0.011 ± 0.011 S-wave fraction 0.408 ± 0.050 ± 0.017 0.694 ± 0.016 ± 0.010 The significance of this magnitude, computed as its value over the sum in quadrature of the statistical and systematic uncertainty, is found to be 10.8 standard deviations. This significance corresponds to the presence of B 0 → K * 0 K * 0 V V decays in the data sample. The S-wave fraction of the decay is equal to 0.408 = 1 − f P B 0 . For the B 0 s → K * 0 K * 0 mode the S-wave fraction is found to be 0.694 ± 0.016 (stat) ± 0.010 (syst).

Systematic uncertainties of the amplitude analysis
Several sources of systematic uncertainty that affect the results of the amplitude analysis are considered and discussed in the following.
Fit method. Biases induced by the fitting method are evaluated with a large ensemble of pseudoexperiments. For each signal decay, samples with the same yield of signal observed in data (see Table 3) are generated according to the PDF of Eq. (8) with inputs set to the results summarised in Table 4. The use of the weights defined in Eq. (9) to account the detector acceptance would require a full simulation and, instead, a parametric efficiency is considered. For each observable, the mean deviation of the result from the input value is assigned as a systematic uncertainty.
Description of the kinematic acceptance. The uncertainty on the signal efficiency relies on the coefficients of Eq. (9) that are estimated with simulation. To evaluate its impact on the amplitude analysis results, the fit to data is repeated several times with alternative coefficients varied according to their covariance matrix. The standard deviation of the distribution of the fit results for each observable is assigned as a systematic uncertainty.
Resolution. The fit performed assumes a perfect resolution on the phase-space variables.
The impact of the detector resolution on these variables is estimated with sets of pseudoexperiments adding per-event random deviations according to the resolution estimated from simulation. For each observable, the mean deviation of the result from the measured value is assigned as a systematic uncertainty.
P-wave mass model. The amplitude analysis is repeated with alternative values of the parameters that define the P-wave mass propagator, detailed in Table 1, randomly sampled from their known values [5]. The standard deviation of the distribution of the amplitude fit results for each observable is assigned as a systematic uncertainty.
S-wave mass model. In addition to the default S-wave propagator, described in Sect. 2, two alternative models are used: the LASS lineshape with the parameters of Table 5, obtained with B 0 → J/ψK + π − decays within the analysis of Ref. [36], and the propagator proposed in Ref. [37]. The amplitude fit is performed with these two alternatives and, for each observable, the largest deviation from the baseline result  Table 5: Alternative parameters of the LASS mass propagator used in the S-wave systematic uncertainty estimation.
is assigned as a systematic uncertainty.
Differences between data and simulation. An iterative method [38], is used to weight the simulated events and improve the description of the track multiplicity and B 0 (s) -meson momentum distributions. The procedure is repeated multiple times and, for each observable, the mean bias of the amplitude fit result is corrected for in the results of Table 4 while its standard deviation is assigned as a systematic uncertainty.
Background subtraction. The data set used in the amplitude analysis is background subtracted using the sPlot method [33,34] that relies in the lineshapes of the fourbody mass fit discussed in Sect. 5. The uncertainty related to the determination of the signal weights is evaluated repeating the amplitude analysis fits with weights obtained fitting the four-body invariant-mass with two alternative models. In the first case, the model describing the signal is defined by the sum of two Crystal Ball functions [31] with a common, free, peak value and the resolution parameter fixed from simulation. In the second case, the model describing the combinatorial background is assumed to be an exponential function. The amplitude fit is performed with the sPlot-weights obtained with the two alternatives and, for each observable, the largest deviation from the baseline result is assigned as a systematic uncertainty. This procedure is also used when addressing the systematic uncertainties in the measured yields of the different subsamples, as discussed in Sect. 8.
Peaking backgrounds. The uncertainty related to the fluctuations in the yields of the Λ 0 b → (pπ − )(K − π + ) and B 0 → ρ 0 K * 0 background contributions are estimated repeating the amplitude-analysis fit with the yield values varied by their uncertainties reported in Sect. 5. For each observable, the largest deviation from the default result is assigned as a systematic uncertainty. This procedure is also used when addressing the systematic uncertainties of the four-body invariant mass yields in Sect. 8.
Time acceptance. The amplitude analysis does not account for possible decay-time dependency of the efficiency, however, the trigger and the offline selections may have an impact on it. This effect only affects B 0 s -meson decays and is accounted for by estimating effective shifts: Γ H = 0.618 → 0.598 ps −1 and Γ L = 0.708 → 0.732 ps −1 , which are obtained with simulation. For each observable, the variation of the result of the fit after introducing these values in the amplitude analysis is considered as a systematic uncertainty.
The resulting systematic uncertainties and the corrected biases, originated in the differences between data and simulation, are detailed in Table 6 for the parameters of the amplitude-analysis fit. The corresponding values for the derived observables are summarised in Table 7. The total systematic uncertainty is computed as the sum in quadrature of the different contributions.

Determination of the ratio of branching fractions
In this analysis, the B 0 → K * 0 K * 0 branching fraction is measured relative to that of B 0 s → K * 0 K * 0 decays. Since both decays are selected in the same data sample and share a common final state most systematic effects cancel. However, some efficiency corrections, eg. those originated from the difference in phase-space distributions of events of the two modes, need to be accounted for. The amplitude fit provides the relevant information to tackle the differences between the two decays.
This branching-fraction ratio is obtained as : Systematic uncertainties for the parameters of the amplitude-analysis fit of the B 0 (s) → (K + π − )(K − π + ) decay. The bias related to differences between data and simulation is included in the results shown in Table 4. Decay mode   Table 2. Also in this case, for the Table 7: Systematic uncertainties for the derived observables of the amplitude-analysis fit of the B 0 (s) → (K + π − )(K − π + ) decay. The bias related to differences between data and simulation is included in the results shown in Table 4.

Decay mode
The detection efficiency is determined from simulation for each channel separately for the different categories discussed in Sect. 6: year of data taking, trigger type and, in addition, the LHCb magnet polarity. An exception is applied to the particle-identification selection whose efficiency is determined from large control samples of D * + → D 0 π + , D 0 → K − π + decays. Differences in kinematics and detector occupancy between the control samples and the signal data are accounted for in this particle-identification efficiency study [39,40].
The different sources of systematic uncertainty in the branching fraction determination are discussed below.
Systematic uncertainties in the factor κ. The uncertainties on the parameters of the amplitude analysis fit described in Sect. 7 affect the determination of the factors κ defined in Eq. (12) as summarised in Table 8. Table 8: Systematic uncertainties in the factor κ defined in Eq. (12) split in categories. The bias originated in differences between data and simulation is corrected for in the κ results shown in Table 9. Systematic uncertainties in the signal yields. As discussed in Sect. 7 uncertainties on the signal yields arise from the model used to fit the four-body invariant mass. The uncertainties from the different proposed alternative signal and background lineshapes are summed in quadrature to compute the final systematic uncertainty.

Decay mode
Systematic uncertainty in the efficiencies. A dedicated data method is employed to estimate the uncertainty in the signal efficiency originated in the PID selection.
The inputs employed for measuring the relative branching fraction are summarised in Table 9. The factor κ is different for the two decay modes because of two main reasons: firstly, the discrepancy between the polarisation assumed in simulation and its measurement is larger for the B 0 s → K * 0 K * 0 than for the B 0 → K * 0 K * 0 decay. Secondly, the different S-wave fraction of the decays. Also, the efficiency ratio of the two modes deviating from one is explained upon the different polarisation of the simulation samples. The LHCb detector is less efficient for values of cos θ 1 (cos θ 2 ) close to unity because of slow pions emitted in K * 0 (K * 0 ) decays and these are more frequent the larger is the longitudinal polarisation.
The final result of the branching-fraction ratio is obtained as the weighted mean of the per-category result obtained with Eq. (11) for the eight categories of Table 9, and found to be Considering that B(B 0 s → K * 0 K * 0 ) = (1.11 ± 0.22 (stat) ± 0.12 (syst)) × 10 −5 , from Ref. [5], the absolute branching fraction for the B 0 → K * 0 K * 0 mode is found to be B(B 0 → K * 0 K * 0 ) = (8.0 ± 0.9 (stat) ± 0.4 (syst)) × 10 −7 . Table 9: Parameters used to determine B(B 0 → K * 0 K * 0 )/B(B 0 s → K * 0 K * 0 ). When two uncertainties are quoted, the first is statistical and the second systematic. The value of f s /f d is taken from Ref. [41]. It is worth noticing that, since the B 0 s → K * 0 K * 0 branching fraction was determined with the B 0 → K * 0 φ decay as a reference [6], the uncertainty on f s /f d , which appears in the ratio of Eq. (13), does not contribute to the absolute branching fraction measurement.

Summary and final considerations
The first study of B 0 → (K + π − )(K − π + ) decays is performed with a data set recorded by the LHCb detector, corresponding to an integrated luminosity of 3.0 fb −1 at centre-ofmass energies of 7 and 8 TeV. The B 0 → K * 0 K * 0 mode is observed with 10.8 standard deviations. An untagged and time-integrated amplitude analysis is performed, taking into account the three helicity angles and the (K + π − ) and (K − π + ) invariant masses in a 150 MeV/c 2 window around the K * 0 and K * 0 masses. Six contributions are included in the fit: three correspond to the B 0 → K * 0 K * 0 P-wave, and three to the S-wave, along with their interferences. A large longitudinal polarisation of the B 0 → K * 0 K * 0 decay, f L = 0.724 ± 0.051 (stat) ± 0.016 (syst), is measured. The S-wave fraction is found to be 0.408 ± 0.050 (stat) ± 0.023 (syst).
A parallel study of the B 0 s → (K + π − )(K − π + ) mode within 150 MeV/c 2 of the K * 0 mass is performed, superseding a previous LHCb analysis [6]. A small longitudinal polarisation, f L = 0.240 ± 0.031 (stat) ± 0.025 (syst) and a large S-wave contribution of 0.694 ± 0.016 (stat) ± 0.012 (syst) are measured for the B 0 s → K * 0 K * 0 decay, confirming the previous LHCb results of the time-dependent analysis of the same data [7].
This result is inconsistent with the prediction of R sd = 16.4 ± 5.2 [13]. Within models such as QCDF or the soft-collinear effective theory, based on the heavy-quark limit the predictions, longitudinal observables, such as the one in Eq. (14), have reduced theoretical uncertainties as compared to parallel and perpendicular ones. The heavy-quark limit also implies the polarisation hierarchy f L f ,⊥ . The measured value for R sd and the f L result of the B 0 s → K * 0 K * 0 decay put in question this hierarchy. The picture is even more intriguing since, contrary to its U-spin partner, the B 0 → K * 0 K * 0 decay is confirmed to be strongly polarised.