Anatomy of a six-parameter ﬁt to the b → s (cid:2) + (cid:2) − anomalies

Discrepancies between measurements of decay modes with an underlying quark level transition b → s (cid:2) + (cid:2) − and standard model (SM) predictions have persisted for several years, particularly for the muon channels. The inade-quacy of the SM becomes more compelling in a global ﬁt. For example, Capdevila et al. (JHEP 01 , 093. arXiv:1704.05340, 2018) described 175 observables by six parameters encoding new physics and quantiﬁed the disagreement with the SM at about the 5 σ level. While certain one and two parameter ﬁts have previously been considered in detail, we establish a framework for the detailed discussion of the full 6d ﬁt. We visualize and quantify the 6d 1 σ region around the best ﬁt point and deﬁne ﬁt uncertainties for both current and future observables. We then deﬁne metrics quantifying the deviations between measurements and both SM and best ﬁt predictions. These metrics relate observables to directions in parameter space, revealing their precise role in the ﬁt, thus providing guidance for future theoretical and experimental work. Some metrics


Introduction
Many measurements have been performed in recent years on decay modes with an underlying quark transition b → s + − , where includes muons and electrons. Not surprisingly, some of these measurements show deviations from the standard model (SM) predictions by a few standard deviations. More interesting is the fact that several of these deviations from the SM appear to be 'in the same direction', in such a way that when quantified by a global fit the discrepancy with the SM is at just over the 5 σ level.
Several global fits that allow for the possibility of new physics (NP) in a model independent fashion described in the framework of effective field theory have been performed in the literature [1][2][3][4][5][6][7][8][9]. The basis for our study will be the fit of Ref. [1] which allows for NP by floating six Wilson coefficients of the effective low energy Hamiltonian responsible for the quark-level transition: The six operators in question are The factorization scale is taken at 4.8 GeV, so that m b (μ b ) is theMS running b-quark mass at that scale, and the SM values of the Wilson coefficients are C SM 7,9,10 = −0.29, 4.07, −4.31 and C SM 7 ,9 ,10 = 0. New physics would be parametrized in a model independent way by deviations in these coefficients from their SM values, C i ≡ C SM i +C NP i (i = 7 ( ) , 9 ( ) , 10 ( ) , = μ), that is, we only treat the muon coefficients as free parameters. We assume these deviations to be real and drop any CP violating observables. We note that this fit does not include (pseudo) scalar or tensor operators.
Reference [1] discussed several scenarios with one or two non-zero NP contributions at a time in detail and analysed the general scenario where NP contributions are allowed to all six muonic operators. Updating their numbers 1 we find the best fit point (BF) to be, 2 C NP 7 = 4.3 × 10 −3 , C NP 7 = 0.019, C NP 9μ = −1.06, C NP 9 μ = 0.37, C NP 10μ = 0.34, C NP 10 μ = −0.04 .
These coefficients lead to a χ 2 lower than that of the SM by 39.9, indicating that the best fit differs from the SM at the level of 5 σ . In this particular fit, the coefficients C NP 9μ , C NP 9 μ , C NP 10μ and C NP 10 μ are so labelled because they are only allowed to differ from the SM for the muons, breaking lepton flavor universality. Bearing in mind that the six coefficients identified explicitly in Eq. 3 are the only ones allowed to differ from their SM value in this study, we will proceed to drop the superscript NP and the muon label for notational simplicity.
From this starting point, our paper provides a comprehensive statistical analysis in the full six-dimensional space, by visualizing aspects of it with the aid of the grand and guided tour [14], a systematic comparison between predictions of the SM and BF point against measurements, and a quantitative analysis relating parameter directions and individual observables, based on the Hessian approximation to the χ 2 function. Throughout this analysis we show for the first time how knowledge of correlated uncertainties is important in judging the relative importance of individual observables to the fit.
Our paper is organized as follows. In Sect. 2 we visualize the neighborhood of the best fit point, we present the Hessian matrix of second derivatives at the minimum of the χ 2 function and use it to construct twelve points characterizing the boundary of the one sigma region. These points will serve as benchmarks for our quantitative studies. We also use the Hessian to quantify the difference between the six-dimensional BF point and simple one and two parameter fits studied previously in the literature. In Sect. 3 we define metrics to compare the predictions of the fit at different parameter sets and apply them to all the observables in Sect. 4. In Sect. 5 we discuss future observables designed to further test lepton flavor universality and assess their potential impact on the global fit. Finally, in Sect. 6 we summarize our results and conclude. 1 The experimental results used here were up to date as of February 2019. 2 The minimisation of the χ 2 function leading to the BF above is performed by means of the Markov-Chain Monte Carlo Metropolis-Hastings algorithm. The small differences between the numbers quoted here and the results of Ref. [1] are a manifestation of the intrinsic error that a numerical minimisation routine always carries.

The best fit point and its neighborhood
The parameter space in our problem corresponds to the set of six Wilson coefficients C i encoding physics beyond the SM as per Eq. 1. Predictions for the experimental observables, as well as their theoretical uncertainties, are functions of these six parameters. The minimization of the χ 2 function has been performed numerically, and in fact the analytic form of the function is not known. It is thus desirable to have an approximate analytic form for this function which we construct in terms of the Hessian matrix of second derivatives of χ 2 at the minimum. This information will form the basis of our analysis in this section.
We obtain this matrix numerically in two different ways. We first use the minimization procedure to provide us with a set of six-dimensional points (4959) near the minimum, S χ 2 (C i ), along with their χ 2 . With this set of points we construct an approximation to the χ 2 function near its minimum and we derive from it the Hessian.Alternatively, once the BF is determined, the Hessian matrix can also be computed by explicit evaluation of the corresponding second partial derivatives of the χ 2 function with numerical routines. We verify agreement between these two numerical Hessian matrices within the uncertainty of the approximations.

