^a^ainstitutetext: University of Hawai‘i at Manoa, Honolulu, HI 96822, USA^b^binstitutetext: Department of Physics, Indian Institute of Technology Gandhinagar, Palaj, Gujarat 382355, India^c^cinstitutetext: Indian Institute of Technology Hyderabad (IITH), Telangana 502285, India

A new Monte Carlo Generator for BSM physics in $B\to K^{*}\ell^{+}\ell^{-}$ decays with an application to lepton non-universality in angular distributions

Alexei Sibidanov a Thomas E. Browder a Shawn Dubey a Shahab Kohani b Rusa Mandal c Saurabh Sandilya a Rahul Sinha a and Sven E. Vahsen

Abstract

Within the widely used EvtGen framework, we have added a new event generator model for $B\to K^{*}\ell^{+}\ell^{-}$ 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 estimate the statistical sensitivity of a simulated experiment 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 the 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\to K^{*}\ell^{+}\ell^{-}$ 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 $\Delta$ -observables significantly reduces uncertainties in the SM expectations due to QCD and resonance effects and is ideally suited for Belle II with the large data sets expected in the next decade. Thanks to the clean experimental environment of an $e^{+}e^{-}$ machine, Belle II should be able to probe BSM physics in the Wilson coefficients $C_{7}$ and $C_{7}^{\prime}$ , which appear at low $q^{2}$ in the di-electron channel.

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\to s$ transitions by multiple flavor physics experiments, LHCb, Belle and BaBar. Here we investigate the prospects for the decay $B\to K^{*}\ell^{+}\ell^{-}$ with $\ell=e,\mu$ . We have implemented a new Monte Carlo (MC) signal generator model inside the widely used event generator framework EvtGen Lange:2001uf , 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}^{\prime},C_{9}^{\prime},C_{10}^{\prime}$ . We allow for imaginary contributions to $\delta C_{i}$ . However, in this paper, we do not consider CP violating observables. Amplitudes for additional BSM physics contributions to $B\to K^{*}\ell^{+}\ell^{-}$ , which we call $\delta C_{i}=C_{i}^{\rm eff}-C_{i}^{\rm SM}$ , and which correspond to a general BSM signal resulting from dimension-six operators, can be chosen by the user. Until recently, non-zero $\delta C_{9}$ and $\delta_{C_{10}}$ were preferred by global theoretical fits to experimental data Altmannshofer:2021qrr ; Alguero:2021anc ; Hurth:2021nsi ; Ciuchini:2021smi ; Chrzaszcz:2018yza . While the significance of lepton-flavor universality violating contributions (LFUV) is reduced and large signals as indicated in [2-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 $\delta C_{9}$ and $\delta C_{10}$ for the di-muon final state are preferred by theoretical fits to current data Ciuchini:2022wbq ; Alguero:2023jeh; Hurth:2023jwr.

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\to K^{*}\ell^{+}\ell^{-}$ decays (Fig. 1): the invariant mass of the lepton pair, $q^{2}$ , the lepton helicity angle from the virtual $Z/\gamma$ -decay, $\cos\theta_{\ell}$ , the $K^{*}$ helicity angle, $\cos\theta_{K}$ , and the angle between the $K^{*}$ and di-lepton decay planes, $\chi$ . 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, $\delta C_{9}$ will lead to correlated signatures in the angular asymmetries $A_{FB}$ versus $q^{2}$ and $S_{5}$ versus $q^{2}$ . Another striking signature of $C_{7}^{\prime}$ is found in the modulation of the angle $\chi$ at low $q^{2}$ in $B\to K^{*}e^{+}e^{-}$ . Belle II should have excellent statistical sensitivity to such BSM physics in $C_{7}$ and $C_{7}^{\prime}$ .

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 $\delta 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\to K^{*}\ell^{+}\ell^{-}$ decays.

In addition to BSM contributions, the observed angular asymmetries could originate from SM QCD or $c\bar{c}$ resonance effects, which appear in the experimental final state $K^{*}\ell^{+}\ell^{-}$ . Given the difficulties in reliably computing all of the possible hadronic effects, the use of $\Delta$ -observables appears to be the only possible approach to unambiguously distinguish between these SM effects and lepton-universality violating BSM physics in $b\to s\mu^{+}\mu^{-}$ . 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 Belle:2016fev 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 $\Delta$ -observables are:

$\displaystyle\Delta A_{\rm FB}(B\to K^{*}\ell^{+}\ell^{-})$	$\displaystyle\equiv A_{\rm FB}(B\to K^{}\mu^{+}\mu^{-})-A_{\rm FB}(B\to K^{}e^{+}e^{-}),$	(1)
$\displaystyle\Delta S_{5}(B\to K^{*}\ell^{+}\ell^{-})$	$\displaystyle\equiv S_{5}(B\to K^{}\mu^{+}\mu^{-})-S_{5}(B\to K^{}e^{+}e^{-}),\ {\rm and}$	(2)
$\displaystyle\Delta C_{9}^{\rm eff}(B\to K^{*}\ell^{+}\ell^{-})$	$\displaystyle\equiv C_{9}^{\rm eff}(B\to K^{}\mu^{+}\mu^{-})-C_{9}^{\rm eff}(B\to K^{}e^{+}e^{-}).$	(3)

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 $\Delta$ -like variables was introduced in Ref. Capdevila:2016ivx . Belle reported a first attempt to measure $\Delta$ -type observables in $B\to K^{*}\ell^{+}\ell^{-}$ using a 0.7 ab^-1 data sample Belle:2016fev . The $\Delta$ -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 B2TIPbook .

Here we introduce a new Monte Carlo generator to enable further searches for BSM contributions in $B\to K^{*}\ell^{+}\ell^{-}$ . 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\to K^{*}e^{+}e^{-}$ and $B\to K^{*}\mu^{+}\mu^{-}$ 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\to K^{*}\ell^{+}\ell^{-}$ modes in LHCb datasets was given in MauriLFV:2018vbg . However, in contrast to MauriLFV:2018vbg , we have developed a MC simulation generator to be used by experimenters to directly explore the effects of BSM Wilson coefficients when large $B\to K^{*}\ell^{-}\ell^{+}$ 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.

2 Overview

2.1 Current experimental anomalies in $b\to 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 LHCb:2021trn ; LHCb:2017avl , in the low di-lepton mass-squared region. However, recent updated results LHCb:2023qnv 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\to K^{*}\mu^{+}\mu^{-}$ LHCb:2020lmf and $B_{s}\to\phi\mu^{+}\mu^{-}$ LHCb:2021zwz decays, which seem to be experimentally robust.

Several theoretical groups fit all $b\to s\ell^{+}\ell^{-}$ observables ( $R_{K}$ type results and angular asymmetries) for BSM Wilson coefficients, $\delta C_{i}$ , in $b\to s\mu^{+}\mu^{-}$ and suggest a coherent picture of BSM physics in $C_{9}$ and $C_{10}$ . BSM physics in $\delta C_{9}$ and $\delta C_{10}$ is found with a range of significances, though this varies between various global analyses, based on different assumptions of hadronic corrections Altmannshofer:2021qrr ; Alguero:2021anc ; Hurth:2021nsi ; Ciuchini:2021smi ; Ciuchini:2022wbq ; Alguero:2023jeh; Hurth:2023jwr.

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.

2.2 Definition of observables and angular asymmetries

For the $B\to K^{*}\ell^{+}\ell^{-}$ decay with $K^{*}\to K\pi$ the differential decay rate can be described in terms of the $K\pi$ invariant mass $m_{K\pi}$ , the di-lepton mass squared $q^{2}$ , and three angles $\theta_{\ell}$ , $\theta_{K}$ , and $\chi$ . The angle $\theta_{\ell}$ 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 $\theta_{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 $\chi$ 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.

Refer to caption — Figure 1: The $B\to K^{*}\ell^{+}\ell^{-}$ decay and the subsequent $K^{*}\to K\pi$ decay kinematic parameters.

In the SM angular asymmetries such as $A_{\rm 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_{\text{FB}}(q^{2})$ is defined as

\displaystyle A_{\text{FB}}(q^{2})=\dfrac{\left[\left(\displaystyle\int_{0}^{1}-\displaystyle\int_{-1}^{0}\right)\mathop{\textnormal{d}\cos{\theta_{\ell}}}\right]\mathop{\textnormal{d}(\Gamma-\bar{\Gamma})}}{\displaystyle\int_{-1}^{1}\mathop{\textnormal{d}\cos{\theta_{\ell}}}\mathop{\textnormal{d}(\Gamma+\bar{\Gamma})}},

(4)

where $\Gamma$ and $\bar{\Gamma}$ denote the decay rate of $\bar{B}^{0}\to\bar{K}^{0*}\ell^{+}\ell^{-}$ and the $CP$ -conjugate channel $B^{0}\to K^{0*}\ell^{+}\ell^{-}$ , respectively (see Eqs. 24 and 44 for full angular distributions). Other important angular asymmetries involve the angle $\chi$ between the $K^{*}$ and di-lepton decay planes Mandal:2014kma :

\displaystyle S_{4}(q^{2})=-\frac{\pi}{2}\frac{\Big{[}\displaystyle\int_{-\pi/2}^{\pi/2}-\displaystyle\int_{\pi/2}^{3\pi/2}\Big{]}\mathop{\textnormal{d}\chi}\Big{[}\displaystyle\int_{0}^{1}-\displaystyle\int_{-1}^{0}\Big{]}\mathop{\textnormal{d}\cos}\theta_{K}\Big{[}\int_{0}^{1}-\int_{-1}^{0}\Big{]}\mathop{\textnormal{d}\cos}\theta_{\ell}\mathop{\textnormal{d}(\Gamma+\bar{\Gamma})}}{\displaystyle\int_{0}^{2\pi}\mathop{\textnormal{d}\chi}\int_{-1}^{1}\mathop{\textnormal{d}\cos\theta_{K}}\int_{-1}^{1}\mathop{\textnormal{d}\cos\theta_{\ell}}\,\mathop{\textnormal{d}(\Gamma+\bar{\Gamma})}}\,,

(5)

and

\displaystyle S_{5}(q^{2})=\frac{4}{3}\frac{\Big{[}\displaystyle\int_{-\pi/2}^{\pi/2}-\displaystyle\int_{\pi/2}^{3\pi/2}\Big{]}\mathop{\textnormal{d}\chi}\Big{[}\int_{0}^{1}-\int_{-1}^{0}\Big{]}\mathop{\textnormal{d}\cos\theta_{K}}\displaystyle\int_{-1}^{1}\mathop{\textnormal{d}\cos\theta_{\ell}}~{}\mathop{\textnormal{d}(\Gamma-\bar{\Gamma})}}{\displaystyle\int_{0}^{2\pi}\mathop{\textnormal{d}\chi}\int_{-1}^{1}\mathop{\textnormal{d}\cos\theta_{K}}\int_{-1}^{1}\mathop{\textnormal{d}\cos\theta_{\ell}}\,\mathop{\textnormal{d}(\Gamma+\bar{\Gamma})}}.

(6)

Note that the angular observable $P_{5}^{\prime}$ , widely used in the literature, is related to $S_{5}$ via $P_{5}^{\prime}\equiv S_{5}/\sqrt{F_{L}(1-F_{L})}$ , where $F_{L}$ is the longitudinal polarization fraction of the $K^{*}$ meson.

2.3 SM and BSM Lorentz structures

The matrix element for the decay $B\to K^{*}\ell^{+}\ell^{-}$ , where $\ell$ is $e$ , $\mu$ or the $\tau$ lepton, is the following:

\begin{split}{\mathcal{M}}\ =\frac{G_{F}\alpha}{\sqrt{2}\pi}V_{tb}V_{ts}^{*}\bigg{\{}\bigg{[}\langle K^{*}|{\bar{s}\gamma^{\mu}({C_{9}^{\text{eff}}P_{L}+C_{9}^{\prime}P_{R}})b}|{\bar{B}}\rangle\\ -\frac{2m_{b}}{q^{2}}\langle{K^{*}}|{\bar{s}i\sigma^{\mu\nu}q_{\nu}(C_{7}^{\text{eff}}P_{R}+C_{7}^{\prime}P_{L})b}|{\bar{B}}\rangle\bigg{]}(\bar{\ell}\gamma_{\mu}\ell)\\ +\langle{K^{*}}|{\bar{s}\gamma^{\mu}({C_{10}P_{L}+C_{10}^{\prime}P_{R}})b}|{\bar{B}}\rangle(\bar{\ell}\gamma_{\mu}\gamma_{5}\ell)\bigg{\}},\end{split}

(7)

where the SM Wilson coefficients are known at NNLL accuracy Altmannshofer:2008dz :

$\displaystyle C_{7}^{\rm eff}$	$\displaystyle=$	$\displaystyle-0.304,$
$\displaystyle C_{9}^{\rm eff}$	$\displaystyle=$	$\displaystyle C_{9}+Y(q^{2})=4.211+Y(q^{2})\,,$	(8)
$\displaystyle C_{10}$	$\displaystyle=$	$\displaystyle-4.103\,,$

and where the function $Y(q^{2})$ is described in Section A. Note that we include the lepton masses ( $m_{\mu}$ and $m_{e}$ ) in all results. We assume BSM states are heavy, far above the $\mu=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}^{\prime}$ is highly suppressed while the $C_{9}^{\prime}$ $C_{10}^{\prime}$ operators are absent in the SM but can arise in the presence of BSM physics. Note that the SM Wilson coefficient $C_{9}^{\rm eff}$ varies with $q^{2}$ ; as shown in Fig. 2 separately for real and imaginary parts.

The $B\to K^{*}$ matrix elements in Eq. 7 can be expressed in terms of seven form factors that depend on the momentum transfer $q^{2}$ between the $B$ and the $K^{*}$ ( $q^{\mu}=p^{\mu}-k^{\mu}$ ):

\begin{split}\langle\bar{K}^{*}(k)|\bar{s}\gamma_{\mu}(1\mp\gamma_{5})b|\bar{B}(p)\rangle=\mp i\epsilon^{*}_{\mu}(m_{B}+m_{K^{*}})A_{1}(q^{2})\\ \pm i(2p-q)_{\mu}(\epsilon^{*}\cdot q)\,\frac{A_{2}(q^{2})}{m_{B}+m_{K^{*}}}\\ \pm iq_{\mu}(\epsilon^{*}\cdot q)\,\frac{2m_{K^{*}}}{q^{2}}\,\left[A_{3}(q^{2})-A_{0}(q^{2})\right]\\ +\epsilon_{\mu\nu\rho\sigma}\epsilon^{*\nu}p^{\rho}k^{\sigma}\,\frac{2V(q^{2})}{m_{B}+m_{K^{*}}},\end{split}

(9)

{\rm with\ }A_{3}(q^{2})=\frac{m_{B}+m_{K^{*}}}{2m_{K^{*}}}\,A_{1}(q^{2})-\frac{m_{B}-m_{K^{*}}}{2m_{K^{*}}}\,A_{2}(q^{2})\mbox{~{}~{}and~{}~{}}A_{0}(0)=A_{3}(0);

(10)

$\displaystyle\langle\bar{K}^{*}(k)\|\bar{s}\sigma_{\mu\nu}q^{\nu}(1\pm\gamma_{5})b\|\bar{B}(p)\rangle$	$\displaystyle=$	$\displaystyle i\epsilon_{\mu\nu\rho\sigma}\epsilon^{*\nu}p^{\rho}k^{\sigma}\,2T_{1}(q^{2})$	(11)
		$\displaystyle\pm T_{2}(q^{2})\left[\epsilon^{}_{\mu}(m_{B}^{2}-m_{K^{}}^{2})-(\epsilon^{*}\cdot q)\,(2p-q)_{\mu}\right]$
		$\displaystyle\pm T_{3}(q^{2})(\epsilon^{}\cdot q)\left[q_{\mu}\!-\displaystyle\frac{q^{2}}{m_{B}^{2}-m_{K^{}}^{2}}\,(2p-q)_{\mu}\right],$

with $T_{1}(0)=T_{2}(0)$ . Here $\epsilon_{\mu}$ is the polarization vector of the $K^{*}$ . We used the following convention for the Levi-Civita tensor $\epsilon_{0123}=+1$ .

We adopt the formalism developed in Ref. Bharucha:2015bzk where the extrapolation of form factors from the calculated LCSR input points ( $q^{2}\lesssim 0\,$ GeV) to larger $q^{2}$ values is performed by a simple pole form with a $z$ -expansion

F_{i}^{B\to K^{*}}(q^{2})\equiv\frac{1}{1-q^{2}/m_{R,i}^{2}}\,\sum_{n=0}^{2}\alpha_{n}^{i}\left[z(q^{2})-z(0)\right]^{n}\,,

(12)

where $z(t)\equiv\displaystyle\frac{\sqrt{t_{+}-t}-\sqrt{t_{+}-t_{0}}}{\sqrt{t_{+}-t}+\sqrt{t_{+}-t_{0}}}$ , $t_{\pm}=(m_{B}\pm m_{K^{*}})^{2}$ and $t_{0}\equiv t_{+}\left(1-\sqrt{1-t_{-}/t_{+}}\right)$ .

The values for the fit parameters $\alpha_{0,1,2}^{i}$ 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. Bharucha:2015bzk .

2.4 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. Ball:2004rg . In our new EvtGen model BTOSLLNP (see Section 2.5), we have now implemented the most recent hadronic form factors from Ref. Bharucha:2015bzk , 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. 12.

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. Bharucha:2015bzk , where the form factors $A_{2}$ and $T_{3}$ are extracted using the following expression:

A_{12}(q^{2})=\frac{(m_{B}+m_{K\pi})^{2}(m_{B}^{2}-m_{K\pi}^{2}-q^{2})A_{1}(q^{2})-\lambda(q^{2})A_{2}(q^{2})}{16m_{B}m_{K\pi}^{2}(m_{B}+m_{K\pi})},

(13)

and

T_{23}(q^{2})=\frac{(m_{B}^{2}-m_{K\pi}^{2})^{2}(m_{B}^{2}+3m_{K\pi}^{2}-q^{2})T_{2}(q^{2})-\lambda(q^{2})T_{3}(q^{2})}{8m_{B}m_{K\pi}^{2}(m_{B}-m_{K\pi})}.

(14)

Here $\lambda(q^{2})$ is the Källén-function defined in Eq. 43 and $m_{K\pi}$ 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.

2.5 BSM Modifications to the EvtGen Monte Carlo Generator

We implement the $B\to K^{*}\ell^{+}\ell^{-}$ 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 Kuhr:2018lps ). In the BTOSLLNP model, the user can specify the following BSM input parameters: $\delta C_{7}$ , $C_{7}^{\prime}$ , $\delta C_{9}$ , $C_{9}^{\prime}$ , $\delta C_{10}$ , $C_{10}^{\prime}$ , $C_{S}-C_{S}^{\prime}$ , and $C_{P}-C_{P}^{\prime}$ . 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 $B\to K^{*}\ell^{+}\ell^{-}$ MC generator:

## the first argument is the Cartesian(0) or polar(1) representation of
## complex BSM coefficients, which are three consecutive numbers
##   {id, Re(C), Im(C)} or {coeff id, |C|, Arg(C)}

## id==0 delta C_7  -- BSM addition to NNLO SM value
## id==1 delta C_9  -- BSM addition to NNLO SM value
## id==2 delta C_10 -- BSM addition to NNLO SM value
## id==3 C’_7  -- BSM right-handed coefficient
## id==4 C’_9  -- BSM right-handed coefficient
## id==5 C’_10 -- BSM right-handed coefficient
## id==6 (C_S - C’_S) -- BSM scalar left- and right-handed coefficient
## id==7 (C_P - C’_P) -- BSM pseudo-scalar left- and right-handed
##                       coefficient


Decay anti-B0
## The SM Case
1.000 anti-K*0 e+ e- BTOSLLNP;
Enddecay

Decay anti-B0
## delta C_9eff = (-0.87, 0.0) all other coefficients correspond to the
## SM values
1.000 anti-K*0 e+ e- BTOSLLNP 0 1 -0.87 0.0 ;
Enddecay

Decay anti-B0
## delta C_9eff = (-0.87, 0.0)  and delta C_10eff = (+0.87, 0.0)
## all other coefficients correspond to the
## SM values
1.000 anti-K*0 e+ e- BTOSLLNP 0 1 -0.87 0.0 2 0.87 0.0 ;
Enddecay

Decay anti-B0
## delta C_7prime eff = ( 0.04, 0.0) all other coefficients correspond to the
## SM values
1.000 anti-K*0 e+ e- BTOSLLNP 0 3 0.04 0.0 ;
Enddecay

2.6 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\to K^{*}\ell^{+}\ell^{-}$ 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. 24, where the matrix element defined in Eq. 7 is already squared, and which explicitly depends only on the four variables $q^{2}$ , $\cos{\theta_{\ell}}$ , $\cos{\theta_{K}}$ , $\chi$ , 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¹¹1Since this is a simplified sensitivity study we do not take into account resonances and thus results in resonance region are not physical..

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}\to K^{*0}(K^{+}\pi^{-})\ell^{+}\ell^{-}$ 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 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:

\delta C_{i}=C_{i}^{\text{BSM}}-C_{i}^{\text{SM}}.

The SM values are taken from NNLO calculations Altmannshofer:2008dz . 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 $\delta C_{9}$ correlated signatures in $A_{\rm FB}$ and $S_{5}$ (and $\mathop{\textnormal{d}\mathcal{B}(B\to K^{*}\ell^{+}\ell^{-})}/\mathop{\textnormal{d}q^{2}}$ )

Figure 6 shows how the observables $S_{5}$ and $A_{\text{FB}}$ are modified when a BSM contribution ( $\delta C_{9}$ ) to $C_{9}$ is present. We assume that $\delta C_{9}=-0.87$ , which corresponds to a representative value from global theory fits to $b\to s$ experimental data in the $\mu\mu$ channel Altmannshofer:2021qrr . There are clear differences between the SM and BSM distributions of $S_{5}$ and $A_{\text{FB}}$ , and these signatures are correlated. However, we do not find any significant effect of $\delta C_{9}$ on the differential distribution $\mathop{\textnormal{d}\mathcal{B}(B\to K^{*}\ell^{+}\ell^{-})}/\mathop{\textnormal{d}q^{2}}$ .

The observables $S_{5}$ and $A_{\text{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}^{\prime}$ , $C_{9}^{\prime}$ , $C_{10}^{\prime}$ and $C_{S}$ are zero) the correlations observed in the results from the Monte Carlo generator between $S_{5}$ and $A_{\text{FB}}$ in Figs. 7, 8, and 9, can be understood from the following analytical expression relating the two observables:

\frac{S_{5}}{A_{\text{FB}}}=-\displaystyle\frac{16\sqrt{2}}{3\sqrt{q^{2}}}\frac{A_{12}}{A_{1}}\frac{m_{B}m_{K^{*}}}{m_{B}+m_{K^{*}}}\displaystyle\frac{{\cal R}\text{e}(C_{9})+\displaystyle\frac{m_{b}}{q^{2}}\text{Re}(C_{7})Z_{1}}{\text{Re}(C_{9})+\displaystyle\frac{m_{b}}{q^{2}}\text{Re}(C_{7})Z_{2}}

(15)

where $Z_{1}$ and $Z_{2}$ are purely functions of form factors and masses and are given by

Z_{1}=\displaystyle\frac{2T_{1}A_{12}(m_{B}+m_{K^{*}})^{2}+T_{23}Vq^{2}}{2VA_{12}(m_{B}+m_{K^{*}})},

(16)

Z_{2}=\displaystyle(m_{B}+m_{K^{*}})\frac{T_{1}}{V}+(m_{B}-m_{K^{*}})\frac{T_{2}}{A_{1}}.

(17)

It is interesting to note that Eq. (15) is independent of $\text{Im}(C_{7})$ , $\text{Im}(C_{9})$ and $C_{10}$ . Eq. (15) implies that the sensitivity to both $\text{Re}(C_{7})$ , and $\text{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}\approx 0$ . The region $q^{2}=(1-8)~{}\text{G\text{e\kern-0.6458ptV}}^{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_{\text{FB}}$ .

3.2 Unbinned four-dimensional maximum likelihood fitting for BSM

Using a high statistics sample of $B\to K^{*}e^{+}e^{-}$ and $B\to K^{*}\mu^{+}\mu^{-}$ MC simulated events, we can compare the distributions in the four kinematic variables for the SM and a BSM scenario with $\delta C_{9}=-0.87$ in three selected regions of $q^{2}$ ( $q^{2}<1~{}\rm{GeV}^{2}$ , $1<q^{2}<6~{}\rm{GeV}^{2}$ and $q^{2}>15~{}\rm{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\theta_{K^{*}}$ and $\cos\theta_{\ell}$ . In all three regions, there are changes in the forward-backward asymmetry of the leptonic system, which are visible in the $\cos\theta_{\ell}$ distributions. At low $q^{2}$ in $B\to K^{*}e^{+}e^{-}$ , the SM pole contribution arising from ${\cal O}_{7}$ dominates and BSM effects are difficult to distinguish. However, in $B\to K^{*}\mu^{+}\mu^{-}$ 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 Altmannshofer:2021qrr ; Alguero:2021anc ; Hurth:2021nsi ; Ciuchini:2021smi that determines $\delta C_{i}$ , we fit directly for $\delta C_{i}$ using 4D unbinned likelihood fits to $q^{2}$ , $\cos{\theta_{\ell}}$ , $\cos{\theta_{K}}$ , and $\chi$ with $\delta C_{i}$ as a free parameter. For example, Fig. 13 shows the distribution of $\delta C_{9}$ resulting from about 2000 pseudo-experiments performing such multidimensional fits repeatedly, assuming for now that $\delta C_{9}=0$ to check the statistical sensitivity. According to the fit results the real and imaginary parts of $\delta 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 $\delta 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.

3.3 $\delta C_{10}$ sensitivity

The statistical sensitivity to $\delta C_{10}$ is estimated from our multidimensional fit. According to the fit results the real and imaginary parts of $\delta C_{10}$ can be constrained to be 0.21 and 0.28, corresponding to 5 and 7 % of $|C_{10}|\approx 4.1$ , respectively. The corresponding distributions are shown in Appendix B.

3.4 $C_{9}^{\prime}$ and $C_{10}^{\prime}$ sensitivity

Next, we investigate the statistical sensitivity to BSM physics in the form of right-handed currents. The right-handed contribution $C_{9}^{\prime}$ can be constrained at level of 3 % of $|C_{9}|$ in the di-muon mode. Similarly, the right-handed contribution $C_{10}^{\prime}$ 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.

3.5 $C_{7}$ , $C_{7}^{\prime}$ signatures at low $q^{2}$ and in the $\chi$ angle distribution

Sensitivity in the di-electron mode at low $q^{2}$ , $q^{2}<2\,\text{G\text{e\kern-0.6458ptV}}^{2}/c^{4}$ , is critical as the new physics terms interfere with the photon pole. Even moderate non-zero $C_{7}^{\prime}$ (right-handed BSM physics) produces a prominent modulation in the $\chi$ distribution, while left-handed BSM physics ( $\delta C_{7}$ ) signatures appear at low $q^{2}$ and in the $\cos\theta_{K}$ distribution. These BSM physics signatures are shown in Fig. 14.

The SM NNLO value of $|C_{7}|$ is $\approx 0.3$ so we expect, according to our study, to constrain the real and imaginary parts of $C_{7}$ to better than 2.5 and 7.5 %, respectively from the di-muon mode alone. The contribution from BSM physics in the form of right-handed currents, which corresponds to $C_{7}^{\prime}$ , is expected to be constrained to better than 6 % of $C_{7}$ from the di-muon mode alone. In the di-electron mode where we have better access to the photon pole we might constrain the real and imaginary parts of $C_{7}$ to better than 1.5 and 6.5 %, respectively and $C_{7}^{\prime}$ will be constrained to better than 3 % of $|C_{7}|$ . The distributions for these fits and the resulting summary table are given in Appendix B.

4 Disentangling QCD and resonance effects from BSM Physics in $B\to K^{*}\ell^{+}\ell^{-}$

The semileptonic decays $B\to K^{*}\ell^{+}\ell^{-}$ 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^{*}\ell^{+}\ell^{-}$ includes contributions from $c\bar{c}$ resonances such as $B\to J/\psi K^{*}$ , $B\to\psi(2S)K^{*}$ as well as other broader $c\bar{c}$ resonances with higher mass, where the $c\bar{c}$ resonances decay to $\ell^{+}\ell^{-}$ . It is well known that these $b\to c\bar{c}s$ modes have non-factorizable contributions Belle:2002otd ; BaBar:2007rbr ; LHCb:2013vga . In addition, complications arise due to non-local contributions from electromagnetic corrections to purely hadronic operators Beneke:2001at , which cannot be calculated from first principles. Only limited success has been achieved in computing the effects of charm-loop contributions Khodjamirian:2010vf . Given the difficulties in reliably computing all the QCD effects, our choice of $\Delta$ observablesB2TIPbook 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 Bharucha:2015bzk . 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 Bharucha:pc.

An example result is shown in Fig. 15 for the $\delta C_{9}$ variable, where a clear correlation between $\delta C_{9}$ extracted in the muon and electron modes is seen. The $\delta 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 $\delta 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. 23 both in the EvtGen event generator and in the likelihood function. Figure 18 shows the $q^{2}$ distribution for $\bar{B}\to\bar{K}^{*}\mu^{+}\mu^{-}$ 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}\in[1,8]\cup[11,12.5]\cup[15,q^{2}_{\text{max}}]$ . For simplicity we use a likelihood function with two resonances only: $J/\psi$ and $\psi(2S)$ . SM MC data is generated without resonances. With the above $q^{2}$ selection applied, the average number of $K^{*}\mu^{+}\mu^{-}$ 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 $\delta 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 $\Delta$ 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 $\Delta$ observables.

The remaining uncertainty then becomes statistical, rather than being limited by hadronic uncertainties. This highlights the importance of the technique utilizing the $\Delta$ variables, and the importance of large experimental statistics for both muon and electron modes.

5 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 $c\bar{c}$ resonances and thus we exclude the region $\pm 0.25\ \text{GeV}^{2}/c^{4}$ around the $J/\psi$ and $\Psi(2S)$ resonances from the fit data. In addition, for $C_{9}$ , $C_{9}^{\prime}$ , $C_{10}$ , and $C_{10}^{\prime}$ coefficients the region $q^{2}<1\,\text{GeV}^{2}/c^{4}$ is excluded to remove the photon pole effects whereas the fits to extract the $\delta C_{7}$ and $C_{7}^{\prime}$ 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 $\Delta$ observable.

6 Conclusions

We have upgraded the widely used MC event generator EvtGen to model $B\to K^{*}\ell^{+}\ell^{-}$ 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\to K^{*}\ell^{+}\ell^{-}$ . 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\to c\bar{c}s$ resonant contributions and up-to-date SM hadronic form factors Bharucha:2015bzk .

In addition to these modifications to EvtGen, we also reviewed and optimized the computational efficiency of EvtGen for $B\to K^{*}\ell^{-}\ell^{+}$ . 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_{FB}(q^{2})$ and $S_{5}(q^{2})$ angular asymmetries. The correlation in Eq. (15) depends on the Wilson coefficients. Therefore, the correlation improves the direct determination of Wilson coefficients. The use of a four-dimensional 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\to K^{*}\ell^{+}\ell^{-}$ . We have demonstrated the experimental feasibility of an approach using $\Delta$ -observables, which allow one to unambiguously distinguish between hadronic effects and BSM physics in the case of LFV BSM physics. The $\Delta$ -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, $\Delta 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}^{\prime}$ , which appear at low $q^{2}$ in the di-electron channel. For example, in the di-electron mode $C_{7}^{\prime}$ should be constrained to 3% of $C_{7}$ . A new MC generator for direct simulation of Wilson coefficients in $B\to K^{*}\ell^{+}\ell^{-}$ 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\to J/\psi K^{*}$ must be simulated and constrained.

Note Added: This paper is a significantly modified and improved version of Ref. Sibidanov:2022gvb , which was submitted to the US Community Summer Study on the Future of Particle Physics (Snowmass 2021). Ref. Sibidanov:2022gvb 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.

7 Acknowledgments

We thank Emi Kou (IJCLab) for useful discussions concerning correlations between angular observables. T.E.B., S.D., S.K., A.S., and S.E.V. thank the DOE Office of High Energy Physics for support through DOE grant DE-SC0010504. R.M. acknowledges the support of SERB Grant SPG/2022/001238. S.S. acknowledges the support of SERB Grant SRG/2022/001608.

Appendix A Description of the theoretical framework

The matrix element introduced in Eq. 7 can further be extended via the inclusion of four dimension six operators namely

\mathcal{O}_{S}^{(\prime)}=(\bar{s}P_{R(L)}b)(\bar{\mu}\mu)~{}~{}{\rm and}~{}~{}\mathcal{O}_{P}^{(\prime)}=(\bar{s}P_{R(L)}b)(\bar{\mu}\gamma_{5}\mu)

with the corresponding Wilson coefficients $C_{S}^{(\prime)}$ and $C_{P}^{(\prime)}$ which can arise in presence of BSM physics. The SM prediction for $C_{i}$ is known to NNLL accuracy and in the amplitude for $B\to K^{*}(\to K\pi)\ell^{+}\ell^{-}$ certain combinations of these Wilson coefficients appear, which are given by Altmannshofer:2008dz

$\displaystyle C_{7}^{\rm eff}$	$\displaystyle=$	$\displaystyle C_{7}-\frac{1}{3}\,C_{3}-\frac{4}{9}\,C_{4}-\frac{20}{3}\,C_{5}\,-\frac{80}{9}\,C_{6}\,,$
$\displaystyle C_{9}^{\rm eff}$	$\displaystyle=$	$\displaystyle C_{9}+Y(q^{2})\,,{\rm with}$
$\displaystyle Y(q^{2})$	$\displaystyle=$	$\displaystyle h(q^{2},m_{c})\left(\frac{4}{3}\,C_{1}+C_{2}+6C_{3}+60C_{5}\right)$	(18)
		$\displaystyle{}-\frac{1}{2}\,h(q^{2},m_{b})\left(7C_{3}+\frac{4}{3}\,C_{4}+76C_{5}+\frac{64}{3}\,C_{6}\right)$
		$\displaystyle{}-\frac{1}{2}\,h(q^{2},0)\left(C_{3}+\frac{4}{3}\,C_{4}+16C_{5}+\frac{64}{3}\,C_{6}\right)$
		$\displaystyle{}+\frac{4}{3}\,C_{3}+\frac{64}{9}\,C_{5}+\frac{64}{27}\,C_{6}\,.$

The function

	$\displaystyle h(q^{2},m_{q})=-\frac{4}{9}\,\left(\ln\,\frac{m_{q}^{2}}{\mu^{2}}-\frac{2}{3}-z\right)\qquad{}$
	$\displaystyle\quad{}-\frac{4}{9}\,(2+z)\sqrt{\|z-1\|}\times\left\{\begin{array}[]{l@{\quad}l}\displaystyle\arctan\,\frac{1}{\sqrt{z-1}}&z>1\\[10.0pt] \displaystyle\ln\,\frac{1+\sqrt{1-z}}{\sqrt{z}}-\frac{i\pi}{2}&z\leq 1\end{array}\right.$		(21)

with $z=4m_{q}^{2}/q^{2}$ , is related to the basic fermion loop. The limiting function $h(q^{2},0)$ is given by

h(q^{2},0)=\frac{8}{27}+\frac{4}{9}\left(\ln\frac{\mu^{2}}{q^{2}}+i\pi\right).

(22)

One possible but simplistic approach to include the resonance effects is to use a Breit-Wigner ansatz where the $h$ function in Eq. 18 for the $c\bar{c}$ contributions should be replaced by Kruger:1996cv

	$\displaystyle h(m_{c},q^{2})\rightarrow h(m_{c},q^{2})$
	$\displaystyle\qquad{}-\frac{3\pi}{\alpha^{2}}\sum_{V=J/\psi,\psi^{\prime},..}\frac{m_{V}\text{Br}(V\rightarrow\ell^{+}\ell^{-})\Gamma^{V}_{\text{total}}}{q^{2}-m_{V}^{2}+im_{V}\Gamma^{V}_{\text{total}}}\,,$		(23)

where $m_{V}$ , $\Gamma^{V}_{\text{total}}$ and $\Gamma^{V}_{\text{had}}$ are the mass, total decay width and hadronic decay width of the vector boson $V$ , respectively. Equation 23 assumes no non-factorizable contributions, and therefore no relative strong phase between the two terms. Using this approach $C_{9}$ is modified by the $c\bar{c}$ contributions as shown in Fig. 24.

The differential decay rate dependence on $q^{2}$ , $\cos\theta_{\ell}$ , $\cos\theta_{K}$ , and $\chi$ assuming zero-width $K^{*}$ can be written as

	$\displaystyle\displaystyle\frac{\mathop{\textnormal{d}}^{4}\Gamma(\bar{B}^{0}\to\bar{K}^{0*}\ell^{+}\ell^{-})}{\mathop{\textnormal{d}q^{2}}\,\mathop{\textnormal{d}\cos{\theta_{\ell}}}\,\mathop{\textnormal{d}\cos{\theta_{K}}}\,\mathop{\textnormal{d}\chi}}=I(q^{2},{\theta_{\ell}},{\theta_{K}},\chi)=$		(24)
	$\displaystyle\frac{9}{32\pi}\Big{[}I_{1}^{s}\sin^{2}{\theta_{K}}+I_{1}^{c}\cos^{2}{\theta_{K}}+(I_{2}^{s}\sin^{2}{\theta_{K}}+I_{2}^{c}\cos^{2}{\theta_{K}})\cos 2{\theta_{\ell}}$
	$\displaystyle+I_{3}\sin^{2}{\theta_{K}}\sin^{2}{\theta_{\ell}}\cos 2\chi+I_{4}\sin 2{\theta_{K}}\sin 2{\theta_{\ell}}\cos\chi+I_{5}\sin 2{\theta_{K}}\sin{\theta_{\ell}}\cos\chi$
	$\displaystyle+I_{6}^{s}\sin^{2}{\theta_{K}}\cos{\theta_{\ell}}+I_{6}^{c}\cos^{2}{\theta_{K}}\cos{\theta_{\ell}}+I_{7}\sin 2{\theta_{K}}\sin{\theta_{\ell}}\sin\chi$
	$\displaystyle+I_{8}\sin 2{\theta_{K}}\sin 2{\theta_{\ell}}\sin\chi+I_{9}\sin^{2}{\theta_{K}}\sin^{2}{\theta_{\ell}}\sin 2\chi\Big{]}\,.$

The angular coefficients can be expressed in terms of transversity amplitudes as

$\displaystyle I_{1}^{s}=$	$\displaystyle\frac{(2+\beta^{2})}{4}\Big{[}\|{{\cal A}_{\perp}^{L}}\|^{2}+\|{{\cal A}_{\parallel}^{L}}\|^{2}+(L\to R)\Big{]}+\frac{4m^{2}}{q^{2}}\Re({{\cal A}_{\perp}^{L}}{{\cal A}_{\perp}^{R}}^{}+{{\cal A}_{\parallel}^{L}}{{\cal A}_{\parallel}^{R}}^{}),$	(25)
$\displaystyle I_{1}^{c}=$	$\displaystyle\|{{\cal A}_{0}^{L}}\|^{2}\!+\!\|{{\cal A}_{0}^{R}}\|^{2}\!+\!\frac{4m^{2}}{q^{2}}\big{[}\left\|{\cal A}_{t}\right\|^{2}\!+\!2\Re({{\cal A}_{0}^{L}}{{\cal A}_{0}^{R}}^{*})\big{]}+\beta^{2}{\|\cal{A}_{S}\|}^{2},$	(26)
$\displaystyle I_{2}^{s}=$	$\displaystyle\frac{\beta^{2}}{4}\Big{[}\|{{\cal A}_{\perp}^{L}}\|^{2}+\|{{\cal A}_{\parallel}^{L}}\|^{2}+(L\to R)\Big{]},$	(27)
$\displaystyle I_{2}^{c}=$	$\displaystyle-\beta^{2}\Big{[}\|{{\cal A}_{0}^{L}}\|^{2}+(L\to R)\Big{]},$	(28)
$\displaystyle I_{3}=$	$\displaystyle\frac{\beta^{2}}{2}\Big{[}\|{{\cal A}_{\perp}^{L}}\|^{2}-\|{{\cal A}_{\parallel}^{L}}\|^{2}+(L\to R)\Big{]},$	(29)
$\displaystyle I_{4}=$	$\displaystyle\frac{\beta^{2}}{\sqrt{2}}\Big{[}\Re({{\cal A}_{0}^{L}}{{\cal A}_{\parallel}^{L}}^{*})+(L\to R)\Big{]},$	(30)
$\displaystyle I_{5}=$	$\displaystyle\sqrt{2}\beta\Big{[}\Re({{\cal A}_{0}^{L}}{{\cal A}_{\perp}^{L}}^{})-(L\to R)-\frac{m}{\sqrt{q^{2}}}\,\Re({{\cal A}_{\parallel}^{L}}{A_{S}^{}}+{{\cal A}_{\parallel}^{R}}{A_{S}^{*}})\Big{]},$	(31)
$\displaystyle I_{6}^{s}=$	$\displaystyle 2\beta\Big{[}\Re({{\cal A}_{\parallel}^{L}}{{\cal A}_{\perp}^{L}}^{*})-(L\to R)\Big{]},$	(32)
$\displaystyle I_{6}^{c}=$	$\displaystyle 4\beta\frac{m}{\sqrt{q^{2}}}\,\Re\left[{{\cal A}_{0}^{L}}{A_{S}^{*}}+(L\to R)\right],$	(33)
$\displaystyle I_{7}=$	$\displaystyle\sqrt{2}\beta\Big{[}\Im({{\cal A}_{0}^{L}}{{\cal A}_{\parallel}^{L}}^{})-(L\to R)+\frac{m}{\sqrt{q^{2}}}\,{\Im}({{\cal A}_{\perp}^{L}}{A_{S}^{}}+{{\cal A}_{\perp}^{R}}{A_{S}^{*}})\Big{]},$	(34)
$\displaystyle I_{8}=$	$\displaystyle\frac{1}{\sqrt{2}}\beta^{2}\Big{[}\Im({{\cal A}_{0}^{L}}{{\cal A}_{\perp}^{L}}^{*})+(L\to R)\Big{]},$	(35)
$\displaystyle I_{9}=$	$\displaystyle\beta^{2}\Big{[}\Im({{\cal A}_{\parallel}^{L}}^{*}{{\cal A}_{\perp}^{L}})+(L\to R)\Big{]},$	(36)

where

\beta=\sqrt{1-\frac{4\,m^{2}}{q^{2}}}

and we have dropped the explicit $q^{2}$ dependence of the transversity amplitudes ${\cal A}_{\perp,\parallel,0}^{L,R}$ and ${\cal A}_{t}$ for notational simplicity. Note that throughout the paper, we retain lepton masses and include non-zero lepton masses.

$\displaystyle A_{\perp L,R}$	$\displaystyle=N\sqrt{2}\lambda^{1/2}\bigg{[}\left[(C_{9}^{\text{eff}}+C_{9}^{{\text{eff}}\prime})\mp(C_{10}^{\text{eff}}+C_{10}^{{\text{eff}}\prime})\right]\frac{V(q^{2})}{m_{B}+m_{K^{\!*}}}$
	$\displaystyle+\frac{2m_{b}}{q^{2}}(C_{7}^{\text{eff}}+C_{7}^{{\text{eff}}\prime})T_{1}(q^{2})\bigg{]},$	(37)
$\displaystyle A_{\parallel L,R}$	$\displaystyle=-N\sqrt{2}(m_{B}^{2}-m_{K^{\!}}^{2})\bigg{[}\left[(C_{9}^{\text{eff}}-C_{9}^{{\text{eff}}\prime})\mp(C_{10}^{\text{eff}}-C_{10}^{{\text{eff}}\prime})\right]\frac{A_{1}(q^{2})}{m_{B}-m_{K^{\!}}}$
	$\displaystyle+\frac{2m_{b}}{q^{2}}(C_{7}^{\text{eff}}-C_{7}^{{\text{eff}}\prime})T_{2}(q^{2})\bigg{]},$	(38)
$\displaystyle A_{0L,R}$	$\displaystyle=-\frac{N}{2m_{K^{\!*}}\sqrt{q^{2}}}\bigg{\{}\left[(C_{9}^{\text{eff}}-C_{9}^{{\text{eff}}\prime})\mp(C_{10}^{\text{eff}}-C_{10}^{{\text{eff}}\prime})\right]$
	$\displaystyle\qquad\times\bigg{[}(m_{B}^{2}-m_{K^{\!}}^{2}-q^{2})(m_{B}+m_{K^{\!}})A_{1}(q^{2})-\lambda\frac{A_{2}(q^{2})}{m_{B}+m_{K^{\!*}}}\bigg{]}$
	$\displaystyle\qquad+{2m_{b}}(C_{7}^{\text{eff}}-C_{7}^{{\text{eff}}\prime})\bigg{[}(m_{B}^{2}+3m_{K^{\!}}^{2}-q^{2})T_{2}(q^{2})-\frac{\lambda}{m_{B}^{2}-m_{K^{\!}}^{2}}T_{3}(q^{2})\bigg{]}\bigg{\}},$	(39)
$\displaystyle A_{t}$	$\displaystyle=\frac{N}{\sqrt{q^{2}}}\lambda^{1/2}\left[2(C_{10}^{\text{eff}}-C_{10}^{{\text{eff}}\prime})+\frac{q^{2}}{m_{\mu}}(C_{P}-C_{P}^{\prime})\right]A_{0}(q^{2}),$	(40)
$\displaystyle A_{S}$	$\displaystyle=-2N\lambda^{1/2}(C_{S}-C_{S}^{\prime})A_{0}(q^{2}),$	(41)

where

N=V_{tb}^{\vphantom{*}}V_{ts}^{*}\left[\frac{G_{F}^{2}\alpha^{2}}{3\cdot 2^{10}\pi^{5}m_{B}^{3}}q^{2}\lambda^{1/2}\beta\right]^{1/2},

(42)

with

\lambda=\lambda(q^{2})=m_{B}^{4}+m_{K^{*}}^{4}+q^{4}-2(m_{B}^{2}m_{K^{*}}^{2}+m_{K^{*}}^{2}q^{2}+m_{B}^{2}q^{2}).

(43)

The differential distribution for the $CP$ -conjugate mode $B^{0}\to K^{0*}\ell^{+}\ell^{-}$ can be written as

\displaystyle\frac{\mathop{\textnormal{d}}^{4}\Gamma(B^{0}\to K^{0*}\ell^{+}\ell^{-})}{\mathop{\textnormal{d}q^{2}}\,\mathop{\textnormal{d}\cos{\theta_{\ell}}}\,\mathop{\textnormal{d}\cos{\theta_{K}}}\,\mathop{\textnormal{d}\chi}}

\displaystyle=\bar{I}(q^{2},{\theta_{\ell}},{\theta_{K}},\chi)\,,

(44)

where $\bar{I}(q^{2},{\theta_{\ell}},{\theta_{K}},\chi)$ is obtained with the following replacements in Eq. 24: $I_{1,2,3,4,7}\to\bar{I}_{1,2,3,4,7}$ and $I_{5,6,8,9}\to-\bar{I}_{5,6,8,9}$ Kruger:1999xa . Here $I_{i}$ is equal to $\bar{I}_{i}$ with all weak $CP$ phases conjugated.

Appendix B Detailed fit result distributions

In Section 3 we studied sensitivity to various signatures of BSM physics. Here we provide the detailed results for the extracted Wilson coefficients with the maximum likelihood fit to MC data sets produced by the generator in SM and various BSM scenarios. A numerical summary of the fit results is shown in Table 1.

Mode	Electron		Muon
Wilson coefficient	Mean	RMS	Mean	RMS
$\Re e~{}\delta C_{9}$	$-0.8642$	$0.0872$	$-0.8623$	$0.1088$
$\Im m~{}\delta C_{9}$	$-0.0389$	$0.3000$	$-0.0446$	$0.3184$
$\Re e~{}\delta C_{10}$	$-$	$-$	$0.0032$	$0.2147$
$\Im m~{}\delta C_{10}$	$-$	$-$	$-0.0556$	$0.2786$
$\Re e~{}\delta C_{9}^{\prime}$	$-$	$-$	$0.0045$	$0.1039$
$\Im m~{}\delta C_{9}^{\prime}$	$-$	$-$	$-0.0003$	$0.1565$
$\Re e~{}\delta C_{10}^{\prime}$	$-$	$-$	$0.0033$	$0.0947$
$\Im m~{}\delta C_{10}^{\prime}$	$-$	$-$	$-0.0037$	$0.1394$
$\Re e~{}\delta C_{7}$	$0.0006$	$0.0040$	$0.0010$	$0.0074$
$\Im m~{}\delta C_{7}$	$0.0001$	$0.0020$	$0.0014$	$0.0230$
$\Re e~{}\delta C_{7}^{\prime}$	$0.0006$	$0.0088$	$0.0000$	$0.0168$
$\Im m~{}\delta C_{7}^{\prime}$	$0.0000$	$0.0091$	$0.0003$	$0.0186$

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

\delta C_{9}=-0.87

. Other values are extracted assuming the SM.

References

(1) D.J. Lange, The EvtGen particle decay simulation package, Nucl. Instrum. Meth. A 462 (2001) 152.
(2) W. Altmannshofer and P. Stangl, New physics in rare B decays after Moriond 2021, Eur. Phys. J. C 81 (2021) 952 [2103.13370].
(3) M. Algueró, B. Capdevila, S. Descotes-Genon, J. Matias and M. Novoa-Brunet, $b\to s\ell\ell$ Global Fits after $R_{K_{S}}$ and $R_{K^{*+}}$ , 4, 2021 [2104.08921].
(4) T. Hurth, F. Mahmoudi, D.M. Santos and S. Neshatpour, More indications for lepton nonuniversality in $b\to s\ell^{+}\ell^{-}$ , Phys. Lett. B 824 (2022) 136838 [2104.10058].
(5) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini and M. Valli, New Physics without bias: Charming Penguins and Lepton Universality Violation in $b\to s\ell^{+}\ell^{-}$ decays, 2110.10126.
(6) M. Chrzaszcz, A. Mauri, N. Serra, R. Silva Coutinho and D. van Dyk, Prospects for disentangling long- and short-distance effects in the decays $B\to K^{*}\mu^{+}\mu^{-}$ , JHEP 10 (2019) 236 [1805.06378].
(7) B. Capdevila, S. Descotes-Genon, J. Matias and J. Virto, Assessing lepton-flavour non-universality from $B\to K^{*}\ell\ell$ angular analyses, JHEP 10 (2016) 075 [1605.03156].
(8) Belle collaboration, Lepton-Flavor-Dependent Angular Analysis of $B\to K^{\ast}\ell^{+}\ell^{-}$ , Phys. Rev. Lett. 118 (2017) 111801 [1612.05014].
(9) Belle-II collaboration, The Belle II Physics Book, PTEP 2019 (2019) 123C01 [1808.10567].
(10) A. Mauri, N. Serra and R. Silva Coutinho, Towards establishing lepton flavor universality violation in $\bar{B}\to\bar{K}^{*}\ell^{+}\ell^{-}$ decays, Phys. Rev. D 99 (2019) 013007 [1805.06401].
(11) LHCb collaboration, Test of lepton universality in beauty-quark decays, 2103.11769.
(12) LHCb collaboration, Test of lepton universality with $B^{0}\rightarrow K^{*0}\ell^{+}\ell^{-}$ decays, JHEP 08 (2017) 055 [1705.05802].
(13) LHCb collaboration, Test of lepton universality in $b\rightarrow s\ell^{+}\ell^{-}$ decays, Phys. Rev. Lett. 131 (2023) 051803 [2212.09152].
(14) LHCb collaboration, Measurement of $CP$ -Averaged Observables in the $B^{0}\rightarrow K^{*0}\mu^{+}\mu^{-}$ Decay, Phys. Rev. Lett. 125 (2020) 011802 [2003.04831].
(15) LHCb collaboration, Branching Fraction Measurements of the Rare $B^{0}_{s}\rightarrow\phi\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow f_{2}^{\prime}(1525)\mu^{+}\mu^{-}$ - Decays, Phys. Rev. Lett. 127 (2021) 151801 [2105.14007].
(16) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini and M. Valli, Constraints on Lepton Universality Violation from Rare $B$ Decays, 2212.10516.
(17) R. Mandal, R. Sinha and D. Das, Testing New Physics Effects in $B\to K^{*}\ell^{+}\ell^{-}$ , Phys. Rev. D 90 (2014) 096006 [1409.3088].
(18) W. Altmannshofer, P. Ball, A. Bharucha, A.J. Buras, D.M. Straub and M. Wick, Symmetries and Asymmetries of $B\to K^{*}\mu^{+}\mu^{-}$ Decays in the Standard Model and Beyond, JHEP 01 (2009) 019 [0811.1214].
(19) A. Bharucha, D.M. Straub and R. Zwicky, $B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules, JHEP 08 (2016) 098 [1503.05534].
(20) P. Ball and R. Zwicky, $B_{d,s}\to\rho,\omega,K^{*},\phi$ decay form-factors from light-cone sum rules revisited, Phys. Rev. D 71 (2005) 014029 [hep-ph/0412079].
(21) R.R. Horgan, Z. Liu, S. Meinel and M. Wingate, Rare $B$ decays using lattice QCD form factors, PoS LATTICE2014 (2015) 372 [1501.00367].
(22) Belle-II Framework Software Group collaboration, The Belle II Core Software, Comput. Softw. Big Sci. 3 (2019) 1 [1809.04299].
(23) Belle collaboration, Measurements of branching fractions and decay amplitudes in $B\to J/\psi K^{*}$ decays, Phys. Lett. B 538 (2002) 11 [hep-ex/0205021].
(24) BaBar collaboration, Measurement of decay amplitudes of $B\to J/\psi K^{*},\psi(2S)K^{*}$ , and $\chi_{c1}K^{*}$ with an angular analysis, Phys. Rev. D 76 (2007) 031102 [0704.0522].
(25) LHCb collaboration, Measurement of the polarization amplitudes in $B^{0}\to J/\psi K^{*}(892)^{0}$ decays, Phys. Rev. D 88 (2013) 052002 [1307.2782].
(26) M. Beneke, T. Feldmann and D. Seidel, Systematic approach to exclusive $B\to V\ell^{+}\ell^{-}$ , $V\gamma$ decays, Nucl. Phys. B 612 (2001) 25 [hep-ph/0106067].
(27) A. Khodjamirian, T. Mannel, A.A. Pivovarov and Y.M. Wang, Charm-loop effect in $B\to K^{(*)}\ell^{+}\ell^{-}$ and $B\to K^{*}\gamma$ , JHEP 09 (2010) 089 [1006.4945].
(28) A. Sibidanov, T.E. Browder, S. Dubey, S. Kohani, R. Mandal, S. Sandilya et al., A New Tool for Detecting BSM Physics in $B\to K^{*}\ell\ell$ Decays, in 2022 Snowmass Summer Study, 3, 2022 [2203.06827v2].
(29) F. Kruger and L.M. Sehgal, Lepton polarization in the decays $B\to X_{s}\mu^{+}\mu^{-}$ and $B\to X_{s}\tau^{+}\tau^{-}$ , Phys. Lett. B 380 (1996) 199 [hep-ph/9603237].
(30) F. Kruger, L.M. Sehgal, N. Sinha and R. Sinha, Angular distribution and CP asymmetries in the decays $\bar{B}\to K^{-}\pi^{+}e^{-}e^{+}$ and $\bar{B}\to\pi^{-}\pi^{+}e^{-}e^{+}$ , Phys. Rev. D 61 (2000) 114028 [hep-ph/9907386].

A new Monte Carlo Generator for BSM physics in B→K∗​ℓ+​ℓ−B\to K^{*}\ell^{+}\ell^{-} decays with an application to lepton non-universality in angular distributions