A new Monte Carlo Generator for BSM physics in B → K ∗ ℓ + ℓ − decays with an application to lepton non-universality in angular distributions

: Within the widely used EvtGen framework, we have added a new event generator model for B → K ∗ ℓ + ℓ − with improved standard model (SM) decay amplitudes and possible BSM physics contributions, which are implemented in the operator product expansion in terms of Wilson coefficients. This event generator can then be used to investigate the experimental sensitivity to the most general BSM signal resulting from dimension-six operators. We describe the advantages and potential of the newly developed ‘Sibidanov Physics Generator’ in improving experimental sensitivity of searches for lepton non-universal BSM physics and clarifying signatures. The new generator can properly simulate BSM scenarios, interference between SM and BSM amplitudes, and correlations between different BSM observables as well as acceptance bias. We show that exploiting such correlations substantially improves experimental sensitivity. As a demonstration of the utility of the MC generator, we examine the prospects for improved measurements of lepton non-universality in angular distributions for B → K ∗ ℓℓ decays from the expected 50 ab − 1 data set of the Belle II experiment, using a four-dimensional unbinned maximum likelihood fit. We describe promising experimental signatures and correlations between observables. The use of lepton-universality violating ∆-observables significantly reduce uncertainties in the SM expectations due to QCD and resonance effects, are ideally suited for Belle II with the large data sets expected in the next decade. Our simulation studies also show that Belle II should have excellent sensitivity to BSM physics in the Wilson coefficients C 7 and C ′ 7 , which appears at low q 2 in the di-electron

1 Introduction Flavor-changing neutral current (FCNC) processes in the weak interaction do not occur at tree-level in the standard model (SM), which makes them sensitive probes of physics beyond the standard model (BSM).Intriguing BSM hints have been reported in b → s transitions by multiple flavor physics experiments, LHCb, Belle and BaBar.Here we investigate the prospects for the decay B → K * ℓ + ℓ − with ℓ = e, µ.We have implemented a new Monte Carlo (MC) signal generator model inside the widely used event generator framework EvtGen [1], with improved SM decay amplitudes implemented in the operator product expansion in terms of the Wilson Coefficients C 7 , C 9 , C 10 and their right-handed counterparts C ′ 7 , C ′ 9 , C ′ 10 .We allow for imaginary contributions to δC i .However, in this paper, we do not consider CP violating observables.Amplitudes for additional BSM physics contributions to B → K * ℓ + ℓ − , which we call δC i = C eff i − C SM i , and which correspond to a general BSM signal resulting from dimension-six operators, can be chosen by the user.Until recently, non-zero δC 9 and δ C 10 were preferred by global theoretical fits to experimental data [2][3][4][5][6].While the significance of lepton-flavor universality violating contributions (LFUV) is reduced and large signals as indicated in [2][3][4][5][6] are not expected, LFUV contributions are still possible and might become significant once again with larger data sets.
Our study considers such a case.Non-zero δC 9 and δC 10 for the di-muon final state are preferred by theoretical fits to current data [16? ? ].
We describe the advantages and the potential of the newly developed 'Sibidanov Physics Generator' in improving statistical sensitivity of BSM physics searches and clarifying experimental signatures.We have examined potential experimental correlations in the four variables that characterize B → K * ℓ + ℓ − decays (Fig. 1): the invariant mass of the lepton pair, q 2 , the lepton helicity angle from the virtual Z/γ-decay, cos θ ℓ , the K * helicity angle, cos θ K , and the angle between the K * and di-lepton decay planes, χ.In addition, there are angular asymmetries that are a function of q 2 .Inclusion of BSM amplitudes in the MC generator allows one to properly simulate acceptance bias in BSM scenarios, interference between SM and BSM amplitudes, and correlations between different BSM observables.For example, δC 9 will lead to correlated signatures in the angular asymmetries A F B versus q 2 and S 5 versus q 2 .Another striking signature of C ′ 7 is found in the modulation of the angle χ at low q 2 in B → K * e + e − .Belle II should have excellent statistical sensitivity to such BSM physics in C 7 and C ′ 7 .We find a number of signatures that were previously not fully appreciated.We show that exploiting correlations between angular observables substantially improves experimental statistical sensitivity.A four-dimensional likelihood function is used to fit for the δC i coefficients.We describe the potential for improved measurements from the expected 50 ab −1 data set of the Belle II experiment with fits to the angles (see Fig. 1) and q 2 in B → K * ℓ + ℓ − decays.
In addition to BSM contributions, the observed angular asymmetries could originate from SM QCD or cc resonance effects, which appear in the experimental final state K * ℓ + ℓ − .Given the difficulties in reliably computing all of the possible hadronic effects, the use of ∆observables appears to be the only possible approach to unambiguously distinguish between these SM effects and lepton-universality violating BSM physics in b → sµ + µ − .Using the BSM Monte Carlo generator, we will demonstrate the feasibility of this experimental approach.We have used the overall detection efficiency from Belle analyses [8] to estimate statistical uncertainties.However, it should be noted that the MC simulation results described here do not yet include backgrounds, final state radiation, detector resolution and efficiency functions.
The three most important ∆-observables are: (1.1) The first two of these observables can be derived from angular asymmetries, while the third is obtained from the four-dimensional, unbinned maximum likelihood fits described above.An example of ∆-like variables was introduced in Ref. [7].Belle reported a first attempt to measure ∆-type observables in B → K * ℓ + ℓ − using a 0.7 ab −1 data sample [8].The ∆-observables appear ideally suited for Belle II (which has comparable sensitivities for di-electrons and di-muons) with the large data sets expected in the next decade [9].
Here we introduce a new Monte Carlo generator to enable further searches for BSM contributions in B → K * ℓ + ℓ − .As an example of the utility of this MC tool, we apply it to generate a number of BSM scenarios and then carry out four-dimensional unbinned maximum likelihood fits to future B → K * e + e − and B → K * µ + µ − Belle II datasets and estimate the statistical sensitivity to lepton flavor violating contributions to the coefficient C 9 .A proposal to perform unbinned fits to angular distributions of dimuon and dielectron B → K * ℓ + ℓ − modes in LHCb datasets was given in [10].However, in contrast to [10], we have developed a MC simulation generator to be used by experimenters to directly explore the effects of BSM Wilson coefficients when large B → K * ℓ − ℓ + data samples are available.This tool allows for straightforward determinations of acceptance and efficiency corrections in BSM scenarios.In particular, using our BSM MC generator, the variations of selection efficiencies and acceptance effects with respect to the underlying physics model can be directly accessed.