Visualization of the BF region
We begin by visualising the six dimensional parameter region using the sample of points S 1σ (points of S χ 2 (C i ) which are at most 1σ away from the BF) and using the tour algorithm to obtain a sequence of 2-d projections that result in an animation of scatter plots (as shown here). 3 To produce the animation we first center all parameter directions, the centered points are then scaled with the standard deviation in each direction within S 1σ , such that all directions have comparable scales. 4 We then search for appropriate projections that illustrate the separation between the BF point and the SM point. Finally, to show these projections we center the view by subtracting the mean position of the points in S 1σ .
In this animation, the corresponding projection matrix for each view is shown as "axes". Notice that the axes are centered with respect to the point cloud and therefore the origin is somewhat shifted with respect to the BF point. Watching the animation builds intuition about the features of S 1σ , for example it illustrates how features observed in lower dimensional studies embed in the full space. Fig. 1 Visualization of the parameter space in six dimensions via a general two dimensional projection. The set S 1σ of points within 1σ of the BF is shown in yellow, the black symbols mark the SM point (box); the BF point (diamond); one dimensional best fits (upwards pointing triangle); and two dimensional best fits (downwards pointing triangle). The S 1σ cloud is seen to be separated from the SM point mostly along the C 9 direction Most notably in this example we observe how the SM point moves away from S 1σ . The animation provides a striking picture of the separation between the BF region and the SM that is mostly along the C 9 direction, a result that is well known in the literature. We further illustrate this with a static picture in Fig. 1, where the projection has been selected for showing the large distance between the SM and BF point. For comparison we illustrate the positions of the SM, BF and selected one and two dimensional best fits from Ref. [1] (listed in the first column of Table 1) as described in the caption. Another interesting feature that can be seen in the animation is that all these lower dimensional best fits are found roughly in the same half of S 1σ .

Quadratic approximation
The Hessian matrix is the standard tool to construct a quadratic approximation to the χ 2 function in the vicinity of the global minimum. The eigenvectors of this matrix correspond to the directions of the principal axes of the sixdimensional confidence level ellipsoids around the BF that occur in this approximation. There are twelve sets of points C i defined by the intersections between the 1σ confidence level hyper-ellipsoid and its principal axes that serve to gauge the behavior of the fit as one moves away from the BF within S 1σ . As these points are constructed through a singular value decomposition (SVD) of the Hessian matrix, we will refer to them as the SVD points in what follows.
The Hessian at the minimum, in the basis (C 7 , C 9 , C 10 , C 7 , C 9 , C 10 ), is given by Table 1 Comparison of certain lower dimensional global χ 2 one or two parameter fits to the global six-dimensional best fit point. We also include the comparison of these scenarios to the SM as presented in Ref. [1]. The listed Pull indicates that all these scenarios are better than the SM by at least 5 σ and that most of them are between 1 σ and 2 σ worse than the 6D BF

Scenario
Best fit [1] P u l l SM [1] Best fit (Quad) Pull

Fig. 2
Three projections of six-dimensional points in rescaled C i space (i.e. each coefficient takes values between 0 and 1) ∈ S 1σ for the exact numerical calculation (yellow) and the quadratic approximation (blue).
The black diamonds are the twelve SVD points listed in Table 9 In diagonal form, H D = diag(6621, 5647, 115.6, 72.6, 44.7, 6.1), exhibits a hierarchy in its eigenvalues, corresponding to some directions being much more constrained than others. In terms of H one has where the deviations of the six Wilson coefficients C i from their BF values are written as the vector (C − C B F ). The SVD points are reproduced in Table 9 in the appendix along with their corresponding χ 2 (with respect to the BF). The first column "EV" labels the SVD directions from 1 to 6 by decreasing eigenvalue, the ± refer to the two possible ways of moving along a particular direction away from the BF point. Within the quadratic approximation to the χ 2 function, all these twelve points have χ 2 = 7.1 and lie precisely 1σ away from the BF. The exact values of the χ 2 for these points are also shown in the table and are an indication of how well the approximation works in this case, they range approximately between 0.85 σ and 1.4 σ . A comparison between S 1σ (yellow) and its quadratic approximation (blue) is shown as an animation (here). To produce this visualization we have rescaled each parameter direction to range between 0 and 1 to facilitate comparisons. Three static projections are shown in Fig. 2, where the black diamonds show the 12 SVD points, and the yellow (blue) regions are the corresponding projections of S 1σ (quadratic approximation to S 1σ ). While the approximation is seen to work reasonably well in most directions, the first two views illustrate some inadequacies. This visualization also reveals new features of the fit, for example, the third view shows the projection onto C 10 vs C 9 illustrating a strong correlation found between these two parameters.
The Hessian approximation also allows us to quantify the distance between the global best fit point, Eq. 3, and certain one or two-parameter scenarios singled out previously in the literature which can be associated with simple NP models. The comparison of these scenarios against the SM was already presented in Ref. [1] in terms of the measure Pull SM , with a large value of this quantity indicating a large deviation from the SM.
This so-called Pull is a statistical measure, presented in units of Gaussian standard deviations (σ ), quantifying the level of agreement between two different parametric hypotheses H 0 and H 1 (H 0 ⊂ H 1 ) in describing a given data set. Let χ 2 min,H 0 and χ 2 min,H 1 be the minimum values of the χ 2 statistic under H 0 and H 1 , respectively (since H 0 is contained in the family of parametrizations described by H 1 , then χ 2 min,H 0 ≥ χ 2 min,H 1 ). For large fits, Wilks Theorem [19] guarantees that χ 2 H 0 H 1 = χ 2 min,H 0 − χ 2 min,H 1 will follow a χ 2 distribution with n dof = n pars,H 1 − n pars,H 0 degrees of freedom, being n pars,H 1 and n pars,H 0 the number of free parameters characterizing the two hypotheses. Then, the Pull comparing H 0 and H 1 reads, where F is the χ 2 cumulative distribution function (CDF).
In Table 1 we present one and two dimensional fits using the quadratic approximation to the χ 2 function. For the sake of comparison, we analyse all the scenarios (without electronic NP contributions) already studied in Ref. [1]. For each scenario, we assess its statistical significance with respect to both the SM and the BF, introducing Pull 6D in analogy with Pull SM to quantify the preference for the six-dimensional BF over the different 1D and 2D hypotheses considered.

