Measurement of the phase difference between short- and long-distance amplitudes in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^+} \!\rightarrow {{{K}} ^+} {\mu ^+\mu ^-} $$\end{document}B+→K+μ+μ- decay

A measurement of the phase difference between the short- and long-distance contributions to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^+} \!\rightarrow {{{K}} ^+} {\mu ^+\mu ^-} $$\end{document}B+→K+μ+μ- decay is performed by analysing the dimuon mass distribution. The analysis is based on pp collision data corresponding to an integrated luminosity of 3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{\,fb}^{-1}$$\end{document}fb-1 collected by the LHCb experiment in 2011 and 2012. The long-distance contribution to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^+} \!\rightarrow {{{K}} ^+} {\mu ^+\mu ^-} $$\end{document}B+→K+μ+μ- decay is modelled as a sum of relativistic Breit–Wigner amplitudes representing different vector meson resonances decaying to muon pairs, each with their own magnitude and phase. The measured phases of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{J}/\psi }$$\end{document}J/ψ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\psi {(2S)}$$\end{document}ψ(2S) resonances are such that the interference with the short-distance component in dimuon mass regions far from their pole masses is small. In addition, constraints are placed on the Wilson coefficients, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}_{9}$$\end{document}C9 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}_{10}$$\end{document}C10, and the branching fraction of the short-distance component is measured.


Introduction
The decay B + → K + μ + μ − receives contributions from short-distance b → s + − flavour-changing neutral-current (FCNC) transitions and long-distance contributions from intermediate hadronic resonances. In the Standard Model (SM), FCNC transitions are forbidden at tree level and must occur via a loop-level process. In many extensions of the SM, new particles can contribute to the amplitude of the b → s + − process changing the rate of the decay or the distribution of the final-state particles. Decays like B + → K + μ + μ − are therefore sensitive probes of physics beyond the SM.
In this paper, a first measurement of the phase difference between the contributions to the short-distance and the narrow-resonance amplitudes in the B + → K + μ + μ − decay is presented. 1 For the first time, the branching fraction of the short-distance component is determined without interpolation across the J/ψ and ψ(2S) regions. The measurement is performed through a fit to the full dimuon mass spectrum, m μμ , using a model describing the vector resonances as a sum of relativistic Breit-Wigner amplitudes. This approach is similar to that of Refs. [13,24], with the difference that the magnitudes and phases of the resonant amplitudes are determined using the LHCb data rather than using the external information on the cross-section for e + e − → hadrons from the BES collaboration [25]. The model includes the ρ, ω, φ, J/ψ and ψ(2S) resonances, as well as broad charmonium states (ψ(3770), ψ(4040), ψ(4160) and ψ(4415)) above the open charm threshold. Evidence for the ψ(4160) resonance in the dimuon spectrum of B + → K + μ + μ − decays has been previously reported by LHCb in Ref. [26]. The continuum of broad states with pole masses above the maximum m μμ value allowed in the decay is neglected.
The measurement presented in this paper is performed using a data set corresponding to 3 fb −1 of integrated luminosity collected by the LHCb experiment in pp collisions during 2011 and 2012 at √ s = 7 TeV and 8 TeV . The paper is organised as follows: Section 2 describes the LHCb detector and the procedure used to generate simulated events; the reconstruction and selection of B + → K + μ + μ − decays are described in Sect. 3; Section 4 describes the m μμ distribution of B + → K + μ + μ − decays, including the model for the various resonances appearing in the dimuon mass spectrum; the fit procedure to the dimuon mass spectrum, including the methods to correct for the detection and selection biases, is discussed in Sect. 5. The results and associated systematic uncertainties are discussed in Sects. 6 and 7. Finally, conclusions are presented in Sect. 8.

