Determination of short-and long-distance contributions in B 0 → K ∗ 0 µ + µ − decays

,


Introduction
Flavour-changing neutral current decays of beauty hadrons, mediated through electroweak loop diagrams, constitute powerful indirect probes of the Standard Model (SM).Several studies of these processes have shown intriguing discrepancies with respect to predictions from the SM, notably in final states with a pair of muons [1][2][3][4][5][6].In particular, angular analyses of B 0 → K * 0 µ + µ − decays1 [7][8][9][10][11][12] reported a tension in the measured value of the angular observable P ′ 5 [13] with respect to its calculations based on the SM [14][15][16].In order to further investigate these anomalies, it is convenient to describe rare b → s decays within an effective field theory in terms of local (short-distance) operators, O i , that only contain light SM fields, and corresponding effective couplings, C i , known as Wilson coefficients.Only a subset of dimension-six operators is relevant for B 0 → K * 0 µ + µ − decays [17]: the radiative operator O 7 , the semileptonic operators O 9 and O 10 that drive the vector and axial leptonic currents, respectively, and their corresponding right-handed counterparts O ′ 7,9,10 .Within this approach, all non-perturbative QCD effects are contained in the local (form factors) and non-local hadronic matrix elements, while physics beyond the SM would appear as a shift in the values of one or more Wilson coefficients.The lack of knowledge on these hadronic matrix elements prevents a clean determination of the Wilson coefficients.Non-local contributions originating from b → scc virtual loops are particularly difficult to assess reliably from first principles as they receive non-perturbative corrections from soft-gluon emissions [14].Over the last two decades, different approaches have been proposed to determine the impact of these contributions, from phenomenological [18][19][20] and data [21,22] analyses to the use of analytic [16,[23][24][25] and dispersion [26] relations, where all these methods rely on the perturbative treatment of four-quark operators in the regime of q 2 ≪ 4m 2 c , with q 2 the dimuon invariant mass squared and m c the mass of the charm quark.While the study and development of all these techniques enabled a clear progress in the understanding of such contributions, a definitive consensus on their impact has not been reached yet [27].A deeper comprehension of these non-local hadronic effects is therefore crucial for the understanding of the anomalies seen in rare b → s hadron decays.
The analysis described in this paper aims at the first determination of such non-local hadronic contributions in B 0 → K * 0 µ + µ − decays from data by performing a q 2 -unbinned amplitude analysis of the decay.Non-local hadronic contributions are determined following two alternative strategies, the first relies purely on experimental data and the second includes state-of-the-art theoretical inputs.The four Wilson coefficients C 9 , C 10 , C ′ 9 and C ′ 10 are determined simultaneously from data together with the local and non-local hadronic matrix elements, while the coefficients C 7 and C ′ 7 are already precisely known from radiative B meson decays [28] and are fixed to their SM values [29,30].This measurement therefore provides a complementary and more in-depth set of information with respect to previous binned analyses of the same decay [1,9].A more concise description of this work is reported in a companion article [31].The employed dataset corresponds to an integrated luminosity of 4.7 fb −1 of pp collisions collected with the LHCb experiment during the years 2011, 2012 and 2016 at centre-of-mass energies of 7, 8 and 13 TeV, respectively.This paper is structured as follows: the formalism employed to describe B 0 → K * 0 µ + µ − decays is presented in Sec.2; the LHCb detector and the procedure used to generate simulated samples is introduced in Sec.3; the selection of signal candidates is described in Sec.4; the amplitude fit method is discussed in Sec.5; the fit to data is presented in Sec.6; sources of systematic uncertainties are discussed in Sec.7; Sec. 8 shows the obtained results; Sec. 9 is dedicated to a final discussion and conclusions are presented in Sec.10.

Formalism
The B 0 → K * 0 µ + µ − decay, where the K * 0 meson is reconstructed through the decay K * 0 → K + π − , can be described by five kinematic variables: the dimuon invariant mass squared, q 2 , the K + π − invariant mass squared, k 2 , and the three decay angles ⃗ Ω = (θ ℓ , θ K , ϕ).Here, θ ℓ is the angle between the µ + (µ − ) and the direction opposite to that of the B 0 (B 0 ) in the rest frame of the dimuon system, θ K is 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, and ϕ is 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.A full description of the angular basis is provided in Ref. [7].The differential decay rate of B 0 → K * 0 (→ K + π − )µ + µ − decays, where the K + π − pair comes from an intermediate resonance of spin-1 (P-wave), can be written as [8,32] I 1s sin 2 θ K + I 1c cos 2 θ K + I 2s sin 2 θ K cos 2θ ℓ + I 2c cos 2 θ K cos 2θ ℓ + I 3 sin 2 θ K sin 2 θ ℓ cos 2ϕ + I 4 sin 2θ K sin 2θ ℓ cos ϕ + I 5 sin 2θ K sin θ ℓ cos ϕ + I 6 sin 2 θ K cos θ ℓ + I 7 sin 2θ K sin θ ℓ sin ϕ + I 8 sin 2θ K sin 2θ ℓ sin ϕ + I 9 sin 2 θ K sin 2 θ ℓ sin 2ϕ , where the I i are the q 2 -dependent angular coefficients defined in Appendix A.1, which can be conveniently expressed in terms of the amplitudes A L,R λ and A t .The former are transversity amplitudes where the symbol λ = 0, ⊥, ∥ refers to the transversity state of the K * 0 meson and the indices L and R denote the chirality of the lepton current.The amplitude A t corresponds to a time-like polarisation of the virtual vector boson decaying into a dimuon pair and longitudinal polarisation of the K * 0 meson [32].These amplitudes describe the decay process and can be expressed as [17,23] where the functions F where µ /q2 , V tb and V * ts are elements of the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix, G F is the Fermi constant and α e is the fine structure constant.The exact definition of the form factors F (T ) λ (q 2 , k 2 ) is given in Appendix A.2, while the definition of the non-local functions H λ (q 2 ) follows what has been proposed in Refs.[16,[23][24][25] where the analytic properties of the hadronic matrix elements are exploited through the mapping [33,34] where t + = 4M 2 D , with M D the D 0 meson mass, and t 0 can be arbitrarily chosen such that z(q 2 = t 0 ) = 0.After this transformation, 2 the non-local hadronic functions can be expressed as where the first and second terms remove the singularities due to the J/ψ and ψ(2S) poles.The Ĥλ (z) are analytic functions which can be further decomposed as where ϕ −1 λ (z) are so-called outer functions that ensure the correct kinematic dependence [24], e.g.H 0 (q 2 = 0) = 0, and α λ,k are the coefficients of a polynomial expansion.
The K + π − system in the final state can also appear in a scalar (S-wave) configuration, which introduces two additional amplitudes [35], with three new form factors, f + , f T and f 0 whose definitions can be found in Appendix A.2.
In the following, contributions from non-local hadronic matrix elements to the scalar amplitudes are ignored.This assumption is studied as a source of systematic uncertainty in Sec. 7.
In order to better separate the contribution of the S-wave amplitudes from those of the P-wave, the k 2 dependence is included in the model.This is achieved by multiplying the decay amplitudes of Eq. 2 by a relativistic Breit-Wigner function for the resonant P-wave [36] and the scalar amplitudes of Eq. 7 by the LASS parameterisation [37] for the slowly varying S-wave, i.e.
where f P and f S encode the P-and S-wave k 2 dependence, respectively.In principle, the proportion between P-and S-waves can be determined from the normalisation of the decay amplitudes, however, an accurate prediction of the two relative contributions involves a complete analysis of the B 0 → K + π − form factors in the full K + π − spectrum [36,38].
Given the limited k 2 range analysed around the K * 0 mass, it is therefore practical to decouple the normalisation of the P-and S-wave amplitudes and introduce an additional complex coefficient to control the relative magnitude and phase between the two.A detailed definition of f P and f S is given in Appendix A.3.Finally, when taking into account the full set of P-wave and S-wave amplitudes, the total B 0 → K + π − µ + µ − differential decay rate reads as 32π 9 + Ĩ4 sin 2θ ℓ + Ĩ5 sin θ ℓ sin θ K cos ϕ where the first term corresponds to the P-wave differential decay rate of Eq. 1 extended to the k 2 dimension via Eq. 8, I S i are pure S-wave angular coefficients and Ĩi denote interference terms between the S-and P-wave amplitudes, as defined in Appendix A.1.

LHCb detector and simulation
The LHCb detector [39,40] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks.The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet.The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c.The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors.Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.The 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 the hardware stage, depending on the data-taking conditions signal candidates are required to have at least one muon with a p T greater than 1 to 2 GeV/c or a pair of muons with the product of their p T above 1 to 4 GeV 2 /c 2 .The software trigger requires a two-, threeor four-track secondary vertex with a significant displacement from any PV.At least one charged particle must have p T greater than 1 GeV/c and be inconsistent with originating from a PV.A multivariate algorithm [41] is used for the identification of secondary vertices consistent with the decay of a b-hadron.
Simulation is required to model the effects of the detector acceptance and the imposed selection requirements.In the simulation, pp collisions are generated using Pythia [42] with a specific LHCb configuration [43].Decays of unstable particles are described by EvtGen [44], in which final-state radiation is generated using Photos [45].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [46] as described in Ref. [47].
Corrections derived from the data are applied to the simulation to account for mismodelling of the charge multiplicity of the event, the B 0 momentum spectrum and the B 0 vertex quality.Similarly, the simulated particle identification (PID) performance is corrected to match that determined from control samples selected from the data [48,49].

Selection of signal candidates
The same dataset analysed in Ref. [9] is considered for this measurement, as well as the same simulated samples and data-simulation corrections.Signal candidates are formed from a pair of oppositely charged tracks that are identified as muons, combined with a K * 0 candidate, formed from two tracks identified as a kaon and a pion, respectively.Signal candidates are required to have a reconstructed K + π − µ + µ − invariant mass, m Kπµµ , in the range [5170, 5700] MeV/c 2 and a reconstructed mass of the K + π − system, m Kπ , in a window of ±100 MeV/c 2 around the known mass of the K * 0 meson [50].In addition, the following requirements are imposed, following the strategy in Ref. [9].The tracks of the four final-state particles must have a significant impact parameter with respect to all PVs in the event and must form a good-quality vertex.The B 0 candidate must be compatible with being produced in one of the PVs, i.e. the IP is small, and its decay vertex must be significantly displaced from the identified PV.The angle between the reconstructed B 0 momentum and the vector connecting the identified PV to the reconstructed B 0 decay vertex is also required to be small.Two types of backgrounds are considered: combinatorial background, where the selected final-state particles do not originate from a single b-hadron decay; and peaking backgrounds, where specific background processes can mimic the signal if their final-state particles are misidentified or misreconstructed.Combinatorial background events are further reduced using a boosted decision tree (BDT) algorithm [51,52].The BDT classifier used for the present analysis is identical to that described in Ref. [9] and the same working point is used.The BDT selection rejects more than 97% of the remaining combinatorial background, while retaining more than 85% of the signal.
Peaking backgrounds arise from decays of B 0 s → ϕ(1020)(→ K + K − )µ + µ − and Λ 0 b → pK − µ + µ − where a kaon or proton is misidentified as a pion, and B 0 → K * 0 µ + µ − in which the kaon and pion from the K * 0 decay are both misidentified.These contributions are rejected by applying dedicated vetoes combining PID information with the invariant mass of the candidates recomputed with the relevant change in the particle mass hypotheses.Cabibbo-suppressed decays of B 0 s → K * 0 µ + µ − , which have the correct final state, contribute at the level of 1% of the signal [53] and are neglected when considering the B 0 → K * 0 µ + µ − decay.Decays of B + → K + µ + µ − can pass the above selection if a low-momentum pion from the rest of 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 fall into the upper-mass sideband.Such decays can therefore distort the estimate of the residual combinatorial background distributions which are assessed from this sideband and are suppressed by removing candidates with m Kµµ within 60 MeV/c 2 (corresponding to approximately three times the detector resolution) around the known B + mass [50].The residual contamination from these peaking backgrounds is studied using simulated events and is estimated to be at the level of 1% or less.Therefore, the peaking backgrounds are sufficiently small to be neglected in the analysis and are considered only as sources of systematic uncertainty.

Amplitude fit method
The q 2 spectrum of B 0 → K * 0 µ + µ − decays is shaped by a variety of contributions, as illustrated in Fig. 1.In this analysis, signal candidates are only considered in two q 2 regions, i.e. [1.1, 8.0] and [11, 12.5] GeV 2 /c 4 .The region below 1.1 GeV 2 /c 4 is not investigated due to the presence of light-quark resonances, e.g.η, ρ 0 (770), ω(792) and ϕ(1020), that would need to be explicitly modelled in the fit.Similarly, the q 2 region above the ψ(2S) resonance is populated by broad charmonium states, e.g.ψ(3770), ψ(4040), that lie beyond the region of validity of the H λ (q 2 ) parameterisation described in Sec. 2. Future theoretical developments on the treatment of non-local hadronic matrix elements are necessary in order to incorporate such contributions in the analysis.On the contrary, narrow charmonium resonances appearing in B 0 → ψ n K * 0 tree-level decays, with ψ n = J/ψ or ψ(2S) and ψ n → µ + µ − , are embedded in the described formalism and appear as poles of the non-local hadronic parameterisation as shown in Eq. 5. Regions with q 2 ∈ [8.0, 11.0] and [12.5, 14.5] GeV 2 /c 4 , which are dominated by these resonant decays, are removed from the signal selection and considered only as control regions for the analysis.
Since B 0 and B 0 decays are treated together in this analysis, no sensitivity to the imaginary part of the Wilson coefficients can be achieved, and thus all these coefficients are assumed to be real without loss of generality.The parameters of interest are determined using an extended unbinned maximum-likelihood fit.The fit is performed simultaneously on the 2011 and 2012 (referred to as Run 1) and 2016 datasets and on the two q 2 regions of the analysis, with all the physics signal parameters shared in the fit.
10 narrow cc broad cc and D D thresholds The dashed line corresponds to the pure rare semileptonic decay, while the solid line includes the impact of different charmonium resonances.The gray bands correspond to the regions dominated by B 0 → J/ψK * 0 and B 0 → ψ(2S)K * 0 tree-level decays, which are selected as control channels.Magnitudes and phases of cc resonant components have been arbitrarily chosen for illustrative purposes.The dominant Wilson coefficients in each region of the spectrum are also highlighted for reference.

Efficiency correction
The reconstruction and selection of signal candidates distort the distributions of the decay angles and q 2 .To account for this acceptance effect, a four-dimensional acceptance function is introduced according to k,l,m,n where L m (x) are Legendre polynomials in x of order m, and q 2 is parameterised in the range between 1 and 14 GeV 2 /c 4 , slightly beyond the q 2 domain of this analysis in order to avoid border effects.The ϕ angle and q 2 are each rescaled to the range [−1, 1] when evaluating the Legendre polynomial.The acceptance is modelled using polynomials of the lowest order that show a good description of the detector response; Legendre polynomials up to order three are used for q 2 ; for the decay angles, polynomials up to order five, four and six are used for cos θ K , cos θ ℓ and the angle ϕ, respectively.Only even orders are used for the angle ϕ.The coefficients c klmn are determined separately for the Run 1 and 2016 datasets by exploiting the orthogonality of the Legendre polynomials on large samples of simulated three-body B 0 → K * 0 µ + µ − phase space decays.Resolution effects are effectively taken into account by evaluating the acceptance function on the reconstructed angles and q 2 after final-state radiation effects.The impact of the resolution on k 2 , θ ℓ , θ K and ϕ is minimal.The acceptance parameterisation is also found to be independent of k 2 in the range considered in this analysis.In addition, as the acceptance is parameterised in terms of all the kinematic variables needed to describe the decay, it does not depend on the model used in the simulation.The validity of the obtained acceptance parameterisation is tested on a large set of simulated signal events generated based on the physics model provided by EvtGen [44,54].The differential decay rate of Eq. 1 is multiplied by the acceptance function of Eq. 10 and the extracted signal parameters are found to be in good agreement with the generated values.

Signal modelling
The parameterisation of the signal is built from the five-dimensional differential decay rate of Eq. 9 multiplied by the acceptance function of Eq. 10.The signal yields are directly extracted from the m Kπµµ distribution, which is included as an additional dimension in the fit to separate signal from background.The probability density function (p.d.f.) of m Kπµµ for the signal is parameterised with the sum of two Gaussian functions with power-law tails on both sides of the peak.The two Gaussian functions share the same mean and tail parameters but are allowed to have different widths.All the parameters of the signal m Kπµµ model are fixed from the simulation with the exception of the mean and widths, which receive a correction from a fit to the m Kπµµ invariant mass distribution of the B 0 → J/ψK * 0 control channel.The resulting six-dimensional signal p.d.f.depends on the four Wilson coefficients C (′) 9 and C (′) 10 , which are varied freely in the fit, and on a large number of parameters that are constrained to different external inputs as described in the next subsections.These are the 4 CKM parameters, 19 (9) local P-wave (S-wave) form factor parameters and between 18 and 30 non-local hadronic parameters depending of the truncation order of H λ (z).External constraints are applied in the fit by multiplying the likelihood function by multi-dimensional Gaussian functions.Finally, the two parameters describing the relative magnitude and phase between the P-and S-waves, as introduced in Appendix A.3, are also free to vary in fit.

External inputs
Several external inputs are used to constrain different parts of the signal decay amplitudes.The CKM elements V tb V * ts are expressed through the Wolfenstein parameterisation, whose parameters {λ, A, ρ, η} are constrained to the values obtained from an SM fit of the unitarity triangle [55], as summarised in Table 1.
Local form factors can be assessed either from first principles through lattice QCD (LQCD) simulations [56,57], or from quark-hadron-duality arguments through QCD light-cone sum rules (LCSR) [58,59].In this analysis, B 0 → K * 0 form factors F λ (q 2 ) are constrained to the combination of LCSR [16,60] and LQCD [61] results provided by the Eos software package [62], with the correlations amongst the form factor parameters fully taken into account.The current knowledge on the B 0 → [Kπ] J=0 scalar form factors is very limited and only recently these effects have been investigated in the context of [38].The resulting q 2 dependence is found to be almost entirely driven by the kinematic factors.In the following, the S-wave form factors f + , f T and f 0 are constrained to those obtained for B + → K + transitions from LQCD [63], with uncertainties inflated by a factor of three to compensate for the different meson species.This assumption is further studied as a source of systematic uncertainty with alternative parametric expressions [64] and it is found to be consistent with the numerical analysis of Ref. [38].
While the determination of the non-local hadronic contributions H λ (q 2 ) is one of the key aspects of this analysis, additional external inputs can be used to constrain these elements.First, experimental measurements of the decays B 0 → ψ n K * 0 can provide information on the magnitude and phase of the non-local matrix elements evaluated at the pole of the corresponding resonance.The amplitude of these decays can be expressed in terms of the residues of the functions H λ (q 2 ) at the ψ n poles as [23] Res where the M ψn and f * ψn are masses and decay constants for the ψ n resonances [14,50,65], respectively, and A ψn λ are the transversity amplitudes of the tree-level decays to charmonia.These amplitudes are experimentally related to the measured polarisation fractions and phase differences [66][67][68][69] and to branching fraction measurements [67,68,70] where M K * 0 is the pole mass of the K * 0 meson and V cb and V * cs are elements of the CKM matrix.Table 2 collects all the external inputs from B 0 → ψ n K * 0 measurements used to constrain the values of H λ (q 2 ) in the analysis.
Theoretical predictions of the ratio H λ /F λ can be obtained at q 2 < 0 GeV 2 /c 4 .In this regime, the effects of four-quark operators can be treated perturbatively and SM theory predictions can be accessed by a light-cone operator-product expansion [14].Precise predictions of the real and imaginary parts of the ratio H λ /F λ are provided by Refs.[16,24] at three q 2 points, i.e. q 2 ∈ {−7, −5, −3} GeV 2 /c 4 .While information at q 2 = −1 GeV 2 /c 4 is also available, in this measurement it is used instead as a test point to check the consistency of the result of the fit.In order to explore the impact and the consistency of these theory predictions, two fit configurations are considered in the analysis: the first includes external constraints to the above-mentioned theoretical prediction on the values of H λ /F λ at negative q 2 , hereafter referred to as "q 2 < 0 constraints"; while the second ignores those constraints and provides a pure data-driven determination of the non-local hadronic matrix elements, referred to as "q 2 > 0 only".
Table 2: Summary of external inputs from B 0 → ψ n K * 0 measurements used in the analysis.When measurements from different experiments are available, these are averaged taking into account correlations.A shift of +π is introduced for the phase differences δ ψn ⊥,∥ to account for the different convention between this analysis and Refs.[66,68,69].

Constraints from branching fraction determination
The estimation of the signal yield can be easily incorporated in the amplitude analysis by performing an extended unbinned maximum likelihood fit.In addition, the observed signal yield in the k 2 and q 2 mass windows can be conveniently related to the signal branching fraction via where τ B is the lifetime of the B 0 meson and the two-dimensional differential decay rate is obtained by integrating the full P-wave decay rate over the three angles.The factor ( 2 3 ) is the squared Clebsch-Gordan coefficient related to K * 0 → Kπ that corresponds to the K * 0 → K + π − decay probability, while the term F P , built from the fraction of the P-wave over the total P-and S-wave decay rates, is introduced to correct for the S-wave contribution observed in the signal region.Similarly, the observed yield in the control channel, N J/ψKπ , which is determined from a fit to the K + π − µ + µ − invariant mass distribution as shown in Fig. 2, includes a combination of P-and S-wave as well as contributions of exotic nature, such as B 0 → Z c (4430) − K + and B 0 → Z c (4200) − K + decays, with Z − c decaying to J/ψπ − .Given that all these contributions share the same final state, they appear as single component in the fit.The observed N J/ψKπ yield is found to be 348 500 ± 600 and 328 500 ± 600 for Run 1 and 2016, respectively, where uncertainties are statistical only, while contamination from B 0 s → J/ψK * 0 decays are found to be of the order of 1%.Since an exact separation between the P-wave, S-wave and exotic components can only be achieved by a dedicated amplitude analysis that comprises the J/ψπ − invariant mass dimension [67], the signal branching fraction of Eq. 15 is normalised to the inclusive branching fraction of B 0 → J/ψK + π − decays.This is taken from Ref. [67]  and is scaled by a correction factor, f B 0 →J/ψKπ ±100MeV , which takes into account the fraction of B 0 → J/ψK + π − decays that fall into the m Kπ window of this analysis.This correction factor is estimated with a series of pseudoexperiments generated from the amplitude model of Ref. [67].Unfortunately, neither the correlations nor the systematic uncertainties associated to the individual amplitudes are provided in Ref. [67] and f B 0 →J/ψKπ ±100MeV can only be evaluated under the conservative assumption of uncorrelated uncertainties.Within this hypothesis, f B 0 →J/ψKπ ±100MeV is found to be 0.644 ± 0.010, whose central value is fixed in the fit and the associated uncertainty is treated as a source of systematic uncertainty.Finally, the charmonium branching fraction, B(J/ψ → µ + µ − ), is taken from Ref. [50] and the relative efficiency between the signal and control modes, R ε , is obtained from simulated samples.The use of simulation to model the efficiency is validated on the B 0 → ψ(2S)K * 0 control channel, where the ratio of branching fractions B(B 0 → ψ(2S)K * 0 )/B(B 0 → J/ψK * 0 ) is measured precisely and found to be consistent with the value of Ref. [50].
It is important to note that the constraint on the branching fraction determination sets the scale of the Wilson coefficients.Without this constraint, |C 9 | and |C 10 | would not be bounded and only ratios of Wilson coefficients could be assessed.

Background modelling
The combinatorial background that passes the selection requirements of Sec. 4 must be modelled in the fit.These candidates are described by an exponential function in m Kπµµ , by Chebyshev polynomials of up-to second order for the decay angles and q 2 , and by the sum of a linear function and a Breit-Wigner amplitude-squared for k 2 .The Breit-Wigner model accounts for background candidates originating from a real K * 0 resonance associated with an unrelated µ + µ − pair.The ratio between this resonant component and the linear term is fixed to the value observed in the upper-mass sideband.In addition, a sculpting effect is introduced by the veto against B + → K + µ + µ − decays, strongly distorting a specific region of the (cos θ K , q 2 , m Kπµµ ) phase space.This effect is studied using the data: events reconstructed in the lepton-flavour-violating channel B 0 → K * 0 µ ± e ∓ are used to study the impact of the veto on the variables of interest in absence of peaking B + → K + µ ± e ∓ backgrounds.A three-dimensional efficiency map is extracted from the ratio of templates obtained with and without the veto selection applied, and is multiplied to the background shape as a correction factor to retrieve the original distributions.Besides this term, all remaining functions are assumed to factorise.This assumption, together with the choice of the order of the polynomials, is investigated as a source of systematic uncertainty in Sec. 7. Finally, the parameterisation of the combinatorial background is treated separately for the Run 1 and 2016 datasets as well as for the two considered q 2 regions [1.1, 8.0] and [11.0, 12.5] GeV 2 /c 4 .All coefficients are allowed to vary in the fit.

Fit to data
The polynomial expansion introduced in Eq. 6 used to parameterise the non-local hadronic matrix elements H λ must be truncated at a certain order z n .Furthermore, the central point of the expansion t 0 is a free parameter of the model and its choice can have an impact on how fast the polynomial expansion converges.In general, a sensible choice is a value of t 0 that minimises |z| on the domain of the expansion.As originally proposed by Ref. [23], the choice of is the one that minimises |z| in the interval −7 GeV 2 /c 4 ≤ q 2 ≤ M 2 ψ(2S) ; this value is taken as the default for the fit configuration with q 2 > 0 information only.However, due to the strong constraints provided by the three q 2 points, t 0 is fixed to zero for the fit model with the q 2 < 0 constraints in order to best accommodate the theoretical inputs in the negative q 2 region.Following this choice, the truncation order z n is determined based on a data-driven procedure: fits are repeated with increasing truncation order for the polynomial sums, i.e. n = 2, 3, 4, 5, and the Akaike information criterion [71] is used to infer the importance of each additional set of coefficients.The improvement between two subsequent orders is considered to be significant if 2∆ log L > 2∆N pars , where N pars is the number of parameters of the model and each additional order z n+1 brings one complex coefficient for each of the three polarisations, for a total of six additional parameters.For the fit model using only inputs at q 2 > 0, it is found that a polynomial expansion truncated at z 2 is sufficient to describe the data.For fits with q 2 < 0 constraints, a significant improvement is found with the inclusion of terms up to z 4 , as shown in Table 3.The quality of the fit is assessed using an unbinned goodness-of-fit test based on point-to-point dissimilarity methods [72] and the p-value is found to be larger than 10%.
Figure 3 shows the distributions of events for the Run 1 and 2016 combined datasets.The result of the fit without the q 2 < 0 constraints are overlaid in the figure, no difference in the projections is visible when including q 2 < 0 constraints.

Systematic uncertainties
There are different categories of systematic uncertainties that affect the extraction of the parameters of interest, from the choice of the baseline amplitude model and the inclusion of external inputs in the fit, to imperfect modelling of experimental effects.The distinct sources of systematic uncertainties are discussed in detail below and are summarised in Table 3: Log-likelihood differences between the fits to data with different truncation orders of the non-local hadronic parameterisation H λ [z n ] for the two considered fit configurations.
2∆ log L q 2 < 0 constraints q 2 > 0 only 8.6 -Table 4. The size of each systematic uncertainty is estimated using pseudoexperiments generated using the observed signal and background yields.The parameters of interest are determined from these pseudoexperiments under the baseline and the systematically varied hypotheses.In most of the cases, the average difference between the two results is taken as an estimation of the systematic uncertainty.Exceptions to this are the systematic uncertainties associated with the use of external inputs and the statistical uncertainty of the efficiency correction, where the standard deviation of the difference of the two results is used instead.
The main sources of systematic uncertainty on the Wilson coefficients C 9 and C 10 arise from the use of the external inputs in the determination of the signal branching fraction of Eq. 14.This primarily concerns the uncertainty on the normalisation branching fraction B(B 0 → J/ψK + π − ) = (1.15 ± 0.01 ± 0.05) × 10 −3 [68] and the fraction of B 0 → J/ψK + π − decays that fall in the m Kπ window of the analysis f B 0 →J/ψKπ ±100MeV = 0.644 ± 0.010.The systematic uncertainties associated to the use of these external inputs are provided separately in view of possible future improvements on these quantities.Contributions from the uncertainty on the branching fraction of the J/ψ → µ + µ − decay, B(J/ψ → µ + µ − ) = (5.96± 0.03) × 10 −2 [50], the uncertainty on the efficiency ratio R ε , and the uncertainty on the observed yield of the control channel N J/ψKπ are also considered and reported together under "others" in Table 4.The uncertainty on R ε is due to the finite size of the simulation samples and assumptions made in the simulation model.The model dependence of the simulation is studied by varying the signal model in multiple ways: hadronic parameters are extensively varied within and beyond the SM prediction, large variations of the Wilson coefficients are artificially inserted, and an S-wave component compatible with what is observed in the fit to data is added.The different sources are found to contribute to the relative uncertainty on R ε at the level of 1-2%, depending on the q 2 region.Finally, the measurement of the yield in the control channel is found to be systematically dominated, with the prime sources of uncertainty associated to the choice of the signal mass model and assumptions about the residual background contribution from Λ 0 b → pK − J/ψ decays.The different sources contribute to the relative uncertainty on N J/ψKπ at a level below 1%.
The main source of systematic uncertainty for C ′ 9 comes from ignoring the nonlocal hadronic contribution in the S-wave component.In absence of any theoretical study on non-local hadronic effects on Kπ-scalar amplitudes, pseudoexperiments are generated assuming a non-local hadronic component which is identical to the one of the longitudinal P-wave amplitude.Other sources of systematic uncertainties associated to  the modelling of the S-wave amplitudes are related to the choice of the S-wave form factors and k 2 parameterisation.The former is assessed by generating pseudoexperiments with the alternative model from Ref. [64], while the latter is assessed by replacing the LASS lineshape with an isobar model built from the sum of the K * 0 (700) and K * 0 (1430) resonances.
For the combinatorial background modelling, three sources of systematic uncertainty are considered.The first is associated with the choice of second-order polynomials to model the background angular and q 2 distributions.Due to the small number of background candidates, the BDT requirement is relaxed and the background candidates selected in the upper-mass sideband are fitted with a fourth-order polynomial in each of the angles and q 2 .This model is used as an alternative model for the generation of pseudoexperiments.The second is associated with the modelling of the k 2 distribution, where the value of the fraction of the resonant component introduced in Sec.5.5 is varied within its uncertainty.The third is associated to the assumption of complete factorisation of the background distributions.This is studied in the upper-mass sideband.A mild non-factorisation between angles ϕ and θ ℓ is observed and an alternative background model that does not assume factorisation in these two variables is used for the generation of pseudoexperiments.
In addition, systematic uncertainties are assessed for the different sources of peaking background that are neglected in the analysis.The distribution of residual peakingbackground events is studied in data, after removing PID information from the BDT and inverting the background vetoes.Events are then drawn from the selected background samples and injected into the pseudoexperiment data.GKvD'18 + LQCD Fit result q 2 > 0 only Fit result q 2 < 0 constr.5 10 Form factor results as a function of q 2 obtained from the amplitude fit in the two fit configurations, compared to the predictions from Refs.[16,60,61] that are used as external constraint in the fit.
Finally, two sources of systematic uncertainties are associated to the determination of the acceptance function.The first is related to the finite size of the simulated samples used to derive the acceptance coefficients and is studied by sampling the obtained coefficients within their covariance matrix.The second is associated to the choice of the order of the Legendre polynomials used and is investigated by considering a higher-order acceptance parameterisation built from polynomials of order six, seven, eight and four for cos θ ℓ , cos θ K , ϕ and q 2 , respectively.

Local form factors
Figure 4 shows the form factor results obtained from the amplitude fit in the two tested configurations.A tendency to slightly adjust the ratio F ⊥,∥ /F 0 towards lower values with respect to the theoretical predictions used as external input to the fit is observed, as shown in Fig. 4 (right).Both the fit configurations with and without the q 2 < 0 constraints manifest this behaviour coherently.

Non-local hadronic contributions
Figure 5 shows the real and imaginary parts of the non-local hadronic contributions obtained for the two fit configurations normalised to the size of the local form factors.The two results are compatible, however some discrepancy is visible in their imaginary parts, especially in Im(H ∥ ).The theoretical predictions at q 2 < 0 impose an extremely strong constraint on the shape of these contributions, which are in fact forced to be approximately constant (and have an imaginary part very close to zero) at negative q 2 .The size of Im(H λ (q 2 )) is then found to rise in the physical region.At finite truncation order, the presence of the constraint at q 2 < 0 limits the flexibility of Im(H λ (q 2 )) in the physical region and overconstrains their contribution towards smaller values.The q 2 > 0 only q 2 < 0 constr.) normalised to the size of the local form factors F λ (q 2 ) obtained for the two fit configurations.The black dots correspond to the predictions at q 2 < 0 from [16,24].The result obtained without the constraints is also extrapolated to the negative q 2 region for comparison.The shaded gray regions correspond to the vetoed J/ψ and ψ(2S) regions.
behaviour of these functions in the transition between the unphysical and physical regions of q 2 is further investigated in Appendix B and the imaginary part of H λ (q 2 ) is found to rise more rapidly than the theoretical predictions.It is interesting to note that, while phase differences between the amplitudes are predicted to be tiny at low q 2 , significant differences are measured between the amplitudes for B 0 → J/ψK * 0 decays.One of the advantages of the parameterisation proposed in Refs.[16,24] is the introduction of a dispersive bound to provide control over the systematic truncation errors on the z-expansion.This states that, under a particular choice of polynomial functions, the sum of the coefficients squared over all b → sℓℓ processes must be less than unity.However, the dispersive bound is found to be irrelevant for this analysis since it is very far from being fulfilled, as the sum of the coefficients squared, after the appropriate basis transformation, is found to be O(10 −3 ), for the fit result without the constraints at negative q 2 .Finally, a good compatibility between the input values and corresponding fit results is observed on all the B 0 → ψ n K * 0 observables.Moreover, in addition to the differences of phases provided by B 0 → ψ n K * 0 external measurements, this analysis introduces another phase difference that can be determined from the model, namely the difference between the phase of A ψn 0 and the local amplitudes.The phase difference of the J/ψ longitudinal amplitude (at the J/ψ mass pole) with respect to the rare mode is found to be −1.55 +0.22   −0.18   for the fit result with the q 2 < 0 constraints and −1.61 +0.22  −0.20 for the one without these constraints, 3 showing a good agreement between the two fit configurations.This result is also compatible with one of the two solutions obtained in the measurement of the phase 3 The fit result with q 2 > 0 only information shows a second solution at about ϕ J/ψ 0 → ϕ J/ψ 0 + π, which is however excluded at more than 3 σ.
difference between B + → K + µ + µ − and B + → J/ψK + decays [21], which are ruled by the same rare-electroweak and tree-level underlying transitions, respectively, but with a different spectator quark.The phase difference of A ψ(2S) 0 with respect to the rare mode shows an almost complete degeneracy and cannot be determined precisely from the fit.

Wilson coefficients
Table 5 reports the values of the Wilson coefficients for the two fit configurations, together with their confidence intervals (C.I.) and compatibility with the Standard Model.For each of the four Wilson coefficients, confidence intervals are built from the one-dimensional profile likelihood scans shown in Fig. 6.The 68% (95%) C.I. range is identified with the interval where the negative log-likelihood difference, ∆NLL, is smaller than 0.5 (2).The difference between the best fit values and the corresponding SM predictions obtained are ∆C 9 = −0.93+0. 53 −0.57(−0.68 for the fit configuration without (with) constraints at negative q 2 , where the SM prediction at the b-quark energy scale is taken to be C SM 9 = 4.27, C SM 10 = −4.17 and C ′ SM 9,10 = 0 [29,30].The coefficient that shows the largest difference with respect to the SM is C 9 , whose compatibility with the SM is found to be at the level of 1.9 and 1.8 standard deviations, for fit models using only q 2 > 0 information and with the q 2 < 0 constraints, respectively.Two-dimensional profile likelihood contours for the Wilson coefficients are shown in Fig. 7, where the 68% (95%) C.I. range is identified with the region where the ∆NLL is smaller than 1.15 (3.09).A shift of approximately 0.2 is observed in the central values of all the Wilson coefficients between the two fit configurations, with the fit result with the q 2 < 0 constraints being closer to the SM.While from a theoretical perspective one could expect that non-local hadronic contributions would only affect C 9 , the experimental determination of the Wilson coefficients is affected by the strong correlations of the system: a modification of the non-local hadronic contributions is found to influence the result on the form factors (as shown in Fig. 4), which in turn have an impact on the Wilson coefficients.This behaviour has been studied with pseudoexperiments, where the same generated dataset is fitted with and without the constraints at negative q 2 replicating the procedure adopted on data, and the variation measured in data is found to be compatible with what is observed in the pseudoexperiments.
Finally, the global compatibility with respect to the SM is evaluated by inspecting the likelihood difference in the four-dimensional space given by the four considered Wilson coefficients.Taking into account the systematic uncertainties, the observed difference in twice the log-likelihood between the best fit and SM point is found to be 2.99 (3.25).Considering the four degrees of freedom of the system, this corresponds to 1.3 (1.4) standard deviations with respect to the SM, for the fit without (with) the negative q 2 constraints.
Table 5: Best fit value, confidence intervals and deviation from the SM predictions [29,30] for the four Wilson coefficients and the two fit configurations.For each Wilson coefficient, the likelihood has been profiled over the other coefficients.The SM predictions at the b-quark energy scale [29,30] are also reported for reference.

Comparison to binned observables
Conventional angular observables accessed by binned angular analyses [7][8][9] can be determined from the fit results by dividing the angular coefficients, I i (q 2 , k 2 ), by the differential decay rate, d 2 Γ P /dq 2 dk 2 , both integrated over k 2 .The determination of these angular observables offers an important perspective for the validation and interpretation of the results.Figures 8 and 9 show the q 2 -dependent angular observables derived from the amplitude fit results.The contributions from non-local effects to the so-called CPaveraged S i [32] and corresponding optimised P i [13] series of observables, ∆S(P ) bscc i , is also illustrated in the plots.In general, the post-fit determination of the angular observables agrees very well with the dedicated measurement of Ref. [9] and the overall impact of non-local hadronic contributions on the angular observables is found to be compatible between the two tested fit configurations.The only exception is observed in the S 7 (and the related P ′ 6 ) observable, which is related to the imaginary part of the product of the longitudinal and parallel amplitudes, where the fit result that includes the theory points at q 2 < 0 does not have enough freedom to fully accommodate the shape observed in the physical region.This is a reflection of the different behaviour of the imaginary part of H λ (q 2 ) between the two fit configurations observed in Sec.8.2.In addition, a closer look at the P ′ 5 observable indicates that non-local hadronic contributions are responsible for a positive shift in P ′ 5 of the order of 0.1 ± 0.1 in the region between 4 and 8 GeV 2 /c 4 .This is found to be true for both the fit configurations with and without the q 2 < 0 constraints, with the latter characterised by a naturally larger uncertainty.
Similarly, the signal branching fraction can be derived from the amplitude fit parameters q 2 > 0 only q 2 < 0 constr.through Eq. 15. Figure 10 shows the unbinned determination of the B 0 → K * 0 µ + µ − differential branching fraction obtained from the amplitude fit compared to the values reported in Ref. [1].Despite the larger sample analysed in this measurement, a visible systematic shift toward lower values is observed in the analysis presented here.This difference is due to the improved understanding of the B 0 → J/ψK + π − normalisation channel from Ref. [67], which provides the external inputs required for the branching fraction determination as in Eq. 14.It is observed that when the same external inputs of Ref. [1] are used, the agreement improves dramatically.Finally, Fig. 10 also shows the obtained fraction of S-wave in the signal region, which is found to be in line with previous estimations.

Discussion
This analysis provides an innovative investigation of B 0 → K * 0 µ + µ − decays that enables the determination of the signal amplitude parameters directly from data.Non-local hadronic contributions are parameterised through a z-expansion, where the truncation of the series unavoidably introduces some model dependence.Two complementary fit configurations are developed with the intention to provide the highest level of information  that can be extracted from the data.The first one provides a direct determination of non-local hadronic contributions which relies entirely on experimental data; the second one includes the best theory knowledge on those matrix elements available to date [16].Both models provide a result which is in line with global fits to b → sµ + µ − mediated decays available in the literature [16,27,[73][74][75][76], favouring a shift in C 9 of about −0.9 or −0.7, depending on the theory assumptions.The angular analysis of Ref. [9], with which this analysis shares the analysed dataset, evaluated the C 9 tension with respect to the SM at 3.3 σ with the Flavio software package [77].This evaluation was, however, based on binned inputs and dimensional estimates of the hadronic uncertainties.In addition, the sole variation of C 9 was considered, i.e. all other Wilson coefficients were fixed to their SM values.On the other hand, Ref. [16], which employs the same parameterisation used in this analysis, considers a simultaneous variation of the Wilson coefficients C 9 and C 10 which results in a similar compatibility with the SM as the one obtained in this analysis.The main difference between this analysis and the work in Ref. [16] is the use of the data on an event-by-event basis.This unbinned approach allows one to explore experimental data in its full complexity and therefore achieves the highest possible precision on the physics parameters of interest.In addition, one of the advantages of the analysis presented in this paper is a more flexible investigation of the non-local contributions, which can be studied in relation to what is observed in the data.Some tension is observed in the imaginary part of those contributions where the theory points seem to overconstrain the behaviour of Im(H λ ) to the physical region.The effect is mostly visible in the S 7 angular observable, where non-zero values can only originate from strong phase differences in the non-local hadronic amplitudes.The fit model with only q 2 > 0 information follows naturally the measurement provided by the binned angular analysis, while the fit with theory cannot depart significantly from zero.
Concerning the P ′ 5 observable, non-local hadronic contributions are found to reduce the tension with respect to the Standard Model.However, despite the high degree of freedom allowed by the inclusion of non-local effects in the model, the fit still prefers to insert a shift in C 9 rather than in the hadronic parameters, assuming the theoretical inputs on local form factors are reliable.The use of predictions at q 2 < 0 does not significantly change the level of agreement between the Wilson coefficients and their SM predictions.
Finally, local form factor predictions are currently the limiting factor for the understanding of the tension observed in the branching fraction measurements of many b → sµµ decay channels.Any further indication on the contribution of these elements to the decay rate is therefore extremely valuable.This analysis favours a small shift in the ratio of the form factors with respect to the theory predictions.

Conclusions
This paper presents the first unbinned amplitude analysis of B 0 → K * 0 µ + µ − decays.This analysis provides a rich and complementary set of information compared to more standard analyses of the same decay that considered only binned observables [1,9].By exploiting the q 2 dependence of the decay amplitudes, an experimental determination of the impact of both short-and long-distance contributions to the decay is obtained for the first time.This provides the most accurate parameterisation of the impact of long-distance effects on the B 0 → K * 0 µ + µ − decay process to date.To favour possible reinterpretation of the analysis, the signal amplitude parameters are made publicly available as supplemental material [78].The fit results indicate a shift in the C 9 coefficient between −0.7 and −0.9, depending on the theory assumptions employed.This shift is consistent with results of global analyses of binned observables.The contribution from non-factorisable corrections alone cannot fully explain the SM discrepancies seen in these decays but does reduce the tension observed in C 9 to the level of 1.8 and 1.9 standard deviations.The global compatibility with the Standard Model is also evaluated by simultaneously considering the Wilson coefficients C 9 , C 10 , C ′ 9 and C ′ 10 and is found to be between 1.3 and 1.4 standard deviations, where the two sets of numbers reflect the choice of theoretical assumptions employed in the analysis.The LHCb result from Ref. [9] is overlaid for comparison, together with the SM predictions from DHMV [14,15] and (for P ′ 5 ) GRvDV [16].q 2 > 0 only q 2 < 0 constr.LHCb Run 1 Figure 10: Branching fraction and S-wave fraction obtained a posteriori from the fit results of the two fit configurations.The published Run 1 measurement from LHCb [1] has been overlaid for comparison.The SM branching fraction prediction from GRvDV [16] is also reported.

Appendices A Formalism
A.1 Angular coefficients In the SM operator basis, the P-wave angular coefficients entering in Eq. 1 are where the negative sign in front of I 4,6s,7,9 is due to the adoption of the LHCb angular convention that differs from the theory convention of Refs.[17,23] by the relations given in Ref. [79] and β l = 1 − 4m 2 l /q 2 , with m l the mass of the lepton.Similarly, the introduction of the S-wave contribution gives origin to the following additional set of angular coefficients and B ′ L is the Blatt-Weisskopf barrier function [82], which depends on the momentum of the decay products and on the size of the decaying particle, known as the meson radius parameter, which is fixed to d = 1.6 GeV −1 c [67].The systematic uncertainty associated with the choice of this value is negligible.In Eq. 20, the first term is a pure kinematic phase space factor while (B ′ L p L ) is the orbital angular momentum barrier factor that accounts for spin-dependent effects in the conservation of the angular momentum for the K * 0 → K + π − decay.The angular momentum between the K * 0 meson and the dimuon system is considered to be zero.
The S-wave component of the signal is modelled using the LASS parametrisation [37], given by where m 0 and Γ 0 are the pole mass and width of the K * 0 (1430) 0 resonance and where a and r are empirical parameters whose values are fixed to a = 1.95 GeV −1 c and r = 1.78 GeV −1 c as from the LASS experiment [37].In order to assess the systematic effect of this choice, these parameters are also fixed to values used in Ref. [1], a = 3.83 GeV −1 c and r = 2.86 GeV −1 c and the resulting systematic uncertainty is found to be negligible.The second term of Eq. 22 is equivalent to a Breit-Wigner function for the K * 0 (1430) resonance.Phase space and orbital angular momentum barrier factors associated to B 0 → K + π − µ + µ − decays employed in Refs.[1,8] have been omitted in Eqs.20 and 22 since these terms are already included in the form factors and amplitude normalisation of Eqs. 2, 3 and 7.
Finally, the P-and S-wave k 2 -dependent lineshapes to be included in the decay amplitudes are defined as where fBW ( fLASS ) is the the Breit-Wigner (LASS) function of Eq. 20 (22) normalised to unity and the coefficients g S and δ S determine respectively the relative magnitude and phase between the two P-and S-wave contributions.
B Non-local contributions at q 2 = −1 GeV 2 /c 4 The baseline fit with q 2 < 0 constraints includes only three of the four q 2 points where SM predictions on the non-local contributions are available [24].Hence, the remaining point at q 2 = −1 GeV 2 /c 4 can be used to test the compatibility of the fit result with the theoretical prediction.Figure 11 shows the obtained H λ (q 2 ) function at low q 2 .While the real part of H λ /F λ is consistent with the theory prediction at q 2 = −1 GeV 2 /c 4 , the imaginary part tends to rise more rapidly than the theoretical predictions.Note that all theory points are strongly correlated, hence the compatibility with the point at q 2 = −1 GeV 2 /c 4 is poor.In fact, to include the theory point at q 2 = −1 GeV 2 /c 4 as part of the constraints to the amplitude fit, it is necessary to further increase the truncation of the polynomial expansion by one additional order.However, this additional degree of freedom uniquely modifies the behaviour of the functions H λ around that point, without providing any significant change to the quality of the fit to data.Hence, no additional information is provided by the inclusion of the point at q 2 = −1 GeV 2 /c 4 and all conclusions remain unchanged.

Supplemental Material -Bootstrapped sets of fit parameters
In order to favour possible reinterpretation of the analysis, bootstrapped samples of signal parameters are uploaded as Supplemental Material [78].These are obtained by sampling the data and repeating the fit for each bootstrapped dataset, where the main sources of systematic uncertainty are included in the fit.The published set of parameters includes the Wilson coefficients as well as P-wave local and non-local hadronic parameters, for the two studied fit configurations.These sets of parameters can be employed to derive confidence intervals on any physics quantity of interest, i.e. checking the compatibility of the result of this analysis with alternative signal parametrisations.

(
form factors) and H λ encode the local and non-local hadronic matrix elements, respectively, and m b and M B correspond to the b quark and B 0 meson masses.The coefficient N is a normalisation factor given by

Figure 3 :
Figure 3: Distribution of the fit variables for the combined Run 1 and 2016 datasets.The distributions of the three angles, q 2 , and k 2 are given for candidates within 50 MeV/c 2 of the known B 0 mass.The total fit projections together with the individual signal and background components are overlaid.

Figure 5 :
Figure5: Real and imaginary part of the non-local contributions H λ (q 2 ) normalised to the size of the local form factors F λ (q 2 ) obtained for the two fit configurations.The black dots correspond to the predictions at q 2 < 0 from[16,24].The result obtained without the constraints is also extrapolated to the negative q 2 region for comparison.The shaded gray regions correspond to the vetoed J/ψ and ψ(2S) regions.

Figure 6 :
Figure 6: One-dimensional profile likelihood scan of the Wilson coefficients.Shaded regions correspond to the 68% and 95% confidence intervals.

Figure 7 :
Figure 7: Two-dimensional profile likelihood scan of the Wilson coefficients.Shaded areas correspond to the 68% C.I. and 95% C.I. contour regions.Dotted contours in the top left plot assume right-handed Wilson coefficients fixed to their SM values, i.e.C ′ 9 = C ′ 10 = 0.

Figure 8 :
Figure8: Angular observables (S-basis) obtained a posteriori from the fit results of the two fit configurations; the subfigures isolate the contribution from non-local effects to the given angular observables.The LHCb result from Ref.[9]  is overlaid for comparison, together with the SM prediction from GRvDV[16].

Figure 9 :
Figure9: Angular observables (P -basis) obtained a posteriori from the fit results of the two fit configurations; the subfigures isolate the contribution from non-local effects to the given angular observables.The LHCb result from Ref.[9]  is overlaid for comparison, together with the SM predictions from DHMV[14, 15]  and (for P ′ 5 ) GRvDV[16].

Figure 11 :
Figure11: Results for the ratio H λ (q 2 )/F λ (q 2 ) obtained from the amplitude fit model with theory.The central values are shown with solid lines, while shaded bands indicate the 68% confidence intervals.The theoretical prediction at q 2 < 0 overlaid for comparison.

Table 1 :
[55]ral values and uncertainties used to constrain the CKM parameters in the fit.All values are taken from the CKMfitter group as of Summer 2019[55].

Table 4 :
Summary of the systematic uncertainties on the Wilson coefficients.The individual sources are described in the text.The subtotal and total values are obtained by adding individual sources in quadrature.