Hessian eigenvectors and the Wilson coefficients
The SVD directions are mostly aligned with one or two of the Wilson coefficients as shown in Table 2. Thus, to a good approximation there is a simple correspondence between a given SVD direction and the C i as given in the table. The direction corresponding to the largest eigenvalue of the Hessian is that in which the χ 2 function changes most rapidly near its minimum and is thus most strongly constrained by the data. In this case the first two directions have similarly large eigenvalues (more than 50 times larger than the rest) and, as can be seen from Table 2, they correspond to the parameters C 7 and C 7 to a very good approximation. The interpretation that these two parameters are essentially fixed and there is very little room to play with them is compatible with statements made in the literature [2]. Interestingly a combination of C 10 and other coefficients is the third most constrained direction, in particular we observe a large correlation with C 9 , see Fig. 2 (right panel), with a correlation coefficient of 0.82. Curiously, this correlation approximates the pattern that would be expected from right-handed currents.
The constraints on C 9 ,10 arise primarily from P 1 and P 4 . As can already be seen in Table 1 of [2], both P 1 , in all its q 2 bins, and P 4 in the high end of the low-q 2 region, are very sensitive to these two coefficients, with the sensitivity to C 10 generally more pronounced. As already observed in that reference, these observables are within 1σ of the SM and thus restrict the range C 9 and C 10 can take. The coefficient C 10 is further constrained by measurement of Br(B s → μ + μ − ), leaving little room for departures from 0.

Metrics for quantitative comparisons
Here we introduce several quantitative metrics to compare different sets of predictions and observations in the next section.

Comparing predictions and observations
In order to quantify the comparison between the different predictions and the experimental results we will use the following metrics 1. The Pull, which measures the difference between the prediction T ( p) for a given set of parameters p and the observation O in terms of the uncertainty constructed by adding experimental and theoretical errors in quadrature and ignoring correlations: 5 We will use the parameter sets, p, corresponding to the SM and the BF point in what follows. 2. The pull difference will be used to compare the BF against the SM, quantified as where σ −1/2 is the square root of the inverse of the full covariance matrix (including both experimental and theoretical errors evaluated at the SM point). Notice that these measures allow an explicit and systematic comparison between Pulls in the SM and BF scenarios. However, the distribution of (Pull) will not follow a χ 2 distribution (but is measuring differences in units of the total uncertainty). 6 The absolute value ensures that a positive number indicates that the BF prediction is in better agreement with the observation, and a negative value signals better agreement of the SM prediction with the observation. (Pull) captures the improvement in matching the observed value for the BF as compared to the SM. Notice that (Pull) and σ (Pull) i have different connotations. While (Pull) measures the absolute preference for the BF over the SM for a given observable, σ (Pull) i corresponds to the difference in conditional Pull, taking into account the correlation with other observables. 5 This Pull is defined for each observable, while the Pull quoted in Table 1 is defined in the parameter space of the six Wilson coefficients. 6 Alternatively, one may wish to define a measure comparing how much each measurement contributes to the total χ 2 in each scenario, i.e. comparing quadratic Pulls rather than linear ones. We find that such a definition gives qualitatively similar results here (an explicit comparison is given in Appendix C). Since our aim is a direct comparison between SM and BF predictions and measured values (rather than statistical interpretation of the results), we will not consider such a definition here.

Examining variations in the fit within S 1σ
We turn our attention to variations in the predictions, ignoring agreement with experiment. The goal is to associate specific observables with particular SVD directions (and therefore specific NP parameters). These metrics are thus constructed to single out large contributions to χ 2 as the parameters move away from their best fit value. Several definitions are possible and we will use the following ones: where the index i labels the observables and the different δs are all calculated for the SVD directions. The different definitions have the simple interpretations: 1. δ is a straightforward comparison between the points on the 1σ ellipsoid and the best fit. The difference is quantified in terms of one standard deviation that ignores correlations. 2. δ is a variation which takes into account the different theoretical errors at different points in parameter space.
Note that δ → δ if the theoretical errors are the same. The comparison between δ and δ thus contains information about the variation in the theoretical errors. 3. δ σ is a generalization of δ to include correlations between observables. Neglecting the dependence of the theoretical errors on the parameters, we define δ σ,i such that We interpret this as the conditional definition of δ. 4. Finally,δ σ is a definition based on to the conditional variance 7 which turns out to be approximately equal to δ σ for the points we discuss below.  Table that illustrates the largest δ's for the different definitions is provided in Appendix B.
To provide some intuition for the interpretation of measures defined in terms of the covariance matrix, we present a short discussion of the two dimensional scenario in the Appendix D.
It is important at this stage to reiterate that the fit procedure of Ref. [1], which forms the basis of our analysis, uses a single covariance matrix with theory errors evaluated at the standard model point.

Correlation between observables
The covariance matrix is an important ingredient to the global fit which, at present, mainly includes correlated theory uncertainties. These theory correlations arise, for example, from definitions in terms of common coefficients, common form factor dependencies among observables or from common hadronic parameters. There are also important experimental correlations which we include when available, but note that not much information on error correlations has been released by the experimental collaborations. The resulting correlations between observables are illustrated in Fig. 3. This figure shows, for instance, the correlation between various B → K μ + μ − branching ratios (IDs 63-73), F L and A F B observables (a table listing all observables and their corresponding ID is presented in Appendix F). This is expected since these three observables are related by d dq 2 , the q 2 distribution of the process [20].
We see also that large q 2 bins are correlated between themselves, but not so much with lower q 2 bins. This appears in the correlation map, for example, observables with IDs between 63 and 73 are highly correlated except for IDs 68 and 73 which are only correlated with each other. For LHCb measurements of angular observables (IDs 15-62) experimental correlations for each bin are strongest for q 2 ∈ [2.5 − 4] (IDs 31-38, for example, stand out in the experimental correlation map). See Appendix E for separate theoretical and experimental correlation maps.

