Maximising the physics potential of B ± → π ± µ + µ − decays

: We present a method that maximises the experimental sensitivity to new physics contributions in B ± → π ± µ + µ − decays. This method relies on performing an unbinned maximum likelihood fit to both the measured dimuon q 2 distribution of B ± → π ± µ + µ − decays, and theory calculations at spacelike q 2 , where QCD predictions are most reliable. We exploit the known analytic properties of the decay amplitude and employ a dispersion relation to describe the non-local hadronic contributions across spacelike and timelike q 2 regions. The fit stability and the sensitivity to new physics couplings and new sources of CP -violation are studied for current and future data-taking scenarios, with the LHCb experiment as an example. The proposed method offers a precise and reliable way to search for new physics in these decays.

Global analyses of these updated measurements point predominantly to anomalous couplings between a left-handed sb current and a vectorial lepton current [14,15].Such a hint is quantitatively supported by, separately, branching ratios and angular b → s data [16], whose systematic uncertainties are generally very different.A more mundane explanation of the experimental measurements involves underestimating hadronic contributions in the SM [15].Such hadronic effects involve non-local matrix elements of four-quark operators that are hard to compute from first principles.However, recent re-appraisals of these hadronic components suggest they are less likely to be the cause of the observed anomalies in b → sµ + µ − decays [17].Such a conclusion could be validated by suitable observables at high q 2 , which share the very same short-distance sensitivity while not suffering from the same long-distance issues [16].These observables include B s → µ + µ − γ [18] and the inclusive B → Xsµ + µ − [19] among the others.
Traditionally, measurements of b → sµ + µ − transitions involve binning the data in regions of the invariant mass of the dimuon system squared (q 2 ) and performing measurements of decay rates and angular observables within each of these bins.Recent developments in theory and experiment have opened up the possibility of fitting the entirety of the differential decay rate of B → K ( * ) µ + µ − transitions to determine new physics couplings and hadronic contributions from the data [20][21][22].
The additional CKM (Cabibbo-Kobayashi-Maskawa matrix) suppression in the SM of b → dℓ + ℓ − relative to b → sℓ + ℓ − processes makes observables of the former even more sensitive probes of new physics [23].In light of the tensions with SM predictions in b → sℓ + ℓ − processes, maximising the experimental sensitivity in B ± → π ± µ + µ − decays is of paramount importance to ascertain a more complete picture of the flavour structure of these tensions, be they due to new physics or hadronic effects.Recently, branching fraction measurements of 26,27] decays have been combined to constrain new physics contributions in b → dℓ + ℓ − processes [28,29].However, such analyses suffer from limited experimental precision and coarse information regarding the q 2 distribution of B ± → π ± µ + µ − processes.The analysis of Ref. [30] uses a dispersive model for the non-local contributions in b → dℓ + ℓ − transitions to predict lepton flavour universality ratios, for which hadronic uncertainties largely cancel.However, in order to ascertain new physics contributions in lepton-flavour-specific final states, it is imperative to separate long-and short-distance effects.This can only be done through an unbinned fit to the dimuon spectrum of B ± → π ± µ + µ − transitions adopting an effective field theory description of the decay amplitudes.Additionally, as will be demonstrated, employing QCD factorisation and light-cone sum rules (LCSR) predictions at negative q 2 to constrain the size of hadronic contributions is essential to maximise sensitivity to new physics in these decays, the incorporation of this information is the primary innovation of this paper.
This paper is organised as follows: Sec. 2 introduces the theoretical background and provides a description of the model used, Sec. 3 describes how the fits to pseudo-datasets are set up, Sec. 4 details our results and finally Sec. 5 provides a conclusion.
The matrix elements arising from these effective operators can be classified as either local form factors or non-local form factors. Local form factors enter the amplitudes through the hadronic matrix elements of a two-parton current, e.g., from the semileptonic operators or the QED radiative operators i = 7, 7 ′ .Non-local form factors enter the amplitudes through the time-ordered product of the electromagnetic current with effective operators: the four-quark current-current or QCD penguin operators; and radiative operators with i = 8, 8 ′ .In Fig. 1 we provide a schematic overview of the two classes of contributions.
In the case of B → π transitions, there exist only three local form factors, which are labelled f + , f 0 and f T .Other form factors must vanish due to Lorentz invariance and parity conservation within the strong interaction.The three form factors are defined via: ) The form factors are scalar-valued functions of the momentum transfer q 2 , which requires some form of parametrization.Here, we use the nominal parametrization and numerical results from Ref. [33].The parametrization used is based on the original BCL parametrization [34].The numerical results are obtained from a combined fit to lattice QCD [35][36][37] and LCSR [33,38,39] inputs.We display these form factor results in Fig. 2. The Lagrangian density Eq. (2.1) gives further rise to non-local contributions, stemming either from the full set of four-quark operators or the radiative operators with i = 8, 8 ′ .In the case of B → π transitions, there exists only a single Lorentz structure for these non-local contributions: The non-local contributions can be recast into a shift to the Wilson coefficient C 9 via: . (2.6) Due to the CP -violating nature of the weak interaction, we must take care to define such a shift separately for the B + and the B − initial state.
Using the above definitions, the differential decay rate for the where This decay rate is defined separately for (2.8)

