Anatomy of a six-parameter fit to the $b\to s \ell^+\ell^-$ anomalies

Discrepancies between measurements of decay modes with an underlying quark level transition $b\to s \ell^+\ell^-$ and the corresponding standard model predictions have persisted for several years, particularly for the muon channels. Some of the largest deviations from the SM occur for the low $q^2$ bins of the observable dubbed $P_5^\prime$ as well as for branching ratios for the channel $B_s\to \Phi \mu\mu$, although the anomalies found in $R_K$ and $R_{K^*}$ are also very attractive since they make a case for lepton flavour universality violation. The apparent inadequacy of the SM becomes more compelling from the point of view of a global fit. For example, the authors of Ref. [1] have described 175 observables in terms of six parameters encoding new physics and have quantified the disagreement with the SM at about the $5\sigma$ level. These authors, and others, have considered in detail certain one and two parameter fits. In this paper we consider the full six-dimensional fit, defining metrics to quantify in detail where the deviations from the SM occur and presenting a visualization of some of our results in the six-dimensional parameter space. This allows us to explore the relation between single observables and the parameter space in the context of the global fit, and in particular also in the context of important correlations between observables.


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]. 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 transition: The six operators in question are The factorisation 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 ( ) ). Ref. [1], following the approach introduced in [2,9], discussed some 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 (SM and chirally-flipped) muonic operators. Here we will focus on this six-dimensional fit and provide a thorough discussion of its features. The best sixdimensional fit (BF) to the 175 observables considered in the latest iteration of the global fit results in the set of coefficients, 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 intrinsec error that a numerical minimisation routine always carries.
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. Our goal in this paper is to study this six-dimensional fit, to visualize aspects of it with the aid of the grand and guided tour [10], and to quantify the importance of the individual observables in the resulting fit.
Our paper is organized as follows. In Section 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 that differ from the BF by one sigma which 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 Section 3 we define metrics to compare the predictions of the fit at different parameter sets and apply them to all the observables. In Section 4 we discuss future observables designed to further test lepton universality and assess their potential impact on the global fit. Finally, in Section 5 we summarize our results and conclude.

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 in Eq. 1. Predictions for the experimental observables, as well as their theoretical uncertainties, are functions of these six parameters. A χ 2 function is constructed as described in Ref. [2] and its minimization leads to the BF point quoted above. To study this fit in detail and to quantify its uncertainties we will explore the neighbourhood of the global minimum of the χ 2 function. The minimization of this function has been performed numerically, and in fact the analytic form of the function is not known. To construct an approximate form of this function, we use the minimization procedure to provide us with a set of six-dimensional points (4959) near the minimum, S χ 2 (C i ), along with the corresponding value of their χ 2 function. This procedure also provides a numerical approximation to the Hessian matrix of second derivatives of the χ 2 at the minimum. With this set of points we construct an approximation to the χ 2 function near its minimum and we derive from it the Hessian. In this way we verify agreement with the Hessian found during minimization within the uncertainty of the approximations. This information will form the basis of our analysis.

