Measurement of ${C\!P}$ violation parameters and polarisation fractions in ${B_s^0\to J/\psi \overline{K}^{*0}}$ decays

The first measurement of ${C\!P}$ asymmetries in the decay ${B_s^0\to J/\psi \overline{K}^{*}(892)^{0}}$ and an updated measurement of its branching fraction and polarisation fractions are presented. The results are obtained using data corresponding to an integrated luminosity of $3.0\,fb^{-1}$ of proton-proton collisions recorded with the LHCb detector at centre-of-mass energies of $7$ and $8\,\mathrm{TeV}$. Together with constraints from ${B^0\to J/\psi \rho^0}$, the results are used to constrain additional contributions due to penguin diagrams in the ${C\!P}$-violating phase ${{\phi}_{s}}$, measured through ${B_s^0}$ decays to charmonium.


Introduction
The CP -violating phase φ s arises in the interference between the amplitudes of B 0 s mesons decaying via b → ccs transitions to CP eigenstates directly and those decaying after oscillation. The phase φ s can be measured using the decay B 0 s → J/ψ φ. Within the Standard Model (SM), and ignoring penguin contributions to the decay, φ s is predicted to be −2β s , with β s ≡ arg(−V cb V * cs /V tb V * ts ), where V ij are elements of the CKM matrix [1]. The phase φ s is a sensitive probe of dynamics beyond the SM (BSM) since it has a very small theoretical uncertainty and BSM processes can contribute to B 0 s -B 0 s mixing [2][3][4][5]. Global fits to experimental data, excluding the direct measurements of φ s , give −2β s = −0.0363 ± 0.0013 rad [6]. The current world average value is φ s = −0.015 ± 0.035 rad [7], dominated by the LHCb measurement reported in Ref. [8]. In the SM expectation of φ s [6], additional contributions to the leading b → ccs tree Feynman diagram, as shown in Fig. 1, are assumed to be negligible. However, the shift in φ s due to these contributions, called hereafter "penguin pollution", is difficult to compute due to the non-perturbative nature of the quantum chromodynamics (QCD) processes involved. This penguin pollution must be measured or limited before using the φ s measurement in searches for BSM effects, since a shift in this phase caused by penguin diagrams is possible. Various methods to address this problem have been proposed [9][10][11][12][13][14], and LHCb has recently published upper limits on the size of the penguin-induced phase shift using B 0 → J/ψ ρ 0 decays [15].
Tree and penguin diagrams contributing to both B 0 s → J/ψ φ and B 0 s → J/ψ K * 0 decays are shown in Fig. 1. In this paper, the penguin pollution in φ s is investigated using B 0 s → J/ψ K * 0 decays 1 , with J/ψ → µ + µ − and K * 0 → K − π + , following the method first proposed in Ref. [9] for the B 0 → J/ψ ρ 0 decay and later also discussed for the B 0 s → J/ψ K * 0 decay in Refs. [11,13]. This approach requires the measurement of the branching fraction, direct CP asymmetries, and polarisation fractions of the B 0 s → J/ψ K * 0 decay. The measurements use data from proton-proton (pp) collisions recorded with the LHCb detector corresponding to 3.0 fb −1 of integrated luminosity, of which 1.0 (2.0) fb −1 was collected in 2011 (2012) at a centre-of-mass energy of 7 (8) TeV. The LHCb collaboration previously reported a measurement of the branching fraction and the polarisation fractions using data corresponding to 0.37 fb −1 of integrated luminosity [16].
The paper is organised as follows: a description of the LHCb detector, reconstruction and simulation software is given in Sect. 2, the selection of the B 0 s → J/ψ K * 0 signal candidates and the B 0 → J/ψ K * 0 control channel are presented in Sect. 3 and the treatment of background in Sect. 4. The J/ψ K − π + invariant mass fit is detailed in Sect. 5. The angular analysis and CP asymmetry measurements, both performed on weighted distributions where the background is statistically subtracted using the sPlot technique [17], are detailed in Sect. 6. The measurement of the branching fraction is explained in Sect. 7. The evaluation of systematic uncertainties is described in Sect. 8 along with the results, and in Sect. 9 constraints on the penguin pollution are evaluated and discussed.

Experimental setup
The LHCb detector [18,19] 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 placed downstream of the magnet. The tracking system provides a measurement of 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, the impact parameter, 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 electromag-netic 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. In this analysis, candidates are first required to pass the hardware trigger, which selects muons with a transverse momentum p T > 1.48 GeV/c in the 7 TeV data or p T > 1.76 GeV/c in the 8 TeV data. In the subsequent software trigger, at least one of the final-state particles is required to have both p T > 0.8 GeV/c and impact parameter larger than 100 µm with respect to all of the primary pp interaction vertices (PVs) in the event. Finally, the tracks of two or more of the final-state particles are required to form a vertex that is significantly displaced from any PV. Further selection requirements are applied offline in order to increase the signal purity.
In the simulation, pp collisions are generated using Pythia [20,21] with a specific LHCb configuration [22]. Decays of hadronic particles are described by EvtGen [23], in which final-state radiation is generated using Photos [24]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [25,26] as described in Ref. [27].

