Prospects for disentangling long- and short-distance effects in the decays B → K∗μ+μ−

Marcin Chrzaszcz, 2, 3 Andrea Mauri, Nicola Serra, Rafael Silva Coutinho, and Danny van Dyk 4 European Organization for Nuclear Research (CERN), Geneva, Switzerland Universität Zürich, Physik Institut, Winterthurer Strasse 190, 8057 Zürich, Switzerland Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of Sciences, Kraków, Poland Physik Department, Technische Universität München, James-Franck-Straße 1, 85748 Garching, Germany


I. INTRODUCTION
The sensitivity of the decay B → K * µ + µ − to effects of beyond the Standard Model (SM) physics is well known (see e.g. [1] for a review). Consequently, this decay is the standard candle in indirect searches for New Physics (NP) effects. A recent analysis of this decay by the LHCb collaboration has first established [2] the so-called P 5 [3] "anomaly", i.e. a deviation in measurements of the eponymous observable from the SM predictions by ∼ 3σ. Subsequent analyses by both LHCb [4] and Belle [5] further increased the tension between the SM predictions and the data. Non-standard measurements of the Lepton-Flavor-Universality (LFU) ratios in b → s processes -such as of R K and R K * [6] by LHCb [7,8] -suggest that a NP explanation of the P 5 anomaly could simultaneously be LFU violating.
In their attempts to understand the anomalies, many phenomenological studies of this decay strive to modelindependently constrain the effects of New Physics. This is usually achieved within the framework of an effective field theory. Within the latter, a subset of the Wilson coefficients C i for the basis of dimension-six operators O i are fitted from data. For the purpose of this letter, we use the effective weak Lagrangian [9], , C 3,..., 6 , α s C 8 , with the current-current operators and radiative/semileptonic operators Our conventions are the same as the ones of ref. [10], which we follow closely.
A detriment to extracting the C i , for i = 7, 9, 10 from data is our lack of knowledge of the hadronic matrix elements of the operators O i . For local interactions, these matrix elements are expressed in terms of hadronic form factors. The latter can be accessed either from first principles through Lattice QCD simulations, or from quark-hadron-duality arguments through QCD Light-Cone Sum Rules (LCSRs). The matrix elements of non-local operators involving insertions of O 1,2 , however, turn out to be more difficult to determine, and have been the focus of much attention over the last two decades. Presently, the largest systematic uncertainty in determinations of the Wilson Coefficient C 9 arises from our lack of understanding of the non-local hadronic matrix elements where J e.m. denotes the electromagnetic current, and represents the four quark operator involving two charm quark fields. In the above k denotes the four-momentum of the final-state hadron, and q describes the momentum transfer to the virtual photon. It is convenient to decompose H µ into scalar-valued Lorentz-invariant quantities H λ (q 2 ) as in [10]. Here λ = 0, ⊥, denotes the polarisation state of the dilepton system.
The objects H λ can be accessed in the limit of large kaon energy in the B rest frame, or equivalently for q 2 a few GeV 2 m 2 b [11,12]. This QCD-improved Factorisation (QCDF) approach has inspired a larger number of phenomenological analyses. However, all of these studies treat the off-shell contributions from the charm pair as perturbative. This treatment is known to receive substantial corrections from soft-gluon emissions [13] off the charm loop, even for the region 1 GeV 2 ≤ q 2 ≤ 6 GeV 2 that is usually used for phenomenological studies. It is therefore not surprising that the QCDF calculations do not agree with the measured observables, for example in the purely hadronic decays B → K * {J/ψ, ψ(2S)}. An alternative to a theoretical determination of H are data-driven analyses 1 [14][15][16][17]. Some of these analyses show promise in fitting a resonance model to the q 2 spectrum. Others determine the non-local contributions below the J/ψ from data. Both approached therefore give access to model-dependent determinations of the WC C 9 . In addition, information on the model parameters are only available a-posteriori, which precludes genuine SM predictions of the observables.
Both drawbacks, the perturbative treatment of the charm quarks at timelike q 2 , and the model assumptions in the parametrisation of the q 2 spectrum have recently been overcome through a parametrisation that is valid for −7 q 2 ≤ M 2 ψ(2S) [10]. In that study, pseudo observables based on the theoretical predictions are generated at q 2 < 0, for which one expects rapid convergence of the Light-Cone OPE [13]. In addition, the residues of the scalar valued correlators H λ (q 2 ) can be constrained from experimental results on the branching ratios and angular distribution of the hadronic decays B → ψ n K * , where ψ n = J/ψ, ψ(2S). Last but not least, expressing the ratio of the H λ over their corresponding (local) form factors F λ in combination with an expansion in the conformal variable z ensures the correct analytic behaviour. Since we follow the results of reference [10] closely, the parametrisation used for the correlators H λ (q 2 ) does not reproduce the physical light-hadron cut starting at q 2 = 4M 2 π . However, it has been argued that this cut is suppressed [10], and we will discuss -in parts -its numerical impact and the possible model bias that the lack of the cut introduces in section III.
Given this recent progress on the non-local matrix elements we now aim to study the possibility of applying the z expansion to future experimental analyses: First, we want to establish to which extent information concerning the non-local matrix elements can be inferred from experimental data of the semileptonic decay. In order to maximise the sensitivity to the parameters, we will focus on an extended unbinned analysis of the data; see [18] for a related study of unbinned analyses with focus on the Kπ S-wave background. Second, we investigate the convergency of the series expansion at different order of z and what is the residual model-bias introduced by truncations. Finally, we want to establish the smallest amount of theoretical inputs necessary to find evidence for new phenomena in quark-flavor physics in a single b → s process and through a single measurement.