Modelling the non-local contributions
In the q 2 < 0 region, it is possible to compute the size of the non-local contributions to B ± → π ± µ + µ − transitions using QCD factorization and LCSR.nonf,spect , where the B ± index is dropped for legibility.The individual components are provided in Sec. 3 of Ref. [41].These components are summed to compute the full non-local contribution as in the following expression, (2.9) To model H (p),B ± q 2 across the full q 2 range, including the physical q 2 > 0 region, we employ once-subtracted dispersion relations1 as in Eq. ( 41) of Ref. [41].Combining with Eq. (2.6) results in the following relation, To ensure the convergence of the dispersive integral for H (p),B ± (q 2 ), we require one subtraction in the dispersion relation.The emerging subtraction terms are matched to the results of the QCD factorisation and LCSR calculations at the subtraction point q 2 0 , as originally proposed in Ref. [41].For this analysis, we choose the subtraction point q 2 0 = −1.5 GeV 2 .Finally, the various Y B ± (q 2 ) terms are the individual components of the non-local contributions that will be introduced in the following paragraphs.
Each resonance (V ) contribution to ∆C B ± 9 (q 2 ) is described with a relativistic Breit-Wigner distribution as follows, Here η V is the resonance magnitude, δ V its phase 2 , and Γ V (q 2 ) the running width, (2.12) The description of the width involves the breakup momentum p both as a function of q 2 and evaluated at q 2 = m 2 V .Our choice of the description of the residues in terms of two magnitude and phases, one each for the B + and B − decay, facilitates the description of CP -violation in the decay.
Open charm continuum We jointly model the combination of the non-resonant continuum of open charm states and the contributions due to further broad vector charmonia following the model suggested in Ref. [22].This model is governed by an overall coupling strength for the modelled two-particle open charm continuum and further includes terms for the S-and P -wave contributions.As for the resonance terms, we choose to describe each coupling in terms of a magnitude η and a phase δ, to facilitate the description of CP -violation in the decay.In contrast to our modelling of the resonances, we choose to use the same coupling strength for both B + and B − decay.The model expression reads: where ρi are hadronic spectral densities defined in Ref. [22] and we use the same subtraction point q 2 0 = −1.5 GeV 2 as before.We fix the magnitudes η D * D , η D * D * and η DD of the modelled contributions to unity and fix the phases δ D * D , δ D * D * and δ DD to zero.In contrast, the "global" parameters η 2P and δ 2P are allowed to vary in fits to pseudo-data.
The joint modelling of the heavy charmonium resonances as one-body intermediate states and the two-particle continuum amplitudes inevitably leads to some double counting and model error.We expect this to be insignificant compared to the statistical uncertainties achievable with the upcoming LHCb datasets.To validate this assumption, we assess the impact of this model choice on the measurement of the Wilson coefficients C 9 and C 10 .We perform fits to pseudo-data generated with the default non-local amplitudes, including the ones above the open charm threshold, and fit back with variations of the non-local amplitude that involve turning off individual open-charm resonant and two-particle amplitudes.The resulting variations on the extracted values of C 9 and C 10 are found to be negligible compared to the statistical precision of any current or future experiment.
Light-quark continuum Finally, we need to consider the non-local contribution from the "light-quark" continuum, i.e., the continuum of ūu, dd and ss states.In a perturbative picture, this contribution arises from weak annihilation and light-quark loop diagrams.This contribution is modelled using the following integral over hadronic spectral densities, where ρ LO (s) and Γ eff are provided in Eq. 38, Eq. 39 and in the text of Ref. [41], respectively.Using a duality threshold s 0 = 1.5 GeV 2 reduces the impact of any potential double counting between the ρ(770) and the ω(782) and the light-quark continuum.The physical quantities that build up this component are known well enough such that Y B ± LQC (q 2 ) is fixed in the fit.

