Exhaustive Model Selection in $b \to s \ell \ell$ Decays: Pitting Cross-Validation against AIC$_c$

TL;DR

This paper compares model selection methods like cross-validation and AICc to identify the best effective operators explaining recent $b o s \, \ell^+\ell^-$ decay data, highlighting the importance of angular observables and specific operator scenarios.

Contribution

It introduces a comprehensive comparison of model selection tools in the context of new physics in $b \to s \ell \ell$ decays, identifying the most favored operators and their parameters.

Findings

A left-handed quark current with vector muon coupling is the most favored scenario.

Angular observables are crucial in the model selection process.

Certain two-operator scenarios can also explain the $R_{K^{(*)}}$ ratios.

Abstract

In the light of recent data, we study the new physics effects in the exclusive $b \to s \ell^+\ell^-$ decays from a model independent perspective. Different combinations of the dimension six effective operators along with their respective Wilson coefficients are chosen for the analysis. To find out the operator or sets of operators that can best explain the available data in this channel, we simultaneously apply popular model selection tools like cross-validation and the information theoretic approach like Akaike Information Criterion (AIC). There are one, two, and three-operator scenarios which survive the test and a left-handed quark current with vector muon coupling is common among them. This is also the only surviving one-operator scenario. Best-fit values and correlations of the new Wilson coefficients are supplied for all the selected scenarios. We find that the angular observables play the dominant role in the model selection procedure. We also note that while a left-handed quark current with axial-vector muon coupling is the only one-operator scenario able to explain the ratios $R_{K^{(*)}}$ ($R_{K^*}$ for $q^2\in [ 0.045, 1.1] {\rm GeV}^2$ in particular), there are also a couple of two operator scenarios that can simultaneously explain the measured $R_{K^{(*)}}$.

Tables

Table 1. Table 1: SM values and experimental world averages of inclusive and exclusive b → s ​ γ → 𝑏 𝑠 𝛾 b\to s\gamma observables used in our analysis.

Observable	SM prediction	Measurement
	$\times {10}^5$	$\times {10}^5$
$\text{BR}{\left(B\to X_s\gamma \right)}_{E_{\gamma }>1.6\text{GeV}}$	$33.6\pm 2.6$ Misiak et al. (2015)	$32.7\pm 1.4$ Misiak and Steinhauser (2017)
$\text{BR}\left(B^{+}\to K^{\ast }\gamma \right)$	$3.51\pm 0.78$	$4.21\pm 0.18$ Amhis et al. (2014)
$\text{BR}\left(B^0\to K^{\ast }\gamma \right)$	$3.49\pm 0.78$	$4.33\pm 0.15$ Amhis et al. (2014)
$\overline{\text{BR}}\left(B_s\to \varphi \gamma \right)$	$4.33\pm 0.77$	$3.5\pm 0.4$Aaij et al. (2013); Dutta et al. (2015)

Table 2. Table 2: Fit-qualities, model selection criteria, parameter estimates and effects on radiative decays for the ‘best’ selected models with the ‘New’ data-set, with the ‘Moments’ estimate of the angular observables. Selected models are obtained from fig. 3a . Last four columns showcase the deviations (in units of ‘ σ 𝜎 \sigma ’) between the experimental value of the radiative decays and the corresponding value obtained with the fit results.

Model	${\chi }_{\text{Min}}^2/$	p-val	${\omega }^{{\Delta \text{AIC}}_c}$	${\text{MSE}}_{X-\text{val}}$	Parameter	Deviation in $\sigma$
Index	DOF	($\%$)	($\%$)		Values	$\text{B}\to X_s\gamma$	$B^{+}\to K^{\ast }\gamma$	${\Delta \text{B}}^0\to K^{\ast }\gamma$	${\Delta \text{B}}_s\to \varphi \gamma$
$18$	250.28/256	$58.9$	$3.4$	$0.917$	$\begin{array}{c}\Delta C_9\to \text{-1.13}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.25}\pm \text{0.17}\end{array}$	$-$	$-$	$-$	$-$
$2$	252.44/257	$56.9$	$3.2$	$0.933$	$\begin{array}{c}\Delta C_9\to \text{-1.12}\pm \text{0.13}\end{array}$	$-$	$-$	$-$	$-$
$76$	249.12/255	$59.2$	$2.2$	$0.918$	$\begin{array}{c}\Delta C_9\to \text{-1.18}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.34}\pm \text{0.19}\\ C_S\to \text{-0.035}\pm \text{0.016}\end{array}$	$-$	$-$	$-$	$-$
$77$	249.16/255	$59.1$	$2.1$	$0.918$	$\begin{array}{c}\Delta C_9\to \text{-1.18}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.34}\pm \text{0.19}\\ C_S^{\prime }\to \text{0.035}\pm \text{0.016}\end{array}$	$-$	$-$	$-$	$-$
$20$	251.52/256	$56.7$	$1.8$	$0.928$	$\begin{array}{c}\Delta C_9\to \text{-1.15}\pm \text{0.14}\\ C_{10}^{\prime }\to \text{-0.1}\pm \text{0.104}\end{array}$	$-$	$-$	$-$	$-$
$10$	251.97/256	$55.9$	$1.4$	$0.932$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.01}\pm \text{0.015}\\ \Delta C_9\to \text{-1.15}\pm \text{0.14}\end{array}$	$0.31$	$-0.87$	$-1.06$	$1.22$
$46$	250.14/255	$57.4$	$1.3$	$0.922$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.0058}\pm \text{0.0155}\\ \Delta C_9\to \text{-1.15}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.24}\pm \text{0.18}\end{array}$	$0.3$	$-0.87$	$-1.06$	$1.22$
$74$	250.16/255	$57.4$	$1.3$	$0.925$	$\begin{array}{c}\Delta C_9\to \text{-1.16}\pm \text{0.15}\\ C_9^{\prime }\to \text{0.26}\pm \text{0.17}\\ \Delta C_{10}\to \text{-0.041}\pm \text{0.118}\end{array}$	$-$	$-$	$-$	$-$
$75$	250.21/255	$57.3$	$1.2$	$0.923$	$\begin{array}{c}\Delta C_9\to \text{-1.12}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.3}\pm \text{0.26}\\ C_{10}^{\prime }\to \text{0.04}\pm \text{0.157}\end{array}$	$-$	$-$	$-$	$-$

Table 3. Table 3: Same as the table 2 , but with the ‘Likelihood’ estimate of the angular observables. Selected models are obtained from fig. 3b .