II. PRELIMINARIES
Assuming on-shell K * dominance, the decay B → K * (→ Kπ)µ + µ − involves four kinematic variables: the dimuon mass square q 2 , as well as two helicity angles within the µ + µ − and Kπ decay planes, respectively, and the azimuthal angle between the planes (see [19,20] and subsequent publications). The Probability Density Function (PDF) for this decay gives rise to a larger number of angular observables, which can be used to extract information on the short-distance physics. Here, we do not use these angular observables directly, but rather use the angular information of the signal PDF for B → K * (→ Kπ)µ + µ − decays in its entirety.
We work with two signal PDFs: PDF 1 and PDF 2 , defined as  Here the four-differential rate is a sesquilinear form of the set of transversity amplitudes A λ (q 2 ), with polarisation states λ =⊥, , 0, t [9]. For later discussion, it is instrumental to understand that the amplitudes with λ =⊥, , 0 can (in the SM) be written as [10] A L,R In the above, the functions F (T ) λ (q 2 ) stand for (linear combinations) of local form factors (FF), while the non-local matrix elements H λ (q 2 ) have been introduced in section I.
The PDFs describe the combined q 2 and full-angular distribution of the decay in two kinematic regions: Note that we impose no constraints on the angular support in either of the PDFs. For a definition of d 4 Γ/(dq 2 d 3 Ω) we refer to [21,22] and references therein.
Toy ensembles are either labelled as "SM", in which case we fix all Wilson Coefficients (WCs) to their SM values; or "Benchmark Point" ("BMP"), in which case the WC C 9 is shifted by −1 from its SM value. Each toy ensemble consists of N ≡ N 1 + N 2 toy events, which are drawn from the combined log-PDF ln PDF ≡ ln PDF 1 + ln PDF 2 of the aforementioned two q 2 regions (see eq. (8)). The total number of events N is varied to explore the sensitivity for present and future experiments.
In all cases under study, we maximise the total log-likelihood ln L tot ≡ ln L 1 + ln L 2 + ln L B with respect to the nuisance parameters {α j } and additionally -for New Physics fits -the WCs C 9 and (in some cases) C 10 . Here L 1,2 are unbinned likelihoods of N i toy events x n,i ∼ PDF i , The last term, L B , incorporates two Poissonian terms for the integrated branching ratios of the decay in the kinematics regions Q 1 and Q 2 , respectively.
We then perform a series of frequentist fits to determine the viability of the approach, and to determine uncertainty intervals for the various parameters. In our nominal fits up to z 2 , we float the full set of 39 nuisance parameters {α j }, with Gaussian constraints as described above. Exception to this are marked appropriately in the text. All quoted 68% confidence level intervals are determined from profile likelihoods. For our studies we are using the same setup as in [10], however independently implemented and cross-checked against the EOS software [25].
The LHCb experiment has already observed 969 and 330 B 0 → K * 0 µ + µ − events in the bins 1.1 GeV 2 ≤ q 2 ≤ 8.0 GeV 2 and 11.0 GeV 2 ≤ q 2 ≤ 12.5 GeV 2 , respectively [4]. Extrapolating in q 2 to the larger bin widths as defined in eq. (8), we fix N LHCb-Run I ≡ 1850. In the same fashion we obtain N Belle = 56. In order to study the sensitivities for future data sets, we extrapolate the number of events by scaling the luminosities and bb production cross section σ bb (s), where s denotes the designed centre-of-mass energy of the b-quark pair. For the LHCb experiment we use σ bb ∝ √ s, while for Belle-II we use σ bb ∝ s. The exact numbers of simulated events for each experiment are listed in table I.
Modelling of both the detector resolution or detection efficiency is hardly possible without access to (non-public) information of the current B physics experiments Belle (II) and LHCb. We therefore assume perfect resolution and efficiency in our studies herefrom out -unless otherwise stated. As a consequence, all of our following results should be understood as upper bounds on the possible sensitivity of any future experimental analysis following our suggestions.
Note that we do not study the hadronic uncertainties in the context of free floating parameters for further WCs 2 beside C 9 and C 10 , since the remaining WCs of semileptonic operators can be disentangled from C 9 as various global analyses of b → s processes have shown; see e.g. [26][27][28][29][30]. Moreover, the hadronic non-local effects can always be attributed to q 2 -dependent shifts of the WCs C 9 and C 9 (the chirality-flipped counterpart to C 9 ). Sizeable shifts to C 9 require NP contributions of non-SM chirality in the operators O (c) 1 (2 ) , and are not further discussed here; see also [31] for a related discussion in the presence of (pseudo)scalar four-quark operators. Consequently, we are convinced that it suffices to demonstrate the separability of hadronic and NP contributions to C 9 as we set out to do in this article.