Current experimental anomalies in b → s neutral current (FCNC) processes
Tensions with the SM in the lepton-flavor-universality (LFU) violating ratios R K and R * K were indicated by initial results from the LHCb collaboration [11,12], in the low di-lepton mass-squared region.However, recent updated results [13] no longer indicate a significant deviation for the R K and R * K ratios.There remain significant deviations from SM expectations observed in angular asymmetries of B → K * µ + µ − [14] and B s → ϕµ + µ − [15] decays, which seem to be experimentally robust.
Several theoretical groups fit all b → sℓ + ℓ − observables (R K type results and angular asymmetries) for BSM Wilson coefficients, δC i , in b → sµ + µ − and suggest a coherent picture of BSM physics in C 9 and C 10 .BSM physics in δC 9 and δC 10 is found with a range of significances, though this varies between various global analyses, based on different assumptions of hadronic corrections [2][3][4][5]16  In order to explain these deviations, the global fits initially assumed that the BSM contribution is restricted to the di-muon channel and does not affect the di-electron channel.This is plausible if the BSM physics is lepton-flavor dependent and will be the basis for the scenarios that are simulated in this paper.Other scenarios with LFU BSM couplings will be simulated in future work.

Definition of observables and angular asymmetries
For the B → K * ℓ + ℓ − decay with K * → Kπ the differential decay rate can be described in terms of the Kπ invariant mass m Kπ , the di-lepton mass squared q 2 , and three angles θ ℓ , θ K , and χ.The angle θ ℓ is defined as the angle between the direction of the positively charged lepton and the direction of B meson in the di-lepton rest frame.The angle θ K is defined as the angle between the direction of the kaon and the direction of the B meson in the K * -meson rest frame.The angle χ is the angle between the plane formed by the di-lepton pair and the plane formed by the K * decay products in the B-meson rest frame.A graphical representation of the angle definitions is shown in Fig. 1.
In the SM angular asymmetries such as A FB arise due to the interference between different decay amplitudes.In the case of BSM physics there will be additional interference terms that are linear in the BSM contribution, which appear in several observables: the well known forward-backward asymmetry A FB (q 2 ) is defined as where Γ and Γ denote the decay rate of B0 → K0 * ℓ + ℓ − and the CP -conjugate channel B 0 → K 0 * ℓ + ℓ − , respectively (see Eqs. A.5 and A.25 for full angular distributions).Other important angular asymmetries involve the angle χ between the K * and di-lepton decay planes [17]: and Note that the angular observable P ′ 5 , widely used in the literature, is related to S 5 via , where F L is the longitudinal polarization fraction of the K * meson.