Visualization of the BF region
We begin by visualising the six dimensional parameter region using the sample of points S 1σ and using the tour algorithm to obtain a sequence of 2-d projections that result in an animation of scatter plots (as shown here). 1 In this animation, the corresponding projection matrix for each view is shown as "axes". Watching the animation allows to build intuition about the features in the high dimensional distribution, e.g. observing if some feature observed in particular low dimensional projections persists when looking at a series of generic projections. 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. 2 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σ .
Most notably in this example we observe how the SM point moves away from the cloud of BF region points along with the change in the importance of C 9 in the projections. 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 Figure 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 the 6-d 1σ region. Figure 1: Visualization of 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.

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 six-dimensional confidence level ellipsoids around the BF that occur in this approximation. These ellipsoids can be used to study variations in the fit defined by parameter displacements along these principal axes away from the BF point. In this case there are twelve, six dimensional, sets of points C i that result in fits that differ from the BF by approximately one standard deviation. The Hessian at the minimum, in the basis (C 7 , C 9 , C 10 , C 7 , C 9 , C 10 ), is given approximately by In diagonal form, H D = diag(6621, 5647, 115.6, 72.6, 44.7, 6.1), exhibiting a hierarchy in its eigenvalues.
The parameter sets obtained with singular value decomposition (SVD) 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 detailed comparison between S 1σ (yellow) and its quadratic approximation (blue) is shown as an animation (here). Three static projections are shown in Figure 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 1σ). While the approximation is seen to work reasonably well in most directions, the first two views illustrate some inadequacies. The third view shows the projection onto C 10 vs C 9 illustrating the strong correlation found between these two parameters. To produce this visualization we have rescaled each parameter direction to range between 0 and 1 to facilitate comparisons.  Figure 2: Three projections of six-dimensional points in rescaled C i space ∈ 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, which sit on the boundary of this region. For context, the grey points are ∈ S 2σ but / ∈ S 1σ . For efficient rendering these grey points are not included in the animation.
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 giving the preference between two different parametric hypotheses H 0 and H 1 (H 0 ⊂ H 1 ) and presenting it in units of Gaussian standard deviation (σ). The bigger the Pull the greater the preference for H 1 over H 0 . 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 ). Using the notation ∆χ 2 H 0 H 1 = χ 2 min,H 0 − χ 2 min,H 1 , the Pull comparing H 0 and H 1 reads, where F is the χ 2 cumulative distribution function (CDF) and n dof = n pars,H 1 − n pars,H 0 is the difference between the number of free parameters characterizing H 0 and H 1 .
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.

Scenario
Best fit [1] Pull  Table 1: Comparison of certain low 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.

Hessian eigenvectors and the Wilson coefficients
The SVD directions are more than 80% 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 6 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 C 9 is the third most constrained direction, compare also Figure 2 (right panel). As can 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. Concerning the strong correlation between the two WCs we note that their contributions to R K and R K have opposite signs, i.e. while a positive C 9 increases the predicted value of R K and decreases the prediction for R K the opposite is true for C 10 . Noting that both observables have been measured to take values lower than the SM prediction, it is expected that some combinations of C 9 and C 10 will be more/less constrained.

Quantitative comparisons between 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: 3 We will use the parameter sets, p, corresponding to the SM, the BF point and the 1σ shifted points in what follows.
The Pulls are the metrics that compare theoretical predictions against experimental measurements and we present results for Pull(SM) and Pull(BF) below. These, of course, quantify the relative position of a measurement with respect to the SM and BF predictions respectively.
2. The pull difference will be used to compare the BF against the SM, quantified as 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) will then highlight the relative contributions from different observables to deviations by the fit from 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.
3. Different metrics will be used to evaluate variations in the fit itself, ignoring agreement with experiment. These allow one to associate specific observables with the uncertainty in the fit along one of the principal axes of the (approximate) one-sigma confidence level ellipsoid. These quantities are thus constructed to single out specific observables with large contributions to ∆χ 2 as the parameters move away from their best fit value. Several definitions are possible and we will compare the following ones: where the index i labels the observables and the different δs are all calculated for the twelve SVD directions. The different definitions have the simple interpretations: (a) δ 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.
(b) δ 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.
8 (c) δ σ is a generalization of δ to include correlations between observables. Neglecting the dependence of the theoretical errors on the parameters, we define δ σ,i such that 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). We interpret this as the conditional definition of δ.
(d) Finally,δ σ is a definition based on to the conditional variance 4 which turns out to be approximately equal to δ σ for the points we discuss below.
A Table that illustrates the largest δ's for the different definitions is provided in Appendix B.
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 stems mainly from correlated theory uncertainties, for example from definitions in terms of shared coefficients, shared form factor dependences among observables or from common sources of hadronic parameters. There are also important experimental correlations which we include when available. The resulting correlations between observables are illustrated in Figure 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 D). This is expected since these three observables are related by dΓ /dq 2 , the q 2 distribution of the process [15]. 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, between IDs 63 and 73 observables 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 to 38, for example, stand out in the experimental correlation map). See Appendix C for separate theoretical and experimental correlation maps.