Event selection
The selection of B 0 s → J/ψ K * 0 candidates consists of two steps: a preselection consisting of discrete cuts, followed by a specific requirement on a boosted decision tree with gradient boosting (BDTG) [28,29] to suppress combinatorial background. All charged particles are required to have a transverse momentum in excess of 0.5 GeV/c 2 and to be positively identified as muons, kaons or pions. The tracks are fitted to a common vertex which is required to be of good quality and significantly displaced from any PV in the event. The flight direction can be described as a vector between the B 0 s production and decay vertices; the cosine of the angle between this vector and the B 0 s momentum vector is required to be greater than 0.999. Reconstructed invariant masses of the J/ψ and K * 0 candidates are required to be in the ranges 2947 < m µ + µ − < 3247 MeV/c 2 and 826 < m K − π + < 966 MeV/c 2 . The B 0 s invariant mass is reconstructed by constraining the J/ψ candidate to its nominal mass [30], and is required to be in the range 5150 < m J/ψ K − π + < 5650 MeV/c 2 .
The training of the BDTG is performed independently for 2011 and 2012 data, using information from the B 0 s candidates: time of flight, transverse momentum, impact parameter with respect to the production vertex and χ 2 of the decay vertex fit. The data sample used to train the BDTG uses less stringent particle identification requirements. When training the BDTG, simulated B 0 s → J/ψ K * 0 events are used to represent the signal, while candidates reconstructed from data events with J/ψ K − π + invariant mass above 5401 MeV/c 2 are used to represent the background. The optimal threshold for the BDTG is chosen independently for 2011 and 2012 data and maximises the effective signal yield.