SM and BSM Lorentz structures
The matrix element for the decay B → K * ℓ + ℓ − , where ℓ is e, µ or the τ lepton, is the following: where the SM Wilson coefficients are known at NNLL accuracy [18]: and where the function Y (q 2 ) is described in Section A. Note that we include the lepton masses (m µ and m e ) in all results.We assume BSM states are heavy, far above the µ = m b = 4.8 GeV/c 2 scale.We therefore treat BSM contributions to the Wilson coefficients as q 2 -independent offsets.The chirality-flipped operator C ′ 7 is highly suppressed while the C ′ 9 C ′ 10 operators are absent in the SM but can arise in the presence of BSM physics.Note that the SM Wilson coefficient C eff 9 varies with q 2 ; as shown in Fig. 2 separately for real and imaginary parts.The B → K * matrix elements in Eq. 2.4 can be expressed in terms of seven form factors that depend on the momentum transfer q 2 between the B and the K * (q µ = p µ − k µ ): with T 1 (0) = T 2 (0).Here ϵ µ is the polarization vector of the K * .We used the following convention for the Levi-Civita tensor ϵ 0123 = +1.We adopt the formalism developed in Ref. [19] where the extrapolation of form factors from the calculated LCSR input points (q 2 ≲ 0 GeV) to larger q 2 values is performed by a simple pole form with a z-expansion where z(t) The values for the fit parameters α i 0,1,2 including the correlation among them, and the masses of resonances m R,i associated with the quantum numbers of the respective form factor F i are obtained in Ref. [19].

Form factor uncertainties
In the old EvtGen model BTOSLLBALL, used by Belle and Belle II up to now, the form factors were taken from Ref. [20].In our new EvtGen model BTOSLLNP (see Section 2.5), we have now implemented the most recent hadronic form factors from Ref. [19], also known as the ABSZ form factor parameterization.ABSZ updates the previous computations by including current hadronic inputs and several theoretical improvements such as inclusion of higher twist distribution amplitudes of mesons, as well as a combined fit to Light-Cone Sum Rule form factors for the low q 2 region and Lattice QCD estimates valid in the high q 2 range.The form factor fit results are expressed as the coefficients of the expansion shown in Eq. 2.9.
There are seven form factors, which are functions of q 2 : V (q 2 ), A 0 (q 2 ), A 1 (q 2 ), A 2 (q 2 ), T 1 (q 2 ), T 2 (q 2 ), and T 3 (q 2 ).The form factors A 12 and T 23 are defined in Ref. [19], where the form factors A 2 and T 3 are extracted using the following expression: and Here λ(q 2 ) is the Källén-function defined in Eq.A.24 and m Kπ is the invariant mass of the kaon and pion to take into account the finite width of K * .If the K * width is omitted, a singularity appears in the physical region.
The resulting form factors we now use in the updated EvtGen generator are shown in Figs. 3 and 4. Note that the theoretical uncertainties in the form factors shown in the figure limit extraction of new physics.Below, we discuss strategies for overcoming this limitation.11 in Ref. [21].These theoretical uncertainties in the form factors can limit extraction of new physics.

BSM Modifications to the EvtGen Monte Carlo Generator
We implement the B → K * ℓ + ℓ − generator with BSM physics contributions in the EvtGen Monte-Carlo simulation framework.We use the BTOSLLBALL model as a reference to create the BTOSLLNP model, which incorporates the desired matrix element with BSM contributions included.Here we use the generator, BTOSLLNP, only in the standalone mode but it has been also successfully integrated into the Belle II Analysis Software Framework (BASF2 [22]).In the BTOSLLNP model, the user can specify the following BSM input parameters: δC 7 , C ′ 7 , δC 9 , C ′ 9 , δC 10 , C ′ 10 , C S −C ′ S , and C P −C ′ P .In the last two cases only such combinations are relevant.Each of these model parameters is complex-valued.The parameters can be specified in any order and in any combination in the user input file.The default value for each parameter is zero and if no parameters are specified in the file the generator gives the SM result.
To generate BSM physics the user inputs several arguments in the user input file.The first argument specifies the coordinate system-Cartesian(0) or polar(1)-to be used for the complex-valued input arguments.The user can then enter each desired BSM parameter as three consecutive numbers.These triplets can be entered in any order.
Below we present several examples of user input files to illustrate the usage of the BSM

Signal generator performance improvements
In addition to implementing BSM amplitudes and the ABSZ form factors in EvtGen, we also significantly improved the performance of EvtGen.The details of these improvements are briefly described below: • The phase space generation procedure has been reviewed and significantly optimized.
• Tensor contraction operations have been reviewed and optimized.
• The importance sampling for three body decays with a pole has been reviewed and an incorrect treatment of the pole has been identified and fixed.
• The maximal amplitude search has been reviewed and several mistakes have been fixed.
• A special importance sampling procedure has been developed to treat narrow resonances for optimal performance.
Combining together all the EvtGen modifications more than two orders of magnitude performance improvement has been attained as shown in Fig. 5.