Detector and simulation
The LHCb detector [27,28] is a single-arm forward spectrometer, covering the pseudorapidity range 2 < η < 5, designed to study the production and decay of particles containing b or c quarks. The detector includes a high-precision tracking system divided into three subsystems: a silicon-strip vertex detector surrounding the pp interaction region, a largearea silicon-strip detector that is 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 situated 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 momentum scale of tracks in the data is calibrated using the B + and J/ψ masses measured in B + → J/ψ K + decays [29]. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/ p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors (RICH). Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger [30], 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.
A large sample of simulated events is used to determine the effect of the detector geometry, trigger, and selection criteria on the dimuon mass distribution of the B + → K + μ + μ − decay. In the simulation, pp collisions are generated using Pythia 8 [31,32] with a specific LHCb configuration [33]. The decay of the B + meson is described by EvtGen [34], which generates final-state radiation using Photos [35]. As described in Ref. [36], the Geant4 toolkit [37,38] is used to implement the interaction of the generated particles with the detector and its response. Data-driven corrections are applied to the simulation following the procedure of Ref. [23]. These corrections account for the small level of mismodelling of the detector occupancy, the B + momentum and vertex quality, and the particle identification (PID) performance. The momentum of every reconstructed track in the simulation is also smeared by a small amount in order to better match the mass resolution of the data.

Selection of signal candidates
In the trigger for the 7 TeV (8 TeV ) data, at least one of the muons is required to have p T > 1.48 GeV/c ( p T > 1.76 GeV/c) and one of the final-state particles is required to have both p T > 1.4 GeV/c ( p T > 1.6 GeV/c) and an IP > 100 µm with respect to all PVs in the event; if this final-state particle is identified as a muon, p T > 1.0 GeV/c is required instead. Finally, the tracks of two or more of the final-state particles are required to form a vertex that is significantly displaced from all PVs.
In the offline selection, signal candidates are built from a pair of oppositely tracks that are identified as muons. The muon pair is then combined with a charged track that is identified as a kaon by the RICH detectors. The signal candidates are required to pass a set of loose preselection requirements that are identical to those described in Ref. [26]. These requirements exploit the decay topology of B + → K + μ + μ − transitions and restrict the data sample to candidates with good-quality vertex and track fits. Candidates are required to have a reconstructed K + μ + μ − mass, m K μμ , in the range 5100 < m K μμ < 6500 MeV/c 2 .
Combinatorial background, where particles from different decays are mistakenly combined, is further suppressed with the use of a Boosted Decision Tree (BDT) [39,40] using kinematic and geometric information. The BDT is identical to that described in Ref. [26] and uses the same working point. The efficiency of the BDT for signal is uniform with respect to m K μμ .
Specific background processes can mimic the signal if their final states are misidentified or partially reconstructed. The requirements described in Ref. [26] reduce the overall contribution of the background from such decay processes to a level of less than 1% of the expected signal yield in the full mass region. The largest remaining specific background contribution comes from B + → π + μ + μ − decays (including B + → J/ψπ + and B + → ψ(2S)π + ), where the pion is mistakenly identified as a kaon.
The K + μ + μ − mass of the selected candidates is shown in Fig. 1. The signal is modelled by the sum of two Gaussian functions and a Gaussian function with power-law tails on  Fig. 1 Reconstructed K + μ + μ − mass of the selected B + → K + μ + μ − candidates. The fit to the data is described in the text both sides of the peak; these all share a common peak position. A Gaussian function is used to describe a small contribution from B + c decays around the known B + c mass [41]. Combinatorial background is described by an exponential function with a negative gradient. At low m K μμ , the background is dominated by partially reconstructed b-hadron decays, e.g. from B {+,0} → K * {+,0} μ + μ − decays in which the pion from the K * {+,0} is not reconstructed. This background component is modelled using the upper tail of a Gaussian function. The shape of the background from B + → π + μ + μ − decays is taken from a sample of simulated events. Integrating the signal component in a ±40 MeV/c 2 window about the known B + mass [41] yields 980 000 B + → K + μ + μ − decays.
When computing m μμ , a kinematic fit is performed to the selected candidates. In the fit, the m K μμ mass is constrained to the known B + mass and the candidate is required to originate from one of the PVs in the event. For simulated B + → J/ψ K + decays, this improves the resolution in m μμ by about a factor of two.