Pull
Here we compare all measurements directly to the SM and the BF respectively in the top two panels of Fig. 4. We have highlighted in red those observables with Pull greater than 2, i.e. those where the theory prediction differs from the measured value by more than 2σ . The Pull(SM) metric simply updates known results from Ref. [2]: four of the largest over- Values larger than 2 for the top panels (1 and 0.84 for the bottom) have been labeled in red and selected for discussion in the text predictions of the SM occur for BR and R K ( ) LHCb measurements, whereas the largest under-predictions occur for P 5 measurements by three different experiments (IDs 44, 52, 108, 128).
The top-right panel of Fig. 4 presents the Pull(BF) metric for all observables, showing overall better agreement between predictions and measured values. We note however that several tensions remain. Two of the P 5 measurements (IDs 108 and 128) are also above the BF prediction whereas an ATLAS P 4 (ID 127) falls below both the SM and BF predictions by more than 2σ .
In the bottom two panels of Fig. 4 we show (Pull) for all observables. Recall that this metric captures the improved agreement between BF prediction and measurement, as compared to the SM prediction, measured in units of the total uncertainty, and negative values signal preference for the SM. The results are shown both ignoring correlations (lower left panel) and including them (lower right panel).
The majority of the points are clustered at small values of (Pull) (or σ (Pull)), indicating insignificant resolving power between the SM and the BF. The distribution shows, however, that even among these points with small (Pull) there is an average preference for positive values. This is of course just the statement that the global fit prefers the BF to the SM, but the distribution shows how much of this overall preference is built from many small differences that go in the same direction.
Points with values of | (Pull)| ≥ 1 (left) and | σ (Pull)| ≥ 0.84 (left) are highlighted in red in Fig. 4. The particular cutoffs for this are rather arbitrary. For individual pulls, we have a statistical interpretation since (Pull) is normalised to the total uncorrelated errors. When including correlations, we follow the argument sketched in Appendix B.
(Pull) simply combines the two upper panels of Fig. 4 to highlight the observables where the BF is a significant improvement over the SM. This exercise shows, for example, that branching ratio measurements and R K which have a large Pull(SM) also have the largest values of (Pull). On the other hand, the low-q 2 bin of the R K * measurement, ID 99 has a large Pull(SM) but not (Pull), indicating that the BF does Amongst the angular observables, P 5 measurements (ID 44, 52) stand out in both Pull(SM) and (Pull), but the latter also highlights a comparable improvement over the SM in P 2 , ID 49.
To include the effect of correlations, we next compare (Pull) to σ (Pull). The points that are singled out as large by both definitions are: Of the observables in this list, (167) is the only one that has a large preference for the SM over the BF.
There are several observables highlighted in the bottomleft panel as showing 1 ≤ | (Pull)| ≤ 2, that no longer stand out when correlations are included (bottom-right panel). This is the case for measurements of P 5 in the last two bins of the low-q 2 region (IDs 44 and 52). The largest value, (Pull) = 1.3, for a P 5 observable is found for the LHCb measurement in the bin [4,6] (ID 44). When correlations are included this drops to σ (Pull) = 0.75 and indeed no P 5 observables are singled out in the bottom right panel. 8 The reason is that the correlations implied by the covariance matrix for this observable are more consistent with the SM than with the BF. This information is captured and penalised in the conditional definition σ (Pull). Similar effects are insignificant for e.g. R K which is found to have negligible correlation with the other observables. This of course does not imply that these observables are no longer relevant to the fit, as values of σ (Pull) are still significantly larger than zero. In other words, when calculating the differences in χ 2 between the SM and BF point after dropping a subset of observables, we find that the smaller set of observables singled out by σ (Pull) already accounts for the largest differences, while a smaller additional reduction in χ 2 is observed when removing the full set singled out by (Pull).
Another observation that results from this comparison is that there is a systematic offset in (Pull) for the BR measurements in observables 1-14, which is no longer present when considering σ (Pull), indicating that the effect is well described by the covariance matrix.

The BF-SM direction at 1σ from the BF
The quadratic approximation of Sect. 2 permits a quick estimate of the point within S 1σ that lies in the BF-SM direction, it has parameters After evaluating the theory predictions at this point, we find it is exactly at χ 2 = 8.0 from the BF. A comparison of the predictions at this point with both the SM and BF points displays the pattern observed as one moves from the SM towards the BF as can be seen in Fig. 5. Fig. 4, we generally note reduced improvement over the SM as expected, except for observables 155 and 167 (where the reduced tension removes it as an outlier from Fig. 5). Recall these are measurements of Br(B → K * μ + μ − ) in the high q 2 region from CMS and ATLAS respectively. This behaviour illus- trates some internal tension in the fit, the measured value of 155 is best accommodated by smaller NP contributions than found in the BF point, while still showing significant deviation from the SM prediction. Considering now the Pull difference with the BF point (right panel), we see that this direction is most constrained by two LHCb measurements of Br(B s → μ + μ − ) (at intermediate q 2 , observables 91 and 92), the LHCb measurement of Br(B 0 → K + * μμ) at high q 2 (observable 73) and of R K (observable 98).

Fit uncertainty
Here we present a direct comparison between theory and experiment for all the observables. In addition to the usual theory and experimental errors, we include an estimate for the uncertainty in the fit. For this purpose we will follow the framework used to discuss uncertainties in global fits of parton distribution functions [21,22]. The results are shown in Fig. 6 which shows the following for all the 175 observables listed in Appendix F:  [15,22] LHCb • The experimentally measured value of the observable with its error (black). • The SM prediction with the estimated theory error (green). • The BF prediction with the estimated theory error (brown). • Our estimate for the error in the fit calculated as the difference between the maximum deviations in theory predictions along one of the eigenvectors of H , within the 1σ region (purple). Figure 6 allows for a quick assessment of the overall picture. It is particularly interesting to study the prediction range within the 1σ BF region (the error in the fit). For example considering observable 99 we see that the discrepancy between the measured value and theory predictions cannot be resolved anywhere within this region, while other measurements may be better explained by points within 1σ of the BF, for example the previously mentioned observable 155.
The observables with values |δ i | > 1 (or |δ σ i | > 0.84) are listed separately in Table 3, showing which parameter directions are associated with large fit uncertainties. Several branching ratios, as well as R K ( * ) show large deviations in directions ±4, 6 which we saw in Table 2 refer mainly to C 10 and C 9 respectively. Note that this is the case for observable 155 which thus points to smaller values of C 10 than found in the BF point, and negative values of C 9 . There are also a few large values along direction ±3 which is dominated by C 10 . Three P 1 measurements at low q 2 show large deviations along directions ±2 which corresponds mostly to C 7 . When correlations are included, two P 2 measurements single out Distribution of δ σ for the P 5 observable at LHCb in different q 2 bins direction 5 which aligns mostly with C 9 . Finally, the case of δ ± 1 shows that the global fit permits a much larger variation of C 7 than the experimental constraint from B → X s γ . It is important to recall at this stage that a large δ means that the one-sigma region around the BF contains variations in the predictions for that particular observable along the specified direction that are much larger than its corresponding experimental or theoretical errors. The corresponding observables are thus strongly constraining the fit in the given direction.
A quick scan of Table 3 reveals that there are no large δs in the P 5 observables. Given the interest in this observable, we examine it separately in Fig. 7 where we compare the Table 4 Ranking of observables by δ along directions 1, 2, 3. In bold the changes when the SVD points are defined by the exact 1σ surface  Table 5 Ranking of observables by δ along directions 4, 5, 6. In bold the changes when the SVD points are defined by the exact 1σ surface   Table 7 Ranking of observables by δ σ in the last three directions. In bold the changes when the SVD points are defined by the exact 1σ surface different q 2 bins as they show different behaviour. There is a somewhat large variation in direction 6 in the first bin, a larger variation with direction 5 for bins [4,6] and [6,8] and finally direction 3 is behind most of the variation for the last bin.