3
Signatures of BSM Physics in B → K * ℓ + ℓ − and improved sensitivity via multi-dimensional likelihood fits We show the expected statistical sensitivity of the unbinned likelihood fit under ideal conditions using only the primary generator output without taking into account detector effects or hadronic uncertainties as well as resonance effects.Hadronic uncertainties are considered in Section 4. A likelihood function has been implemented according to Eq. A.5, where the matrix element defined in Eq. 2.4 is already squared, and which explicitly depends only on the four variables q 2 , cos θ ℓ , cos θ K , χ, and on the Wilson coefficients as fit parameters.The likelihood function has been cross-checked against EvtGen-generated distribution and asymmetries.The quark masses m b = 4.8 GeV/c 2 , m c = 1.3 GeV/c 2 were used.Good agreement between the EvtGen predictions and integrals of the likelihood function can be seen in Fig. 6 1 .
The expected experimental data sample depends on integrated luminosity and on the selection criteria required to sufficiently suppress backgrounds.In Belle II, we currently expect about 25% selection efficiency for the B 0 → K * 0 (K + π − )ℓ + ℓ − mode and 40 to 50 ab −1 total integrated luminosity.Here we study only the statistical sensitivity of the fit so all internal parameters in the EvtGen decay generator and the likelihood function are the same.In addition, for simplicity we do not include resonance effects in the generator and the likelihood function and use the entire q 2 range without veto windows.The total decay rate is not fixed in the test.Depending on assumptions we expect from 5000 to 10000 events in the di-muon mode and about 25% more events at very low q 2 for the di-electron 1 Since this is a simplified sensitivity study we do not take into account resonances and thus results in resonance region are not physical.mode.In the sensitivity tests below, the number of events lies within this range.
In each fit we extract the deviation of a single Wilson coefficient from the SM value: The SM values are taken from NNLO calculations [18].Since the Wilson coefficients are complex numbers, both the real and imaginary parts are extracted in these fits, and shown in the figures that follow.
3.1 δC 9 correlated signatures in A FB and S 5 (and dB(B ) )  Figure 6 shows how the observables S 5 and A FB are modified when a BSM contribution (δC 9 ) to C 9 is present.We assume that δC 9 = −0.87,which corresponds to a representative value from global theory fits to b → s experimental data in the µµ channel [2].There are clear differences between the SM and BSM distributions of S 5 and A FB , and these signatures are correlated.However, we do not find any significant effect of δC 9 on the differential distribution dB(B → K * ℓ + ℓ − ) / dq 2 .
The observables S 5 and A FB with the Wilson coefficients varied in a more wide range are shown in Fig. 7, 8, and 9.These plots suggest that BSM physics contributions to the different Wilson coefficients affect different q 2 regions in the observables and the correlations between them have to be selected in the proper q 2 regions to maximize statistical sensitivity to BSM.
In the absence of contributions from right-handed and scalar currents, (i.e.assuming C ′ 7 , C ′ 9 , C ′ 10 and C S are zero) the correlations observed in the results from the Monte Carlo generator between S 5 and A FB in Figs. 7, 8, and 9, can be understood from the following analytical expression relating the two observables: where Z 1 and Z 2 are purely functions of form factors and masses and are given by It is interesting to note that Eq. (3.1) is independent of Im(C 7 ), Im(C 9 ) and C 10 .Eq. (3.1) implies that the sensitivity to both Re(C 7 ), and Re(C 9 ) should disappear at large q 2 .Since all asymmetries must vanish at q 2 = 0, no sensitivity is observed in Figs. 7, 8, and 9, at q 2 ≈ 0. The region q 2 = (1 − 8) GeV 2 is therefore highly sensitive to C 7 and C 9 .It may be noted that only ratios of form factors contribute to S 5 /A FB .