Analysis Setup
We generate pseudo-datasets using the decay rate in Eq. (2.7) and keep the Wilson coefficients set to their SM values.The parameters that describe the local form factors are assigned to the central values obtained in Ref. [33].The parameters in the description of the non-local form factors are instead obtained from a χ 2 fit of our model described in Sec.2.1 to the theoretical pseudo data points at q 2 < 0 computed in Ref. [41].Our results for the latter parameters are compatible with those of Ref. [41].We present our results for the non-local contributions expressed in terms of the quantity ∆C 9 in Fig. 3.
To ascertain a realistic expected precision on the parameters of interest from the fit to the q 2 spectrum of B ± → π ± µ + µ − decays, we need to take into account the expected experimental q 2 resolution R(q 2 reco , q 2 ) and the reconstruction efficiency ε(q 2 ).We use the experimental q 2 resolution used in the LHCb analysis of B ± → K ± µ + µ − decays in Ref. [42].Our choice is motivated by the expectation that this resolution is close, if not identical, to the LHCb resolution for B ± → π ± µ + µ − decays.For the reconstruction efficiency, we take the q 2 shape of the efficiency reported in Ref. [42] and extrapolate it linearly to the larger phase space of B ± → π ± µ + µ − decays.The final signal q 2 model is given by the convolution which is obtained through a fast Fourier transform.The signal yield is obtained using the expression where α is a factor that represents all relative efficiency effects such that the calculated signal yield is compatible with that of the measured yields in different q 2 bins in Ref. [24].The factor N B ± →J/ψK ± is the reconstructed yield of B ± → J/ψ(→ µ + µ − )K ± candidates taken from Ref. [42] that used 3 fb −1 of LHCb Run1 data.The ratio of CKM matrix elements converts N B ± →J/ψK ± to the number of expected B ± → J/ψK ± decays.The factor L scales up the yields from 3 fb −1 of LHCb Run1 data to future projections of LHCb integrated luminosities, including the increase in the B-hadron production cross-sections coming from the centre-of-mass energy changes of the LHC.The factor dΓ dq 2 dq 2 |Y ± J/ψ (q 2 )| 2 dq 2 uses our model of the B ± → π ± µ + µ − decay rate to transform the B ± → J/ψπ ± yield into one across the entire q 2 phase space.Finally, the signal purity is estimated from Fig. 3 of Ref. [24], and the q 2 model of the background is taken from Fig. 3 of Ref. [42] and modelled using a kernel density estimator.Pseudo-datasets are generated with a sample size corresponding to that expected in the current LHCb Run1+2 dataset (9 fb −1 ) and future LHCb Upgrades of 23 fb −1 , 45 fb −1 and 300 fb −1 .
We perform unbinned maximum likelihood fits to these pseudo-datasets where the magnitude parameters η B ± J/ψ and the parameters of the local form factors are fixed in the fit.Fixing these parameters incurs a systematic uncertainty, the size of which we assess in Sec.4.2.In this fit configuration, we measure all the phases, including the phase of C 9 , relative to that of C 7 , which is fixed in the fit.
Examples of pseudo-datasets representing 9 fb −1 and 300 fb −1 are presented in Fig. 4 along with the model employed in the pseudo-dataset generation.The non-local, penguin, and interference components of the model are shown separately.