Ranking observables
Absolute values of δ, or rather δ 2 , tell us how much each observable contributes to constraining a given direction in parameter space. Therefore this ranking gives an indication of how the eigendirections get constrained in the global fit, and the hierarchy in values of δ 2 allows us to judge the relative importance of each observable. Without correlations we list the largest five values of δ 2 for each direction in Tables 4, 5. Direction one (mostly C 7 ) for example, is mostly constrained by Br(B → X s γ ) (ID 171), as can be seen both from Table 3 or Table 4. A somewhat different picture is obtained when taking into account correlations, see Tables 6, 7, and direct comparison of the two gives indications of which observables are most sensitive to such effects. Direction two is constrained predominantly by low q 2 bin observations of P 1 , and direction four is dominated by the single observable 98 (LHCb measurement of R K ), especially when taking into account correlation effects. A very different picture is observed in direction three, which does not exhibit a large hierarchy in δ 2 . This indicates that it is really the combination of multiple observables that constrains this direction. Observable 98 is found to be relevant in constraining direction six, for which we find that the list of most sensitive observables is similar to that found for direction three. Direction five (mostly C 9 ) is not especially constrained by a single observable, as indicated by the absence of a particularly large δ 5± . The largest δ 2 in this case occurs for 57 or P 2 (rather than P 5 as one might have expected). Moreover, observable 57 is more constraining in direction 5+ than in 5−.
The most striking difference when including correlations occurs for observable 68 (a large q 2 bin measurement of Br(B 0 → K 0 μ + μ − ) by LHCb). While it appears in the first few positions in the rankings without covariance in several directions (3,4,5,6), it drops below our cutoff when including covariance. As the correlation maps show, 68 is part of a group of highly correlated observables. It is in particular strongly correlated with observables 73 and 155, which also drop out of the rankings when including correlations.

Variance in the fit
The rankings of the previous section provide information about observables in the parameter space defined by the eigendirections. We have already seen that several observables are important in constraining multiple directions. An alternative way of looking at this information is to study which parameter combinations result in the largest variance in theory predictions. One approach is therefore to perform a principal component analysis (PCA) on the set of delta vectors. PCA is an orthogonal linear transformation onto a coordinate system such that the first basis direction is aligned with the maximum variance in the data, the second basis is the direction of maximum variation orthogonal to the first coordinate, and the remaining bases are sequentially computed analogously. It can be used for dimension reduction as the first few principal components (PCs) capture most of the information.
For this we consider each observable as one data point with 12 parameters, the values of δ in the 12 shifted points. The first two principal components, for example, provide the directions with largest variations, and plotting the data points in these projections shows which observables dominate. Different information is captured by looking at each observable in isolation (using δ) or in the context of correlations within the global fit (using δ σ ), and we therefore reproduce this analysis for both cases.
For the PCA analysis the data should first be centered, i.e. the mean in each direction has to be subtracted. In our case, the mean values are close to zero so the effect of centering is not very large. We find very symmetric behavior: the main difference between plus/minus directions is just the sign of δ. This means that we can fully describe the 12 dimensional distribution in the space of the first six PCs. These six remaining PCs are found to contain considerable variance in the distribution: whereas the first PC explains 31% of the variance, the sixth one explains 8% when correlations are ignored. When correlations are kept the first PC explains 20% of the variance and the sixth one explains 13%. This suggests that all six dimensions (i.e. WCs) still allow for considerable variance in the predictions of the observables (recall that this is measured relative to the errors). The full rotation matrix transforming from delta space onto the first six PCs is given explicitly in Table 8. Notice the differences between the two rotation matrices, for example δ 3 (mostly C 10 ) is an important contribution to PC2 based on δ, but not relevant in the first two PCs when considering δ σ (we can already observe the large reduction of variance in that direction by comparing the rankings in δ 2 3 of Tables 4  and 6). We find that δ 6 is the only direction which exhibits strongly asymmetric behavior: for certain observables there are differences between the change in prediction in plus/minus directions, see Fig. 9 (left). This figure compares the values in the two directions of δ 6 σ (a similar but more crowded picture is found plotting δ 6 ), and shows as an extreme example observable 100, R K , for which the theory prediction varies significantly along one direction but not the opposite. δ 6− is the only one of the twelve points with a large negative C 9 , and to a lesser extent C 10 .
We now focus on the first two PCs to study which directions and observables are responsible for the largest variation. To get an overview of the distribution of δs we show the projection of observables onto the first two principal components in Fig. 8 in the form of so-called biplots. These show the projected data points, as well as a visualisation of the projection in the form of labeled arrows pointing outwards from the center. This format makes it easy to relate directions on the projection to the original parameters.
When considering each observable in isolation (left view), clear trends can be observed. For example, observables aligned with direction 6−, 5− and anti-aligned with direction 5+ are mainly branching ratio observations in bins of large q 2 (e.g. IDs 68, 93, …). There are differences in branching ratio observables depending on the final state: notably most observables with negative PC1 but positive PC2 correspond to decays into K * , while decays into K + and K 0 appear to take negative values in PC2 (e.g. IDs 98, 13,14). Angular observables on the other hand show a very different behavior. For example observables 28, 41 and 44 are found to have the largest values of PC1. Large q 2 bins are different, e.g. IDs 56, 57, 60, take small positive values in PC1 but large absolute values in PC2.
The picture changes drastically when considering correlations (right view), where the relevance of large q 2 bins of branching ratio observables is no longer dominant. Note also the different effect that including the correlations has on the angular observables: for some of them the differences appear to be enhanced (e.g. IDs 16 and 49); yet others no longer stand out in this picture such as P 5 (IDs 44 and 60).
It is further instructive to assess the impact of covariance in a more direct fashion, to this end we introduce a difference in terms of the absolute values, With the aid of this metric, most of the differences can be explained in two dimensions as the first two principal components capture about 80% of overall variation. We present the projection onto these first two PCs in Fig. 9 (right). This difference shows the same behaviour of BR observables already observed by comparing the results in