Treatment of peaking backgrounds
After the suppression of most background with particle identification criteria, simulations show residual contributions from the backgrounds The invariant mass distributions of misidentified B 0 → J/ψ π + π − and B 0 s → J/ψ π + π − events peak near the B 0 s → J/ψ K − π + signal peak due to the effect of a wrong-mass hypothesis, and the misidentified B 0 s → J/ψ K + K − candidates are located in the vicinity of the B 0 → J/ψ K + π − signal peak. It is therefore not possible to separate such background from signal using information based solely on the invariant mass of the J/ψ K − π + system. Moreover the shape of the reflected invariant mass distribution is sensitive to the daughter particles momenta. Due to these correlations it is difficult to add the b-hadron to J/ψ h + h − (where h is either a pion, a kaon or a proton) misidentified backgrounds as extra modes to the fit to the invariant mass distribution. Instead, simulated events are added to the data sample with negative weights in order to cancel the contribution from those peaking backgrounds, as done previously in Ref. [8]. Simulated b-hadron to J/ψ h + h − events are generated using a phase-space model, and then weighted on an event-by-event basis using the latest amplitude analyses of the decays , and B 0 → J/ψ π + π − [34]. The sum of weights of each decay mode is normalised such that the injected simulated events cancel out the expected yield in data of the specific background decay mode.
In addition to Λ 0 b → J/ψ pK − and B → J/ψ h + h − decays, background from Λ 0 b → J/ψ pπ − is also expected. However, in Ref. [35] a full amplitude analysis was not performed. For this reason, as well as the fact that the Λ 0 b decays have broad mass distributions, the contribution is explicitly included in the mass fit described in the next section. Expected yields for both B → J/ψ h + h − and Λ 0 b → J/ψ ph − background decays are given in Table 1.
These four modes are statistically disentangled through a fit to the J/ψ K − π + invariant mass. The combinatorial background is described by an exponential distribution, the Λ 0 b → J/ψ pπ − decay by the Amoroso distribution [36] and the B 0 and B 0 s signals by the double-sided Hypatia distribution [37], where K ν (z) is the modified Bessel function of the second kind, K λ (ζ) , and A, B, C, D are obtained by imposing continuity and differentiability. This function is chosen because the event-by-event uncertainty on the mass has a dependence on the particle momenta. The estimate of the number of B 0 → J/ψ K + π − decays lying under the B 0 s peak is very sensitive to the modelling of the tails of the B 0 peak. The fitted fraction is in good agreement with the estimate from simulation.
In the fit to data, the mean and resolution parameters of both the B 0 s and B 0 Hypatia functions are free to vary. All the remaining parameters, namely λ, a 1 , n 1 , a 2 and n 2 , are fixed to values determined from fits to B 0 s and B 0 simulated events. All the Λ 0 b → J/ψ pπ − shape parameters are fixed to values obtained from fits to simulated Λ 0 b → J/ψ pπ − events, while the exponent of the combinatorial background is free to vary.
Due to the small expected yield of Λ 0 b → J/ψ pπ − decays compared to those of the other modes determined in the fit to data, and to the broad distribution of Λ 0 b → J/ψ pπ − decays across the J/ψ K − π + invariant mass spectrum, its yield is included in the fit as a Gaussian constraint using the expected number of events and its uncertainties, as shown in Table 1.
From studies of simulated (MC) samples, it is found that the resolution of B 0 s and B 0 mass peaks depends on both m K − π + and cos(θ µ ), where θ µ is one of the helicity angles used in the angular analysis as defined in Sect. 6. The fit to the J/ψ K − π + invariant mass spectrum, including the evaluation of the sWeights, is performed separately in twenty bins, corresponding to four m K − π + bins of 35 MeV/c 2 width, and five equal bins in cos(θ µ ). The overall B 0 s and B 0 yields are obtained from the sum of yields in the twenty bins, giving where the statistical uncertainties are obtained from the quadratic sum of the uncertainties determined in each of the individual fits. Systematic uncertainties are discussed in Sect. 8. The correlation between the B 0 and B 0 s yields in each bin are found to be smaller than 4%. The ratio of the B 0 s and B 0 yields is found to be N B 0 s /N B 0 = (8.66±0.24(stat) +0.18 −0.16 (syst))× 10 −3 . Figure 2 shows the sum of the fit results for each bin, overlaid on the J/ψ K − π + mass spectrum for the selected data sample.  This analysis uses the decay angles defined in the helicity basis. The helicity angles are denoted by (θ K , θ µ , ϕ h ), as shown in Fig. 3. The polar angle θ K (θ µ ) is the angle between the kaon (µ + ) momentum and the direction opposite to the B 0 s momentum in the K − π + (µ + µ − ) centre-of-mass system. The azimuthal angle between the K − π + and µ + µ − decay planes is ϕ h . The definitions are the same for B 0 s or B 0 s decays. They are also the same for B 0 → J/ψ K * 0 decays.
The shape of the angular distribution of B 0 s → J/ψ K * 0 decays is given by Ref. [38], where λ = 0, ±1 is the J/ψ helicity, α µ = ±1 is the helicity difference between the muons, J is the spin of the K − π + system, H are the helicity amplitudes, and d are the small Wigner matrices. The helicity amplitudes are rotated into transversity amplitudes, which correspond to final P eigenstates, The distribution in Eq. 4 can be written as the sum of ten angular terms, four corresponding to the square of the transversity amplitude of each final state polarisation, and six corresponding to the cross terms describing interference among the final polarisations. The modulus of a given transversity amplitude, A x , is written as |A x |, and its phase as δ x . The convention δ 0 = 0 is used in this paper. The P-wave polarisation fractions are , with i = 0, , ⊥ and the S-wave fraction is defined as . The distribution of the CP -conjugate decay is obtained by flipping the sign of the interference terms which contain |A ⊥ |. For the CP -conjugate case, the amplitudes are denoted as A i . Each A i and the corresponding A i are related through the CP asymmetries, as described in Sect. 6.3.

Partial-wave interference factors
In the general case, the transversity amplitudes of the angular model depend on the K − π + mass (m K − π + ). This variable is limited to be inside a window of ±70 MeV/c 2 around the K * 0 mass. Figure 4 shows the efficiency-corrected m K − π + spectra for B 0 s and B 0 using the nominal sets of sWeights.
In order to account for the m K − π + dependence while keeping the framework of an angular-only analysis, a fit is performed simultaneously in the same four m K − π + bins defined in Sect. 5. Different values of the parameters |A S | 2 and δ S are allowed for each bin, but the angular distribution still contains mass-dependent terms associated with the interference between partial-waves. If only the S-wave and P-wave are considered, such interference terms correspond to the following complex integrals, where m

L(H)
Kπ is the lower (higher) limit of the bin, ε m (m Kπ ) is the acceptance for a K − π + candidate with mass m Kπ (see Appendix A for a discussion on the angular acceptance), Φ stands for the phase space, and P (S) is the P-wave (S-wave) propagator. The phase space term is computed as where p denotes the K * 0 momentum in the B 0 s rest frame and q refers to the K − momentum in the K * 0 rest frame.
The phase θ SP is included in the definition of δ S but the C SP factors, corresponding to real numbers in the interval [0, 1], have to be computed and input to the angular fit. The contribution of D-wave (J = 2) in the m K − π + range considered is expected to be negligible. Therefore the nominal model only includes S-wave and P-wave. To determine the systematic uncertainty due to possible D-wave contributions, C SD and C PD factors are also computed, using analogous expressions to that given in Eq. 9. The C ij factors are calculated by evaluating numerically the integrals using the propagators outlined below, and are included as fixed parameters in the fit. A systematic uncertainty associated to the different possible choices of the propagator models is afterwards evaluated.  Figure 4: Efficiency corrected m K − π + distribution for B 0 s shown in squares (red) and B 0 shown in circles (black) using sWeights computed from the maximum likelihood fit to the J/ψ K − π + invariant mass spectrum.
The S-wave propagator is constructed using the LASS parametrisation [39], consisting of a linear combination of the K * 0 (1430) 0 resonance with a non-resonant term, coming from elastic scattering. The P-wave is described by a combination of the K * (892) 0 and Table 2: The C SP , C SD and C PD factors calculated in each of the four m K − π + bins around the K * 0 peak. , and the D-wave is assumed to come from the K * 2 (1430) 0 contribution. Relativistic Breit-Wigner functions, multiplied by angular momentum barrier factors, are used to parametrise the different resonances. Table 2 contains the computed C SP , C SD and C PD factors.