Differential decay rate
Following the notation of Ref. [42], the CP-averaged differential decay rate of B + → K + μ + μ − decays as a function of the dimuon mass squared, q 2 ≡ m 2 μμ , is given by where |k| is the kaon momentum in the B + meson rest frame.
Here m K and m B are the masses of the K + and B + mesons while m s and m b refer to the s and b quark masses as defined in Ref. [42], m μ is the muon mass and β 2 = 1−4m 2 μ /q 2 . The constants G F , α, and V tq are the Fermi constant, the QED fine structure constant, and CKM matrix elements, respectively. The parameters f 0,+,T denote the scalar, vector and tensor B → K form factors. The C i are the Wilson coefficients in an effective field theory description of the decay. The coefficient C 9 corresponds to the coupling strength of the vector current operator, C 10 to the axial-vector current operator and C 7 to the electromagnetic dipole operator. The operator definitions and the numerical values of the Wilson coefficients in the SM can be found in Ref. [43]. Right-handed Wilson coefficients, conventionally denoted C i , are suppressed in the SM and are ignored in this analysis. The Wilson coefficients C 9 and C 10 are assumed to be real. This implicitly assumes that there is no weak phase associated with the short-distance contribution. In general, CP-violating effects are expected to be small across the m μμ range with the exception of the region around the ρ and ω resonances, which enter with different strong and weak phases [44]. The small size of the CP asymmetry between B − and B + decays is confirmed in Ref. [45]. In the present analysis, there is no sensitivity to CP-violating effects at low masses and therefore the phases of the resonances are taken to be the same for B + and B − decays throughout.
Vector resonances, which produce dimuon pairs via a virtual photon, mimic a contribution to C 9 . These long-distance hadronic contributions to the B + → K + μ + μ − decay are taken into account by introducing an effective Wilson coefficient in place of C 9 in Eq. 1, where the term Y (q 2 ) describes the sum of resonant and continuum hadronic states appearing in the dimuon mass spectrum. In this analysis Y (q 2 ) is replaced by the sum of vector meson resonances j such that where η j is the magnitude of the resonance amplitude and δ j its phase relative to C 9 . These phase differences are one of the main results of this paper. The q 2 dependence of the magnitude and phase of the resonance is parameterised by A res j (q 2 ). The resonances included are the ω, ρ 0 , φ, J/ψ, ψ(2S), ψ(3770), ψ(4040), ψ(4160) and ψ(4415). Contributions from other broad resonances and hadronic continuum states are ignored, as are contributions from weak annihilation [46][47][48]. No systematic uncertainties are attributed to these assumptions, which are part of the model that defines the analysis framework of this paper.
The function A res j (q 2 ) is taken to have the form of a relativistic Breit-Wigner function for the ω, ρ 0 , φ, J/ψ, ψ(2S) and ψ(4040), ψ(4160) and ψ(4415) resonances, where m 0 j is the pole mass of the jth resonance and 0 j its natural width. The running width j (q 2 ) is given by where p is the momentum of the muons in the rest frame of the dimuon system evaluated at q, and p 0 j is the momentum evaluated at the mass of the resonance. To account for the open charm threshold, the lineshape of the ψ(3770) resonance is described by a Flatté function [49] with a width defined as where m D is the mass of the D 0 meson and q 2 0 is the q 2 value at the pole mass of the ψ(3770). The coefficients 1 = 0.3 MeV/c 2 and 2 = 27 MeV/c 2 are taken from Ref. [41] and correspond to the sum of the partial widths of the ψ(3770) to states below and above the open charm threshold. For q 2 < 4m 2 D , the phase-space factor accompanying 2 in Eq. 6 becomes complex.
The form factors are parameterised according to Ref.
[50] as with, for this analysis, N = 3. Here m B * s (m B * s0 ) is the mass of the lowest-lying excited B s meson with J P = 1 − (0 + ). The coefficients b + i are allowed to vary in the fit to the data subject to constraints from Ref. [42], whereas the coefficients b 0 i and b T i are fixed to their central values. The function z is defined by the mapping and