Constraining the non-local contribution
We relate the model of the non-local contribution ∆C B ± 9 , as in Eq. (2.10), to the sum of the various QCD factorisation and LCSR predictions at negative q 2 as in Eq. (2.9).This relationship is visualised in Fig. 3 where the red points denote the QCD factorisation and LCSR predictions, and the black line is our model of the non-local contributions to We can exploit this relation when fitting our model for the B ± → π ± µ + µ − decay rate, Eq. (2.7), to data in the physical q 2 > 0 region through the introduction of a multivariate Gaussian factor to the likelihood function.This factor relates our dispersive non-local model to the theory reference values computed at different q 2 < 0 values, indicated by the red points in Fig. 3.The dimensionality of this multivariate Gaussian constraint is given by the number of negative q 2 points considered multiplied by four3 .
The uncertainties of the reference values and the correlations between these uncertainties need to be taken into account in the multivariate constraint.As Ref. [41] does not provide these correlations, we make the conservative assumption that all the uncertainties used to compute the theory terms are uncorrelated between q 2 points and between real/imaginary B + /B − contributions.Assuming the uncertainties between the q 2 points are uncorrelated reduces the statistical power of the constraint.
Constraints are placed on the magnitude parameters of the resonances V = ρ(770), ω(782) and ψ(2S) using measured central values and uncertainties for the CP -averaged branching fractions of B(B → V (→ µ + µ − )π) and of A V CP [43].These are essential for reliable fit convergence and are employed in all the fits discussed in this paper.

Choosing a q 2 range
The region of q 2 above the open-charm threshold is particularly problematic due to the presence of multiple broad overlapping resonances that interfere with non-resonant contributions.With the number of signal decays expected in the existing LHCb Run1+2 dataset, it is unfeasible to float all the parameters associated with non-local contributions arising from open-charm states.Their impact, however, is sub-dominant for q 2 ≲ 14 GeV 2 .This leads us to fix these parameters and to restrict the phase space region for our analysis.
We use the results from the B + → K + µ + µ − measurement of Ref. [42] scaled by |V cd /V cs | to fix the residues of the open-charm states.We further limit the phase space to q 2 reco < 14.06 GeV 2 .This cut is motivated by the fact that, taking into account resolution effects, contributions above the ψ(3770) are negligible.
In future datasets, such as those expected by LHCb's planned upgrade, the signal yield will be sufficient to fit the entire q 2 phase space with these non-local parameters floating.
Therefore, the open charm region is included in the fits to 300 fb −1 of pseudo-data, as presented in the bottom panels of Fig. 4.

Size of the dataset Relative size
Fit success (%) w/o q 2 < 0 with q 2 < 0 9 fb −1  1. Stability of the fits to pseudo-data.The last column separates fits that do not use theoretical inputs at negative q 2 from those that do.

Contamination from
The decay B ± → K ± µ + µ − with a K ± → π ± misidentification is a potentially dangerous background to measurements of B ± → π ± µ + µ − as it is less CKM suppressed than the B ± → π ± µ + µ − process.However, the binned measurement of the B ± → π ± µ + µ − decay rate presented in Ref. [24] demonstrated that the B ± → K ± µ + µ − background can be brought under control through the use of particle identification information from the ringimaging Cherenkov (RICH) systems of LHCb.In this study, we assume the signal purity of a window of ±40 MeV around the B ± mass as given in Ref. [24].However, the B ± → π ± µ + µ − analysis of Ref. [24] vetoed the regions associated with resonant dimuon contributions from B ± → K ± ψ(→ µ + µ − ) decays, where ψ is J/ψ or ψ(2S).Therefore, our assumed purity of B ± → π ± µ + µ − decays in the q 2 regions near the large charmonia resonances is not valid.In principle, an experimental analysis that attempts to fit the entire q 2 spectrum of B ± → π ± µ + µ − decays would have to adopt stricter particle identification criteria to reduce the background from B ± → K ± ψ(→ µ + µ − ) decays down to a controllable level at the expense of signal efficiency.An experimental analysis may need to undertake some optimisation of the selection, including a background component for B ± → K ± ψ(→ µ + µ − ) backgrounds in the fitted model and studying the impact on the signal precision.Therefore, dealing with this background is beyond the scope of our study.

Experimental precision and prospects
To estimate the expected sensitivity to new physics and understand the impact of the q 2 < 0 constraint, we fit generated pseudo-datasets with and without the theoretical constraint at q 2 < 0 included in the likelihood.Each fit is initialised from multiple starting positions to avoid localised turning points in the likelihood space.The fit result with the largest likelihood is recorded.