CP asymmetries
The direct CP violation asymmetry in the B 0 (s) decay rate to the final state where A (s) i are the transversity amplitudes defined in Sect. 6.1 and the additional index s is used to distinguish the B 0 s and the B 0 -meson. The index i refers to the polarisation of the final state (i = 0, , ⊥, S) and is dropped in the rest of this section, for clarity.
The raw CP asymmetry is expressed in terms of the number of observed candidates by Both asymmetries in Eq. 11 and Eq. 12 are related by [41] where A D (f ) is the detection asymmetry, defined as in Eq. (16), production asymmetry, defined as in Eq. (15), ζ (s) = +1(−1) and κ (s) accounts for the dilution due to B 0 (s) −B 0 (s) oscillations [42]. The κ (s) factor is evaluated by where ε(t) is the time-dependent acceptance function, assumed to be identical for the B 0 s → J/ψ K * 0 and B 0 → J/ψ K * 0 decays. The symbols Γ (s) and ∆m (s) denote the decay width and mass differences between the B 0 (s) mass eigenstates.
The B 0 (s) −B 0 (s) production asymmetry is defined as where σ is the B 0 (s) production cross-section within the LHCb acceptance. The production asymmetries reported in Ref. [43] are reweighted in bins of B 0 (s) transverse momentum to obtain The κ (s) factor in Eq. 14 is determined by fixing ∆Γ (s) , ∆m (s) and Γ (s) to their world average values [30] and by fitting the decay time acceptance ε(t) to the nominal data sample after applying the B 0 sWeights, in a similar way to Ref. [44]. It is equal to 0.06% for B 0 s decays, and 41% for B 0 . This reduces the effect of the production asymmetries to the level of 10 −5 for B 0 s → J/ψ K * 0 and 10 −3 for B 0 → J/ψ K * 0 decays. Other sources of asymmetries arise from the different final-state particle interactions with the detector, event reconstruction and detector acceptance. The detection asymmetry, A D (f ), is defined in terms of the detection efficiencies of the final states, ε det , as The detection asymmetry, measured in bins of the K + momentum in Ref.
[45], is weighted with the momentum distribution of the kaon from the B 0 (s) → J/ψ K * 0 (K * 0 ) decays to give The branching fraction B(B 0 s → J/ψ K * 0 ) is obtained by normalising to two different channels, B 0 s → J/ψ φ and B 0 → J/ψ K * 0 , and then averaging the results. The expression is used for the normalisation to a given B q → J/ψ X decay, where N refers to the yield of the given decay, ε corresponds to the total (reconstruction, trigger and selection) efficiency, The event selection of B 0 s → J/ψ φ candidates consists of the same requirements as those for B 0 s → J/ψ K * 0 candidates (see Sect. 3), with the exception that φ candidates are reconstructed in the K + K − state so there are no pions among the final state particles. In addition to the other requirements, reconstructed φ candidates are required to have mass in the range 1000 < m K − K + < 1040 MeV/c 2 and to have a transverse momentum in excess of 1 GeV/c 2 .

Efficiencies obtained in simulation
A first estimate of the efficiency ratios is taken from simulated events, where the particle identification variables are calibrated using D * ± decays. The efficiency ratios estimated from simulation, for 2011 (2012) data, are