Unbinned four-dimensional maximum likelihood fitting for BSM
Using a high statistics sample of B → K * e + e − and B → K * µ + µ − MC simulated events, we can compare the distributions in the four kinematic variables for the SM and a BSM scenario with δC 9 = −0.87 in three selected regions of q 2 (q 2 < 1 GeV 2 , 1 < q 2 < 6 GeV 2 and q 2 > 15 GeV 2 ).The resulting distributions are shown in Figs. 10, 11 and 12.These three figures show features of BSM physics that were not previously appreciated.For example, in each q 2 region, there are differences in shape between the SM and BSM scenario in cos θ K * and cos θ ℓ .In all three regions, there are changes in the forward-backward asymmetry of the leptonic system, which are visible in the cos θ ℓ distributions.At low q 2 in B → K * e + e − , the SM pole contribution arising from O 7 dominates and BSM effects are difficult to distinguish.However, in B → K * µ + µ − the BSM effects are more prominent in the angular distributions.At higher q 2 , there are differences in the shape of the q 2 distribution, which are correlated to changes in the angular distributions.All such effects can be accommodated with multi-dimensional maximum likelihood fits.
Rather than extracting angular asymmetries in bins of q 2 , and then providing these as inputs to a theory-based fitting package [2][3][4][5] that determines δC i , we fit directly for δC i using 4D unbinned likelihood fits to q 2 , cos θ ℓ , cos θ K , and χ with δC i as a free parameter.For example, Fig. 13 shows the distribution of δC 9 resulting from about 2000 pseudoexperiments performing such multidimensional fits repeatedly, assuming for now that δC 9 = 0 to check the statistical sensitivity.According to the fit results the real and imaginary parts of δC 9 can be constrained with a standard deviation of 0.15 and 0.35, corresponding to 3 and 7 % of SM |C 9 |, respectively.This is a substantial improvement over the traditional analysis method.The pull distributions demonstrate that the obtained fit uncertainties meet expectations.The mitigation of form factor uncertainties is discussed below.
To further test the fit procedure we generate samples with the BSM physics contribution δC 9 = −0.87 and repeat the multidimensional fits.The fit extracts the real part of the BSM contribution correctly with 8 and 10 standard deviation significance for the di-muon and di-electron modes, respectively.The imaginary part is consistent with zero and its uncertainty roughly matches to the SM fit results, albeit a small number of samples show relatively large deviation from zero.The fitted Wilson coefficient distributions obtained in the pseudo-experiments are shown in Appendix B.  -17 - Distribution of δC 9 values obtained from unbinned likelihood fits to simulated B → K * µ + µ − SM decays in the full q 2 range, without resonances included.The real and imaginary parts of δC 9 extracted by the fits are shown in the top row.The bottom row shows the corresponding pull distributions, which demonstrate that the fitter correctly estimates fit uncertainties.About 2000 pseudo-experiments are performed, where in each 40 ab −1 luminosity and 25% efficiency are assumed.

δC 10 sensitivity
The statistical sensitivity to δC 10 is estimated from our multidimensional fit.According to the fit results the real and imaginary parts of δC 10 can be constrained to be 0.21 and 0.28, corresponding to 5 and 7 % of |C 10 | ≈ 4.1, respectively.The corresponding distributions are shown in Appendix B. 3.5 C 7 , C ′ 7 signatures at low q 2 and in the χ angle distribution Sensitivity in the di-electron mode at low q 2 , q 2 < 2 GeV 2 /c 4 , is critical as the new physics terms interfere with the photon pole.Even moderate non-zero C ′ 7 (right-handed BSM physics) produces a prominent modulation in the χ distribution, while left-handed BSM physics (δC 7 ) signatures appear at low q 2 and in the cos θ K distribution.These BSM physics signatures are shown in Fig. 14.The effect of C ′ 7 and δC 7 on angular distributions in the di-electron decay mode for q 2 < 2 GeV 2 /c 4 .There is a striking modulation in the χ angle distribution due to the right-handed BSM physics contribution as well as clear left-handed BSM signatures at low q 2 and in cos θ K .
4 Disentangling QCD and resonance effects from BSM Physics in B → K * ℓ + ℓ − The semileptonic decays B → K * ℓ + ℓ − are expected to have smaller theoretical uncertainties than purely hadronic B decays.However, these decay modes still require QCD corrections, which are not easily calculable.The experimental final state K * ℓ + ℓ − includes contributions from cc resonances such as B → J/ψK * , B → ψ(2S)K * as well as other broader cc resonances with higher mass, where the cc resonances decay to ℓ + ℓ − .It is well known that these b → ccs modes have non-factorizable contributions [23][24][25].In addition, complications arise due to non-local contributions from electromagnetic corrections to purely hadronic operators [26], which cannot be calculated from first principles.Only limited success has been achieved in computing the effects of charm-loop contributions [27].Given the difficulties in reliably computing all the QCD effects, our choice of ∆ observables [9] offers the only existing approach to obtain sensitivity to BSM physics.The global fits assume that BSM effects depend on lepton flavor and are present in the muon mode but are negligible in the electron mode.This is required to explain the deviation in R K and R K * from the SM.The SM hadronic uncertainties due to the strong interaction are the same in the both modes, and thus should cancel to first order in the difference so that only BSM contributions remain.
We test the assumption that hadronic uncertainties affect the di-muon and di-electron modes equally, and cancel in the difference.We fit to pseudo-experiments with MC statistics roughly matched to the expectation for future Belle II experimental data.The MC data are generated using nominal ABSZ form factors [19].In each fit the hadronic form factors in the unbinned likelihood function are varied, but taken to be the same for the di-muon and di-electron modes.Specifically, we randomly pick the parameters of the z-expansion in the ABSZ form factors within the uncertainties defined by the correct covariance matrices and thus automatically take into account correlations between parameters [? ].
An example result is shown in Fig. 15 for the δC 9 variable, where a clear correlation between δC 9 extracted in the muon and electron modes is seen.The δC 9 distribution for the di-muon mode as well as the difference with the di-electron mode is shown in Fig. 16.The average value of the difference is close to zero and its standard deviation is about 40% less than for the di-muon mode alone, thus demonstrating the cancellation of hadronic uncertainties.
To illustrate the magnitude of possible bias on fitted Wilson coefficients from ABSZ hadronic form factor uncertainties, we again perform pseudo-experiments with form factor parameters in the likelihood function that differ from the nominal parameters used to generate the fitted MC data.This time, however, we use the same parameters for each fit, so that the bias does not average out over multiple experiments.In Fig. 17 the δC 9 coefficient extraction is shown as an example.We find that the systematic bias in this coefficient is nearly twice as large as the statistical uncertainty of the fit.With the expected amount of Belle II data, SM theoretical uncertainties might mimic BSM physics effects.
Resonance effects are taken into account by modifying the Wilson coefficient C 9 according to Eq. A.4 both in the EvtGen event generator and in the likelihood function.
Result of 4-D unbinned likelihood fits for the δC 9 variable using both B → K * e + e − and B → K * µ + µ − MC data.Hadronic form factors are varied within their uncertainties.The fits are performed in the full q 2 range without including resonances.Note that the bias on the fitted δC 9 values, resulting from form factor uncertainties, is strongly correlated for the electron and muon modes.Figure 18 shows the q 2 distribution for B → K * µ + µ − decays.We find that resonances affect q 2 regions that extend well beyond the resonance widths, and thus rather large veto windows are required to mitigate their presence.
To illustrate how the resonances might affect the C 9 extraction even with large veto windows, we perform likelihood fits for the region q 2 ∈ [1, 8] ∪ [11, 12.5] ∪ [15, q 2 max ].For simplicity we use a likelihood function with two resonances only: J/ψ and ψ(2S).SM MC data is generated without resonances.With the above q 2 selection applied, the average  number of K * µ + µ − and K * e + e − decays in each fit is 4850.The total branching fraction is constrained with 5 % precision to the SM prediction.The results of the fits are shown in Fig. 19 (top row).The effect of the resonances is significant, biasing the fitted δC 9 values away from zero.This effect is particularly large for the imaginary part, where the bias is substantially larger than the SM value shown in Fig. 2. The ∆ variables shown Fig. 19 (bottom row), however, behave well and are centered close to zero with uncertainties much smaller than the original bias due to the resonances.This again suggests that possible BSM contributions to the muon mode can be extracted even in the presence of large theoretical uncertainties by using ∆ observables.
The remaining uncertainty then becomes statistical, rather than being limited by hadronic uncertainties.This highlights the importance of the technique utilizing the ∆ variables, and the importance of large experimental statistics for both muon and electron modes.