Model	${\chi }_{\text{Min}}^2/$	p-val	${\omega }^{{\Delta \text{AIC}}_c}$	${\text{MSE}}_{X-\text{val}}$	Parameter	Deviation in $\sigma$
Index	DOF	($\%$)	($\%$)		Values	$\text{B}\to X_s\gamma$	$B^{+}\to K^{\ast }\gamma$	${\Delta \text{B}}^0\to K^{\ast }\gamma$	${\Delta \text{B}}_s\to \varphi \gamma$
$132$	217.02/210	$35.5$	$3.7$	$0.985$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.04}\pm \text{0.015}\\ \Delta C_9\to \text{-1.39}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.45}\pm \text{0.2}\\ C_S\to \text{-0.042}\pm \text{0.013}\end{array}$	$0.44$	$-0.82$	$-1.02$	$1.27$
$133$	217.07/210	$35.4$	$3.6$	$0.986$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.04}\pm \text{0.015}\\ \Delta C_9\to \text{-1.39}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.45}\pm \text{0.2}\\ C_S^{\prime }\to \text{0.042}\pm \text{0.013}\end{array}$	$0.44$	$-0.82$	$-1.02$	$1.27$
$130$	217.58/210	$34.5$	$2.8$	$0.976$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.044}\pm \text{0.015}\\ \Delta C_9\to \text{-1.42}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.32}\pm \text{0.19}\\ \Delta C_{10}\to \text{-0.16}\pm \text{0.11}\end{array}$	$0.47$	$-0.81$	$-1.01$	$1.28$
$46$	219.66/211	$32.7$	$2.8$	$0.988$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.04}\pm \text{0.015}\\ \Delta C_9\to \text{-1.34}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.33}\pm \text{0.2}\end{array}$	$0.44$	$-0.82$	$-1.02$	$1.27$
$47$	220.36/211	$31.5$	$2.$	$0.995$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.048}\pm \text{0.015}\\ \Delta C_9\to \text{-1.43}\pm \text{0.15}\\ \Delta C_{10}\to \text{-0.16}\pm \text{0.11}\end{array}$	$0.5$	$-0.8$	$-1.$	$1.29$
$10$	222.46/212	$29.7$	$1.9$	$1.001$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.043}\pm \text{0.015}\\ \Delta C_9\to \text{-1.33}\pm \text{0.13}\end{array}$	$0.46$	$-0.81$	$-1.01$	$1.28$
$257$	216.47/209	$34.7$	$1.7$	$0.99$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.042}\pm \text{0.015}\\ \Delta C_9\to \text{-1.42}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.41}\pm \text{0.21}\\ \Delta C_{10}\to \text{-0.091}\pm \text{0.123}\\ C_S\to \text{-0.036}\pm \text{0.017}\end{array}$	$0.46$	$-0.82$	$-1.01$	$1.28$
$258$	216.51/209	$34.6$	$1.7$	$0.99$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.042}\pm \text{0.015}\\ \Delta C_9\to \text{-1.42}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.41}\pm \text{0.21}\\ \Delta C_{10}\to \text{-0.092}\pm \text{0.123}\\ C_S^{\prime }\to \text{0.036}\pm \text{0.017}\end{array}$	$0.46$	$-0.82$	$-1.01$	$1.28$
$131$	218.98/210	$32.1$	$1.4$	$0.993$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.04}\pm \text{0.015}\\ \Delta C_9\to \text{-1.32}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.47}\pm \text{0.25}\\ C_{10}^{\prime }\to \text{0.11}\pm \text{0.14}\end{array}$	$0.44$	$-0.82$	$-1.02$	$1.27$
$265$	216.92/209	$33.9$	$1.4$	$0.986$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.04}\pm \text{0.015}\\ \Delta C_9\to \text{-1.39}\pm \text{0.13}\\ C_9^{\prime }\to \text{0.44}\pm \text{0.2}\\ C_S\to \text{-0.24}\pm \text{0.6}\\ C_S^{\prime }\to \text{-0.19}\pm \text{0.6}\end{array}$	$0.44$	$-0.82$	$-1.02$	$1.27$

Table 4. Table 4: Fit-qualities, model selection criteria, parameter estimates and effects on radiative decays for the ‘best’ selected models with the ‘Old’ data-set, with the ‘Moments’ estimate of the angular observables. Selected models are obtained from fig. 4a . Last four columns showcase the deviations (in units of ‘ σ 𝜎 \sigma ’) between the experimental value of the radiative decays and the corresponding value obtained with the fit results.

Model	${\chi }_{\text{Min}}^2/$	p-val	${\omega }^{{\Delta \text{AIC}}_c}$	${\text{MSE}}_{X-\text{val}}$	Parameter	Deviation in $\sigma$
Index	DOF	($\%$)	($\%$)		Values	$\text{B}\to X_s\gamma$	$B^{+}\to K^{\ast }\gamma$	${\Delta \text{B}}^0\to K^{\ast }\gamma$	${\Delta \text{B}}_s\to \varphi \gamma$
$2$	245.67/254	$63.5$	$5.$	$0.918$	$\begin{array}{c}\Delta C_9\to \text{-1.26}\pm \text{0.14}\end{array}$	$-$	$-$	$-$	$-$
$10$	244.92/253	$63.1$	$2.6$	$0.916$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.013}\pm \text{0.015}\\ \Delta C_9\to \text{-1.3}\pm \text{0.15}\end{array}$	$0.32$	$-0.87$	$-1.06$	$1.22$
$19$	245.42/253	$62.2$	$2.$	$0.926$	$\begin{array}{c}\Delta C_9\to \text{-1.22}\pm \text{0.16}\\ \Delta C_{10}\to \text{0.061}\pm \text{0.123}\end{array}$	$-$	$-$	$-$	$-$
$21$	245.48/253	$62.1$	$2.$	$0.923$	$\begin{array}{c}\Delta C_9\to \text{-1.27}\pm \text{0.15}\\ C_S\to \text{-0.021}\pm \text{0.026}\end{array}$	$-$	$-$	$-$	$-$
$22$	245.51/253	$62.$	$1.9$	$0.923$	$\begin{array}{c}\Delta C_9\to \text{-1.27}\pm \text{0.15}\\ C_S^{\prime }\to \text{0.02}\pm \text{0.026}\end{array}$	$-$	$-$	$-$	$-$
$18$	245.55/253	$62.$	$1.9$	$0.915$	$\begin{array}{c}\Delta C_9\to \text{-1.25}\pm \text{0.14}\\ C_9^{\prime }\to \text{0.067}\pm \text{0.195}\end{array}$	$-$	$-$	$-$	$-$
$20$	245.59/253	$61.9$	$1.9$	$0.92$	$\begin{array}{c}\Delta C_9\to \text{-1.26}\pm \text{0.14}\\ C_{10}^{\prime }\to \text{-0.03}\pm \text{0.109}\end{array}$	$-$	$-$	$-$	$-$

Table 5. Table 5: Same as table 4 , but with the ‘Likelihood’ estimate of the angular observables. Selected models are obtained from fig. 4b .

Model	${\chi }_{\text{Min}}^2/$	p-val	${\omega }^{{\Delta \text{AIC}}_c}$	${\text{MSE}}_{X-\text{val}}$	Parameter	Deviation in $\sigma$
Index	DOF	($\%$)	($\%$)		Values	$\text{B}\to X_s\gamma$	$B^{+}\to K^{\ast }\gamma$	${\Delta \text{B}}^0\to K^{\ast }\gamma$	${\Delta \text{B}}_s\to \varphi \gamma$
$10$	213.78/209	$39.6$	$5.3$	$0.973$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.028}\pm \text{0.015}\\ \Delta C_9\to \text{-1.37}\pm \text{0.14}\end{array}$	$0.37$	$-0.85$	$-1.04$	$1.24$
$2$	217.19/210	$35.2$	$2.7$	$0.989$	$\begin{array}{c}\Delta C_9\to \text{-1.28}\pm \text{0.13}\end{array}$	$-$	$-$	$-$	$-$
$49$	213.2/208	$38.8$	$2.5$	$0.981$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.029}\pm \text{0.015}\\ \Delta C_9\to \text{-1.4}\pm \text{0.14}\\ C_S\to \text{-0.028}\pm \text{0.019}\end{array}$	$0.38$	$-0.85$	$-1.04$	$1.25$
$50$	213.23/208	$38.7$	$2.5$	$0.981$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.029}\pm \text{0.015}\\ \Delta C_9\to \text{-1.4}\pm \text{0.14}\\ C_S^{\prime }\to \text{0.028}\pm \text{0.019}\end{array}$	$0.38$	$-0.85$	$-1.04$	$1.25$
$47$	213.65/208	$37.9$	$2.$	$0.976$	$\begin{array}{c}{C}_7^{\prime }\to \text{0.029}\pm \text{0.015}\\ \Delta C_9\to \text{-1.39}\pm \text{0.15}\\ \Delta C_{10}\to \text{-0.042}\pm \text{0.117}\end{array}$	$0.38$	$-0.85$	$-1.04$	$1.25$

Equations

$$R_H=\frac{{\int }_{q_{min}^2}^{q_{max}^2}\frac{d\Gamma (B\to H{\mu }^{+}{\mu }^{-})}{d{q}^2}}{{\int }_{q_{min}^2}^{q_{max}^2}\frac{d\Gamma (B\to H{e}^{+}{e}^{-})}{d{q}^2}}$$