Limitations of the Hessian approach
As discussed above and in Appendix A, the quadratic approximation is found to be a good description to the full fit, but does not capture all details. 9 Most notably, asymmetries present in the exact χ 2 function are not captured in the approximation, and this leads to deviations of up to about 30% in (χ 2 ) for some directions. As a consequence several SVD points are not exactly 1σ away from the BF point when measured with the exact χ 2 function. An alternative approach would be to use the Hessian for the identification of the SVD directions, but to define the SVD points by the intersection of said directions with the exact 1σ surface. We have explicitly verified that we reach the same conclusions with both approaches, even though detailed quantitative results will of course be slightly different. A summary of these comparisons is given below.
The observables with δ values above the cutoff are listed in Tables 3 and 11, where underlined (bold) indicates they are only above the cutoff when evaluating the SVD points in the approximation (exact calculation). Several differences are found, most notably along direction 3 + , as may be expected 9 Note that this is expected by construction, the χ 2 function of [1] is defined in the linearized or gaussian regime, see [2]. Asymmetries in the χ 2 function are therefore induced by higher order terms in the expression of the observables as functions of the Wilson coefficients, which are expected to be small for |C NP i | 1.
from the results in Table 9. In those cases values of δ are close to the cutoff in both calculations, with differences typically smaller that 10%. The change also affects the rankings shown in Tables 4, 5, 6 and 7, where differences in ordering when using the exact calculation to find the SVD points are given in bold. While we find the exact values of δ in directions with the largest asymmetries to change notably (as expected), the effect in terms of hierarchies and ordering is very limited, and reordering only happens in a few instances where the absolute values of δ are very similar in both approaches.
Finally we have also re-evaluated the PCA. We find that the projections onto the first two principal components are similar in terms of how each of the δ directions contributes to the projection. Differences become more relevant beyond the first two principal components, since these directions carry less of the overall variance and thus are more sensitive to details. For the same reason these higher PC's are not important in the discussion.
These results confirm that the quadratic approximation is appropriate for the description of the fit to the level performed in this paper.

Additional observables proposed to test lepton flavor universality
Several additional observables have been proposed in the literature to test lepton flavor universality by directly comparing the distributions in modes with muons to those in modes with electrons [23]. In this section we assess their likely future impact on the global fit, pinpointing which of these best constrain each eigendirection. A list of the 48 new observables with their corresponding ID is provided in Appendix F. We begin with a direct comparison of the theoretical predictions for the SM along with their uncertainty (green), a theoretical prediction for the six dimensional BF (brown) and the uncertainty in the fit calculated as before (purple) in Fig. 10. Estimates of the experimental sensitivities expected for measurements of Q observables have been presented by the Belle II collaboration, see [24], Table 67.
Theoretical errors are larger for predictions for the BF point compared to SM predictions. This is because long distance non-perturbative effects, one of the main sources of hadronic uncertainties, cancel in the SM. In the presence of NP that distinguishes between muons and electrons, these uncertainties get reintroduced proportionally to the size of NP. A small set of observables have a pole in their q 2spectrum, which causes predictions within bins containing the pole to become unstable and their errors to diverge. These particular observables, which show large uncertainties are shown separately in Fig. 10.
These two figures already indicate that most of these proposed observables will increase the constraining power of the global fit, as their uncertainties are much smaller than the current uncertainty in the fit. Of course, this will require measurements with experimental errors comparable in size (or smaller) than the corresponding theory errors. We can be a bit more specific by associating with each of these observables the corresponding directions which they can constrain. This is shown in Fig. 11, quantified with the metric δ i . The most sensitive observables to each of the 12 directions are marked in red and labelled in the figure.
There are a handful of observables that stand out regarding their potential to further constrain the parameter space: these are Q 1 , B 5 and B 6s . As a function of Wilson coefficients, Q 1 only depends on C 7 ,9 ,10 , and is one of the most stringent tests for the search of right-handed currents. Since the current fit is poorly constrained in the C 9 direction, this results in large values of δ 6 for Q 1 (IDs 1, 7, 13, 19, 25) in Fig. 11. The case of B 5,6s is rather different. These two observables, particularly in the very low-q 2 region, provide direct access to C 10 and this shows up as large δ 4 (IDs 5, 6,17,18,47,48) in Fig. 11. Surprisingly, some large-q 2 bins of these observables also show large δ 4 , at the same level of some of their low-q 2 counterparts, suggesting high sensitivity to C 10 also in this region.
This illustrates the promising opportunities for setting more precise constraints on the Wilson coefficients (as also stated in [25]). As a caveat we note that our observations could change when correlations are included.