Fit to the m µµ distribution
In order to determine the magnitudes and phases of the different resonant contributions, a maximum likelihood fit in 538 bins is performed to the distribution of the reconstructed dimuon mass, m rec μμ , of candidates with m K μμ in a ±40 MeV/c 2 window about the known B + mass. The m rec i.e. by Eq. 1, multiplied by the detector efficiency, ε, as a function of the true dimuon mass, m μμ , and convolved with the experimental mass resolution R discussed in Sect. 5.2.

Signal model
The magnitudes and phases of the resonances are allowed to vary in the fit, as are the Wilson coefficients C 9 and C 10 . As the contribution of C 7 to the total decay rate is small, it is fixed to its SM value of C SM 7 = −0.304 ± 0.006 [43]. The form factor f + (q 2 ) is constrained in the fit according to its value and uncertainty from Ref. [42]. The form factors f 0 (q 2 ) and f T (q 2 ) have a limited impact on the normalisation and shape of Eq. 1, and are fixed to their values from Ref. [42]. The masses and widths of the broad resonances above the open charm threshold are constrained according to their values in Ref. [51]. The masses and widths of the ρ, ω and φ mesons and the widths of the J/ψ and ψ(2S) mesons are fixed to their known values [41]. The large magnitude of the J/ψ and ψ(2S) amplitudes makes the fit very sensitive to the position of the pole mass of these resonances. Due to some residual uncertainty on the momentum scale in the data, the pole masses of the J/ψ and ψ(2S) mesons are allowed to vary in the fit.
The short-distance component is normalised to the branching fraction of B + → J/ψ K + measured by the B-factory experiments [41]. After correcting for isospin asymmetries in the production of the B + mesons at the ϒ(4S), the branching fraction is B(B + → J/ψ K + ) = (9.95 ± 0.32) × 10 −4 [52]. This is further multiplied by B(J/ψ → μ + μ − ) = (5.96 ± 0.03)×10 −2 [41] to account for the decay of the J/ψ meson. Table 1 Resolution parameters of the different convolution regions in units of MeV/c 2 . The α l and α u parameters are shared between the J/ψ and ψ(2S) regions. The parameters without uncertainties are fixed from fits to the simulated events The branching fraction of the decay B + → K + μ + μ − via an intermediate resonance j is computed from the fit as where τ B is the lifetime of the B + meson. The branching fractions of B + → ρ K + , B + → ωK + , B + → φ K + and B + → ψ(3770)K + are also constrained assuming factorisation between the B decay and the subsequent decay of the intermediate resonance to a muon pair. These branching fractions are taken from Ref. [41].