Sensitivity prospects with Belle II data
To estimate possible statistical sensitivity to BSM physics, we perform fits to generated data corresponding to several Belle II luminosity milestones.The generated data include contributions of cc resonances and thus we exclude the region ±0.25 GeV 2 /c 4 around the J/ψ and Ψ(2S) resonances from the fit data.In addition, for C 9 , C ′ 9 , C 10 , and C ′ 10 coefficients the region q 2 < 1 GeV 2 /c 4 is excluded to remove the photon pole effects whereas the fits to extract the δC 7 and C ′ 7 coefficients still use the photon pole data.The reconstruction efficiency is assumed to be about 25% based on the Belle results.The result of the fits are graphically shown in Fig. 20, 21, 22, and 23.We find that some fits are biased with low statistics albeit the bias is still significantly smaller than the statistical uncertainty and rapidly decreases with the integrated luminosity.As expected, the bias is the same both in the di-electron and di-muon modes and thus cancels in the ∆ observable.

Conclusions
We have upgraded the widely used MC event generator EvtGen to model B → K * ℓ + ℓ − with improved SM decay amplitudes and amplitudes for possible BSM physics contributions, implemented in the operator product expansion in terms of Wilson coefficients.This newly developed event generator model, the 'Sibidanov Physics Generator', implemented in the EvtGen framework, has been used to investigate the potential experimental signatures of lepton-flavor-violating BSM physics signals in B → K * ℓ + ℓ − .We describe the advantages and illustrate the potential of the new physics generator model, which implements the most general BSM contributions.This new generator can improve sensitivity of BSM experimental searches, and clarify BSM physics signatures.Advantages for experimental work include properly simulating BSM scenarios, interference between SM and BSM amplitudes, and correlations between various BSM signatures as well as acceptance bias.To address SM uncertainties, due to the strong interaction, which might mimic BSM effects, we have implemented b → ccs resonant contributions and up-to-date SM hadronic form factors [19].
In addition to these modifications to EvtGen, we also reviewed and optimized the computational efficiency of EvtGen for B → K * ℓ − ℓ + .We achieved a cumulative performance gain of two orders of magnitude.
We have shown that directly using the correlation between angular observables improves experimental statistical sensitivity.For example, in the case of BSM C 7 , C 9 , or C 10 contributions, there must be correlated signatures between A F B (q 2 ) and S 5 (q 2 ) angular asymmetries.The correlation in Eq. (3.1) depends on the Wilson coefficients.Therefore, the correlation improves the direct determination of Wilson coefficients.The use of a fourdimensional unbinned maximum likelihood fit to data naturally takes this correlation into account and improves the experimental sensitivity to BSM physics compared to independent binned fits to angular asymmetries.Note that theoretical fits to binned data reported in literature also take into account this correlation.
QCD uncertainties due to form factors and resonances currently limit the sensitivity of fits for BSM physics in B → K * ℓ + ℓ − .We have demonstrated the experimental feasibility of an approach using ∆-observables, which allow one to unambiguously distinguish between hadronic effects and BSM physics in the case of LFV BSM physics.The ∆-observables appear ideally suited for experimental analysis in Belle II with the large, 50 ab −1 , data sets expected in the next decade.For instance, ∆C 9 will be constrained to better than 3%.We find that Belle II also should have excellent statistical sensitivity to BSM physics in C 7 and C ′ 7 , which appear at low q 2 in the di-electron channel.For example, in the di-electron mode C ′ 7 should be constrained to 3% of C 7 .A new MC generator for direct simulation of Wilson coefficients in B → K * ℓ + ℓ − will be a useful tool for experimenters measuring these decays and determining the acceptance corrections for BSM physics.We demonstrated an application for lepton-flavor violating BSM physics.Further work will use this MC generator for experimental feasibility studies of ML (Machine Learning) based approaches to C i determination as well as to more difficult BSM scenarios with LFU (Lepton Flavor Universal) couplings in which the effects of interference with B → J/ψK * must be simulated and constrained.
Note Added: This paper is a significantly modified and improved version of Ref. [28], which was submitted to the US Community Summer Study on the Future of Particle Physics (Snowmass 2021).Ref. [28] will not, however, appear in the Snowmass proceedings.Improvements in this paper include, but are not limited to, calculations of correlations between several experimental observables and discussions of prospects for BSM-sensitive observables with several benchmark values of Belle II integrated luminosity.with z = 4m 2 q /q 2 , is related to the basic fermion loop.The limiting function h(q 2 , 0) is given by h(q 2 , 0) = 8 27 + 4 9 ln µ 2 q 2 + iπ .(A.3) One possible but simplistic approach to include the resonance effects is to use a Breit-Wigner ansatz where the h function in Eq.A.1 for the cc contributions should be replaced by [29] h where m V , Γ V total and Γ V had are the mass, total decay width and hadronic decay width of the vector boson V , respectively.Equation A.4 assumes no non-factorizable contributions, and therefore no relative strong phase between the two terms.Using this approach C 9 is modified by the cc contributions as shown in Fig. 24.The differential decay rate dependence on q 2 , cos θ ℓ , cos θ K , and χ assuming zero-width  -36 -