Summary and conclusions
In this paper we establish a framework to quantify the importance of individual observables in the results of a global multidimensional fit. We apply our framework to the global fit to b → s + − mediated observables with the six free parameters C NP 7 , C NP 7 , C NP 9μ , C NP 9 μ , C NP 10μ and C NP 10 μ of Ref. [1] (and for notational simplicity have dropped the indices NP and μ).
We began with a direct visualisation of the one sigma region around the BF point and its position relative to the SM in six dimensions. We then provided a quadratic approximation to the global fit in parameter space based on the Hessian matrix of second derivatives at the BF point. This construction was used to find twelve points characterizing the 1σ contour. These 12 points were used to assess the fit uncertainty for each observable in the fit. In addition they represent 1σ shifts along well defined directions in parameter space and are thus representative of parameter directions. In Sect. 3 we defined quantitative metrics to evaluate the relative importance of each of the 175 observables to the global fit. Section 4 presented a systematic study of these measures,  Table 14 for the 12 SVD directions. The largest ones for each case are highlighted in red and labeled comparing both the SM to the BF point and the BF point to the set of 12 representative points, thereby illustrating the interplay between observable and parameter space. Throughout this discussion we have emphasised the role of correlated errors, which we found to be important in the evaluation of single observables. Finally Sect. 5 applied the same framework to assess the likely impact of future measurements on the global fit.
The coefficient C NP 9μ is of particular interest, as it has been singled out by lower dimensional fits, always finding a large deviation from the SM in this coefficient. Indeed among the six parameters considered in the global fit, only a negative C NP 9μ presents the correct patterns for explaining some of the most striking anomalies, increasing the predictions for the P 5 observable while reducing the predicted values of R K and R K * , and it is therefore expected to play a major role in the global fit. Our visualisation of the six dimensional BF region confirms that the cloud of points within 1σ of the BF is clearly separated from the SM mostly along the C NP 9μ direction. At one sigma from the BF, in the direction of the SM, C NP 9μ is still 60% as large and still the largest of the NP parameters. The importance of R K and R K * can be appreciated in Fig. 5, which shows that after moving 1σ from the BF towards the SM, these observables (98, 100) still stand out: they have a large (Pull) preference for the NP point and (especially 98) already shows a large Pull against this move from the BF.
Correlated uncertainties (which at present include mostly the ones in the theory), play an important role in the discus-sion. They reduce the preference for the BF over the SM for angular observables when considering conditional measures such as σ (Pull) in place of the isolated metric Pull(SM), used exclusively in previous works. While Pull(SM) shows the largest discrepancies with the SM for BR, R K ( ) and P 5 observables, we find that the preference of the BF over the SM for some of the discrepancies is less apparent with the inclusion of correlations. This is because patterns in deviations between predictions and measurements may be more/less consistent with expectations from correlated uncertainties depending on the parameter point, and this information is needed to understand the impact of each observable on the total χ 2 . While correlation effects are negligible for the considered R K ( ) observables, important correlations are found for the uncertainties of angular observables. In particular for the P 5 they reduce the Pull difference between BF and SM point when measured by σ (Pull). We further note that some of the P 5 observables also disagree with the BF at more than the 2σ level. As pointed out in [1], the anomaly in P 5 is best described by a larger negative C NP 9μ ( −1.8) than the one obtained in a global fit.
More generally our work illustrates some of the internal tensions in the fit, for example through a global discussion of the Pull for all observables. We also show for the first time a comprehensive picture of all observables entering the fit, comparing SM and BF predictions to the measured value, while also showing experimental, theoretical as well as fit uncertainties, see Fig. 6.
A main focus of this study was then to consider each of the eigendirections of the Hessian matrix to asssociate the most sensitive observables to corresponding directions in parameter space. For this discussion directions are labeled from most to least constrained, with each direction being dominated by one of the six considered Wilson coefficients, see Table 2.
Our key findings for each direction can be summarized as follows (for ease of comparison we list ID numbers used in the figures with each observable): • Along the most constrained direction (direction one, corresponding to C NP 7 ), the most sensitive observable is by far B → X s γ (171), which has a much smaller error than the variation allowed by the global fit as can be seen in Fig. 6. This is confirmed with large values of the δ measures defined in Eq. 9.
• Direction two (corresponding to C NP 7 ) is also strongly constrained by a single observable, the low q 2 measurements of P 1 (16). Again, this is a feature in Fig. 6 where the fit allows a much larger uncertainty for 16 than its experimental (or theoretical) error, and confirmed by large values of δ.
• A different picture is found for direction three (mostly C NP 10 μ ), for which several observables show comparable sensitivities, demonstrating that the constraint accumu-lates from many small contributions in the same direction. Note that in this case correlations become more important, as can be seen by comparing the rankings in δ and δ σ , and neglecting correlated contributions may result in overestimating the sensitivity of particular observables, most notably in this case for measurements of large q 2 bins of Br(B → K μμ) (IDs 68 and 155).
• Direction four (mostly C NP 10μ ) is particularly sensitive to the LFUV ratio R K (98) and can be further probed by the proposed observables B 5,6s . Note that within the current six dimensional fit the tension with the SM and the fit projection along this direction is about 1σ , making such future probes especially interesting.
• Direction five (which mostly corresponds to C NP 9μ ) exhibits a similar behavior of not being especially constrained by a single observable. Interestingly P 2 (41, 49, 57) observables are found to be most sensitive to shifts in this direction, while neither R K (98) nor R K (100) appear in our rankings. Measurements of P 5 (44, 52) are found to be somewhat less sensitive than those of P 2 in our comparison, indicating that the normalisation to the uncertainties is important here.
• The least constrained direction (direction six) is mostly aligned with C NP 9 μ . In this case the parameter points of the +(−) directions are quite different, with C NP 9 μ = 1.72(−0.72). As a consequence we observe very different sensitivities between the ± directions, which is illustrated e.g. in the left panel of Fig. 9. Considering the fit uncertainties for future observables we find large potential of constraining this direction by measurements of Q 1 .
We have thus cataloged sensitivities of (current and future) observables and related them to specific direction in the six dimensional parameter space considered here, which can be used to inform future theoretical and experimental studies. Note also that the Hessian approximation given in this paper allows for quick estimates of χ 2 close to the BF point, and the set of representative points listed in Table 9 can be used to estimate current fit uncertainties on additional observables not considered here. Table 9 Parameter sets obtained with SVD along with their correspondent exact χ 2 for displacements (both ways +/−) away from the BF along the six eigenvector directions ( χ 2 = 7.1 for these points in the quadratic approximation). The eigenvectors are ordered by decreasing value of the corresponding eigenvalues of H EV C 7 C 9 C 10 C 7 C 9 C 10 χ 2 exact A Intersection of the one-sigma ellipsoid with its principal axes SVD produces the twelve points shown in Table 9. The eigenvalues are numbered in decreasing order. For comparison we list in the table the χ 2 difference between these points and the best fit calculated numerically and using the quadratic approximation where χ 2 = 7.1 exactly (this is the condition used to find the twelve points). The column χ 2 exact lists the number extracted from the code of Ref. [1]. The last column is used below to quantify a cut-off for 'large' δ σ . The small differences between the χ 2 exact and 7.1 gives us an indication of uncertainty in our discussion, and we note how this is slightly different along the different SVD directions. For comparison purposes we show in Table 10 the same information with points chosen at the intersections of the SVD directions with the exact 1σ surface.
Of course, the SVD directions are only defined in terms of the Hessian matrix, and the shape of the full 1σ surface may differ significantly from the assumed hyper-ellipsoid. While in this case the true shape is adequately described by the approximation (see Fig. 2), some deviations occur, as illustrated in the left plot of Fig. 12 which shows direction 5 in the C 9 − C 10 plane. In this case the asymmetry of the true 1σ contour results in large differences (as quantified in Table 9) in χ 2 between the quadratic approximation and the  Figure 12 (right) shows the projection onto C 9 − C 10 and illustrates the importance of the SVD directions. The correlations between parameter directions in the fit encode how combinations of parameters are constrained by the measurements. Important correlations are found between C 9 and C 10 , and the shape is indeed captured by the SVD points shown in black. On the other hand, selecting the intersection of WC directions with the 1σ surface will result in loss of information, the selected points would not be representative of the surface envelope.