Correction factors for yields and efficiencies
The signal and normalisation channel yields obtained from a mass fit are affected by the presence of a non-resonant S-wave background as well as interference between S-wave and P-wave components. Such interference would integrate to zero for a flat angular acceptance, but not for experimental data that are subject to an angle-dependent acceptance. In addition, the efficiencies determined in simulation correspond to events generated with an angular distribution different from that in data; therefore the angular integrated efficiency also needs to be modified with respect to simulation estimates. These effects are taken into account using a correction factor ω, which is the product of the correction factor to the angular-integrated efficiency and the correction factor to the P-wave yield: where N B 0 s →J/ψ K * 0 , N Bq→J/ψ X are the yields obtained from the mass fits, ε MC Bq→J/ψ X , ε MC B 0 s →J/ψ K * 0 are the efficiencies obtained in simulation, and ω is calculated as where F X Bq→J/ψ X is the fraction of the P-wave X resonance in a given B q → J/ψ X decay (related to the presence of S-wave and its interference with the P-wave), and c Bq→J/ψ X is a correction to ε MC Bq→J/ψ X due to the fact that the simulated values of the decay parameters differ slightly from those measured. The values obtained for the ω correction factors are The study of penguin pollution requires the calculation of ratios of absolute amplitudes between B 0 s → J/ψ K * 0 and B 0 is very useful. This normalisation is given by where B(K * 0 → K − π + ) = 2/3 and B(φ → K + K − ) = (49.5 ± 0.5)% [30]. Using N B 0 s →J/ψ K − π + as given in Eq. 3, and N B 0 s →J/ψ K + K − = 58 091 ± 243 (stat) ± 319 (syst) as obtained from a fit to the invariant mass of selected B 0 s → J/ψ φ candidates, where the signal is described by a double-sided Hypatia distribution and the combinatorial background is described by an exponential distribution, a value of The normalisation to B 0 → J/ψ K * 0 is given by where N B 0 →J/ψ K + π − and N B 0 s →J/ψ K − π + are given in Eq. 2 and Eq. 3, respectively, and resulting in a value of where the third uncertainty comes from the hadronisation fraction ratio f d /f s = 3.86 ± 0.22 [7]. 8 Results and systematic uncertainties Section 8.1 presents the results of the angular fit as well as the procedure used to estimate the systematic uncertainties, while in Sect. 8.2 the results of the branching fraction measurements and the corresponding estimated systematic uncertainties are discussed.

Angular parameters and CP asymmetries
The results obtained from the angular fit to the B 0 s → J/ψ K * 0 events are given in Table 3  and Table 4 for the P-wave and S-wave parameters, respectively. For comparison, the previous LHCb measurements [16] of f 0 and f were 0.50 ± 0.08 ± 0.02 and 0.19 +0.10 −0.08 ± 0.02, respectively. The angular distribution of the signal and the projection of the fitted distribution are shown in Fig. 5. The statistical-only correlation matrix as obtained from the fit to data is given in Appendix B. The polarisation-dependent CP asymmetries are compatible with zero, as expected in the SM. The polarisation fractions are in good agreement with the previous measurements [16] performed on the same decay mode by the LHCb collaboration using data corresponding to an integrated luminosity of 0.37 fb −1 .
Various sources of systematic uncertainties on the parameters of the angular fit are studied, as summarised in Table 3 and Table 4 for the P-wave and S-wave parameters. Two classes of systematic uncertainties are defined, one from the angular fit model and another from the mass fit model. Since the angular fit is performed on the data weighted using the signal sWeights calculated from the fit to the J/ψ K − π + invariant mass, biases on the mass fit results may be propagated to the sWeights and thus to the angular parameters. Overall, two sources of systematic uncertainties dominate: the angular acceptance and the correlation between the J/ψ K − π + invariant mass and θ µ .

Systematic uncertainties related to the mass fit model
To determine the systematic uncertainty arising from the fixed parameters in the description of the J/ψ K − π + invariant mass, these parameters are varied inside their uncertainties, as determined from fits to simulated events. The fit is then repeated and the widths of the B 0 s and B 0 yield distributions are taken as systematic uncertainties on the branching fractions. Correlations among the parameters obtained from simulation are taken into account in this procedure. For each new fit to the J/ψ K − π + invariant mass, the corresponding set of sWeights is calculated and the fit to the weighted angular distributions is repeated. The widths of the distributions are taken as systematic uncertainties on the angular parameters. In addition, a systematic uncertainty is added to account for imperfections in the modelling  Table 3: Summary of the measured B 0 s → J/ψ K * 0 P-wave properties and their statistical and systematic uncertainties. When no value is given, it means an uncertainty below 5 × 10 −4 , except for the two phases, δ (rad) and δ ⊥ (rad), in which case the uncertainty is below 5 × 10 −3 . of the upper tail of the B 0 and B 0 s peaks. Indeed, in the Hypatia distribution model, the parameters a 2 and n 2 take into account effects such as decays in flight of the hadron, that affect the lineshape of the upper tail and could modify the B 0 leakage into the B 0 s peak. The estimate of this leakage is recalculated for extreme values of those parameters, and the maximum spread is conservatively added as a systematic uncertainty.
Systematic uncertainties due to the fixed yields of the B 0 s → J/ψ K + K − , B 0 s → J/ψ π + π − , B 0 → J/ψ π + π − , and Λ 0 b → J/ψ pK − peaking backgrounds, 3 are evalu- Table 4: Summary of the measured B 0 s → J/ψ K * 0 S-wave properties and their statistical and systematic uncertainties. When no value is given, it means an uncertainty below 5 × 10 −4 , except for the four phases related to the S-wave component, δ S (rad), in which case the uncertainty is below 5 × 10 −3 . The m K − π + binning definition is identical to the one given in Table 2. Correlations between the J/ψ K − π + invariant mass and the cosine of the helicity angle θ µ are taken into account in the nominal fit model, where the mass fit is performed in five bins of cos(θ µ ). In order to evaluate systematic uncertainties due to these correlations, the mass fit is repeated with the full range of cos(θ µ ) divided into four or six equal bins.
used to subtract them is constant in the nominal fit.
For each new mass fit, the angular fit is repeated using the corresponding set of sWeights. The deviations from the nominal result for each of the variations are summed quadratically and taken as the systematic uncertainty.

