Angular analysis of the $B^{0}\rightarrow K^{*0}\mu^{+}\mu^{-}$ decay using $3\,\mbox{fb}^{-1}$ of integrated luminosity

An angular analysis of the $B^{0}\rightarrow K^{*0}(\rightarrow K^{+}\pi^{-})\mu^{+}\mu^{-}$ decay is presented. The dataset corresponds to an integrated luminosity of $3.0\,{\mbox{fb}^{-1}}$ of $pp$ collision data collected at the LHCb experiment. The complete angular information from the decay is used to determine $C\!P$-averaged observables and $C\!P$ asymmetries, taking account of possible contamination from decays with the $K^{+}\pi^{-}$ system in an S-wave configuration. The angular observables and their correlations are reported in bins of $q^2$, the invariant mass squared of the dimuon system. The observables are determined both from an unbinned maximum likelihood fit and by using the principal moments of the angular distribution. In addition, by fitting for $q^2$-dependent decay amplitudes in the region $1.1<q^{2}<6.0\mathrm{\,Ge\kern -0.1em V}^{2}/c^{4}$, the zero-crossing points of several angular observables are computed. A global fit is performed to the complete set of $C\!P$-averaged observables obtained from the maximum likelihood fit. This fit indicates differences with predictions based on the Standard Model at the level of 3.4 standard deviations. These differences could be explained by contributions from physics beyond the Standard Model, or by an unexpectedly large hadronic effect that is not accounted for in the Standard Model predictions.


Introduction
The decay B 0 → K * 0 µ + µ − proceeds via a b → s quark flavour-changing neutral current (FCNC) transition. In the Standard Model (SM) the decay is therefore forbidden at tree level and occurs, at lowest order, via electroweak penguin and box processes. In extensions of the SM, new particles may enter in competing processes and can significantly change the branching fraction of the decay and the angular distribution of the final-state particles. Angular observables are of particular interest, since theoretical predictions of such observables tend to be less affected by the hadronic uncertainties associated with the B 0 → K * 0 transition. Throughout this paper K * 0 is used to refer to the K * (892) 0 resonance.
The LHCb collaboration previously determined a set of angular observables in the B 0 → K * 0 µ + µ − decay, 1 using data collected during 2011, corresponding to an integrated luminosity of 1.0 fb −1 [1]. Different subsets of these angular observables have also been measured by the BaBar, Belle, CDF and CMS collaborations [2-7] and all of these measurements are in good agreement with SM predictions. The LHCb collaboration has also used the 2011 dataset to determine an alternative set of angular observables that have reduced theoretical uncertainties [8]. In contrast to the previous analyses, these observables cannot be extracted from single angle distributions. This second LHCb analysis found a local deviation with respect to the SM prediction in one observable, P 5 , with a significance corresponding to 3.7 standard deviations. Possible interpretations of this discrepancy and the consistency of all of the measurements of b → s transitions have been widely discussed in the literature [9][10][11][12][13][14][15][16][17][18][19][20][21].
The present paper describes an updated angular analysis of the B 0 → K * 0 µ + µ − decay, using the LHCb Run 1 data sample, corresponding to an integrated luminosity of 3.0 fb −1 . The data were recorded in pp collisions at centre-of-mass energies of 7 and 8 TeV during 2011 and 2012, respectively. All previous analyses of the B 0 → K * 0 µ + µ − decay have extracted only part of the information available, by fitting simplified forms of the angular distribution. This paper presents a complete set of observables for the first time, based on the full angular distribution. The simultaneous determination of these observables allows correlations between the measured quantities to be computed, enabling the use of the results in global fits to theoretical models. This is critical to understand whether SM dynamics are sufficient to explain the above discrepancy, or if extensions to the SM are necessary.
The structure of this paper is as follows. In Sec. 2, the angular distribution and observables for the B 0 → K * 0 µ + µ − decay are presented. Section 3 describes the experimental setup. The reconstruction and selection of the B 0 → K * 0 µ + µ − candidates and sources of background are presented in Sec. 4. The method used to correct the angular distribution for experimental effects is detailed in Sec. 5 and the parameterisation of the mass distribution is described in Sec. 6. The determination of the angular observables is detailed in Sec. 7, and Sec. 8 discusses sources of systematic uncertainty. Results are given in Sec. 9 and the compatibility with predictions based on the Standard Model is 2 Angular distribution and observables The final state of the decay B 0 → K * 0 µ + µ − can be described by q 2 , the invariant mass squared of the dimuon system, and three decay angles Ω = (cos θ l , cos θ K , φ). The angle between the µ + (µ − ) and the direction opposite to that of the B 0 (B 0 ) in the rest frame of the dimuon system is denoted θ l . In this analysis, the K * 0 meson is reconstructed through the decay K * 0 → K + π − . The angle between the direction of the K + (K − ) and the B 0 (B 0 ) in the rest frame of the K * 0 (K * 0 ) system is denoted θ K . The angle between the plane defined by the dimuon pair and the plane defined by the kaon and pion in the B 0 (B 0 ) rest frame is denoted φ. More details of the angular basis adopted in this analysis are given in Appendix A of Ref. [1].
The differential decay rates of B 0 → K * 0 µ + µ − and B 0 → K * 0 µ + µ − decays, in terms of q 2 and the three angles, are given by where Γ (Γ) refers to decays involving a b (b) quark and hence a B 0 (B 0 ) meson, the terms f i ( Ω) are formed from combinations of spherical harmonics and the I i (Ī i ) are q 2 -dependent angular observables. The I i can be expressed as bilinear combinations of six complex decay amplitudes, A L,R 0, ,⊥ , which correspond to the different transversity states of the K * 0 meson and the different (left-and right-handed) chiralities of the dimuon system. An additional suffix s or c is conventionally added to some of the I i terms to indicate that they have a sin 2 θ K or cos 2 θ K dependence. When q 2 is sufficiently large (q 2 > ∼ 1 GeV 2 /c 4 ), the muons can be considered massless. The list of the angular terms and observables that remain in this massless limit is given in Table 1.
Following the notation of Ref. [22], q 2 -dependent CP averages, S i , and CP asymmetries, A i , can be defined as (2) In the massless limit, the CP -averaged observables S 1(s,c) and S 2(s,c) obey the relations S 1s = 3S 2s , S 1c = −S 2c and 3 4 (2S 1s + S 1c ) − 1 4 (2S 2s + S 2c ) = 1 (see for example Ref. [22]). These relationships reduce the number of independent CP -averaged observables from eleven to eight. The relations between the observables also hold to a good approximation for q 2 < 1 GeV 2 /c 4 and are therefore adopted for the full q 2 range. The S 1c observable corresponds to the fraction of longitudinal polarisation of the K * 0 meson and is therefore more commonly referred to as F L , with It is also conventional to replace S 6s by the forward-backward asymmetry of the dimuon system A FB , with A FB = 3 4 S 6s . The CP -averaged angular distribution of the B 0 → K * 0 µ + µ − decay can then be written as 1 d(Γ +Γ)/dq 2 d 4 (Γ +Γ) dq 2 d Ω = 9 32π 3 4 (1 − F L ) sin 2 θ K + F L cos 2 θ K + 1 4 (1 − F L ) sin 2 θ K cos 2θ l −F L cos 2 θ K cos 2θ l + S 3 sin 2 θ K sin 2 θ l cos 2φ +S 4 sin 2θ K sin 2θ l cos φ + S 5 sin 2θ K sin θ l cos φ + 4 3 A FB sin 2 θ K cos θ l + S 7 sin 2θ K sin θ l sin φ +S 8 sin 2θ K sin 2θ l sin φ + S 9 sin 2 θ K sin 2 θ l sin 2φ . (4) Additional sets of observables, for which the leading B 0 → K * 0 form-factor uncertainties cancel, can be built from F L and S 3 -S 9 . Examples of such optimised observables include the transverse asymmetry A (2) T [23], where A (2) T = 2S 3 /(1 − F L ), and the P ( ) i series of observables [24]. In this paper the notation used is T , The definition of the P i observables differs from that of Ref. [24], but is consistent with the notation used in the LHCb analysis of Ref. [8].
In addition to the resonant P-wave K * 0 contribution to the K + π − µ + µ − final state, the K + π − system can also be in an S-wave configuration. The addition of an S-wave component introduces two new complex amplitudes, A L,R S , and results in the six additional angular terms that are given in the lower part of Table 1. In the analyses described in Refs. [1,8] the S-wave contribution, which is expected to be approximately 5%, was treated as a systematic uncertainty. The presence of a K + π − system in an S-wave configuration modifies the angular distribution to 1 d(Γ +Γ)/dq 2 d 4 (Γ +Γ) F S sin 2 θ l + 9 32π (S 11 + S 13 cos 2θ l ) cos θ K + 9 32π (S 14 sin 2θ l + S 15 sin θ l ) sin θ K cos φ + 9 32π (S 16 sin θ l + S 17 sin 2θ l ) sin θ K sin φ , where F S denotes the S-wave fraction, and the terms S 11 , S 13 -S 17 arise from interference between the S-and P-wave amplitudes. Note that F S replaces the terms S 10 and S 12 , with F S = 3S 10 = −3S 12 . Throughout this paper, F S and the interference terms between the S-and P-wave are treated as nuisance parameters.
Due to the flavour specific final state of the decay, the CP asymmetries A i can be determined from differences in the angular distributions between B 0 and B 0 decays.
In this analysis, three separate techniques are used to study the angular distribution: 1. An unbinned maximum likelihood fit is used to determine the CP -averaged observables F L , A FB , and S 3 -S 9 , as well as the CP asymmetries A 3 -A 9 , averaged over bins of q 2 . In addition, the P ( ) i observables are determined by reparameterising the likelihood fit. The data are analysed in q 2 bins of approximately 2 GeV 2 /c 4 width and also in wider 1.1 < q 2 < 6.0 GeV 2 /c 4 and 15.0 < q 2 < 19.0 GeV 2 /c 4 bins for which there are particularly precise theoretical predictions. The unbinned maximum likelihood fit is described in Sec. 7.1.
2. The same observables are also determined using principal angular moments. This so-called method of moments gives an approximately 15% less precise determination of the observables than the likelihood fit but is particularly robust for low signal yields and does not require a complex angular fit [25]. This allows the observables to be determined in approximately 1 GeV 2 /c 4 wide q 2 bins, which gives additional shape information that is useful in regions where the observables vary rapidly with q 2 . The method is described in Sec. 7.2. Table 1: Angular observables I j and their corresponding angular terms for dimuon masses that are much larger than twice the muon mass. The terms in the lower part of the table arise from the K + π − S-wave contribution to the K + π − µ + µ − final state. TheĪ i coefficients are obtained by making the substitution A →Ā, i.e. by complex conjugation of the weak phases in the amplitudes. i sin 2θ K sin θ l sin φ  3 Im(A L S A L * ⊥ + A R S A R * ⊥ ) sin θ K sin 2θ l sin φ 3. Finally, the observables S 4 , S 5 and A FB vary as a function of q 2 and are known to change sign in the SM. By fitting for the decay amplitudes as a function of q 2 , the q 2 values at which these observables cross zero can be determined. At leading order these zero-crossing points are free from B 0 → K * 0 form-factor uncertainties and consequently provide a precision test of the SM [26,27]. The method is applied in the range 1.1 < q 2 < 6.0 GeV 2 /c 4 and is described in Sec. 7.3.
The three methods are complementary, but their results are correlated and cannot be combined. Method 1 is the most precise and is therefore used to compare to the SM predictions. The q 2 bins used for the likelihood fit of the angular observables and the method of moments are given in Tables 4 and 7 of Appendix A, respectively.