B Samples of largest deltas
The specific value that makes a given δ large, and thus interesting, is arbitrary. Here we compare the different definithe directions are physically meaningful in how they enter predictions for the considered observables.
tions introduced in Eq. 9. The definition of δ and δ makes it reasonable to set the interest cutoff when they are equal to one, as this number corresponds to one standard deviation as measured by the total uncorrelated error. On the other hand, the definition δ σ can be related to the calculation of χ 2 as shown in the last column of Table 9 with the quantity i δ 2 σ . Since we are comparing points that lie within χ 2 = 7.1 (in the quadratic or Hessian approximation), we define the cutoff for a 'large' δ σ in this case as √ 0.71 ∼ 0.84. This choice singles out those observables that by themselves contribute 10% or more of the shift in χ 2 . Note that this approximation is not as good for directions 5 + , 6 ± , where a better cut-off might be 1. The subjectivity of this cut-off is mitigated with the introduction of the section with rankings in the text. A list of δs above these cutoffs are presented in Table 11.
The columns for δ and δ are nearly identical because the estimate of theoretical errors varies little between the BF and the twelve SVD points. The sole exception shown in the Table for δ 6+ > 1 occurs for observable 20, for which δ 6+ = 0.57. The difference arises from the estimate for the theoretical error, at the BF it is just over 60% larger than the corresponding estimate at the point 6+. Our numerical calculations imply that with the current level of precision in the parametrization of hadronic uncertainties (mainly form factors) the theoretical error estimates have a non-trivial dependence on the NP contributions to the Wilson coefficients. These errors are computed via a multivariate gaussian scan over all the nuisance parameters, i.e. form factors, decay constants, CKM matrix elements, etc. As explained in [2], a rescaling of the errors on the nuisance parameters is needed to ensure gaussian behavior. The rescaling factor is most relevant for some observables such as the Br(B → K * μμ) computed for q 2 bins exceeding 8 GeV 2 (ID 164), while it has almost no impact on observables that are less sensitive to form factors. We have checked explicitly that once the factor ensuring gaussian behavior is found, the results are independent of the exact numerical choice for this rescaling. Finally, the reasonable agreement between the columns δ σ andδ σ confirms that the latter is a fair approximation to the former.

C Quadratic differences between SM and BF Pull
In Sect. 4 we use the linear (Pull) to compare the SM and BF point Pull for each observable. As noted, we may consider a quadratic definition which more accurately captures how each observable contributes differently to the χ 2 function (in the absence of correlated uncertainties). Figure 13 shows 2 (Pull) for each observable. Comparing to the corresponding results in Fig. 4 we note that the qualitative picture remains the same, but larger differences are emphasised when considering 2 (Pull). Notably small systematic effects seen in Fig. 4 are no longer visible in Fig. 13. In addition, observable 14 no longer stands out when considering 2 (Pull), which is because the Pull is close to one for the SM prediction, and

D Interpretation of measures defined in terms of the covariance matrix
The definition of e.g. σ (Pull) and similar measures is not immediately intuitive. Here we present the explicit expression in a two dimensional scenario to aid with the interpretation.
In general the variance-covariance matrix can be decomposed as as σ i j = ik ρ kl l j (14) where ik = diag(σ 1 , σ 2 , . . . , σ N ) and ρ the correlation matrix, i.e. with entries ρ i j equal to 1 if i = j and ρ i j = ρ ji ∈ [−1, 1]. The inverse of the variance-covariance matrix is therefore given as where −1 = diag(1/σ 1 , 1/σ 2 , . . . , 1/σ N ) and in general the inverse correlation matrix is given via the adjugate matrix as In two dimensions the correlation matrix and its inverse can be written as Inverting the correlation matrix introduced a relative minus sign for the off-diagonal terms entering the χ 2 function, which now takes the form and we have introduced the shorthand notation D i = T i − O i for the difference between theory prediction and observed value for observable i. We can recover the uncorrelated definition by setting = 0. Alternatively, in the presence of correlations, the overall factor 1 1− 2 captures the additional information content available in that case. We see that the third term in Eq. 18 will depend on the signs of the D i relative to ρ. If both theory predictions for the considered parameter point differ from the observed value in the same direction, i.e. the D i are either both positive or both negative, the sign of the third term will be determined by that of the correlation coefficient , leading to two possible scenarios: • > 0: This indicates that the patterns in D i are consistent with the correlation obtained by varying the nuisance parameters, as a result the overall value of the χ 2 function is reduced. • < 0: Negative correlation is not consistent with our assumptions about the D i , this situation will therefore result in an overall increase in the value of the χ 2 function.
Analogous considerations hold if the D i have opposite sign, where < 0 will result in reduced values of χ 2 and > 0 will increase the value.
We now want to capture this information on observableby-observable level, e.g. in our definition of the Pull. To avoid working with the square root of the matrix let us consider the following definition for the Pull of observable i (corresponding to usingδ σ ): For our 2-d example, for the first observable we obtain Indeed this results in the same behaviour observed for the χ 2 function, i.e. for = 0 we recover the uncorrelated definition, elsewise the total absolute value of the Pull will be increased/decreased if the pattern in the D i is inconsistent/consistent with that of the correlation matrix ρ, with the additional contribution being weighted according to the correlation obtained via the nuisance parameters.
Notice that while the factor used in the definition of the Pull is somewhat arbitrary, here we recover the χ 2 function as i.e. in analogy to our definition without correlation, wherẽ ii the conditional variance. We recall that numerically the results of δ σ defined in terms of the square root of the covariance matrix, and those found forδ σ are similar, though the exact expression is more complicated already for the two dimensional scenario.
Finally we note that when going beyond two dimensions additional effects need to be taken into account, such as simultaneous correlation between groups of observables. While this yields somewhat more complicated expressions (obtained via the inversion of higher dimensional matrices), the overall picture is similar to that shown by the two dimensional example.

E Experimental and theoretical correlation matrices
We reproduce here the information in Sect. 3.3 but showing separately the theoretical and experimental correlations between observables, see Fig. 14.