Pull
Here we compare all observables directly to the SM and the BF respectively in the top two panels of Figure 4. We have highlighted in red those observables with Pull greater than 2. The Pull(SM) metric simply updates known results from Ref. [2]: four of the largest over-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 Figure 4 presents the Pull(BF) metric for all observables. We observe that three of the standout points against the SM also stand out against the BF. Two of the P 5 measurements (IDs 52 and 108) 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 Figure 4 we show ∆(Pull) for all observables. The idea of this metric is to highlight those observables in better agreement with the BF point than with 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 Figure 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 App. B. The largest values without correlations are found for branching ratio measurements and for R K . The points that are singled out as large by both definitions are:  • 98: 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 other observables highlighted in the bottom-left panel as showing a pull difference between 1 and 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 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.
More generally, we observe that including correlations reduces the preference in pull for the BF over the SM for some angular observables. This is most apparent for P 5 because it is the observable where the biggest individual tensions have been found among all the angular observables. The reason is that the patterns in the residuals (i.e. the Pull) for these observables are more consistent with the covariance matrix from the point of view of the SM compared to the BF point. 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.
∆ σ (Pull) is a relative metric, comparing the BF and the SM given the correlations between observables. For individual observables the statistical preference between the BF and the SM is still encoded in ∆(Pull).