Fit stability
With the signal yields expected from the 9 fb −1 LHCb Run1+2 dataset, we find that the best-fit point of a significant fraction of pseudo-datasets lies in an unphysical region.The decay rate of Eq. (2.7) is not differentiable with respect to C 10 in the point C 10 = 0 due to the |C 10 | 2 dependence in the decay rate.As our likelihood minimisation relies on gradient increase of signal yield expected from future runs of the LHC.We present an illustrative example to highlight this point.We overlay intervals obtained using smaller form factor uncertainties in the lower panels of Fig. 6.Here, we assume improved calculations could produce uncertainties three times smaller.This would be in line with the improvements achieved for B → K ( * ) in Ref. [46].The improvement in the intervals is significant and brings the result much closer to the statistical-only intervals in the top panels of Fig. 6.
Given that the flavour anomalies could be indicating the presence of large lepton flavour universality violating (LFUV) contributions to C τ 9 , the study of C τ 9 through b → d{e + e − , µ + µ − } transitions is an increasingly interesting subject [47][48][49][50][51][52][53][54].As demonstrated in Ref. [22], large non-local contributions from C τ 9 can be imprinted into the q 2 spectrum of B ± → K ± µ + µ − decays.Larger futures datasets of B ± → π ± µ + µ − decays could be used to study the q 2 distribution of B ± → π ± µ + µ − decays by including a C τ 9 contribution for τ τ re-scattering to µµ.Additionally, with larger datasets, it would be possible to lift the model-dependence of the open charm continuum model by floating individual components of the Y B ± 2P,cc (q 2 ) model or by allowing for CP -violation.Increasing the complexity of the non-local model will only increase the relevance of the q 2 < 0 constraint.Finally, in the future, it will be possible to fit the B ± → π ± µ + µ − decay rate for the presence of new physics with scalar and tensor Wilson coefficients.This would require a 2D fit of q 2 and the lepton helicity angle cos(θ ℓ ) using the double-differential decay rate [55,56].Employing the q 2 < 0 information will be essential to maximise sensitivity to new physics in all these studies.

Summary and Conclusions
This paper presents an approach that maximises the sensitivity of new physics searches in B ± → π ± µ + µ − transitions.We employ a dispersive model to perform unbinned maximum likelihood fits to both the measured dimuon q 2 spectrum of B ± → π ± µ + µ − decays and to theoretical constraints on the non-local contributions at q 2 < 0. Our approach ensures that the size and the q 2 dependence of non-local contributions to B ± → π ± µ + µ − transitions in the q 2 < 0 region align with predictions.We perform fits to pseudo-datasets and demonstrate the expected sensitivity to CP -violating and CP -conserving contributions for a variety of upcoming datasets.We observe that including the theoretical constraints markedly increases the fit stability and improves the sensitivity to non-local parameters and, consequently, to the Wilson coefficients.Variations in the modelling of the non-local amplitude above the open-charm threshold were found to have a negligible impact on the extracted values of the Wilson coefficients compared to their statistical precision.We conclude that without increased model dependence, an unbinned analysis of the Run1+2 LHCb dataset would be challenging due to poor fit stability.Instead, we present the expected sensitivity for the future scenarios of 45 fb −1 and 300 fb −1 of LHCb data.We include systematic effects arising from our incomplete knowledge of the B ± → J/ψ(→ µ + µ − )π ± branching fractions and local form-factors.We find that uncertainties due to the local form factor knowledge currently form the dominant systematic uncertainty.This highlights that improving the precision of local form factors will be an essential step to fully exploit the physics potential of future datasets.

. 4 )Figure 1 .
Figure 1.Schematic overview of the two classes of contributions to B ± → π ± µ + µ − decays.The local contributions are presented in the left and central sketches, and one example of the non-local contributions is presented in the right sketch.

Figure 2 .
Figure 2. B → π local form factors obtained in Ref. [33] by a combined fit to lattice QCD and light-cone sum rule estimates.The bands correspond to the 68% interval.

13 ) 2
Contrary to what is done in the description of exclusive b → sµ + µ − decays, the phases in our hadronic model for the non-local contributions are not strong phases; instead, they are superpositions of two strong phases arising from the two terms and the relative weak phase in Eq. (2.1).