Systematic uncertainties related to the angular fit model
In order to account for systematic uncertainties due to the angular acceptance, two distinct effects are considered, as in Ref. [8]. The first is due to the limited size of the simulation sample used in the acceptance estimation. It is estimated by varying the normalisation weights 200 times following a Gaussian distribution within a five standard deviation range taking into account their correlations. For each of these sets of normalisation weights, the angular fit is repeated, resulting in a distribution for each fitted parameter. The width of the resulting parameter distribution is taken as the systematic uncertainty. Note that in this procedure, the normalisation weights are varied independently in each m K − π + bin. The second effect, labelled as data-simulation corrections in the tables, accounts for differences between the data and the simulation, using normalisation weights that are determined assuming the amplitudes measured in Ref. [47]. The difference with respect to the nominal fit is assigned as a systematic uncertainty. The uncertainties due to the choice of model for the C SP factors are evaluated as the maximum differences observed in the measured parameters when computing the C SP factors with all of the alternative models, as discussed below. Instead of the nominal propagator for the S-wave, a combination of the K * 0 (800) 0 and K * 0 (1430) 0 resonances with a non-resonant term using the isobar model is considered, as well as a K-matrix [48] version. A pure phase space term is also used, in order to account for the simplest possible parametrisation. For the P-wave, the alternative propagators considered are the K * (892) 0 alone and a combination of this contribution with the K * 1 (1410) 0 and the K * 1 (1430) 0 using the isobar model. In order to account for the absence of D-wave terms in the nominal fit model a new fit is performed, including a D-wave component, where the related parameters are fixed to the values measured in the K * 2 (1430) 0 region. The differences in the measured parameters between the results obtained with and without a D-wave component are taken as the corresponding systematic uncertainty.
The presence of biases in the fit model itself is studied using parametric simulation. For this study, 1000 pseudoexperiments were generated and fitted using the nominal shapes, where the generated parameter values correspond to the ones obtained in the fit to data. The difference between the generated value and the mean of the distribution of fitted parameter values are treated as a source of systematic uncertainty.
Finally, the systematic uncertainties due to the fixed values of the detection and production asymmetries are estimated by varying their values by ±1 standard deviation and repeating the fit.

Branching fraction
Several sources of systematic uncertainties on the branching fraction measurements are studied, summarised along with the results in Table 5: systematic uncertainties due to the external parameter f d /f s and due to the branching fraction B(φ → K + K − ); systematic uncertainties due to the ratio of efficiencies obtained from simulation and due to the angular parameters, propagated into the ω factors (see Sect. 8.1); and systematic uncertainties affecting the B 0 s → J/ψ K * 0 and B 0 → J/ψ K * 0 yields, which are determined from the fit to the J/ψ K + π − invariant mass and described in Sect. 8.1. Finally, a systematic uncertainty due to the B 0 s → J/ψ φ yield determined from the fit to the J/ψ K + K − invariant mass distribution, described in Sect. 7.3, is also taken into account, where only the effect due to the modelling of the upper tail of the B 0 s peak is considered (see Sect. 8.1.1). For the computation of the absolute branching fraction B(B 0 s → J/ψ K * 0 ) (see Sect. 7.5), two additional systematic sources are taken into account, the uncertainties in the external parameters B(B 0 → J/ψ K * 0 ) and B(B 0 s → J/ψ φ).  Following the strategy proposed in Refs. [9,11,13], the measured branching fraction, polarisation fractions and CP asymmetries can be used to quantify the contributions originating from the penguin topologies in B 0 s → J/ψ K * 0 . To that end, the transition amplitude for the B 0 s → J/ψ K * 0 decay is written in the general form