The BF-SM direction at 1σ from the BF
The quadratic approximation of Section 2 permits a quick estimate of the point within S 1σ that is closest to the SM, it lies in the BF-SM direction and 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 Figure 5. The left panel contains seven different outliers from the |∆(Pull)| plot, two are P 5 measurements and the rest are BR. One of them, ID 155, has opposite behavior to the others, its Pull from this point being larger than its Pull from the BF. Finally, the right panel shows that only two LHCb measurements of Br(B s → Φµ + µ − ) (at intermediate q 2 ) differ from the BF by more than 1σ.

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 [16,17]. The results are shown in Figure 6 which shows the following for all the 175 observables listed in Appendix D: • 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). The observables with values δ i > 1 (or δ σi > 0.84) are listed separately in Table 3. Several   13 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. 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 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 Figure 7 where we compare the 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   Figure 7: Distribution of δ σ for the P 5 observable at LHCb in different q 2 bins.
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.   Table 5: Ranking of observables by δ along directions 4,5,6.

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  Table 7: Ranking of observables by δ σ in the last three directions.  Table 8: Rotation matrix between spaces of δs and principal components for definition 1 (δ left) and definition 3 (δ σ right) 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 considered observables. The full rotation matrix transforming from delta space to the first six PC space is given explicitly in Table 8. 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 Figure 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 .
Another observation is that δ 3 (mostly C 10 ) is an important contribution to PC2 based on δ, but not relevant in the first two PCs when for PCA based on δ σ . This suggests that including correlations reduces the variance in that direction.
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 Figure 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. ID 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. 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,49); yet others no longer stand out in this picture such as P 5 (IDs 44, 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 Figure 9 (right). Branching ratios stand out the most and again, we observe a difference between K and K final states.  4 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 [18]. 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 forty eight new observables with their corresponding ID is provided in Appendix D.
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 Figure 10. Estimates of the experimental sensitivities expected for measurements of Q observables have been presented by the Belle II collaboration, see [19], Table 67.
It is again noteworthy the difference in the size of the theoretical error between the SM and BF cases. This is an effect stemming from the characteristic structure of these observables. Long distance non-perturbative effects, one of the main sources of hadronic uncertainties, cancel for them in the SM, since hadronic effects couple with identical strength to both muons and electrons. However, 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 the mentioned observables have a pole in their q 2 -spectrum, 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 Figure 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 Figure 11, quantified with the metric δ i .  Table 13.
The most sensitive observables to each of the 12 directions are marked in red and labelled in the figure.
We observe that there are a handful of observables that stand out regarding their potential to further constrain the parameter space: these are the so-called Q 1 , B 5 and B 6s . This can be understood by looking at their particular sensitivities to certain Wilson coefficients. As a function of Wilson coefficients, Q 1 only depends on C 7 ,9 ,10 , being one of the most stringent tests for the search of right-handed currents. The combination of this sensitivity of Q 1 and the fact that the current fit is least constraining on C 9 ,10 , highlights the importance of Q 1 for δ 6 (IDs 1, 7, 13, 19, 25) in Figure 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 Figure 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.
As already mentioned in [20] and stated here in a different manner, these observables provide excellent opportunities for setting more precise constraints on the Wilson coefficients encoding possible NP effects. As a caveat we note that our observations could change when correlations are included.  Table 13 for the 12 SVD directions. The largest ones for each case are highlighted in red and labeled.

Summary and Conclusions
We have studied the global fit to observables in the b → s + − system 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 N P and µ from them. 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 lying approximately one sigma away from the BF along the six axes of the one-sigma ellipsoid. We evaluated the theory predictions at these twelve points and used that information do discuss the uncertainty in the global fit in observable space. In Section 3 we defined quantitative metrics to evaluate the relative importance of each of the 175 observables to the global fit. These measures enabled a study of the interplay between observable and parameter space, in particular also in the context of correlated errors. Moreover, we can use 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, 21 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 Figure 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 discussion. 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. 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 9 ( −1.8) than the one obtained in a global fit.
A main focus of this study was then to consider each of the eigendirections of the Hessian matrix to assess the sensitivity of observables to shifts of 1σ away from the BF point. Our main findings can be summarized as follows: • Along direction one we see that the largest δ by far corresponds to B → X s γ (171) and that this is not changed by correlations. The parameter C NP 7 is constrained by the B → X s γ observable, which has a much smaller error than the variation allowed by the global fit as can be seen in Figure 6.
• A similar behavior, but to a lesser extent is shown by observable 16 which has the largest δ along direction two, corresponding to C NP 7 . Generally this direction is tightly constrained by the low q 2 measurements of P 1 . Again, this is also a feature in Figure 6 where the fit allows a much larger uncertainty for 16 than its experimental (or theoretical) error.
• Direction three illustrates clearly that the preference for the BF over the SM accumulates from many small contributions in the same direction rather than from a handful of observables. This is also seen in that the majority of the observations are clustered at small values of ∆(Pull) (or ∆ σ (Pull)), but their distribution shows an average preference for positive values.
• Direction four is sensitive to the LFV ratios R K ( ) and can be further probed by the proposed observables B 5,6s .
• 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 observables are found to be most sensitive to shifts in this direction.
• Direction six behaves differently, with observable 98 (R K ) having the most impact both with and without correlations for direction 6+ and observable 100 (R K ) dominating along direction 6−. IDs 68, 155 (large q 2 bins of Br(B → K µµ)) and 172 (Br(B s → µµ)) also play an important role in this direction.
• We catalogued the effects of the new LFV observables proposed in the literature. We corroborate the importance of Q 1 for constraining on C NP 9 µ,10 µ and of B 5,6s for constraining C NP 10µ Acknowledgments This work was supported in part by the Australian Government through the Australian Research Council. BC acknowledges financial support from the grant FPA 2017-86989-P and Centro de Excelencia Severo Ochoa SEV-2012-0234. We thank Dianne Cook for help with R and with visualization and Joaquim Matias for useful conversations. G.V. thanks the CERN theory group for their hospitality and partial support while this work was completed.
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' δ σ .  Table 9: Parameter sets obtained with SVD along with their correspondent ∆χ 2 for displacements (both ways +/−) away from the BF along the six eigenvector directions. The eigenvectors are ordered by decreasing value of the corresponding eigenvalues of H.
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.

B Samples of largest deltas
The specific value that makes a given δ large, and thus interesting, is arbitrary. Here we compare the different definitions introduced in Eq. 8. 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 Table of δs above these cutoffs are presented here. 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 ID 164 computed for q 2 bins exceeding 8 GeV 2 , 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. The columns δ σ and δ σ show the same large δs which follows from Eq. 8 as we always use the same covariance matrix (with theory errors estimated at the SM point). Finally, the reasonable agreement between the columns δ σ andδ σ confirms that the latter is a fair approximation to the former.

C Experimental and theoretical correlation matrices
We reproduce here the information in Section 3.1 but showing separately the theoretical and experimental correlations between observables, see Figure 12.

D Observables used in the fit
In Tables 11, 12 we list the 175 observables included in the fit along with their corresponding ID used throughout our paper. The additional observables that have been proposed in the literature are listed in Table 13.

ID Observable
Exp.