Figure 2 .
Figure 2. The SM Wilson coefficient C 9 versus q 2 now used in EvtGen (see Eq. 2.5) without cc resonances.

Figure 3 .
Figure 3. Result of the combined fit to LCSR and LQCD calculations from Ref. [19] for B → K * hadronic form factors.The shaded bands show the combined one-standard-deviation uncertainty of the fit according the provided covariance matrix.The points with error bars show LQCD results from Table11in Ref.[21].These theoretical uncertainties in the form factors can limit extraction of new physics.

Figure 4 .
Figure 4.The B → K * hadronic form factors A 2 and T 3 obtained using Eq.2.10 and 2.11.The shaded bands show the propagated one-standard-deviation uncertainty.These theoretical uncertainties in the form factors can limit extraction of new physics.

Figure 5 .
Figure 5.The measured time to generate 100000 B → K * e + e − events with (bottom bars) and without (top bar) EvtGen modifications.

Figure 6 .
Figure 6.Comparison of S 5 (left plot) and A FB (right plot) observables with BSM δC 9 = −0.87 and SM δC 9 = 0 in the di-muon mode.The points are generated with the BSM EvtGen simulation while the curves are the results of integrating the four-dimensional likelihood function.Note that no cc resonance effects are included.

Figure 7 .
Figure 7.Comparison of S 5 (left plot) and A FB (right plot) observables with several values of δC 7 in the di-muon mode obtained by integrating the four-dimensional likelihood function.Note that no cc resonance effects are included.

Figure 8 .
Figure 8.Comparison of S 5 (left plot) and A FB (right plot) observables with several values of δC 9 in the di-muon mode obtained by integrating the four-dimensional likelihood function.Note that no cc resonance effects are included.The expectation that δC 9 dependence at high q 2 should cancel out arises as a consequence of Eq.(3.1).The dependence at high q 2 only occurs for a very large value of δC 9 = −1.74.

Figure 9 .
Figure 9.Comparison of S 5 (left plot) and A FB (right plot) observables with several values of δC 10 in the di-muon mode obtained by integrating the four-dimensional likelihood function.Note that no cc resonance effects are included.