III. INITIAL STUDY
Our main questions concerning the prospects of future analyses pertain to disentangling the hadronic effects of the hadronic correlators from NP effects in the Wilson coefficient C 9 .
These questions are: A Can the parameters {α  Prior theoretical knowledge can be used in the fits through Gaussian constraints on the nuisance parameters α (k) λ describing the non-local correlator. However, these constraints can only be produced for parametrisations with a truncation at z 2 [10]. Here we investigate if, in principle, one could abstain from using Gaussian constraints on these parameters, and instead determine them only from data of the semileptonic decay.
In order to answer this question, we perform our analyses in two ways. For our first analysis, we use the constraints provided in [10], which are based on theory calculations in the negative q 2 region, as well as experimental measurements of the angular distributions of the decays B → K * J/ψ and B → K * ψ(2S). These constraints are included in form of a multivariate Gaussian distribution with correlation information taken into account. Our second analysis does not use either source of information as a constraint on the parameters {α (k) λ }, and floats them instead within the range α A series of 500 toy data sets have been produced with the BMP scenario, and then fitted with and without usage of the constraints. From our toy analyes we conclude the following: 1. The analysis with Gaussian constraints on the parameters {α (k) λ } is able to extract additional information on the hadronic correlators from the fit, i.e., the obtained uncertainties are smaller than the corresponding Gaussian constraints. The uncertainties on the hadronic parameters scale by a factor between 0.5 and 0.8 for the expected statistics of the LHCb Run II.
2. When removing the constraints, we find that the fit still converges, and that we are able to disentangle hadronic effects from NP in C 9 . Our estimator for C 9 is unbiased for a large number of events ∼ 30k. For the expected statistics of the LHCb Run II the uncertainties obtained from the fit on the parameters {α (k) λ } are found to be comparable to the ones from the prior predictions [10].

B. Finite-width effects
The constraints for the parameters {α k } are based [10] on theoretical results (at negative q 2 ) and the residues of the ψ = {J/ψ, ψ(2S)} poles of the correlation functions H λ (q 2 ) in the complex q 2 plane. Information on the residues can be obtained from experimental results for the angular distribution of the decays B → K * (→ Kπ)ψ(→ µ + µ − ) in small mass windows around the respective ψ masses. In reference [10], the narrow charmonia have been assumed to be stable to simplify the discussion. As a consequence, the poles of the proposed parametrisation are located on the real q 2 axis. However, the physical poles are shifted below the real axis by finite-width effects. Moreover, the shift is directly connected to the q 2 shape of the resonances.
In order to study possible bias introduced by neglecting the width of the narrow charmonia, we study an ensemble of 1k toy analyses corresponding to N LHCb Upgrade = 62k events each. We produce the toy events for each fit by using the SM scenario. For each toy analysis we perform two fits to the toy data: one with the nominal PDF, and one for which we modify the PDF such that the poles corresponding to the two narrow charmonium states are shifted below the real q 2 -axis by iM ψ Γ ψ . Note that since this shift is not q 2 dependent, the induced imaginary part does not vanish for q 2 < 4M 2 π . Consequently, our fit PDF does not respect unitarity. However, since it is only used to model possible fit bias at q 2 ≥ 1 GeV 2 4M 2 π , this does not pose a problem here. 3 We carry out both fits, with the nominal and the modified PDF without the usage of any Gaussian constraints on the parameters {α k }. We find that the results of the fits to toys with and without the width effects are indistinguishable even for an ensemble corresponding to the LHCb Upgrade. As such, the bias introduced by neglecting the finite width in the semileptonic regions Q 1 and Q 2 will not play a relevant role for any of the upcoming data sets.

C. Model bias of and sensitivity to higher orders in z
In the analysis [10] the truncation order was chosen as K = 2, in order to ensure that a-priori predictions (i.e., genuine SM predictions without using B → K * µ + µ − data) can be made. In the applications discussed here we are not bound to using the priors presented in [10] as Gaussian constraints; see section III A. Thus, we can explore the sensitivity to coefficients that enter with z 3 or even higher powers in z. Including these in the experimental analysis has the potential to reduce the model-dependence of the results on C 9 due to the z truncation.
We choose to probe the sensitivity to coefficients of order z 3 as follows: We first produce 4k toy analyses with data sets corresponding to the LHCb Run II luminosity. For each study, the pseudo events follow from the BMP scenario. In addition, we introduce the coefficients {α (λ) 3 } for the higher order terms in the correlators H λ (q 2 ). We produce toys for α Consequently, we find no sensitivity to coefficients which are smaller in magnitude than 5 · 10 −3 .
We also investigate if the value of C 9 obtained from a fit to order z 3 is fully compatible with the fit to order z 2 . Our analysis yields C 9 z 3 fit − C 9 z 2 fit = 0.17 (11) and σ(C 9 ) z 2 fit = 0. 19 and σ(C 9 ) z 3 fit = 0.69 .
Our findings can be summarised as follows: • The impact of z 3 terms on the extraction of C 9 amounts to a shift by (formally) less than one standard deviation of the fit to order z 2 , when considering data set sizes up to the LHCb Run-II size.
• The model-dependence of the fit to order z 3 is large in the absence of any theory constraints on the parameters {α (k) λ }. We conclude that a small model-dependence can only be achieved by using more information than only the experimental data of the semileptonic decay B → K * µ + µ − . However, when using order z 3 or higher, the current set of theoretical and experimental constraints does not allow for a staged approach [10], in which the posterior of a C 9 -agnostic fit is used as a prior for the NP fit. We see three possibilities to overcome this problem: 1. through a staged analysis that uses theoretical predictions beyond what has been discussed in [10] (e.g. lattice QCD calculations of the hadronic correlator on or in between the narrow charmonium resonances); 2. through a simultaneous analysis of the semimuonic and semielectronic decays, in which the parameters for the non-local matrix elements are shared (as shown in [32]); 3. through a combined analysis of the theory predictions at negative q 2 , the measurements of the hadronic decays B → K * {J/ψ, ψ(2S)}, and the measurements of the semileptonic decays. This approach is investigated in the next section.

IV. COMBINED UNBINNED ANALYSIS
In this section we extend our previous analysis by performing a combined unbinned fit of the theory predictions at negative q 2 , the measurements of the hadronic decays B → K * {J/ψ, ψ(2S)}, and the measurements of the semileptonic decays. We investigate the following points: A Can the residual model-bias seen in section III C be reduced by including theoretical predictions and measurements of the hadronic decays B → K * {J/ψ, ψ(2S)}?
B What is the individual impact of the theoretical predictions and the hadronic decays inputs?
C In this combined analysis, what model bias is introduced by the truncation assumption at the production level?
D What are the prospects for a simultaneous fit to the WCs C NP 9 and C NP 10 ? E What is the gain in precision when determining the shape of the q 2 -differential angular observables compared to binned analyses?
A. Combined fit to theory points at q 2 < 0, hadronic and semileptonic decays In order to perform a combined fit to all the available (theoretical and experimental) information, we extend the current framework to include the predictions on the hadronic correlator calculated for points at negative q 2 and the two sets of pseudo-observables (three magnitudes and two relative phases for each of the J/ψ and ψ(2S) resonances), as in [10]. Both contributions are included in the fit as multivariate Gaussians of the relevant pseudo-observables. For the production of the ensembles the corresponding central values are shifted to match their predictions given a certain set of parameters {α   observable obtained from z 2 , z 3 , z 4 and z 5 fits for the BMP scenario with Re C NP 9 = −1.
We explore the model-bias introduced by the truncation of the series by repeating the fit with different truncation orders for the z expansion. Two sets of toy analyses with 500 ensembles each are produced, with data sets corresponding to the LHCb Run II and LHCb Upgrade [50 fb −1 ] luminosity, respectively.
In Table II we report the resulting sensitivity to the C NP 9 obtained for z 2 , z 3 , z 4 and z 5 fits. Our findings can be summarised as follows: • The uncertainty on C NP 9 roughly doubles moving from z 2 fits to z 3 fits, for both statistics under consideration.
• For the dataset corresponding to the expected statistics at the LHCb Upgrade [50 fb −1 ], the uncertainty on C NP 9 slightly improves for orders higher than z 3 .
• For the dataset corresponding to the expected statistics at the LHCb Run II we observe the saturation of the uncertainty at the orders z 3 and z 4 . However, for the fit with z 5 the uncertainty starts to increase, pointing to a statistical limitation.
The observed distributions of the hadronic correlator H λ (q 2 )/F λ (q 2 ) for different orders in z are shown in Fig.1, corresponding to the statistics expected at the LHCb Upgrade [50 fb −1 ]. In all cases, the uncertainty drastically increases for higher orders in the ψ(2S) window. For the regions Q 1 and Q 2 , however, we find that the behaviour observed for C NP 9 is reflected in the real part of the hadronic correlator, i.e. the uncertainty saturates for orders higher than z 3 .

B. Exploring the impact of the inputs from theory and hadronic decays
We further investigate the impact of the additional pseudo-observables introduced in the combined fit, to distinguish the benefits obtained from either of the two inputs. In particular, we investigate the following questions: • Is it possible to perform a purely experimental analysis, i.e., excluding the theoretical points at negative q 2 and relying only on semileptonic and hadronic decays?
• The pseudo-observables obtained for the hadronic decays currently constrain only two relative phase between the three polarisations. What is the impact of a hypothetical theory determination of the absolute phase of the hadronic decays or an increased precision of the relative phases?
To address the first point we repeat the analysis removing the constraints introduced by the theoretical calculation at negative q 2 , and we examine the stability of the fit scanning different orders of the expansion. We consider the BMP scenario corresponding to the expected statistics at LHCb Upgrade [50 fb −1 ]. Fig. 2 shows the result of the fit assuming z 4 truncation of the expansion performed with and without the input from the theory points. We find a strong model-bias, similarly to what is presented in section III C. We conclude that a purely experimental analysis that combines information from the semileptonic decay and the hadronic B → K * {J/ψ, ψ(2S)} decays is not currently possible. The desired disentangling of the hadronic effect from possible NP contributions crucially relies on the theory inputs from the points at negative q 2 .
We investigate the benefits from hypothetical improvements to the constraints based on the hadronic decays. First of all we note that, as shown in Fig. 1, the uncertainty on the hadronic correlator evaluated at the J/ψ is extremely small. This is due to the fact that the region of the J/ψ is highly constrained by the interference of theory information at negative q 2 and the events of the semileptonic decay. In fact we find that, already for datasets corresponding to the expected statistics at LHCb Run II, the impact of the pseudo-observables of the J/ψ is negligible. Furthermore, the fit is able to select the absolute phase of the J/ψ with the same precision as the two relative phases. As a consequence, in the following we focus on the impact of the ψ(2S) pseudo-observables on the combined fit and we test whether it would be beneficial to have a measurement of the absolute phase of the ψ(2S) and/or assuming future improvement in the measurement of the pseudo-observables of the ψ(2S) (currently the two relative phases are weakly constrained [10]). We proceed with a hypothetical constraint on the absolute phase of the ψ(2S) and repeat the analysis with two configurations: first, we assume the relative uncertainty of the absolute phase of the ψ(2S) to be similar to the relative phases' uncertainties. Second, we reduce the uncertainties of all phases of the ψ(2S) to reflect the uncertainties of the J/ψ. In both cases the central value of the absolute phase of the ψ(2S) is set to the prediction obtained from the default set of parameters {α (k) λ } used for the production of the ensembles. We find that, even assuming the best case, the improvements on the determination of the WC C NP 9 are negligible. Fig. 3 shows the comparison of the results obtained for the hadronic correlator for the analysis carried out with all the currently available information and the analysis that assumes the future improvements on the ψ(2S) described above. Both analyses are performed with datasets corresponding to the LHCb Upgrade [50 fb −1 ] expected statistics and assume the z 4 truncation in the series expansion. We find that the benefits produced by the assumed improvements on the ψ(2S) pseudo-observables are limited to the region of the ψ(2S).
C. On the truncation of the series at z K All the studies presented so far assumed a fixed set of initial values for the parameters {α (k) λ } in the production of the ensembles (obtained from [10]), in this section we investigate the effect on the fit results of a different choice for the initial values of the hadronic parameters. We investigate two options. First, we produce ensembles with non-zero coefficients for order of the expansion up to z 3 (i.e. α (λ) i≥4 = 0). Second, we produce ensembles with non-zero coefficients for order of the expansion up to z 4 (i.e. α (λ) i≥5 = 0). The choice of the above-mentioned non-zero coefficients is based on the following criteria: they must be realistic (i.e. compatible the theory predictions at negative q 2 and the pseudo-observables from the hadronic decays) and reduce the tension with the P 5 anomaly (i.e. hadronic effects mimic the behaviour of NP). The resulting set of parameters is shown in Fig. 4 together with the value of the P 5 angular observable obtained in the different cases.
For each of the two generated configurations we perform the analysis as described in section IV A, repeating the fit by varying the truncation of the expansion from z 2 to z 5 . Results are shown in Tables III and IV, respectively. Our conclusion can be summarised as follows: • fitting with the expansion truncated at z 2 (i.e. lower order than what is used for the production of the ensembles) introduces a strong bias in the estimator for C NP 9 ; • when the order of the truncation in the fitting procedure catches up the one used for the production of the ensembles the estimator for C NP 9 is unbiased; • the uncertainty on C NP 9 varying the order of the fit follows the pattern observed in section IV A.
• our lack of knowledge on the real description of the hadronic effects in nature can be investigated by scanning the order of the truncation of the series until the central value of C NP 9 stabilises. Obviously, this procedure is bounded to the limit of the available statistics. As mentioned in section II, C NP 10 does not suffer from pollution from hadronic non-local effects. Nevertheless, it is interesting to extend the explored WCs parameter space and study the effect of floating C NP 9 and C NP 10 simultaneously observable obtained from z 2 , z 3 , z 4 and z 5 fits for the BMP scenario when produced with non-zero z 3 coefficients as described in the text.

LHCb Run2
LHCb observable obtained from z 2 , z 3 , z 4 and z 5 fits for the BMP scenario when produced with non-zero z 4 coefficients as described in the text.
in the fit. We repeat the analysis as in section IV A, producing 500 ensembles assuming the BMP scenario (with NP inserted only in C NP 9 , i.e. C NP 10 = 0) with z 2 and performing the fit with truncation at the order z 2 , z 3 , z 4 and z 5 . The 2D pulls are shown in Fig. 5 while the single projections are reported in Tables V and VI. Figure 6 shows the same result for different datasets corresponding to the expected statistics at LHCb Run II, LHCb Upgrade [50 fb −1 -300 fb −1 ] and Belle II [50 ab −1 ]. For simplicity, only the results obtained from the z 3 analysis are shown.
We note that: • The uncertainty on C NP 9 varying the order of the fit follows the pattern observed in section IV A; • due to the correlation between C NP 9 and C NP 10 the projection of the uncertainty on the single WC C NP 9 is larger compared to the case of section IV A when C NP 10 was fixed in the fit; • besides the non-local hadronic effects, a precise determination of the WCs is limited by the uncertainties on the form factors, in particular, the precision on C NP 10 already saturates with the statistics expected to be collected at LHCb Run II (see Fig. 6). The precision on C NP 10 can be substantially improved by including constraints from the decay B s → µ + µ − .

E. Unbinned determination of angular observables
One of the benefits of our proposed approach is that it takes advantage of the full unbinned description of the decay and, additionally, the amplitude fit allows to reproduce confidence intervals for the commonly used angular observables. In the following we investigate the statistical uncertainty expected for the obtained angular observables. We perform the analysis on 500 ensembles generated with the expected statistics at LHCb Run II and we repeat the fit with different truncations at the z 2 , z 3 , z 4 and z 5 order. We find that the uncertainty on all the angular observables is independent on the assumption on the truncation of the series expansion, leading to clean results free from systematic uncertainties on the hadronic parametrisation. It is interesting to compare the statistical uncertainty on the angular observables obtained by the unbinned amplitude fit with respect to the binned approach. We perform a binned fit on the same ensembles generated above, splitting the datasets in 1 GeV 2 q 2 bins, and we fit the obtained angular distributions with the signal PDF d 3 Γ/d 3 Ω as described in [4]. Figure 7 shows the improvement on the statistical uncertainty on the angular observables obtained by the unbinned amplitude fit compared to a q 2 binned angular analysis. and C NP 10 for different non-local hadronic parametrisation models. The contours correspond to 3 σ statistical-only uncertainty bands evaluated with the expected statistics after LHCb Run II (left) and LHCb Upgrade [50 fb −1 ] (right).

V. CONCLUSION
Measurements of angular observables in the decay B → K * µ + µ − have shown discrepancies with respect to Standard Model predictions, mainly in the angular observable known as P 5 . This anomaly has been widely discussed in the literature, in particular since non-local charm contributions are challenging to be predicted from a theory point of view. Here we carried out a sensitivity study of the decay B → K * µ + µ − , taking care to account for hadronic matrix elements of both local and non-local operators in a model-independent fashion. This is done by performing an extended unbinned maximum likelihood fit, which allows to use the full information of the data and extract simultaneously the hadronic parameters and the Wilson Coefficients. Other studies have proposed methods to disentangle NP from non-local hadronic effects in (un)binned likelihood fits. Following the parametrisation of Ref. [10] we studied the properties of the fit for a large number of different scenarios, by using simulated pseudoexperiments. Fitting with the order of the z-expansion used in the production of the ensembles leads to an unbiased value of C NP 9 . Increasing the  order of the z-expansion, the central value of C NP 9 stays unbiased while the uncertainty increases. It is observed that theory constraints strongly mitigate this problem. The increase in the uncertainty on C NP 9 is roughly one order of magnitude smaller than the statistical uncertainty obtained adding one order to the z-expansion in the fit. The fact that the uncertainty on C NP 9 steadily increases with the order of the polynomial of the z-expansion does not allow us to rigorously assign a systematic uncertainty due to the truncation of the z-expansion a-priori. We also found that the unbinned fit allows to extract additional information on the non-local matrix elements from semi-muonic decay events alone. Our study goes beyond previous works, and assesses in a quantitative way the model-dependency due to our ansatz for the non-local hadronic contributions. Our approach allows systematic improvements (through increasing the truncation order in z) and estimation of systematic uncertainties (through varying the truncation order even when the data is described well). In addition, the unbinned fit can be used for the determination of the usual angular observables with precision beyond what can be expected with the standard binned approach. We find that the angular observables obtained with this method do not exhibit any sizeable model bias due to the truncation of the z-expansion. It should be emphasised that the gain in sensitivity cannot be directly read from Fig. 7 since the points of the binned likelihood fit are all uncorrelated. To fully access the comparison of the two approaches, a binned fit to the angular observables with the same model should be performed.
While not strictly necessary, improving our knowledge of the B → K * ψ(2S) amplitudes can give us greater confidence of the obtained fit results. We therefore encourage revisiting their analysis at the present B-physics experiments. The application of the unbinned fit for the decays B → Kµ + µ − and Λ b → Λµ + µ − has not been studied here. It is unclear how the lack of phase information on the J/ψ and ψ(2S) poles will affect the efficiency of the unbinned fit. We therefore encourage dedicated sensitivity studies for these decays as a natural extension of our present work.