Mass resolution
The convolution of the resolution function with the signal model is implemented using a fast Fourier transform technique [53,54]. The fit to the data is performed in three separate regions of dimuon mass: 300 ≤ m rec μμ ≤ 1800 MeV/c 2 , 1800 < m rec μμ ≤ 3400 MeV/c 2 and 3400 < m rec μμ ≤ 4700 MeV/c 2 .
To increase the speed of the fit, the resolution is treated as constant within these regions using the resolution at the φ, J/ψ and ψ(2S) pole masses. The impact of this assumption on the measured phases of the J/ψ and ψ(2S) resonances has been tested using pseudoexperiments and found to be negligible. This is to be expected as the spectra in all other regions vary slowly in comparison to the resolution function. The resolution is modelled using the sum of a Gaussian function, G, and a Gaussian function with power-law tails on the lower and upper side of the peak, C, The component with power-law tails is defined as and is normalised to unity. The parameters describing the resolution model for the J/ψ and ψ(2S) regions ( f , σ C , σ G , n l , n u , α l , α u ) are allowed to vary in the fit to the data. The parameters α l , α u and f are shared between the J/ψ and ψ(2S) regions. The resolution parameters for the φ region can not be determined from the data in this way and are instead fixed to their values in the simulation. The resulting values of the resolution parameters are summarised in Table 1. As a cross-check, a second fit to the m rec μμ distribution is performed using the full m μμ dependence of the resolution model in Eq. 12 and a numerical implementation of the convolution. In this fit to the data, the parameters of the resolution model are taken from simulated B + → K + μ + μ − events and fixed up to an overall scaling of the width of the resolution function. The two fits to m rec μμ yield compatible results.

Efficiency correction
The measured dimuon mass distribution is biased by the trigger, selection and detector geometry. The dominant sources of bias are the geometrical acceptance of the detector, the impact parameter requirements on the muons and the kaon and the p T dependence of the trigger. Figure 2 shows the efficiency to trigger, reconstruct and select candidates as a function of m μμ in a sample of simulated B + → K + μ + μ − candidates. The rise in efficiency with increasing dimuon mass originates from the requirement that one of the muons has p T > 1.
The values of the parameters ε i are fixed from simulated events and are given in Table 2.

Background model
The reconstructed dimuon mass distribution of the combinatorial background candidates is taken from the m K μμ upper mass sideband, 5620 < m K μμ < 5700 MeV/c 2 . When evaluating m rec μμ , m K μμ is constrained to the centre of the sideband rather than to the known B + mass. Combinatorial background comprising a genuine J/ψ or ψ(2S) meson is described by the sum of two Gaussian functions. After applying the mass constraint, the means of the Gaussians do not correspond exactly to the known J/ψ and ψ(2S) masses. Combinatorial background comprising a dimuon pair that does not originate from a J/ψ or ψ(2S) meson is described by an ARGUS function [55]. The lineshape of the background from B + → π + μ + μ − decays, where the pion is mistakenly identified as a kaon, is taken from simulated events.