$$R_{K^{\ast }}=\{\begin{array}{l}0.920\pm 0.007,\text{q2∈[0.045,1.1] GeV2},\\ 0.996\pm 0.002,\text{q2∈[1.1,6] GeV2}.\end{array}$$

$$R_K=0.846_{-0.054-0.014}^{+0.060+0.016},\text{q2∈[1.1,6] GeV2},$$

$$R_{K^{\ast }}=\{\begin{array}{l}0.660_{-0.070}^{+0.110}\pm 0.024,\text{q2∈[0.045,1.1] GeV2},\\ 0.685_{-0.069}^{+0.113}\pm 0.047,\text{q2∈[1.1,6] GeV2}.\end{array}$$

$$R_{K^{\ast }}=\{\begin{array}{l}0.52_{-0.26}^{+0.36}\pm 0.05,\text{q2∈[0.045,1.1] GeV2},\\ 0.96_{-0.29}^{+0.45}\pm 0.11,\text{q2∈[1.1,6] GeV2}.\end{array}$$

$$R_K=0.98_{-0.23}^{+0.27}\pm 0.06,\text{q2∈[1,6] GeV2},$$

$$Br(B_s\to \mu \mu )=\left(3.1\pm 0.6\right)\times 10^{-9}.$$

$$f_{B_s}=0.2284\pm 0.0037\mathrm{GeV}.$$

$${\mathcal{H}}_{eff}=-\frac{4G_F}{\sqrt{2}}\left({\lambda }_t{\mathcal{H}}_{eff}^{(t)}+{\lambda }_u{\mathcal{H}}_{eff}^{(u)}\right)$$

$${\mathcal{H}}_{eff}^{(t)}$$

$$\sum _{i=7,8,9,10,P,S}(C_i{\mathcal{O}}_i+C_i^{\prime }{\mathcal{O}}_i^{\prime }),$$

$${\mathcal{H}}_{eff}^{(u)}$$

$${\mathcal{O}}_7$$

$${\mathcal{O}}_9$$

$${\mathcal{O}}_{10}$$

$${\mathcal{O}}_S$$

$${\mathcal{O}}_P$$

$${\mathrm{AIC}}_c={\chi }_{min}^2+2K+\frac{2K(K+1)}{n-K-1}$$

$$w_i^{\Delta {\text{AIC}}_c}=\frac{e^{(-{\Delta }_i^{AIC}/2)}}{{\sum }_{r=1}^R{e}^{(-{\Delta }_r^{AIC}/2)}},$$

Full text

I Introduction

Decays involving $b\to s\ell\ell$ transitions are suppressed in the standard model (SM). These decay modes are potentially sensitive to new physics effects. Whether the contributions appear at the tree or the loop level depends on the type of the new physics (NP). A lot of attention, both experimental and theoretical, have been given to $B\to K^{(*)}\mu^{+}\mu^{-}$ decays in the last couple of years. There are several angular observables associated with these decays, which are potentially sensitive to the NP effects and are measured by LHCb Aaij et al. (2016a); Aaboud et al. (2018). A couple of them have shown discrepancies with their respective SM predictions Descotes-Genon et al. (2013a, b); Horgan et al. (2014); Bharucha et al. (2016). However, these angular observables are not free from hadronic uncertainties and it is fairly possible that the observed discrepancies are due to poorly known hadronic effects, e.g., see Ciuchini et al. (2016) for details. Furthermore, these modes offer theoretically clean observables like

[TABLE]

where $H$ is either $K$ or $K^{*}$ meson and $q^{2}$ is the dilepton squared mass. These ratios are useful to test lepton flavor universality violation (LFUV) and within appropriately chosen ranges of $q^{2}$, these observables can be predicted very precisely in the SM; see Hiller and Kruger (2004); Bordone et al. (2016) for details. The SM predictions are, respectively, $R(K)=1.0004(8)$, and

[TABLE]

The LHCb collaboration has measured Aaij et al. (2017, 2019)

[TABLE]

and

[TABLE]

We will use the notation $R^{Low}_{K^{*}}$ and $R^{Central}_{K^{*}}$ from now on to represent $R_{K^{*}}$ for values of $q^{2}$ in $[0.045,1.1]{\rm GeV}^{2}$ and $[1.1,6]{\rm GeV}^{2}$, respectively. Very recently Belle has measured the observables $R_{K^{(*)}}$ and the measured values are given by Abdesselam et al. (2019a)

[TABLE]

and Abdesselam et al. (2019b)

[TABLE]

These new measurements have larger uncertainties compared to those from LHCb, but the results are consistent with each other. Belle has also measured separate ratios like $R_{K^{*}}^{0}$ and $R_{K^{*}}^{+}$, but the associated uncertainties are quite large at the moment. On the whole, the deviation between data and SM predictions stand at the level of 2.5 to 3 $\sigma$. Future measurements of these ratios with enough statistical significance would have the potential to discover NP unambiguously.

The observed discrepancies can be explained in various NP models. Different types of new physics interactions (like vector, scalar etc.) with different Lorentz structures may contribute to these decays and explain the data. A lot of work has already been done and it is a difficult task to quote all of them. We are more interested, in the present work, in a model independent analysis. There are a few related analyses available in the literature, which mainly focus on considering one or two operators at a time Aebischer et al. (2019); Alok et al. (2019); Capdevila et al. (2018); Arbey et al. (2018); Descotes-Genon et al. (2013c); Altmannshofer and Straub (2013); Beaujean et al. (2014); Hurth and Mahmoudi (2014); Altmannshofer and Straub (2015); Altmannshofer et al. (2009); Egede et al. (2010); Das and Sinha (2012); Hurth et al. (2014); Ciuchini et al. (2019); Bobeth et al. (2012); Geng et al. (2017); Ciuchini et al. (2017); Hiller and Kruger (2004); Kowalska et al. (2019).

In this article, we have done a model independent analysis of the NP affecting the $b\to s\ell^{+}\ell^{-}$ decay modes. The operator basis is exactly the same as that given in Ref. Altmannshofer et al. (2009). We have considered all possible combinations of these operators and categorized them as independent ‘model’s (scenarios). There are several models capable of describing the observed data and one is thus confronted with the problem of model selection.

In short, the problem of model selection is as follows: any model, used to represent certain observation, will almost never be exact; chances are, that some information will be lost due to the choice of that particular model Box (1976); Box et al. (2005). Choosing a simplistic model with too few parameters can involve making unrealistically simple assumptions and lead to high bias, poor prediction, and consequent missed opportunities for insight. While simplistic models are not flexible enough to describe the sample or the population well, a model with a larger number of parameters can fit the observed data very well. Does this make it a better model? With too many parameters, we face the possibility of just fitting the noise in the data and losing sight of the important trends. The most general problem in model selection is thus the optimization of the parameters required to explain certain observation Burnham and Anderson (2002); Geisser and Eddy (1979). The motivation: there must be a happy medium somewhere.

The most generally applicable, powerful, and reliable method for model comparison (also computationally expensive) is ‘cross-validation’Andrae et al. (2010), which, in addition to testing the predictive power of the model, minimizes the bias and variance together by minimizing the mean-squared-error (MSE). On the other hand, penalized-likelihood information criteria, such as the Akaike Information Criterion (AIC) Akaike (1998), and the Bayesian Information Criterion (BIC; more aptly named as Schwarz Information Criterion) are widely used for model selection. AIC estimates the relative amount of information lost by a given model: the less the information lost by a model, the higher the quality of that model. For a detailed discussion on various model selection procedures and their relative performance with respect to cross-validation, ref. arlot2010. In our earlier publications Bhattacharya et al. (2017, 2019), we have used these criteria in the context of NP model selections in $b\to c\tau\nu_{\tau}$ decays.

Very recently, in a Bayesian analysis of $b\to s\ell^{+}\ell^{-}$ decays Ciuchini et al. (2019), an information criterion has been used. They have shown the use of a criterion closely related to DIC (Deviance Information Criterion; the definition by Gelman et al Gelman et al. (2013)) and BPIC (Bayesian Predictive Information Criterion Ando (2007)) for model selection, which is not only ideal for samples from a Markov Chain Monte Carlo but is also asymptotically equivalent Claeskens and Hjort (2008) to natural model-robust version of AIC.

AIC is easy to calculate in a frequentist analysis, which is not the case for Bayesian analyses. The main difference between that analysis and ours is that they created the model hierarchy by defining $\Delta IC=IC_{SM}-IC_{NP}$. As a result, the quality of a model is determined with respect to the SM, whereas in our case the best model is picked up first and the hierarchy is defined with respect to that. Still, similar to our findings, they have found that the case $\Delta C_{9}$ ($C^{NP}_{9,~{}\mu}$ in their paper) provides the optimal outcome for $B\to K^{*}\ell\ell$ transitions if we consider only the ‘Moments’ data for the angular observables, in addition to the new LFUV data.

In the present analysis, we use both AIC and cross-validation to pin down the best possible model(s), and find out how one can use both procedures in tandem to glean the most out of the data at hand.

The article is organized as follows: in section II, we discuss the present experimental and theoretical status of the observables used in this analysis. Section IV discusses the detailed methodology of the statistical analysis, as well as model selection. We present our results in section V and in section VI we summarize.

II Backdrop

II.1 Experimental

We list the experimental results used in this analysis and the corresponding references below:

• •

Binned data on the angular observables related to the $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ decays have been obtained from refs. Aaij et al. (2016a) (LHCb) and Aaboud et al. (2018) (ATLAS)111We refrain from using the very old (2012) CDF data available from the public note Collaboration et al. (2012) on the angular observables in the $B\to K^{(*)}\mu^{+}\mu^{-}$ decays..

• •

Binned data for the differential branching fraction for $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ have been obtained from ref. Aaij et al. (2016b) and that for $B^{+}\to K^{*+}\mu^{+}\mu^{-}$ from Aaij et al. (2014a). Both of these are LHCb references.

• •

Binned data on the angular observables for $B^{+}\to K^{+}\mu^{+}\mu^{-}$ ($A_{FB}$ and $F_{H}$) have been taken from ref. Sirunyan et al. (2018) (CMS).

• •

Binned data on the differential branching fraction for $B^{+}\to K^{+}\mu^{+}\mu^{-}$ and $B^{0}\to K^{0}\mu^{+}\mu^{-}$ reported by LHCb have been taken from ref. Aaij et al. (2014a).

• •

Binned data on the angular observables for $B_{s}\to\phi\mu^{+}\mu^{-}$ (LHCb) have been taken from ref. Aaij et al. (2015).

• •

The lepton flavor universality violating (LFUV) observables $R_{K^{*}}$, both for the low and central bin, have been obtained from ref. Aaij et al. (2017) (LHCb). We also include the recent measurements on these observables (for the same bins) from Belle Abdesselam et al. (2019a). The old $R_{K}$ data from LHCb has been taken from ref. Aaij et al. (2014b). The updated result on the same has also been included Aaij et al. (2019).

• •

The experimental result for the branching ratio (BR) corresponding to $B_{s}\to\mu^{+}\mu^{-}$ has been taken from hfl which is the average of the measured values by CMS, ATLAS and LHCb. The value is given by

[TABLE]

The decay constant is taken from ref. Charles et al. (2014); Aoki et al. (2019)

[TABLE]

All numerical uncertainties quoted in this analysis, unless otherwise specified, denotes the 1-$\sigma$ (68% c.l.) range.

A few words regarding the data on the angular observables due to LHCb taken from ref. Aaij et al. (2016a) is in order at this point. LHCb has provided the data corresponding to the angular observables in bins of $q^{2}$ ($q=p_{\mu^{+}}+p_{\mu^{-}}$, $p_{\mu}$ being the four-momentum of muon) by performing two separate analyses. The more commonly used dataset in the community is that due to the “Method of Moments”. The angular observables in this case are determined by using a principal moment analysis of the angular distribution without carrying out any angular fit to the data Beaujean et al. (2015); Gratrex et al. (2016). These moments are continuous functions of $q^{2}$. The statistical uncertainties for these angular moments are estimated using a bootstrapping technique Efron (1979) and confidence intervals are defined such that they include the $16^{th}$–$84^{th}$ percentiles of the bootstrap distribution of the observables. The other method termed the “Maximum likelihood fit” involves an unbinned maximum likelihood fit to the invariant mass $m(\mu^{+}\mu^{-}(K^{*}\to)K^{+}\pi^{-})$ and the three decay angles $\cos\theta_{l}$, $\cos\theta_{K}$ and $\phi$ in each $q^{2}$ bin, where:

• •

$\theta_{l}$ is the angle between the $\mu^{+}$ ($\mu^{-}$) and the direction opposite to that of the $B^{0}$ ($\bar{B}^{0}$) in the rest frame of the dimuon system,

• •

$\theta_{K}$ is the angle between the direction of the $K^{+}$ ($K^{-}$) and the $B^{0}$ ($\bar{B}^{0}$) in the rest frame of the $K^{*0}$ ($\bar{K}^{*0}$) where the $K^{*0}$ meson is reconstructed through the decay $K^{*0}\to K^{+}\pi^{-}$, and

• •

$\phi$ is the angle between the plane defined by the dimuon pair and the plane defined by the $K$ and the $\pi$ in the $B^{0}$ ($\bar{B}^{0}$) rest frame.

The bin sizes corresponding to the maximum likelihood analysis are larger than those for the method of moments. This is done since there is a dearth of statistics, and an increase in the bin-size renders the precision comparable with that for the method of moments. With increased number of events in the future, an unbinned likelihood analysis will become the norm, but at the present precision level, the ‘Moments’ data is at least equally dependable, if not more. To examine and point to any fundamental difference between these two data-sets in presence of NP models, we use both these sets as separate cases in our analysis. To the best of our knowledge, this is the first global $b\to s\ell\ell$ analysis that takes both of these datasets into account.

Apart from classifying the data according to whether it corresponds to the “Likelihood” or the “Moment” method for the angular observables, we have also prepared separate datasets which we call:

• •

“Old” dataset, containing the (previous) estimates for the LFUV $R_{K}$ and $R_{K^{*}}$ ratios from refs. Aaij et al. (2014b, 2017) respectively, and

• •

“New” dataset, where the old estimate for $R(K)$  Aaij et al. (2014b) by LHCb is replaced by the new one Aaij et al. (2019), while both the previous Aaij et al. (2017) and the current Abdesselam et al. (2019a) estimates for the $R_{K^{*}}$ ratio have been included. We also include the most recent measurement for $R_{K}$ due to Belle Abdesselam et al. (2019b).

We should mention here that we have only taken the “low bins” ($q^{2}\leq 6$ GeV2) corresponding to the experimental data referred to above. This is done so that we can avoid the region around the $J/\psi$ resonance (the “broad charmonium” region) since a trustworthy theoretical estimate for this region is challenging. Hence, we do not include the $R_{K^{*}}$ data from the low-recoil region provided by the recent Belle measurement from Abdesselam et al. (2019a). We take care of the systematic and statistical correlations separately in the data as and when they have been reported.

II.2 Theoretical

The effective Hamiltonian and the operator basis for exclusive $b\to s\mu^{+}\mu^{-}$ decays are taken from Altmannshofer et al. (2009); Altmannshofer and Straub (2015) and is written as:

[TABLE]

with the CKM combination $\lambda_{i}=V_{ib}V_{is}^{*}$ and

[TABLE]

We consider NP effects in the following operators:

[TABLE]

The NP contributions to operators $\mathcal{O}_{9,10}$ is given by $\Delta C_{9,10}$. In these decays, when the final state contains a vector meson, one can construct various helicity amplitudes. These helicity amplitudes are used to form angular coefficients which are relevant in defining the CP-symmetric and asymmetric observables measured by the different experimental collaborations. The details about various transversity amplitudes and the respective angular coefficients can be obtained from Altmannshofer et al. (2009). The two major components that go into the formation of the helicity amplitudes are the Wilson coefficients (WC) of different operators and the form factors which are defined as the hadronic matrix elements of various operators. We follow ref. Bharucha et al. (2016) for the form factors in $B\to K^{*}$ and $B_{s}\to\phi$ decays 222Although the latest LCSR estimates for the $B\to K^{*}$ decays is calculated in ref. Gubernari et al. (2019), it does not include the corresponding $B_{s}\to\phi$ matrix elements. Hence we refrain from using these results..

For the $B\to K$ sector we closely follow the methodology communicated in ref. Bobeth et al. (2007). This includes expressing the differential decay distribution in terms of a polynomial in $\cos\theta$, where $\theta$ denotes the angle between the direction of motion of the parent $B$ meson and the positively charged lepton in the dilepton center of mass frame. The coefficients of these terms can then be expressed as combinations of the corresponding WC and form factors. For the form factors, we use the results from ref. Altmannshofer and Straub (2015), where the authors perform a combined fit to the lattice computation in ref. Bouchard et al. (2013) as well as LCSR predictions at $q^{2}=0$ Ball and Zwicky (2005); Bartsch et al. (2009), using the parametrization and conventions of Bouchard et al. (2013). The method is described in details in the appendix of ref. Buras et al. (2015).

We also take care of the correlations among these form factor elements as reported in these references, in order to propagate them to form the theoretical correlations and errors for the corresponding observables.

Since our aim is to perform a global model selection based on the plethora of available $b\to s\ell\ell$ data discussed in sec. II.1, there is a possibility that amongst the selected models the operator with $C_{7}^{\prime}$ as coefficient may appear as a plausible solution. Such an operator is also relevant for the radiative decays like inclusive and exclusive $b\to s\gamma$. For such scenarios, we have checked whether parameter spaces which are allowed by $b\to sll$ data are also allowed by the inclusive $B\to X_{s}\gamma$ measurement, alongwith the branching ratios for the three exclusive radiative modes $B^{+}\to K^{*+}\gamma$, $B^{0}\to K^{*0}\gamma$ and the time integrated $\overline{\text{BR}}(B_{s}\to\phi\gamma)$ 333We refrain from using measurements for the CP asymmetries since our NP Wilson coefficients are taken to be real, thus excluding the possibility of CP violation in NP. The definitions and formulae for these modes are taken from ref. Paul and Straub (2017). We provide the experimental values and the SM estimates used in our analysis in Table. 1. The corresponding theoretical (for $B\to X_{s}\gamma$) and experimental references are provided therein. Our SM values are consistent with the estimates of ref. Paul and Straub (2017), within $1~{}\sigma$.

III NP and current data on $R_{K}$ and $R_{K^{*}}$

Before pursuing a detailed discussion on model selection, let us look for the NP effects in $b\to s\mu\mu$ decays, only in the light of recently updated measurements on $R_{K}$ and $R_{K^{*}}$ in this section, focusing on the measurements of $R_{K^{*}}$ in the low $q^{2}$ bins. There is some discrepancy between this particular data and the corresponding predicted value in the SM. However, one needs to remember that the angular observables are not free from hadronic uncertainties. In this part of the study, we do not consider any of the angular observables, neither do we carry out any fit to data. We simply check the dependencies of $R^{Low}_{K^{*}}$, $R^{Central}_{K^{*}}$, and $R_{K}$ on various WCs in one and two operator scenarios. We do not include $C_{S}$, $C_{P}$, $C^{\prime}_{S}$, and $C^{\prime}_{P}$, since the corresponding operators by themselves, or combinations including such operators, are unable to explain the observed data in $R^{Low}_{K^{*}}$. These WCs also suffer from tight constraints due to $B_{s}\to\mu\mu$ decays Alonso et al. (2014). Also, the new electromagnetic dipole operator ${\cal O}^{\prime}_{7}$ alone would not be able to explain the observed data on $R_{K}^{(*)}$ and the branching fractions in the above mentioned radiative decays simultaneously. Hence, we have not considered the effects of this operator in this part of the analysis.

The results of our analysis are presented in figures 1 and 2. In single operator scenarios, the correlations between all the above-mentioned observables are shown in figures 1a and 1b. It would be difficult to explain the observed data for $R^{Low}_{K^{*}}$ within their $1\sigma$ ranges will be difficult in the single-operator scenarios. Although the allowed region is tightly constrained, ${\cal O}_{10}$ (with WC $\Delta C_{10}$) is the only operator that can simultaneously explain all the data on $R_{K^{(*)}}$ except $R^{Low}_{K^{*}}$ from LHCb. The required value of $\Delta C_{10}$ lies in between 0.5 and 1.5, which is consistent with the measured value of $Br(B_{s}\to\mu\mu)$ within its 2-$\sigma$ range for detail see figure 2k. However, there are several candidates in the two operator scenarios that could explain all the data simultaneously. Among various possible combinations, the highly probable scenarios are the operators with the WCs $[\Delta C_{9},C^{\prime}_{10}]$, $[\Delta C_{10},C^{\prime}_{9}]$, $[\Delta C_{10},C^{\prime}_{10}]$, and $[\Delta C_{10},\Delta C_{9}]$. The other possible scenario $[\Delta C_{9},C^{\prime}_{9}]$ is less favored but allowed by the data. Also, the allowed values of $\Delta C_{10}$ and/or $C^{\prime}_{10}$ can explain $Br(B_{s}\to\mu\mu)$ within its 1-$\sigma$ range, wherever applicable; see figure 2k for details. In fig. 2k, we have shown the variations of $Br(B_{s}\to\mu\mu)$ w.r.t. the parameters $[\Delta C_{10},C^{\prime}_{10}]$. However, there are scenarios where only $\Delta C_{10}$ or $C^{\prime}_{10}$ appears. In such cases, depending on the scenario, one needs to look at the plot with either $\Delta C_{10}=0$ or $C^{\prime}_{10}=0$. To conclude this section, we would like to mention that the three or more operator scenarios could also be relevant to explain the present data on $R_{K^{(*)}}$ simultaneously. The take home message is that simultaneous contributions from various operators are required for a simultaneous explanation of the $R_{K^{(*)}}$ data alone. The results of this section will be useful for a better understanding of the results in the following section.

IV Methodology

IV.1 Parameter Estimation

The methodology adopted in this paper for the model selection is as follows:

Define Models:

Considering the NP Wilson coefficients real, we take all possible combinations (511 in total) of the coefficients forming a predefined global set of different scenarios. Each scenario with a specific combination of coefficients thus constitutes a potential ‘model’ to explain the experimental results.

Numerical Optimization:

Next, for each such ‘model’ $k$, as mentioned above, we perform a Frequentist statistical analysis optimizing a $\chi^{2}$ statistic which is a function of the Wilson coefficients. Whenever applicable, statistical (systematic) covariance matrices $V^{stat(syst)}$, are constructed by taking separate correlations. Theoretical uncertainties are propagated separately and are introduced in the $\chi^{2}$ in terms of a ‘theoretical’ covariance matrix $V^{th}$. The effect of the interplay of the SM uncertainties and the NP parameters come in the fit at a higher order and are neglected without any loss of generality. Following section II.1, we perform 4 types of fit for each ‘model’:

1. (a)

New data with Likelihood data for angular observables, a total of 214 observables. 2. (b)

New data with Moments data for angular observables, a total of 258 observables. 3. (c)

Old data with Likelihood data for angular observables, a total of 211 observables. 4. (d)

Old data with Moments data for angular observables, a total of 255 observables.

All fits are done in batch using Mathematica© in the form of a package Patra (2019). The chosen optimization method is ‘Differential Evolution’, a stochastic parallel direct search evolution strategy Storn and Price (1997) 444Capable of handling non-differentiable, nonlinear and multi-modal objective functions. Considerably faster than generic genetic algorithms and extremely able to find the global minimum..

Post-process:

In the post-process for each fit, we obtain fit-quality using $p$-value and find outliers by constructing a ‘Pull’ (related to Studentized residuals; for our purpose, the difference between the fitted and experimental results, normalized by the uncertainty of the data, including theory uncertainties Gross (2017); Demortier and Lyons (2008)) for each data-point. We also check the normality of the ‘Pull’-distribution (i.e. consistency with a Gaussian of $\mu=0$ and $\sigma=1$) to ensure the applicability of the $\chi^{2}$ as the fit-statistic. We use the “Cramér-von Mises” criterion Harald (1928) for the normality check. Scenarios not satisfying the normality criterion are dropped from the analysis.

Parameter-space:

Parameter uncertainties are obtained both from the Fisher matrix555In case of approximately Gaussian parameter-profile likelihoods, it is possible to obtain the ‘HESSE’ errors James and Roos (1975), which are, obviously, symmetric. and the profile-likelihood curve666Range of the $1\sigma$ confidence level (CL) of the profile likelihoods of the said parameter. One and two dimensional profile likelihoods in this analysis will be depicted as 1-CL plots, closely following the PROB method followed in ref. Aaij et al. (2016c)

With the remaining scenarios, we perform a model-selection procedure for each data-set. In the following sub-section, we elaborate the methods used to do the multi-model selection procedure.

IV.2 Model Selection

Following the ‘concept of parsimony’ Box and Jenkins (1970), we need to optimize the dimension (measure of the degree of structure) of the model explaining our data. All model selection methods, to some extent, depend on the principle of parsimony Breiman (1992). In statistical terms, this is expressed as a bias versus variance trade-off. In general, bias decreases and variance increases as the model-dimension increases.

IV.2.1 Cross-Validation:

As we have mentioned in the introduction, ‘cross-validation’ is the most generally applicable, powerful, reliable, and computationally expensive method for model comparison. The most straightforward and the most expensive flavor of cross-validation is “leave-one-out cross-validation” (LOOCV). In LOOCV, one of the data points is left out and the rest of the sample (“training set”) is optimized for a particular model. Then that result is used to find the predicted residual for the left out data point. This process is repeated for all data points and a mean-squared-error (MSE) is obtained using all those residuals. This process is repeated for all models. The models with the least MSE are the best ones.

IV.2.2 Criteria from Information Theory:

In addition to the extreme computational cost demanded by cross-validation methods, especially LOOCV, its applicability is questionable to very small sample sizes Beleites et al. (2013); Varoquaux (2017). Due to this reason, in our earlier works Bhattacharya et al. (2017, 2019), we have shown the importance and use of the information-theoretic criterion ‘AIC’Akaike (1998) and its second order variant ‘AICc’Sugiura (1978). It has been shown that minimizing AIC is asymptotically equivalent to cross-validationShibata (1981). For a detailed discussion on AICc, we point the reader to those papers and references therein. Here, let us reiterate the main important aspects of AICc in with respect to model selection in the present work:

AICc:

If the full reality or truth is noted as $f$ and an approximating model in terms of probability distribution is $g$, then we can define a model selection criterion in terms of the $\chi^{2}_{min}$ (equivalent to the maximum point of the empirical log-likelihood function) in the parameter space:

[TABLE]

where $n$ is the number of data points and $K$ is the number of estimable parameters777A more preferable way of calculating ‘number of estimable parameters’ is to calculate the $p$-value of the fit from toy Monte-Carlo (MC) method. Under the assumption that the fit-statistic follows a $\chi^{2}$ distribution, this can give us the number of degrees of freedom, and thus the number of estimable parameters. Still, as we need the differences between the AICc values instead of the absolute ones, the naive way of parameter counting works fine..

$w^{\Delta\text{AIC}_{c}}_{i}$:

The model which is the ‘closest’ to the unknown reality generating the data should have the smallest value of AICc among the considered models. Simple differences of them ($\Delta^{AIC}_{i}={\rm AIC}^{i}_{c}-{\rm AIC}^{min}_{c}$) estimate the relative expected information loss between $f$ and $g_{i}$ allowing comparison and ranking of candidate models in increasing order of $\Delta^{AIC}_{i}$. Generally, the level of empirical support in favor of $g_{i}$ is considered substantial when $\Delta$AICc is between 0 and 2 ($\Delta$AIC${}_{c}\leq 4$ is considered to be a conservative and loose bound). We can also quantify the weight of evidence in favor of model $i$ by defining a set of positive “Akaike weights”:

[TABLE]

adding up to $1$ Burnham and Anderson (2002). As these depend on the entire set, adding or dropping a model during an analysis requires re-computation for all models in the new set.

In the present analysis, we have a unique opportunity to not only test the relative capability of MSE from cross-validation and $w^{\Delta\text{AIC}_{c}}_{i}$, but also the validity of the empirical rule of selecting models with $\Delta\text{AIC}_{c}\leq 2$. To that end, we first select a large number of competing models by using the conservative limit of $\Delta\text{AIC}_{c}\leq 4$, and then distribute them in the plane of MSE${}_{\text{X-val}}$ vs. $w^{\Delta\text{AIC}_{c}}$ and check how they are clustered. Models with a low value of MSE${}_{\text{X-val}}$ and a high value of $w^{\Delta\text{AIC}_{c}}$ are the undoubtedly the best ones to explain the data.

V Results

V.1 Model Selection

As explained in section IV, we perform the fit for four different sets. After applying the normality check (as explained in sec. IV.1) on all the 511 models (thus ensuring only valid fits remain in our data-set), we pick out the large set of models with $\Delta\text{AIC}_{c}\leq 4$. The list of models, thus selected, are shown in figures 3 and 4, which are based on the analysis of all the available data sets given in sec II.1. Each point in these figures represents a ‘selected’ model (i.e. a model for which $\Delta\text{AIC}_{c}\leq 4$). Among these, the indices (labels for points) for models selected by $\Delta\text{AIC}_{c}\leq 2$ are framed. As is evident from the figure (and explained earlier in section IV.2.1), the lower the MSE${}_{\text{X-val}}$, the better the model. One can clearly see three separate clusters depending only on MSE${}_{\text{X-val}}$, and we can safely label the lowest one as the cluster of the best models (from MSE${}_{\text{X-val}}$) and discard the rest. Similarly, there are three clusters in the $w^{\Delta\text{AIC}_{c}}$ direction as well, where the cluster with the largest $w^{\Delta\text{AIC}_{c}}$ value contains only one model. All models with $\Delta\text{AIC}_{c}\leq 2$ lie in the two rightmost clusters. Following section IV.2.2, we know that a model with a larger value of $w^{\Delta\text{AIC}_{c}}$ is more probable to explain the data. So, we put a commensurate bound on $w^{\Delta\text{AIC}_{c}}$. We note that out of the various possible combinations, only a few are ‘selected’ by the combined criteria of MSE${}_{\text{X-val}}$ and $\text{AIC}_{c}$.

Figures 3a and 3b respectively compare the selected models in the analyses of moments and likelihood data on angular observables along with the new data on $R_{K^{(*)}}$. Similar comparisons are done with the old data on $R_{K^{(*)}}$ in figures 4a and 4b, respectively. We note that apart from the two-operator scenario $[{\cal O}^{\prime}_{7},{\cal O}_{9}]$, the likelihood and moments data pick up very different combinations of operators after model selection. The likelihood data prefer scenarios with more operators than those that are required for moments data. This could be since most of the observables which are determined from an unbinned maximum likelihood fit are consistent with their respective SM predictions. A large value of a single WC may lead to a discrepancy between the measured values and the corresponding SM predictions. It is thus preferable to have simultaneous contributions from different operators. The same observation can be made from the comparison of the quality of fits with moments and likelihood data which are shown in table 2 and 3, respectively. The values of the reduced $\chi^{2}$ and corresponding $p$-values, indicating the quality of the fit and the relative quality of the various models (scenarios) for a given set of data are provided in each corresponding table. We note that in general, the quality of the fits is better in the analysis with moments data compared to that with likelihood data. It can be seen that the ranking of the models depends not only on the quality of fit but also on the penalty function defined in eq. 12. This is one of the major advantage of AICc over data-fitting. Results of a similar analysis after dropping $R_{K^{(*)}}$ are presented in figure 9 in Appendix B. We note that after dropping the LFUV observables from the fits for the respective cases, the multi-operator scenarios can not survive the competition from those scenarios with relatively less number of operators.

As we know, the data from Belle has significant errors compared to those given by LHCb. To check whether this data has any influence over the selected scenarios, we have performed an analysis where we drop the $R_{K^{(*)}}$ data from Belle only. We note that the selected scenarios are the same as that given in figure 3. Also, the quality of fits or the respective $p$-values do not change. The reason is straightforward: the Belle data being consistent with LHCb, and the extracted $R(K)$ and $R(K^{*})$ in the selected models in figure 3 are consistent with both LHCb and Belle. Hence, the conclusions remain unchanged even after dropping Belle data.

There are a few other relevant observations here as well. For example, both in the case of the ‘Old’ and ‘New’ data sets, the majority of the selected scenarios are with two, three, or four operators and all of them contain the operator corresponding to the WC $\Delta C_{9}$. Also, the only scenario with a single operator selected by the ‘New’ data set (in particular the ‘Moments’ dataset) has the WC $\Delta C_{9}$ (model 2). Model 2 is clearly the better option for all the datasets. However, for the ‘New Moments’ dataset, scenarios with two operators corresponding to the WCs [$\Delta C_{9}$, $C^{\prime}_{9}$] (model 18), [$\Delta C_{9}$, $C^{\prime}_{7}$] (model 10), and [$\Delta C_{9}$, $C^{\prime}_{10}$] (model 20) are other possible choices. An important feature of the data is that the ‘Likelihood’ and ‘Moments’ datasets pick completely different types of models. For likelihood, an explanation of the observed data with a single operator is less plausible.

We note that some models have low MSE scores but are rejected by AICc. Although leave-one-out-cross-validation (LOO-CV) is asymptotically equivalent to AIC, there are differences between the two. Theoretical considerations aside, AIC is just likelihood penalized by the degrees of freedom. Evidently, AIC accounts for uncertainty in the data (-2Log(L)) and assumes that more parameters lead to a higher risk of overfitting (2k). Cross-validation (CV) just looks at the test set performance of the model, with no further assumptions. There is no explicit measure of model complexity, unlike AIC. Clearly, AIC penalizes model complexity more than CV. The accepted practice in the literature is that if one cares mostly about making predictions and assumes that the test set(s) to be reasonably similar to the validation sets, one should go for CV (only with a large number of data) arlot2010.

Following the reasoning in the paragraph above, it seems clear a priori that models with low MSE but not selected by AICc may have more parameters. In our analysis, the models are organized in such a way that more complex models with a larger number of parameters have a higher index. A quick examination of the plots in figures 3 and 4 validates this assumption, as the left sides of the plots are littered with models of higher indices.

On another note, as the title suggests, this study clearly finds that under the simultaneous application of both AIC and CV, the models cluster in such a way that lets us carry out a tighter selection of models than would have been possible by the use of any one of these criteria. So the fact that some of the selected models by CV are further discarded by AIC due to relative complexity is actually quite a non-trivial finding.

The best fit values of the new WCs with the corresponding errors for the selected models in Figure 3 are given in tables 2 and 3. We note that for all the fit results, the allowed values of $\Delta C_{9}$ are negative and greater than one. Other WCs appearing along with $\Delta C_{9}$ as probable solutions have values $<<1\ \text{or}\ -1$. The selected models with $C_{7}^{\prime}$ as one of the WCs will impact radiative decays. Hence, along with the allowed values for the new WCs, we have shown that in all the relevant models the branching fractions for $B\to K^{*}\gamma$, $B\to\phi\gamma$ and $B\to X_{s}\gamma$ are consistent with their respective measurements within 1$\sigma$ confidence interval in the relevant tables. The fitted values of the selected WCs remain almost unchanged in the analysis obtained after dropping $R_{K^{(*)}}$. The allowed parameter spaces of the respective WC’s are similar in the analyses with old data-set which are given in the Appendix in tables 4, 5. However, note that the fit qualities are relatively better in the analysis with old data which is an indicator of relatively poor NP-description of the new LFUV observables along with the angular observables.

We compare our fit results for various angular observables with their respective measured values and the SM predictions in figures 5 and 6 for a set of selected models. We see that the data in a few bins are inconsistent with their corresponding SM predictions, particularly in the data set for the ‘Moment’ analysis. For the ‘Likelihood’ analysis, most of the measured values in different bins are consistent with their SM predictions, albeit with exceptions. Assuming the observed discrepancies are due to NP, most of them can be resolved by our selected models. For the likelihood data, our analysis shows that the angular observables obtained from the selected models are fully consistent with their SM predictions as well as measured values. On the contrary, the fitted angular observables in the selected models in the moment-data analysis are shifted from their respective SM predictions in a few bins. Moreover, in some of those, the fit results shift from their measured values as well. One can also notice the correlations among various angular observables in the presence of different new operators from these figures. For example, $A_{FB}$ and $F_{L}$ are positively correlated, $F_{L}$ and $S_{5}$ are negatively correlated, etc. For some of the bins of $F_{L}$, the values predicted by the selected models are shifted from respective measured values. Similar plots obtained in our analysis after dropping $R_{K^{(*)}}$ is provided in fig. 10 in the appendix.

In figure 7, the predicted values of different observables and their respective correlations are shown for a few of the models (only one or two operator scenarios) selected in our analysis. We note that $R^{Central}_{K^{*}}$ and $R_{K}$ are correlated differently in different two operator scenarios. However, all of them satisfy the current experimental bounds on these observables. Interestingly, the single operator scenario with $\Delta C_{9}$ and the two operator scenario with $[\Delta C_{9},C^{\prime}_{7}]$ are unable to satisfy the current experimental bounds on $R^{Low}_{K^{*}}$. In these scenarios, satisfying the experimental bounds on $R_{K}$, it is hard to get a value of $R^{Low}_{K^{*}}$ below $0.85$. However, the scenarios with $[\Delta C_{9},C^{\prime}_{9}]$ and $[\Delta C_{9},C^{\prime}_{10}]$ as WCs can explain the observed data on $R_{K}$ and $R_{K^{*}}$. From the respective correlations, one can also see that $R^{Low}_{K^{*}}$ and $R_{K}$ are negatively correlated and values of $R^{Low}_{K^{*}}$ lower than $0.88$ prefers $R_{K}>0.8$. In figures 7c, 7d, 7e, and 7f, we provide the predicted values of the branching fractions of different radiative decays and their correlations with $R^{Low}_{K^{*}}$ in scenarios corresponding to the WCs $[\Delta C_{9},C^{\prime}_{7}]$. We note that there are no noticeable correlations between these branching fractions and $R^{Low}_{K^{*}}$, or for that matter with $R_{K}$ and $R^{Central}_{K^{*}}$.

The $q^{2}$ distributions and the zero crossing of the angular observables $A_{FB}$, $S_{4}$ and $S_{5}$ corresponding to the SM and the selected models are shown in fig. 8. We have noted discrepancies between the $q^{2}$ distributions for the SM and the selected models in the cases of $A_{FB}$ and $S_{5}$ while $S_{4}(q^{2})$ in the selected models are fully consistent with the SM. The $q^{2}$ distributions of these angular observables in different selected models overlap with each other. Hence it is hard to discriminate models from these distributions. In the Appendix in figure 12, we compare the values of $q^{2}$ at the zero crossing ($q_{0}^{2}$) between SM, our selected models and the measured values for the above mentioned observables.

The uncertainties of the observables, in terms of the parameters, are obtained by two different techniques for figures 7 and 8. For figure 7, which contains spaces for any 2 observables, we take, from our fit-results, all the information (best-fit values, uncertainties, and correlations) of only the NP parameters occurring in those observables, create a multivariate distribution out of those, and sample a large number of points ($\sim 5000$) from that distribution using Monte-Carlo. For each of those points, we get sets of values of the observables, which in turn lets us draw the $1~{}(39.35\%)$ and $2~{}(63.21\%)~{}\sigma$ contours (for 2-observable plots) from these data-sets.

For figure 8, we are dealing with only one observable at a time and that too, for a specific value of $q^{2}$. Taking the $1\sigma~{}(68\%)$ confidence levels for the marginal likelihoods around the central (for SM/nuisance parameters) or best-fit values (for NP) of the parameters as uncertainties and the corresponding correlations between them, we propagate the uncertainties to get the central values and uncertainties of the observables at some $q^{2}$. Doing this for many points over the allowed $q^{2}$ range and interpolating between them gives us the plots of $q^{2}$ distribution.

We should also mention here that the scenario $\Delta C_{9}=-\Delta C_{10}$ that arises, for example, in some Leptoquark models (among other model-dependent origins) does not pass the criteria of $\Delta\text{AIC}_{c}\leq 4$. We hence do not display or discuss this scenario.

VI Summary

In this article, we have analyzed the semileptonic $b\to s\ell\ell$ decays in a model independent framework with the relevant dimension six effective operators invariant under the strong and electromagnetic gauge groups. Our chosen set of operators does not include the four quark operators, chromomagnetic operators and tensor operators. Different possible combinations of all the effective operators have been considered, and following the statistical tools like cross-validation and the small-sample-corrected Akaike Information Criterion ($AIC_{c}$), we have found out the combinations which best explain the available data. We have provided separate analyses for the data on angular observables obtained from an unbinned maximum likelihood fit and that due to the principal moments of the angular distribution in $B\to K^{*}\mu^{+}\mu^{-}$ decay.

Among all the possible combinations, a relatively small number of one, two and three-operator scenarios satisfy the criterion of a selected ‘best’ model. All the selected scenarios contain a left-handed quark current with vector muon coupling as an operator (${\cal O}_{9}$). This is also the only one-operator scenario that survives the exclusion test in our search for the ‘best’ model(s). We have noted differences between the selected models in the analysis with angular data from likelihood fit and those from the principal moments analysis. The $R_{K^{(*)}}$, along with the angular observables associated with $B\to K^{*}\mu^{+}\mu^{-}$ decays, have played an important role in this selection. In the analysis with the new data on $R_{K^{(*)}}$, the scenarios with three, four and five operators are selected. This could be due to the fact that the tension between the updated measured values of $R_{K}$ and $R_{K^{*}}$ and their respective SM predictions have reduced in comparison to that for their old experimental measurements. Therefore, in order to explain all the data simultaneously, simultaneous contributions from different operators are required. We have noticed changes in the selected scenarios when we drop $R_{K^{(*)}}$ from the list of inputs. We have performed the analysis with and without the 2019 updates on $R_{K^{(*)}}$ from Belle and have compared them. We have noticed changes in the allowed parameter spaces for the Wilson coefficients of the selected scenarios.

We have compared our fit results for the angular observables with the corresponding SM predictions and the measured values in different bins for a few selected models. While our fit results are fully consistent with both the measured values and the respective SM predictions in the analysis of likelihood data, in the analysis of moment data, there are discrepancies between our fitted results and the respective SM predictions and the measured values in some bins. For some of the selected scenarios, we have studied the correlations between different observables, which show that the operator ${\cal O}_{9}$ and the combination of ${\cal O}_{9}$ and ${\cal O}^{\prime}_{7}$ (flavor changing electromagnetic dipole operator) can not explain all the available data on $R_{K^{(*)}}$ simultaneously. In particular, they have difficulty in explaining the observed results of $R_{K^{*}}$ in the low $q^{2}$ bin ($q^{2}\in\ [0.045,1.1]\ {\it GeV}^{2}$).

We have studied the NP effects in $R_{K^{(*)}}$ only, and noticed that the operator with a left-handed quark current with an axial-vector muon coupling (${\cal O}_{10}$) is the only one-operator scenario that can explain the data. Also, the parameter space for the corresponding Wilson-coefficient $\Delta C_{10}$, allowed by $R_{K^{(*)}}$, is tightly constrained by the measured values of $Br(B_{s}\to\mu^{+}\mu^{-})$. However, there are a few two-operator scenarios which have the potential to explain the current observation. Those operators are obtained from possible combinations of ${\cal O}_{9}$, ${\cal O}_{10}$ and operators like right-handed quark current with vector or axial-vector muon couplings (${\cal O}^{\prime}_{9}$, ${\cal O}^{\prime}_{10}$). In the two-operator scenarios, the allowed parameter spaces for $\Delta C_{10}$ or/and $C^{\prime}_{10}$ can comfortably explain the observed data on $Br(B_{s}\to\mu^{+}\mu^{-})$.

Note:The complete list of all models used in this analysis is too long to include in this draft. We have added an ancillary file named “models.json” along with this draft (to be found within the ‘arXiv’ source file). This file contains all combinations of WCs, relating them with their corresponding indices in our analysis.

Acknowledgements.

This project, S.N., and S.K.P. are supported by the Science and Engineering Research Board, Govt. of India, under the grant CRG/2018/001260.

Appendix A Best fit values of the new WCs in the analysis with old dataset

The selected models in the analysis with old datasets are shown in figure 4. The best fit values of the corresponding WCs along with their respective errors are given in tables 4 and 5, respectively.

Appendix B Results for model-selection with datasets other than ”New-moments” and ”New-likelihood”

A comparison of the selected models in figures 3 and 9 is helpful in understanding the impact of the new $R_{K^{(*)}}$ measurements by LHCb and Belle on the process of model selection. With more precise measurements on $R_{K}$ and $R_{K^{*}}$, the consistency between the data and the corresponding SM measurements may increase further. Therefore, it is important to gain insights on the probable NP effects in the angular observables. In figure 9, we provide the results for the model selection corresponding to the ”New-moments” and the ”New-likelihood” datasets after dropping the LFUV observables $R_{K}$ and $R_{K^{*}}$. The selected scenarios can be compared with the one given in figure 3. We note that in both the likelihood and moments data under the given selection setup, the number of selected scenarios reduce after we drop $R(K)$ and $R(K^{*})$ from the inputs. A comparison of between figures 9a and 3a indicates that the three-operator scenarios become less favourable. Similarly, from figures 9b and 3b we see that the four and five-operators scenarios are less favourable in the analysis without $R(K)$ and $R(K^{*})$. As explained in the main text, the explanation of new data prefers NP scenarios with more than two operators like three, four or five-operator scenarios. In particular, we have noted that the fit qualities improve once we drop these LFUV observables from the fits in general, which is at par with our expectations. For comparison with the ‘New’ data-set, figure 10 lists the angular observables from SM, experiment, and our fit results.

Figure 11 depicts the allowed parameter spaces of the most commonly occurring one and two parameter scenarios selected from different types of fits and data-sets. For one operator scenario (${\cal O}_{9}$), the allowed parameter space of the corresponding WC $\Delta C_{9}$, is shown in figure 11a. For the two operator scenarios, figures 11b, 11c and 11d shows the correlations between the WCs. We note that the allowed values of $C^{\prime}_{7}$ and $\Delta C_{9}$ have reasonably small ranges and they are negatively correlated. The corresponding value is $-0.316$. Large (negative) values of $\Delta C_{9}$ prefers large positive values of $C^{\prime}_{7}$. The other two plots show the allowed parameter spaces and the correlations of $\Delta C_{9}$ with $C^{\prime}_{9}$ with $C^{\prime}_{10}$, respectively. The fitted values of $C^{\prime}_{9}$ and $C^{\prime}_{10}$ have large errors. Guessing the exact correlations between them from the figures alone is therefore difficult. However, one can see that in the analysis with new data, the value of the correlation between $\Delta C_{9}$ and $C^{\prime}_{10}$ ( + 0.24) is greater than that between $\Delta C_{9}$ and $C^{\prime}_{9}$ which is $-0.11$.

In figure 12 we compare the values of $q^{2}$ at the zero crossing ($q_{0}^{2}$) between SM, our selected models and the measured values for the above mentioned observables. We note that in our selected models, the $q_{0}^{2}$ for all these three observables are in good agreement with the corresponding measured values.