3. 4
C ′ 9 and C ′ 10 sensitivity Next, we investigate the statistical sensitivity to BSM physics in the form of right-handed currents.The right-handed contribution C ′ 9 can be constrained at level of 3 % of |C 9 | in the di-muon mode.Similarly, the right-handed contribution C ′ 10 can be constrained at level of 3 % of |C 10 | in the di-muon mode.Details of the fit results are shown in Appendix B.

Figure 14 .
Figure 14.The effect of C ′ 7 and δC 7 on angular distributions in the di-electron decay mode for q 2 < 2 GeV 2 /c 4 .There is a striking modulation in the χ angle distribution due to the right-handed BSM physics contribution as well as clear left-handed BSM signatures at low q 2 and in cos θ K .

Figure 16 .
Figure16.Distribution of the real part of the δC 9 coefficient obtained from fits to B → K * µ + µ − events (left plot) and the differences between those values and corresponding results for the B → K * e + e − mode (right plot) for randomly chosen hadronic form factors within the uncertainties as shown in Fig.15.Note that uncertainties are reduced by subtracting the fit results for the two modes.The remaining uncertainty is due to limited statistics in each fit.

Figure 17 .
Figure 17.Distribution of δC 9 obtained from unbinned likelihood fits to simulated B → K * µ + µ − SM decays with a randomly selected but fixed set of 7 hadronic form factors, consistent with the SM within theoretical uncertainties.The shift of the mean (0.27 in ℜ δC 9 ) from zero cannot be distinguished from BSM physics using this B → K * µ + µ − final state alone.

Figure 18 .
Figure18.The q 2 distribution of B → K * µ + µ − decay in the presence of cc resonances.The histogram is the result from the EvtGen generator, the green curve shows the result of the likelihood integration without resonances, and the red(magenta) curve is the result of the likelihood integration with the strong phase 0(π) when resonances are included.The contribution of these resonances (and non-factorizable effects) will be a limiting uncertainty in the extraction of BSM Wilson coefficients from B → K * µ + µ − .

Figure 19 .
Figure 19.A test of the resonance effects in the fit for the δC 9 parameter.The top row shows correlation plots for the real and imaginary parts of δC 9 .The bottom row shows the difference between modes.

Figure 20 .
Figure 20.The expected constraint on the δC 9 Wilson coefficient based on the unbinned maximum likelihood fit to MC data corresponding to various Belle II luminosity milestones.The mean value and its uncertainty (error bars) are shown in the left plot for the real part of the coefficient.The plot on the right shows the total uncertainty in the ∆C 9 variable.

Figure 21 .
Figure 21.The expected constraint on the δC 10 Wilson coefficient based on the unbinned maximum likelihood fit to MC data corresponding to the Belle II luminosity milestones.The mean value and its uncertainty (error bars) are shown in the left plot for the real part of the coefficient.The plot on the right shows the total uncertainty in the ∆C 10 variable.

Figure 22 .
Figure 22.The expected constraint on the δC 7 Wilson coefficient based on unbinned maximum likelihood fits to MC data corresponding to various Belle II luminosity milestones.The mean value and its uncertainty (error bars) are shown in the left plot for the real part of the Wilson coefficient.The plot on the right shows the total uncertainty in the ∆C 7 variable.

Figure 23 .
Figure 23.The expected constraint in the C ′ 7 Wilson coefficient based on unbinned maximum likelihood fits to MC data corresponding to various Belle II luminosity milestones.The mean value and its uncertainty (error bars) are shown in the left plot for the real part of the coefficient.The plot on the right shows the total uncertainty in the ∆C ′ 7 variable.

Figure 24 .
Figure 24.The SM Wilson coefficient C 9 versus q 2 (see Eq. 2.5) with cc resonances included.The top row shows the full range and the bottom row shows a zoomed-in version.The contribution of these resonances will be a limiting uncertainty in the extraction of BSM Wilson coefficients from B → K * µ + µ − .

Figure 32 .Figure 33 .
Figure 32.Distribution of δC 7 values obtained from unbinned likelihood fits to simulated B → K * e + e − SM decays in the full q 2 range, without including resonances.
Signatures of BSM Physics in B → K * ℓ + ℓ − and improved sensitivity via multi-dimensional likelihood fits 11 3.1 δC 9 correlated signatures in A FB and S 5 (and dB(B → K * ℓ + ℓ − ) / dq 2 ) 7 , C ′ 7 signatures at low q 2 and in the χ angle distribution 204 Disentangling QCD and resonance effects from BSM Physics in B → K * ℓ + ℓ

Table 1 .
Summary of the fit results for 1000 samples of 50 ab −1 assuming 25% experimental reconstruction efficiency.For extraction of the C 9 coefficient the samples are simulated with the NP case δC 9 = −0.87.Other values are extracted assuming the SM.