Results
The dimuon mass distributions and the projections of the fit to the data are shown in Fig. 3. Four solutions are obtained with almost equal likelihood values, which correspond to ambiguities in the signs of the J/ψ and ψ(2S) phases. The values of the phases and branching fractions of the vector meson resonances are listed in Table 3. The posterior values for the f + form factor are reported in Table 4. A χ 2 test between the data and the model, with the binning scheme used in Fig. 3, results in a χ 2 of 110 with 78 degrees of freedom. The largest disagreements between the data and the model are localised in the m μμ region close to the J/ψ pole mass and around 1.8 GeV/c 2 . The latter is discussed in Sect. 7.
The branching fraction of the short-distance component of the B + → K + μ + μ − decay can be calculated by integrating Eq. 1 after setting the amplitudes of the resonances to zero. This gives where the statistical uncertainty includes the uncertainty on the form-factor predictions. The systematic uncertainty on the branching fraction is discussed in Sect. 7. This measurement is compatible with the branching fraction reported in Ref. [22]. The two results are based on the same data and therefore should not be used together in global fits. The branching fraction reported in Ref. [22] is based on a binned measurement in q 2 regions away from the narrow resonances (φ, J/ψ and ψ(2S)) and then extrapolated to the full q 2 range. The contribution from the broad resonances was thus included in that result.  Table 3 Branching fractions and phases for each resonance in the fit for the four solutions of the J/ψ and ψ(2S) phases. Both statistical and systematic contributions are included in the uncertainties. There is a common systematic uncertainty of 4.5%, dominated by the uncertainty on the B + → J/ψ K + branching fraction, which provides the normalisation for all measurements  LHCb Fig. 4 Two-dimensional likelihood profile for the Wilson coefficients C 9 and C 10 . The SM point is indicated by the blue marker. The intervals correspond to χ 2 probabilities with two degrees of freedom A two-dimensional likelihood profile of C 9 and C 10 is also obtained as shown in Fig. 4. The intervals correspond to χ 2 probabilities assuming two degrees of freedom. Only the quadrant with C 9 and C 10 values around the SM prediction is shown. The other quadrants can be obtained by mirroring in the axes. The branching fraction of the short-distance component provides a good constraint on the sum of |C 9 | 2 and |C 10 | 2 (see Eq. 1). This gives rise to the annular shape in the likelihood profile in Fig. 4. In addition, there is a modest ability for the fit to differentiate between C 9 and C 10 through the interference of the C 9 component with the resonances. The visible interference pattern excludes very small values of |C 9 |. Overall, the correlation between C 9 and C 10 is approximately 90%. The best-fit point for the Wilson coefficients (in a given quadrant of the C 9 and C 10 plane) and the corresponding B + → K + μ + μ − branching fraction are the same for the four combinations of the J/ψ and ψ(2S) phases. Including statistical and systematic uncertainties, the fit results deviate from the SM prediction at the level of 3.0 standard deviations. The uncertainty is dominated by the precision of the form factors. The best-fit point prefers a value of |C 10 | that is smaller than |C SM 10 | and a value of |C 9 | that is larger than |C SM 9 |. However, if C 10 is fixed to its SM value, the fit prefers |C 9 | < |C SM 9 |. This is consistent with the results of global fits to b → s + − processes. Given the model assumptions in this paper, the interference with the J/ψ meson is not able to explain the low value of the branching fraction of the B + → K + μ + μ − decay while keeping the values of C 9 and C 10 at their SM predictions.

Systematic uncertainties
Sources of systematic uncertainty are considered separately for the phase and branching fraction measurements. In both cases, the largest systematic uncertainties are accounted for in the statistical uncertainty as they are included as nuisance parameters in the fit. For smaller sources of uncertainty, the fit is repeated with variations of the inputs and the difference is assigned as a systematic uncertainty. A summary of the remaining systematic uncertainties can be found in Table 5.
The parameters governing the behaviour of the tails of the resolution function are particularly correlated with the phases. The systematic uncertainty on the resolution model is included in the statistical uncertainty by allowing the resolution parameter values to vary in the fit. If the tail parameters are fixed to their central values, the statistical uncertainties on the phase measurements decrease by approximately 20%. The choice of parameterisation for the resolution model is validated using a large sample of simulated events and no additional uncertainty is assigned for the choice of model. For the branching fraction measurement, the uncertainty arising from the resolution model is negligible compared to other sources of systematic uncertainty.
Similarly to the resolution model, the systematic uncertainty associated with the knowledge of the f + (q 2 ) form factor is included in the statistical uncertainty. If the formfactor parameters are fixed to their best-fit values, the statistical uncertainties on the phases decrease by 4% (1%) for the J/ψ (ψ(2S)) measurements. For the branching fraction, the uncertainty is 2%, which is of similar size as the statistical uncertainty. At around m μμ = 1.8 GeV/c 2 there is a small discrepancy between the data and the model (see Fig. 3). This is interpreted as a possible contribution from excited ρ, ω or φ resonances. Given the limited knowledge of the masses and widths of the states in this region, these broad states are neglected in the nominal fit. They are, however, visible in e + e − → hadrons vacuum polarisation data [41]. To test the effect of such states on the phases of the J/ψ and ψ(2S) mesons, an additional relativistic Breit-Wigner amplitude is included with a width and mass that are allowed to vary in the fit. The inclusion of this Breit-Wigner amplitude marginally improves the fit quality around m μμ = 1.8 GeV/c 2 and changes the J/ψ (ψ(2S)) phase by 40% (20%) of its statistical uncertainty, which is added as a systematic effect. The magnitude of the amplitude is not statistically significant and its mean and width do not correspond to a known state. The phases of the other resonances in the fit have larger statistical uncertainties and the inclusion of this additional amplitude has a negligible effect on their fit values. Given that the contribution of this amplitude is small compared to the shortdistance component, its effect on the branching fraction is only around 1%.
Other, smaller systematic uncertainties include modelling of the combinatorial background, calculation of the efficiency as a function of q 2 and the uncertainty on the B + → J/ψ K + branching fraction. The latter affects the branching fraction measurement and is obtained from Ref. [52], which results in a 4% uncertainty.