Detector and simulation
The LHCb detector [28, 29] 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 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 pp interaction vertex (PV), 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 (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 [30]. Simulated signal events are used to determine the impact of the detector geometry, trigger, reconstruction and candidate selection on the angular distribution of the signal. In addition, simulated samples are used to estimate the contribution of possible background processes. In the simulation, pp collisions are generated using Pythia [31] with a specific LHCb configuration [32]. Decays of hadronic particles are described by EvtGen [33], in which final-state radiation is generated using Photos [34]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [35] as described in Ref. [36]. Data-driven corrections are applied to the simulation to account for a small level of mismodelling of the detector occupancy, B 0 momentum and B 0 vertex quality. Similarly, the simulated particle identification (PID) performance is corrected to match that determined from control samples selected from the data.

Selection of signal candidates
The B 0 → K * 0 µ + µ − signal candidates are required to pass a hardware trigger, which selects events containing at least one muon with 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 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.
Signal candidates are formed from a pair of oppositely charged tracks that are identified as muons, combined with a K * 0 meson candidate. The K * 0 candidate is formed from two charged tracks that are identified as a kaon and a pion, respectively. The four tracks of the final-state particles are required to have a significant impact parameter with respect to all PVs in the event. The tracks are then fitted to a common vertex, which is required to be of good quality. The impact parameter of the B 0 candidate with respect to one of the PVs is required to be small and the vertex of the B 0 candidate is required to be significantly displaced from the same PV. The angle θ DIRA between the reconstructed B 0 momentum and the vector connecting the PV to the reconstructed B 0 decay vertex is required to be small. Candidates are required to have reconstructed B 0 invariant mass, m(K + π − µ + µ − ), in the range 5170 < m(K + π − µ + µ − ) < 5700 MeV/c 2 . Finally, the reconstructed mass of the K + π − system, m(K + π − ), is required to be in the range 796 < m(K + π − ) < 996 MeV/c 2 .
Background formed by combining particles from different b-and c-hadron decays (referred to as combinatorial background) is further reduced using a boosted decision tree (BDT) [37,38], which is trained using data. As a proxy for the signal decay, B 0 → J/ψ K * 0 decays are used to train the BDT, where the J/ψ is reconstructed through its decay into µ + µ − . Candidates from the upper mass sideband 5350 < m(K + π − µ + µ − ) < 7000 MeV/c 2 are used as a proxy for the background. As input variables, the BDT uses the reconstructed B 0 lifetime and vertex fit quality, the momentum and transverse momentum of the B 0 candidate, cos θ DIRA , particle identification information from the RICH detectors and the muon system, as well as variables describing the isolation of the final state tracks [39]. To best exploit the data available for training, the k-folding technique [40] is employed with k = 10. At the chosen working point, the BDT has a background rejection of 97% and a signal efficiency of 85%. The signal efficiency and background rejection of the BDT is uniform in m(K + π − µ + µ − ) and m(K + π − ). The distortion induced in q 2 and the angular distributions is discussed in Sec. 5.
The K + π − µ + µ − invariant mass versus q 2 for candidates that pass the full selection is shown in Fig. 1. The B 0 → K * 0 µ + µ − signal candidates are clearly visible as a vertical band. The contributions from the decays B 0 → J/ψ K * 0 and B 0 → ψ(2S)K * 0 , which proceed through tree-level b → ccs transitions, have a dimuon mass consistent with the known J/ψ or ψ(2S) meson mass and m(K + π − µ + µ − ) consistent with that of the known B 0 meson mass. The horizontal bands are formed from combinatorial background comprising a genuine J/ψ or ψ(2S) meson and a K * 0 candidate selected from elsewhere in the event. : Invariant mass of the K + π − µ + µ − system versus q 2 . The decay B 0 → K * 0 µ + µ − is clearly visible inside the dashed vertical lines. The horizontal lines denote the charmonium regions, where the tree-level decays B 0 → J/ψ K * 0 and B 0 → ψ(2S)K * 0 dominate. These candidates are excluded from the analysis.

Background composition
In addition to combinatorial background, there are several sources of background that accumulate in m(K + π − µ + µ − ) and can potentially mimic the signal decay if they are misreconstructed in the detector. These are referred to as peaking backgrounds. Contamination from peaking backgrounds is estimated using samples of simulated events. The tree-level decays B 0 → J/ψ K * 0 and B 0 → ψ(2S)K * 0 dominate in the regions 8.0 < q 2 < 11.0 GeV 2 /c 4 and 12.5 < q 2 < 15.0 GeV 2 /c 4 , respectively, and these q 2 regions are therefore excluded from the analysis of the signal decay. However, these decays can still form a source of background if the µ − (µ + ) is misidentified as a π − (K + ) and the π − (K + ) is misidentified as a µ − (µ + ). To remove this background, candidates are rejected if the π − (K + ) satisfies the muon identification criteria and the mass of the π − µ + (K + µ − ) system, when the π − (K + ) is assigned the muon mass, is consistent with that of a J/ψ or ψ(2S) meson. Possible pollution from B 0 → K * 0 φ(→ µ + µ − ) decays is removed by excluding from the analysis the q 2 region 0.98 < q 2 < 1.10 GeV 2 /c 4 .
The decay Λ 0 b → pK − µ + µ − , which can proceed via e.g. the Λ(1520) resonance, can be a source of peaking background if the proton is misidentified as a pion. This background is suppressed by rejecting candidates where the pion is not unambiguously identified by the RICH detectors and which have a mass close to the known Λ 0 b mass, when the pion is assigned the proton mass. Similarly, Λ 0 b → pK − µ + µ − backgrounds with double misidentification of the hadrons, i.e. where the proton is misidentified as a kaon and the kaon is misidentified as a pion, are suppressed using PID information.
The decay B 0 s → φ(→ K + K − )µ + µ − can mimic the signal decay if one of the kaons is misidentified as a pion. This background is suppressed by requiring stringent PID criteria if, after assigning the kaon mass to the pion candidate, the reconstructed invariant masses of the B 0 and K * 0 candidates are consistent with the known B 0 s and φ masses. The decay B + → K + µ + µ − can form a background if a low momentum pion from elsewhere in the event is added to form a four-particle final state. The resulting invariant mass m(K + π − µ + µ − ) will be larger than the known B 0 mass but can contribute to the upper mass-sideband. Such decays can therefore distort the estimate of the angular distribution of the residual background, which is assessed from this sideband. This background is suppressed by removing candidates with 5220 < m(K + µ + µ − ) < 5340 MeV/c 2 . It is also possible to have backgrounds from B 0,+ → K * 0,+ µ + µ − decays, where the pion from the K * meson is replaced by another pion from the rest of the event. This background does not peak in the signal region and is considered as part of the combinatorial background. Finally, B 0 → K * 0 µ + µ − decays can form a background to B 0 → K * 0 µ + µ − decays (and vice versa) if the K + (K − ) is misidentified as the π + (π − ) and the π − (π + ) is misidentified as the K − (K + ). These misidentified decays are referred to as signal swaps and are suppressed using PID information.
After all vetoes are applied, the largest peaking background contribution is from The residual background from these decays is expected to be at a level of (1.0 ± 0.4)% of the signal yield. The next largest backgrounds are B 0 → K * 0 µ + µ − signal swaps at (0.64 ± 0.06)%, misidentified B 0 s → φµ + µ − events at (0.33 ± 0.12)% and B 0 → J/ψ K * 0 decays with double misidentification at (0.05 ± 0.05)% of the signal yield. All of the sources of peaking background are sufficiently small such that they are neglected in the angular analysis but are considered further as sources of systematic uncertainty. The background from b-hadron decays where two hadrons are misidentified as muons is negligible. The largest residual background is combinatorial in nature and varies smoothly with m(K + π − µ + µ − ), m(K + π − ) and the decay angles.

Angular acceptance
The triggering, reconstruction and selection criteria distort the distributions of the decay angles θ l , θ K and φ, as well as the q 2 distribution, giving rise to so-called acceptance effects. The dominant acceptance effects come from momentum and impact parameter requirements. In particular, the implicit momentum threshold that is required for tracks to traverse the magnetic spectrometer removes low momentum particles. In contrast to the previous LHCb analyses [1,8], the acceptance is not assumed to factorise in the three decay angles. Instead, the efficiency is parameterised in four dimensions, according to where the terms L h (x) denote Legendre polynomials of order h and the observables are rescaled to the range −1 < x < +1 when evaluating the polynomial. For cos θ l , cos θ K and φ, the sum in Eq. (8) encompasses L h (x) up to fourth, fifth and sixth order, respectively. The q 2 parameterisation comprises L h (x) up to fifth order. The coefficients c ijmn are determined using a principal moment analysis of simulated three-body B 0 → K * 0 µ + µ − phase-space decays. As the efficiency is parameterised in terms of all of the relevant kinematic variables needed to describe the decay, it does not depend on the model used in the simulation. The angular acceptance in cos θ l , cos θ K and φ is shown for 0.10 < q 2 < 0.98 GeV 2 /c 4 and 18.0 < q 2 < 19.0 GeV 2 /c 4 in Fig. 2. The acceptance varies smoothly as a function of q 2 between these extremes. The acceptance as a function of q 2 , after integrating over the decay angles, is also shown. The description of the angular acceptance is cross-checked, for q 2 = m 2 (J/ψ ), using the decay B 0 → J/ψ K * 0 . This decay can be selected in the data with background contamination below 1% and the angular structure has been determined by measurements made by the BaBar, Belle and LHCb collaborations [41][42][43]. With the acceptance correction derived using the above method, the B 0 → J/ψ K * 0 angular observables obtained from the LHCb data are in good agreement with these previous measurements. The angular fit of the B 0 → J/ψ K * 0 data is shown in Fig. 15 of Appendix B.
6 The K + π − µ + µ − and K + π − mass distributions The K + π − µ + µ − invariant mass is used to discriminate between signal and background. The distribution of the signal candidates is modelled using the sum of two Gaussian functions with a common mean, each with a power-law tail on the low-mass side. The parameters describing the signal mass-shape are determined from a fit to the B 0 → J/ψ K * 0 decay in the data, as shown in Fig. 3, and are subsequently fixed when fitting the B 0 → K * 0 µ + µ − candidates. In samples of simulated B 0 → K * 0 µ + µ − decays, the mass resolution is observed to vary with q 2 by 2-8%. A scale factor is therefore taken from the simulation and is used to correct the width of the Gaussian functions in the different q 2 bins. A component is included in the fit to account for B 0 s → J/ψ K * 0 decays, which are at a level of 0.8% of the B 0 → J/ψ K * 0 signal yield [44]. However, the B 0 s → K * 0 µ + µ − decay is neglected when fitting B 0 → K * 0 µ + µ − candidates. Combinatorial background is described well by a single exponential distribution in m(K + π − µ + µ − ). The B 0 → K * 0 µ + µ − signal yield integrated over the q 2 ranges 0.1 < q 2 < 8.0 GeV 2 /c 4 , 11.0 < q 2 < 12.5 GeV 2 /c 4 and 15.0 < q 2 < 19.0 GeV 2 /c 4 is determined to be 2398 ± 57. The signal yield in the range 1.1 < q 2 < 6.0 GeV 2 /c 4 is 624 ± 30.
As detailed in Secs. 7.1 and 7.2, the likelihood fit and the method of moments use additional information from the m(K + π − ) distribution to constrain the fraction of K + π − S-wave present in the data. To describe this distribution, the K * 0 signal component is modelled using a relativistic Breit-Wigner function for the P-wave component and the LASS parameterisation [45] for the S-wave component. The combinatorial background is described by a linear function in m(K + π − ). There is no evidence for a K * 0 component in the background m(K + π − ) distribution. When fitting for K * 0 decay amplitudes, m(K + π − ) is integrated over, as detailed in Sec. 7.3.  Figure 2: Relative efficiency in cos θ l , cos θ K , φ and q 2 , as determined from a principal moment analysis of simulated three-body B 0 → K * 0 µ + µ − phase-space decays. The efficiency as a function of cos θ l , cos θ K and φ is shown for the regions 0.1 < q 2 < 0.98 GeV 2 /c 4 (black solid line) and 18.0 < q 2 < 19.0 GeV 2 /c 4 (red dashed line). The efficiency as a function of q 2 is shown after integrating over the decay angles. The histograms indicate the distribution of the simulated three-body B 0 → K * 0 µ + µ − phase-space decays used to determine the acceptance.

Angular analysis of the decay
The three methods used to determine the CP -averaged angular observables, CP asymmetries and the zero-crossing points of S 4 , S 5 and A FB are detailed below. Section 7.1 describes the determination of the observables in bins of q 2 using a maximum likelihood fit. Section 7.2 discusses the determination of the same set of observables using a principal moment analysis. Finally, Sec. 7.3 describes a fit to the angular and q 2 distribution of the decay, parameterised in terms of the decay amplitudes rather than the observables. This fit is used to determine the zero-crossing points of S 4 , S 5 and A FB .

Determination of angular observables with a likelihood fit
In each q 2 bin, an unbinned maximum likelihood fit to m(K + π − µ + µ − ) and the three decay angles cos θ l , cos θ K and φ is used to determine the angular observables introduced in Sec. 2. The angular distribution of the signal is described using Eq. (6). The background angular distribution is modelled with second order polynomials in cos θ l , cos θ K and φ, the parameters of which are left free in the fit. The angular distribution is assumed to factorise in the three decay angles. This assumption has been validated in the upper mass sideband.
In order to describe the signal angular distribution, the angular acceptance discussed in Sec. 5 must be accounted for. The acceptance is treated in one of two ways, depending on the q 2 range being fitted. In the narrow q 2 bins, the acceptance is treated as being constant across each bin and is included in the fit by multiplying Eq. (6) by the acceptance function evaluated at the bin centre. In the wider 1.1 < q 2 < 6.0 GeV 2 /c 4 and 15.0 < q 2 < 19.0 GeV 2 /c 4 bins, the shape of the acceptance can vary significantly across the bin. In this case, the candidates are weighted in the likelihood fit by the inverse of their efficiency. The event weights are scaled such that this pseudo-likelihood fit has confidence intervals with the correct coverage.
The K + π − µ + µ − invariant mass is included in the fit to separate signal from background. The signal and background mass distributions are parameterised as described in Sec. 6. In order to better constrain the S-wave fraction, a simultaneous fit of the m(K + π − ) distribution is performed using the parameterisation described in Sec. 6. The signal fraction and F S are common parameters in the simultaneous fits to the m(K + π − ) distribution and to the angular and m(K + π − µ + µ − ) distributions. Figure 4 shows the projections of the fitted probability density function on the angular and mass distributions for the 1.1 < q 2 < 6.0 GeV 2 /c 4 q 2 bin. Good agreement of the fitted function with the data is observed. Projections for the other q 2 bins are provided in Appendix B.
The P ( ) i observables introduced in Sec. 2 are determined by reparameterising Eq. (4) using a basis comprising F L , P 1,2,3 and P 4,5,6,8 . The CP asymmetries are determined by modifying the angular convention, introducing a relative sign between the angular terms : Angular and mass distributions for 1.1 < q 2 < 6.0 GeV 2 /c 4 . The distributions of m(K + π − ) and the three decay angles are given for candidates in the signal mass window ±50 MeV/c 2 around the known B 0 mass. The candidates have been weighted to account for the acceptance. Overlaid are the projections of the total fitted distribution (black line) and its different components. The signal is shown by the blue shaded area and the background by the red hatched area.
f 3 ( Ω)-f 9 ( Ω) for B 0 and B 0 decays, such that Eq. (4) is given in terms of F L and the CP asymmetries A 3 -A 9 . The B 0 or B 0 flavour is determined from the charge of the final-state kaon.
To ensure correct coverage for the uncertainties of the angular observables, the Feldman-Cousins method [46] is used with nuisance parameters treated according to the plug-in method [47]. Angular observables are considered one at a time, with the other angular observables treated as nuisance parameters. The nuisance parameters also include the signal fraction, the background parameters, F S and the angular terms that arise from interference between the S-and P-wave.

Determination of angular observables using the method of moments
The angular observables are also determined using a principal moment analysis of the angular distribution, without making any angular fit to the data [25,48]. As a continuous function of q 2 , the moments are given by The average M i (q 2 ) in a bin of q 2 is estimated by replacing the integral in Eq. 9 with a sum over the candidates in the dataset. The angular acceptance is accounted for by weighting the candidates in the sum, The sum is evaluated for candidates within ±50 MeV/c 2 of the B 0 mass. The weight, w e , is the reciprocal of the candidate's efficiency and is computed as described in Sec. 5. The mass window contains more than 95% of the signal candidates. This sum is also computed for candidates with 5350 < m(K + π − µ + µ − ) < 5700 MeV/c 2 and the resulting value of M i is used to subtract the background contribution from the ±50 MeV/c 2 window. The functions f i are given in Table 1. Due to their dependence on the spherical harmonics, most of the angular terms are orthogonal. For f i=3-9 , such that the moments give the CP -averaged observables S 3 to S 9 with a coefficient, λ i , that takes into account the normalisation. In the limit of massless muons, the moments are related to the observables by the expressions where, as noted previously, A FB = 3 4 S 6s . The relevant signal and background yields and the S-wave fraction F S are determined from a two-dimensional extended unbinned maximum likelihood fit to the m(K + π − µ + µ − ) and m(K + π − ) distributions, using the shapes described in Sec. 6.
The statistical uncertainties of the angular moments are estimated using a bootstrapping technique [49]. Confidence intervals are defined such that they include the 16 th -84 th percentiles of the bootstrap distribution of the observables. When computing the P ( ) i observables, bootstrapped data with unphysical F L (F L < 0 or F L > 1) are added at ±∞ to ensure that the resulting intervals do not undercover. As in the likelihood method, the CP asymmetries are determined by flipping the sign of the relevant B 0 angular terms. The resulting moments are then used to determine the CP asymmetries by substituting In the moment analysis, an additional angular observable that is not present in the massless limit is determined. This observable is sensitive to large new scalar or tensor contributions to the decay and is associated with a new forward-backward asymmetry of the dimuon system, f 6c ( Ω) = cos 2 θ K cos θ l [22,50]. The corresponding observable is highly correlated to A FB but can be determined from the moments M 6c and M 6s , using S 6c = 8M 6c − 2M 6s .

Determination of zero-crossing points using the decay amplitudes
In the 1.1 < q 2 < 6.0 GeV 2 /c 4 region, it is also possible to determine the amplitudes A L,R 0, ,⊥ , appearing in Table 1, using a smoothly varying q 2 -dependent parameterisation. For q 2 > ∼ 6.0 GeV 2 /c 4 , resonant cc states make a simple parameterisation of the q 2 -dependence impossible. A similar problem exists below 1.1 GeV 2 /c 4 due to the presence of light resonances.
The amplitudes A L,R 0, ,⊥ are complex functions of q 2 and therefore, at each point in q 2 , the decay B 0 → K * 0 µ + µ − is described by twelve real degrees of freedom. Several symmetries leave the angular distribution of the final-state particles unchanged [51]. These symmetries allow four components of the amplitudes to be set to zero. This simplification results in eight independent degrees of freedom which completely describe the B 0 → K * 0 µ + µ − decay. In this paper, the choice is made and the P-wave amplitudes are expressed using the form where α L,R i , β L,R In the q 2 range considered, the S-wave amplitudes are expected to vary slowly with q 2 [53]. To simplify the fit, these amplitudes are therefore assumed to be constant in q 2 and are described with a single complex parameter. The systematic uncertainty related to this approximation is negligible. After applying the symmetry constraints, the B 0 and B 0 decays are each described by 24 real parameters for the P-wave amplitudes and four real parameters for the S-wave amplitudes. With the 3 fb −1 dataset, it is not possible to determine the parameters describing both the B 0 and B 0 decays separately. It is therefore assumed that CP symmetry holds in the decay such that the amplitudes describing the B 0 and B 0 decays are identical.
An unbinned maximum likelihood fit to the distributions of m(K + π − µ + µ − ), cos θ l , cos θ K , φ and q 2 is used to determine the amplitude parameters. The integral of the angular distribution is required to be consistent with the number of signal candidates in the fit. For simplicity, m(K + π − ) is not included in the fit. The variation of the amplitudes with m(K + π − ) is accounted for by replacing products of amplitudes where g i (m(K + π − )) describes the variation of the amplitude A L,R i with m(K + π − ). The same models are used for the S-and P-wave lineshapes as in Secs. 7.1 and 7.2. The acceptance as a function of cos θ l , cos θ K , φ and q 2 (as described in Sec. 5) is included in the amplitude fit. The combinatorial background is parameterised by a linear function in q 2 . The background angular distribution is assumed to be independent of q 2 and is described by the product of three second-order polynomials. The background model and its factorisation in the decay angles and q 2 is checked using candidates in the m(K + π − µ + µ − ) sideband. Figure 5 shows the projections of the fitted probability density functions on the angular and q 2 distributions of the candidates. In contrast to Fig. 4, the effect of the selection efficiency on the angles and q 2 is included in the signal distribution. The figure therefore shows the distribution of the candidates rather than the candidates weighted by the inverse of their selection efficiency. Good agreement of the fitted function and the data is observed.
The amplitude parameters are used to construct observables as continuous functions of q 2 . The observables S 4 , S 5 and A FB have zero-crossing points and these are determined by solving a quartic equation. The different solutions of this equation are separable based on the sign of the slope of the observable in the vicinity of the zero-crossing point. Only zero-crossing points in the range 1.1 < q 2 < 6.0 GeV 2 /c 4 with a local slope consistent with the data above 6.0 GeV 2 /c 4 are retained.
The large number of parameters floating in the fit, coupled with the limited number of signal candidates present in the dataset, results in a non-parabolic likelihood surface. Therefore, as in the determination of the angular moments, the statistical uncertainties of the q 2 -dependent observables and their corresponding zero-crossing points are determined using a bootstrapping technique [49]. The statistical coverage of the resulting intervals is checked using simulated events and is found to be correct for the observables S 4 , S 5 and A FB . Despite the coverage being correct, approximately 10% of the bootstrapped datasets result in no zero-crossing point with the correct slope in the q 2 range 1.1 < q 2 < 6.0 GeV 2 /c 4 . In these cases, the zero-crossing point is added to the bootstrap distribution at ±∞ to ensure that the method does not undercover. The determination of the q 2 -dependent amplitudes in principle allows the full observable basis to be determined. However, pseudoexperiments indicate that a larger dataset is required in order to guarantee the correct coverage of the uncertainties on the observables other than S 4 , S 5 and A FB .

Systematic uncertainties
Effects that can alter the mass or angular distribution of either the signal or background candidates are sources of systematic uncertainty. The various sources of systematic uncertainty are discussed in detail below and are summarised in Table 2. In general, the systematic uncertainties are significantly smaller than the statistical uncertainties.
The size of each systematic uncertainty is estimated using pseudoexperiments in which one or more parameters are varied. The angular observables are determined from these pseudoexperiments using the nominal model and the systematically varied model. For each Table 2: Summary of the different sources of systematic uncertainty on the angular observables. Upper limits or typical ranges are quoted for the different groups of observables. The column labelled q 2 0 corresponds to the zero-crossing points of S 4 , S 5 and A FB .
.01 Det. and prod. asymmetries --< 0.01 < 0.02 observable, in each q 2 region, the systematic uncertainty is then taken as the average of the difference between the two models. The pseudoexperiments are generated with signal yields many times larger than that of the data, in order to render statistical fluctuations negligible.
The main systematic effects associated with the signal modelling arise from the estimate of the angular acceptance. Four separate sources of systematic uncertainty are considered: the statistical uncertainty on the acceptance correction resulting from the limited size of the simulation sample from which it is determined; an uncertainty associated with the parameterisation that is used to describe the acceptance function; an uncertainty arising from residual data-simulation differences; and, for the likelihood fit of the angular observables in narrow q 2 bins, an uncertainty associated with evaluating the acceptance at a fixed point in q 2 .
The statistical uncertainty on the acceptance function is evaluated using pseudoexperiments that are generated by coherently fluctuating the acceptance parameters according to the covariance matrix for the angular moments of the acceptance function. To evaluate the uncertainty associated with the particular choice of order for the polynomials used to describe the acceptance function, pseudoexperiments are produced in which the polynomial order is simultaneously increased by two in q 2 and in each of the angles.
After the B 0 momentum spectrum, detector occupancy and PID performance of the simulation are corrected to match the data, there is very good agreement between the properties of simulated and genuine B 0 → J/ψ K * 0 decays. There are, however, some small remaining differences in the momentum and transverse momentum spectra of the reconstructed pion that can affect the determination of the acceptance correction. A new acceptance correction is derived after re-weighting the simulated phase-space sample to account for the observed differences. A more conservative variation has also been considered in which an acceptance correction is derived without any of the data-simulation corrections applied. The larger of the variations observed is added as a systematic uncertainty.
When determining the angular observables in the narrow q 2 bins with the maximum likelihood fit, the acceptance is evaluated using the q 2 value of the bin centre. Pseudoexperiments are generated to assess the bias generated by this choice, using instead the value of q 2 of the left-or right-hand bin boundary.
Possible contributions from the tails of higher mass K * states in the 796 < m(K + π − ) < 996 MeV/c 2 window are also considered. Simulation studies indicate that any bias arising from these states is negligible compared to the statistical uncertainty on the angular observables.
For the background modelling, two sources of systematic uncertainty are considered. The first source is associated with the choice of second-order polynomials to model the background angular distribution in the fits of the angular observables and the q 2 -dependent decay amplitudes. It is not possible to fit a more complex model to the data because of the small number of background candidates. Therefore, to test the model, the BDT requirement is relaxed and the background candidates are fitted with a fourth-order polynomial in each of the three angles. This shape is used when generating the pseudoexperiments. The second source is associated with the fit to the q 2 -dependent decay amplitudes. In this case, the q 2 dependence of the background model is modified from a linear function to a third-order polynomial.
Systematic uncertainties are assessed for the different sources of peaking background that are neglected in the analysis. As detailed in Sec. 4.1, the most important backgrounds are those from Λ 0 b → pK − µ + µ − and B 0 s → φµ + µ − decays, where a kaon or proton is misidentified as pion; and B 0 → K * 0 µ + µ − decays, where the kaon and pion are both misidentified. Taking the angular distribution of the background from simulated events, pseudoexperiments are generated with these backgrounds included, and the angular observables determined as if the background were not present. Pseudoexperiments are also generated in which the angular distribution of the B 0 s → φµ + µ − and Λ 0 b → pK − µ + µ − decays are taken from data. These decays are selected by removing PID information from the BDT and inverting the background vetoes.
Systematic uncertainties are also assessed for the signal mass-modelling in m(K + π − µ + µ − ) and m(K + π − ). To assess the model of m(K + π − µ + µ − ), a fit is performed to B 0 → J/ψ K * 0 data using the sum of two Gaussian distributions without the power law tails. To assess the modelling of m(K + π − ), pseudoexperiments are produced by systematically varying the S-and P-wave line-shape parameters. For the S-wave, the LASS line-shape is also exchanged for the sum of resonant K * 0 (800) 0 (sometimes referred to as the κ resonance) and K * 0 (1430) 0 contributions. For the fit to the q 2 -dependent decay amplitudes, an additional uncertainty is assigned for the choice of the q 2 parameterisation of the S-wave components. As described in Sec. 7.3, the S-wave amplitudes are taken to be constant in q 2 . Motivated by Ref.
[53], a systematic variation is considered by assuming that the S-wave amplitudes A L,R S have the same q 2 dependence as the longitudinal P-wave amplitudes A L,R 0 .
The measured CP asymmetries can be biased due to detection and production asymmetries. The B 0 production asymmetry is measured to be less than 1% [54,55]. The effect of this asymmetry is further suppressed due to B 0 -B 0 mixing. The kaon detection asymmetry was measured in Ref. [56]. In contrast to the other sources of systematic uncertainty, the shift due to the detection and production asymmetries is calculated directly without generating pseudoexperiments. The systematic uncertainty on the angular observables A i due to production and detection asymmetries is found to be less than 0.01. The effect of these asymmetries on the CP -averaged observables is negligible.
In the q 2 -bin 0.10 < q 2 < 0.98 GeV 2 /c 4 , the muon mass-squared is comparable to q 2 and the relations between S 1(s,c) , S 2(s,c) and F L (see Sec. 2) are only approximate. The assumption that these relations hold has no impact on the measured values of S 3 -S 9 or A 3 -A 9 but results in a biased estimate of F L and hence of the P ( ) i observables. In pseudoexperiments based on the SM, this bias is typically at the level of 0.02. This can be accounted for in the SM predictions for this q 2 -bin and hence is not considered as a source of systematic uncertainty.
For F L and A FB , the largest source of systematic uncertainty comes from the datasimulation comparison of the pion momenta. The systematic uncertainty assigned to this effect is at the level of 0.01 − 0.02, depending on the q 2 bin. This uncertainty constitutes up to 30% of the statistical uncertainty on F L and 20% of the statistical uncertainty on A FB . For S 5 and A 5 , the largest source of systematic uncertainty comes from the choice of polynomial order for the angular acceptance. If polynomials two orders higher are used, a variation of ∼ 0.01 is observed. For the remaining CP -averaged and CPasymmetric observables, the uncertainties arising from the data-simulation comparison and the acceptance are small. However, there are three other non-negligible sources of systematic uncertainty. Throughout the full q 2 range, peaking backgrounds introduce a systematic uncertainty at the level of 0.01 or less. For the likelihood fit of the angular observables, in the first two q 2 bins (where the acceptance changes most rapidly), the uncertainty arising from using the bin centre, as opposed to a bin edge, is at the level of 0.01 or less. Finally, at high q 2 , the statistical precision on the acceptance correction leads to a systematic uncertainty at the level of 0.01 or less. For the P ( ) i observables, the situation is more complex and the systematic uncertainty is shared more evenly between the different sources (see Table 2). The dominant sources of systematic uncertainty can all be reduced in future analyses with larger datasets.
Propagating the above sources of systematic uncertainty to the zero-crossing points yields uncertainties at the level of 0.07 GeV 2 /c 4 for S 4 , 0.02 GeV 2 /c 4 for S 5 and 0.03 GeV 2 /c 4 for A FB . These uncertainties are negligible compared to the statistical uncertainties.

Results
The CP -averaged observables that are obtained from the likelihood fits are shown together with the SM predictions in Fig. 6. The CP asymmetries are shown in Fig. 7. The SM predictions are based on the prescription of Ref. [19]. In contrast to the alternative SM predictions that are available in Refs. [20,[57][58][59][60][61][62], these predictions update the calculations from Ref. [63] to account for the known correlation between the different form factors [64]. Light-cone sum rule predictions, which are valid in the low-q 2 region, are also combined with lattice determinations at high q 2 [65,66] to yield more precise determinations of the form factors over the full q 2 range. The predictions are made in the regions 0.1 < q 2 < 6.0 GeV 2 /c 4 and 15.0 < q 2 < 19.0 GeV 2 /c 4 . No predictions are made for the region close to the narrow cc resonances, the J/ψ and ψ(2S), where many of the assumptions that go into the SM predictions are thought to be invalid. Reference [19] does not make predictions for the S 7,8,9 and the A i observables. These observables are all expected to be close to zero in the SM.
The results of the fits for the optimised angular observables are shown together with their SM predictions in Fig. 8. For the P ( ) i observables, predictions from Ref. [14] are shown using form factors from Ref. [67]. The SM predictions are restricted to the q 2 range q 2 < 8.0 GeV 2 /c 4 . The variation in the size of the uncertainties on the P ( ) i observables arises from their dependence on F L . When F L is large, the (1 − F L ) term in the definition of the observables gives a large uncertainty that can be significantly asymmetric.
The results of the likelihood fit for the angular observables are given in Tables 3-6 of Appendix A, with the statistical and systematic uncertainties separated. In general, the correlations between the observables are small. The most notable exceptions are the correlations between A FB and F L , which can be as large as 60%, and the correlations between the different P ( ) i observables in the range 2.5 < q 2 < 4.0 GeV 2 /c 4 . The correlations between A FB and F L arise from the requirement that the differential decay rate in Eq. (4) be positive across the entire phase space. The correlations between the P ( ) i observables originate from their common dependence on F L . The correlation matrices for all of the q 2 bins are available in Appendices C, D and E. The values of F S obtained from the fits are consistent with the S-wave contribution of approximately 5% observed in B 0 → J/ψ K * 0 data [41-43]. Considering the observables individually, the results appear largely in agreement with the SM predictions. The exception to this is the observable S 5 and the related observable P 5 . Small differences can also be seen in the measured A FB distribution, where the data lie systematically below the SM predictions in the range 1.1 < q 2 < 6.0 GeV 2 /c 4 . No significant CP asymmetry is seen.
The discrepancy in P 5 confirms the result of the previous LHCb analysis [8], where a difference was seen between the data and the SM predictions in the q 2 range 4.30 < q 2 < 8.68 GeV 2 /c 4 . In the present analysis, a deviation from the SM prediction is observed in each of the 4.0 < q 2 < 6.0 GeV 2 /c 4 and 6.0 < q 2 < 8.0 GeV 2 /c 4 bins at a level of 2.8 and 3.0 standard deviations, respectively. The SM predictions for the optimised observables that are used in this analysis are taken from Ref. [14]. The predictions are an update of the SM calculation from Ref. [68], which was used to compare the previous LHCb P The results of the moment analysis are shown in Figs. 9, 10 and 11 and given in Tables 7, 8 and 9 of Appendix A. The same behaviour is seen as in the likelihood fit, where some differences are observed between the SM predictions and the data in S 5 (and P 5 ) at low values of q 2 . The observable S 6c is also included in Table 7 and shown in Fig. 12.  Figure 6: The CP -averaged observables in bins of q 2 , determined from a maximum likelihood fit to the data. The shaded boxes show the SM predictions based on the prescription of Ref. [19].     This observable is consistent with zero, as expected in the SM. As a cross-check, the observables have also been determined by a moment analysis in the approximately 2 GeV 2 /c 4 q 2 bins used in the likelihood fit of the angular observables. The differences between the central values of the two methods are compatible with those expected from pseudoexperiments. The correlation matrices for all of the q 2 bins are available in Appendices F-H.    Figure 13: The observables S 4 , S 5 and A FB determined by fitting for the q 2 dependent decay amplitudes. The line indicates the best-fit to the dataset. The band indicates the 68% interval on the bootstraps at each point in q 2 . Note that, the correlation between points in the bands means it is not possible to extract the uncertainty on the zero-crossing points from these figures.

Compatibility with the Standard Model
The EOS software package [57] is used to determine the level of compatibility of the data with the SM. It provides predictions for the observables integrated over the q 2 bins used in the analysis. A χ 2 fit is performed to the CP -averaged angular observables F L , A FB and S 3 -S 9 obtained from the likelihood fit to the data. The χ 2 fit uses observables in the range q 2 < 8.0 GeV 2 /c 4 and a wide q 2 bin covering the range 15.0 < q 2 < 19.0 GeV 2 /c 4 . Previous analyses [9-12, [19][20][21] have shown that the existing measurements of decays involving a b → s quark transition, including the previous LHCb P 5 result from Ref. [8], can be accounted for by modifying only the real part of the vector coupling strength of the decays, conventionally denoted Re(C 9 ). An analysis considering additional effective couplings would require a global fit to all of the measurements of b → s quark transitions and is beyond the scope of this paper. Note that modifying just the axial-vector coupling strength, C 10 , would lead to a branching fraction for the B 0 s → µ + µ − decay that is excluded by existing measurements [72].
In the χ 2 fit, the correlations between the different observables are taken into account. The floating parameters are Re(C 9 ) and a number of nuisance parameters associated with the form factors, CKM elements and possible sub-leading corrections to the amplitudes. The sub-leading corrections to the amplitudes are expected to be suppressed by the size of the b-quark mass relative to the typical energy scale of QCD. The nuisance parameters are treated according to the prescription of Ref. [11] and are included in the fit with Gaussian constraints. In the χ 2 minimisation procedure, the value of each observable (as derived from a particular choice of the theory parameters) is compared to the measured value. Depending on the sign of the difference between these values, either the lower or upper (asymmetric) uncertainty on the measurement is used to compute the χ 2 .
The minimum χ 2 corresponds to a value of Re(C 9 ) shifted by ∆Re(C 9 ) = −1.04 ± 0.25 from the SM central value of Re(C 9 ) = 4.27 [11] (see Fig. 14). From the difference in χ 2 between the SM point and this best-fit point, the significance of this shift corresponds to 3.4 standard deviations. As discussed in the literature [9-12, 14-21], a shift in C 9 could be caused by a contribution from a new vector particle or could result from an unexpectedly large hadronic effect.
If a fit is instead performed to the CP -averaged observables from the moment analysis in the same q 2 ranges, then ∆Re(C 9 ) = −0.68 ± 0.35 is obtained. As expected, the uncertainty on ∆Re(C 9 ) is larger than that from the likelihood fit. Taking into account the correlations between the two methods, the values of ∆Re(C 9 ) are statistically compatible. LHCb SM Figure 14: The ∆χ 2 distribution for the real part of the generalised vector-coupling strength, C 9 . This is determined from a fit to the results of the maximum likelihood fit of the CP -averaged observables. The SM central value is Re(C SM 9 ) = 4.27 [11]. The best fit point is found to be at ∆Re(C 9 ) = −1.04 ± 0.25.

Conclusions
This paper presents the first analysis of the full angular distribution of the B 0 → K * 0 µ + µ − decay. The analysis uses the complete LHCb Run 1 dataset and supersedes the results presented in Refs. [1,8]. In addition to CP -averaged observables, a complete set of CP asymmetries of the angular distribution are measured for the first time. Correlations between the different observables are computed to allow the results to be included in global fits of b → s data.
Three separate techniques are used to analyse the data. An unbinned maximum likelihood fit to the full angular distribution is made in approximately 2 GeV 2 /c 4 wide q 2 bins. Observables are also determined by computing moments of the angular distribution in q 2 bins approximately 1 GeV 2 /c 4 wide. In addition, for the first time, a q 2 -dependent fit is performed to the angular distribution in order to determine the six complex decay amplitudes that describe the decay. The position in q 2 at which several observables cross zero is determined using these amplitudes.

Appendices A Tables of results
The results of the likelihood fits described in Sec. 7.1 are given in Tables 3-6 below. The results of the method of moments described in Sec. 7.2 are given in Tables 7-9 below. Table 3: CP -averaged angular observables evaluated by the unbinned maximum likelihood fit, in the range 1.1 < q 2 < 6.0 GeV 2 /c 4 and 15.0 < q 2 < 19.0 GeV 2 /c 4 . The first uncertainties are statistical and the second systematic.

C Correlation matrices for the CP -averaged observables from the maximum likelihood fit
Correlation matrices between the CP -averaged observables in the different q 2 bins are provided in Tables 10-19 for the likelihood fit. 1.00 1.00 1.00 Table 13: Correlation matrix for the CP -averaged observables from the maximum likelihood fit in the bin 4.0 < q 2 < 6.0 GeV 2 /c 4 . 1.00 Table 14: Correlation matrix for the CP -averaged observables from the maximum likelihood fit in the bin 6.0 < q 2 < 8.0 GeV 2 /c 4 . 1.00 1.00 1.00 1.00 Table 18: Correlation matrix for the CP -averaged observables from the maximum likelihood fit in the bin 1.1 < q 2 < 6.0 GeV 2 /c 4 . 1.00 1.00 D Correlation matrices for the CP -asymmetric observables from the maximum likelihood fit Correlation matrices between F L and the CP -asymmetric observables in the different q 2 bins are provided in Tables 20-29 for the likelihood fit. 1.00 1.00 Table 22: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 2.5 < q 2 < 4.0 GeV 2 /c 4 . 1.00 Table 23: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 4.0 < q 2 < 6.0 GeV 2 /c 4 . 1.00 Table 24: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 6.0 < q 2 < 8.0 GeV 2 /c 4 . 1.00 Table 25: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 11.0 < q 2 < 12.5 GeV 2 /c 4 . 1.00  1.00 Table 27: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 17.0 < q 2 < 19.0 GeV 2 /c 4 . 1.00 Table 28: Correlation matrix for the CP -asymmetric observables from the maximum likelihood fit in the bin 1.1 < q 2 < 6.0 GeV 2 /c 4 . 1.00 1.00 E Correlation matrices for the optimised angular observables from the maximum likelihood fit Correlation matrices between F L and the optimised P ( ) i basis of observables in the different q 2 bins are provided in Tables 30-39 for the likelihood fit.  1.00 −0.10 P 8 1.00 1.00 0.08 P 8 1.00  1.00 0.10 P 8 1.00      1.00  1.00 1.00  1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 G Correlation matrices for the CP -asymmetric observables from the method of moments Correlation matrices between the CP asymmetries in the different q 2 bins are provided in Tables 55-69 for the moment analysis. The correlations are determined by a bootstrapping technique. 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 −0.08 A 9 1.00 1.00 1.00 1.00 1.00 −0.08 A 9 1.00 1.00 1.00