Relative branching fraction
where λ = |V us | = 0.22548 +0.00068 −0.00034 [6] and i labels the different polarisation states. In the above expression, A i is a CP -conserving hadronic matrix element that represents the tree topology, and a i parametrises the relative contribution from the penguin topologies. The CP -conserving phase difference between the two terms is parametrised by θ i , whereas their weak phase difference is given by the angle γ of the Unitarity Triangle.
Both the branching fraction and the CP asymmetries depend on the penguin parameters a i and θ i . The dependence of A CP i is given by [9] To use the branching fraction information an observable is constructed [9]: and Φ(x, y) = (1 − (x − y) 2 )(1 − (x + y) 2 ) is the standard two-body phase-space function. The primed quantities refer to the B 0 s → J/ψ φ channel, while the non-primed ones refer to B 0 s → J/ψ K * 0 . The penguin parameters a i and θ i are defined in analogy to a i and θ i , and parametrise the transition amplitude of the B 0 s → J/ψ φ decay as Assuming SU(3) flavour symmetry, and neglecting contributions from exchange and penguin-annihilation topologies, 4 which are present in B 0 s → J/ψ φ but have no counterpart in B 0 s → J/ψ K * 0 , we can identify The contributions from the additional decay topologies in B 0 s → J/ψ φ can be probed using the decay B 0 → J/ψ φ [13]. The current upper limit on its branching fraction is B(B 0 → J/ψ φ) < 1.9 × 10 −7 at 90% confidence level (C.L.) [50], which implies that the size of these additional contributions is small compared to those associated with the penguin topologies.
The H i observables are constructed in terms of the theoretical branching fractions defined at zero decay time, which differ from the measured time-integrated branching fractions [51] due to the non-zero decay-width difference ∆Γ s of the B 0 s meson system [7]. The conversion factor between the two branching fraction definitions [51] is taken to be where η i is the CP eigenvalue of the final state, and y s = ∆Γ s /2Γ s . Taking values for Γ s , ∆Γ s and φ SM s from Refs. [6,7], the conversion factor is 1.0608 ± 0.0045 (0.9392 ± 0.0045) for the CP -even (-odd) states. For the flavour-specific B 0 s → J/ψ K * 0 decay η i = 0, resulting in a conversion factor of 0.9963 ± 0.0006. The ratios of hadronic amplitudes |A i /A i | are calculated in Ref.
[52] following the method described in Ref. [53] and using the latest results on form factors from Light Cone QCD Sum Rules (LCSR) [54]. This leads to Assuming Eq. 28 and external input on the Unitarity Triangle angle γ = 73.2 +6.3 For the longitudinal polarisation state the phase θ is unconstrained. Correlations between the experimental inputs are ignored, but the effect of including them is small. The two-dimensional confidence level contours are given in Fig. 6. This figure also shows, as different (coloured) bands, the constraints on the penguin parameters derived from the individual observables entering the χ 2 fit. The thick inner darker line represents the contour associated with the central value of the input quantity, while the outer darker lines represent the contours associated with the one standard deviation changes. For the parallel polarisation the central value of the H observable does not lead to physical solutions in the θ -a plane, and the thick inner line is thus absent. When decomposed into its different sources, the angle φ s takes the form where −2β s is the SM contribution, φ BSM which is in good agreement with the values measured in Ref. [15], and with the predictions given in Refs. [12][13][14].
The above results are obtained assuming SU(3) flavour symmetry and neglecting contributions from additional decay topologies. Because a i e iθ i represents a ratio of hadronic amplitudes, the leading factorisable SU(3)-breaking effects cancel, and the relation between a i e iθ i and a i e iθ i is only affected by non-factorisable SU(3)-breaking. This can be parametrised using two SU(3)-breaking parameters ξ and δ as The above quoted results assume ξ = 1 and δ = 0. The dependence of the uncertainty on ∆φ J/ψ φ s,i on the uncertainty on ξ is illustrated in Fig. 7, while the dependence on the uncertainty on δ is negligible for the solutions obtained for {a i , θ i }.