Conclusions
This paper presents the first measurement of the phase difference between the short-and long-distance contributions to the B + → K + μ + μ − decay. The measurement is performed using a binned maximum likelihood fit to the dimuon mass distribution of the decays. The long-distance contributions are modelled as the sum of relativistic Breit-Wigner amplitudes representing different vector meson resonances decaying to muon pairs, each with their own magnitude and phase. The short-distance contribution is expressed in terms of an effective field theory description of the decay with the Wilson coefficients C 9 and C 10 , which are taken to be real. These are left free in the fit and all other components set to their corresponding SM values. The B → K hadronic form factors are constrained in the fit to the predictions from Ref. [42].
The fit results in four approximately degenerate solutions corresponding to ambiguities in the signs of the J/ψ and ψ(2S) phases. The values of the J/ψ phases are compatible with ± π 2 , which means that the interference with the shortdistance component in dimuon mass regions far from their pole masses is small. The negative solution of the J/ψ phase agrees qualitatively with the prediction of Ref. [47], where long-distance contributions are calculated at negative q 2 and extrapolated to the q 2 region below the J/ψ pole-mass using a hadronic dispersion relation. The fit model, which includes the conventional J PC = 1 −− cc resonances, is found to describe the data well, with no significant evidence for the decays B + → ψ(4040)K + or B + → ψ(4415)K + . The values of the ψ(3770) and ψ(4160) phases are compatible with those reported in Ref. [13].
The measurement of the Wilson coefficients prefers a value of |C 10 | < |C SM 10 | and a value of |C 9 | > |C SM 9 |. If the value of C 10 is set to that of C SM 10 , the measurement favours the region |C 9 | < |C SM 9 |. These results are similar to those reported previously in global analyses. The interference between the short-and long-distance contributions in the regions around the ρ, ω and the φ, and in the region q 2 > m 2 ψ(2S) , results in the exclusion of the hypothesis that C 9 = 0 at more than 5 standard deviations. The dominant uncertainty on the measurements of C 9 and C 10 arises from the knowledge of the B → K hadronic form factors. The current data set allows the uncertainties on these hadronic parameters to be reduced. Improved inputs on the form factors from lattice QCD calculations and the larger data set that will be available at the end of the LHC Run 2 are needed to further improve the measurement of the Wilson coefficients.
A similar strategy to the one applied in this paper can be extended to other b → s + − decay processes to understand the influence of hadronic resonances on global fits for C 9 and C 10 . However, the situation is more complicated in decays where the strange hadron is not a pseudoscalar meson as the amplitudes corresponding to different helicity states of the hadron can have different relative phases.
Finally, a measurement of the branching fraction of the short-distance component of B + → K + μ + μ − decays is also reported and is found to be B(B + → K + μ + μ − ) = (4.37 ± 0.15 (stat) ± 0.23 (syst)) × 10 −7 , where the first uncertainty is statistical and second is systematic. In contrast to previous analyses, the measurement is performed across the full q 2 region accounting for the interference with the long-distance contributions and without any veto of resonance-dominated regions of the phase space. The value of the branching fraction is found to be compatible with previous measurements [22], but smaller than the SM prediction [42].