Combination with
The information on the penguin parameters obtained from B 0 s → J/ψ K * 0 can be complemented with similar information from the SU(3)-related mode B 0 → J/ψ ρ 0 [15]. Both modes describe ab →ccd transition, and are related by exchanging the spectator s ↔ d quarks. The decay amplitude of B 0 → J/ψ ρ 0 is also parametrised as which is the equivalent of Eq. 23. In contrast to B 0 s → J/ψ K * 0 , however,ã i andθ i also include contributions from exchange and penguin-annihilation topologies, which are present in B 0 → J/ψ ρ 0 but have no counterpart in B 0 s → J/ψ K * 0 . Assuming SU(3) symmetry, and neglecting the contributions from the additional decay topologies in B 0 s → J/ψ φ and B 0 → J/ψ ρ 0 , the relation in Eq. 28 can be extended to which allows a combined fit to be performed to the CP asymmetries and branching fraction information in B 0 s → J/ψ K * 0 and B 0 → J/ψ ρ 0 .  The B 0 → J/ψ ρ 0 decay exhibits decay-time-dependent CP violation, which is described by two parameters, the direct CP asymmetry C i , which in the SU(3) limit is related to A CP i as C i = −A CP i , and the mixing-induced CP asymmetry S i . Their dependence on the penguin parametersã i andθ i is given by where η i is the polarisation-dependent CP eigenvalue of the B 0 → J/ψ ρ 0 decay, and φ d is a CP -violating phase arising from the interference between B 0 -B 0 mixing and the subsequent B 0 decay. The use of S i to constrain the penguin parameters a i and θ i requires external information on the CP phase φ d . The most precise value of φ d is determined from B 0 → J/ψ K 0 decays, but this determination is also affected by penguin pollution. A recent study of the penguin effects in B + → J/ψ π + , B + → J/ψ K + , B 0 → J/ψ π 0 and B 0 → J/ψ K 0 S decays is performed in Ref.
In addition, a second set of H i observables can be constructed by replacing B 0 s → J/ψ K * 0 by B 0 → J/ψ ρ 0 in Eq. 25. To minimise the theoretical uncertainties associated with the use of these H i observables, the strategy proposed in Ref.
[13] is adopted. That is, the relation between the hadronic amplitudes in B 0 s → J/ψ K * 0 and B 0 → J/ψ ρ 0 is assumed, and therefore relying on theoretical input from LCSR is no longer needed. Instead, the ratio |A /A| can be determined directly from the fit, providing experimental information on this quantity. Effectively, the three CP asymmetry parameters A CP i , C i and S i determine the penguin parameters a i and θ i . Thus, this result for a i and θ i predicts the values of the two observables H i (B 0 s → J/ψ K * 0 ) and H i (B 0 → J/ψ ρ 0 ). By comparing these two quantities with the branching fraction and polarisation information on B 0 s → J/ψ K * 0 , B 0 → J/ψ ρ 0 and B 0 s → J/ψ φ, the hadronic amplitude ratios |A i /A i | can be determined. The impact of the H i observables on the penguin parameters a i and θ i is negligible in the combined fit.
For the combined analysis of B 0 s → J/ψ K * 0 and B 0 → J/ψ ρ 0 a modified least-squares fit is performed. External inputs on γ = 73.2 +6.3  [16], with precision improved by a factor of 2 − 3. The shift on φ s due to penguin pollution is estimated from a combination with the B 0 → J/ψ ρ 0 channel [15], and is found be to compatible with the result from the earlier analysis.

Appendices A Angular acceptance
To take into account angular acceptance effects, ten normalisation weights, ξ ij , are computed and embedded in the normalization integral of the angular distribution given in Eq. 4 following the procedure described in Ref. [57]. Using the transversity amplitude basis, the fitting PDF can be written as where the real or imaginary angular functions F ij (θ K , θ µ , ϕ h ) are obtained when combining Eq. 4 and Eq. 5-8, and where Ω (θ K , θ µ , ϕ h ) denotes the angular acceptance. The normalization weights correspond to the integrals In the absence of acceptance effects, the normalisation weights related to the interference terms are equal to zero by definition, whereas those related to each polarisation amplitude squared are equal to unity. Eight sets of normalisation weights are calculated separately, one for each m K − π + bin and kaon charge. In order to correct both for imperfections in the detector simulation and for the absence of any S-wave component in the simulation sample, the weights are refined using an iterative procedure where the angular acceptance is re-evaluated recursively until it does not change significantly. Table 6 gives one set of normalisation weights after the Table 6: Corrected angular acceptance weights for K − π + events lying in the first m K − π + bin. The ξ ij weights are normalised with respect to the ξ 00 weight.
ij ξ ij /ξ 00 1 (00) 1.000 2 ( ) +1.379 ± 0.029 3 (⊥⊥) +1.388 ± 0.030 4 ( ⊥) +0.035 ± 0.019 5 (0 ) −0.003 ± 0.012 6 (0⊥) +0.010 ± 0.011 7 (SS) +1.190 ± 0.019 8 (S ) −0.042 ± 0.017 9 (S⊥) +0.029 ± 0.016 10 (S0) −0.929 ± 0.024 iterative procedure. The effect of this correction is below one standard deviation for all the normalisation weights except for the (S0) weight. This is expected due to the rapid efficiency drop close to cos θ K = 1 which directly impacts the (S0) weight. At each step of this procedure the simulation sample is corrected both for the absence of an S-wave component and for the imperfections in the detector simulation. For the first correction, the angular fit result to data is used, whereas for the second the kaon and muon track momentum distributions of data are used. In both cases the correction is implemented by assigning weights to each event of the simulation sample.

B Correlation matrix
The statistical-only correlation matrix of the angular parameters obtained from the fit to data, as described in Sect. 8.1, is given in Table 7. Here, the superscript l = 0, 1, 2, 3 in F l S and δ l S represent the number of the m K − π + bin as defined in Table 2.