Flavour anomalies in supersymmetric scenarios with non-minimal flavour violation

Motivated by tensions between experimental measurements and SM predictions in $b\to s \ell^+\ell^-$ transitions, we present the first study of non-minimal flavour-violating Minimal Supersymmetric Standard Model (MSSM) scenarios contributing to the relevant Wilson coefficients to address the observed anomalies using SuperIso and MARTY. We calculate the full one-loop analytical contributions of the general MSSM to Wilson coefficients relevant for flavour anomalies, together with the anomalous muon magnetic dipole moment $(g-2)_\mu$. We show that, after imposing theoretical constraints on the flavour-violating parameters we can find scenarios in agreement with the experimental measurements that can address at the same time the tensions in flavour observables and in $(g-2)_\mu$.


Introduction
In recent years, impressive progress has been achieved in studying and measuring semileptonic B decays. In particular, neutral currents with b → s transitions offer a plethora of clean observables that have been under scrutiny as they present tensions with the Standard Model (SM) predictions. The first tension, at the level of 3σ, was reported in 2013 in the measurement of angular observables related to B → K * µ + µ − decay [1]. Since then, similar tensions have been observed in several decays, such as B → Kµ + µ − , B s → φµ + µ − and Λ b → Λµ + µ − [2][3][4][5]. In addition, LHCb measured lepton-flavour universality violating (LFUV) ratios R K ( * ) = BR(B → K ( * ) µµ)/BR(B → K ( * ) ee), that are predicted very precisely in the SM, and confirmed the tension with the SM with about 3σ significance for low dilepton mass squared (q 2 ) [6,7]. Interestingly, all these deviations point to a coherent and consistent pattern, and can find a common explanation from new physics (NP) contributing to the Wilson coefficients C 9 (as was shown in e.g. Refs. [8][9][10][11]).
While LFUV observables have theoretical uncertainties at the percent level (or below) due to the cancellation of hadronic uncertainties in the ratios, the rest of the b → s observables are subject to assumptions made for the non-local hadronic effects and generally suffer from larger theoretical uncertainties [12]. In this analysis, we consider the Minimal Supersymmetric Standard Model (MSSM) [13,14], which predicts a superpartner particle (sparticle) to each SM field, together with an additional Higgs doublet. As supersymmetry (SUSY) is not observed at low energy scales, it needs to be a broken symmetry of nature. To preserve some of the nice features of supersymmetry, it should be "softly" broken, namely by introducing a SUSY-violating effective Lagrangian L SOF T , that contains all necessary couplings and masses, adding up to 105 new free parameters. Until very recently, due to obvious computational challenges, the whole MSSM has been little studied. Indeed, more constrained SUSY models were devised to allow for doable calculations and computations, through well-motivated and seemingly reasonable, but not physically founded assumptions. Such models with simplifications at the GUT scale consider a handful of parameters, like the constrained MSSM (cMSSM) [15]. More recently, the phenomenological MSSM (pMSSM), which considers CP conservation and Minimal Flavour Violation (MFV) simplifications [16], entered within computational reach with its 19 free parameters [17][18][19][20][21]. These models fail to provide a SUSY scenario fully compatible with the aforementioned flavour anomalies if R-parity is conserved [22] (for R-parity violating models, see e.g. Refs. [23,24]).
In this work, following our preliminary results [25], we will go one step further and consider for the first time a more general setup, based on the assumptions of the pMSSM but including in addition Non Minimal Flavour Violation (NMFV) in the squark sector, as a candidate for the explanation of the flavour anomalies in the b → sll transitions. NMFV allows for sizeable Flavour Changing Neutral Currents (FCNC) effects coming directly from the squark mass matrices at the weak scale, whose off-diagonal entries are then considered as new free parameters with respect to MFV scenarios. We will first consider NMFV contributions to Wilson coefficients through the Mass Insertion Approximation (MIA) and show that the new FCNCs can highly affect the value of C 9 in some scenarios, while still being compatible with the rest of the b → s constraints. Then, we will present the first analytical calculation of the general contributions to C 7 , C 9 , C 10 and also (g − 2) µ in the full MSSM with 105 parameters, and their evaluation for particular NMFV scenarios with 42 parameters.
The paper is organized as follows: Section 2 describes the theoretical context of our analysis. In Section 3 the flavour-violating parameters are introduced in the mass insertion approximation, and their new contributions to Wilson coefficients are defined. In Section 4, the numerical setup for our scans is presented. Section 5 shows and discusses how the NMFV models may fit the flavour anomalies. In Section 6, we present the first full analytical evaluation of the Wilson coefficients and (g−2) µ , using MARTY [26,27], in the MSSM and their evaluation in NMFV scenarios, confronting the results to the expected experimental values. Finally, the conclusions are given in Section 7.

Theoretical Context
In the MSSM, the most general soft supersymmetry-breaking Lagrangian can be written as: −L SOF T = −L gaugino − L sfermions − L Higgs − L tril. , where the different terms are [16,28] • Mass terms for the gluinos, winos and binos: whereB,W andG are the bino, wino and gluino fields, respectively.
• Mass terms for the scalar fermions: 2) whereQ i andL i are the left-handed squarks and sleptons, respectively, with their right-handed counterparts :Ũ ,D andẼ (no right-handed (s)neutrinos are assumed). The indices i, j run over generation, and all scalar squared mass matrices are Hermitian.
• Mass and bilinear terms for the Higgs bosons: where µ is the supersymmetric Higgs mass parameter.
• Trilinear couplings between sfermions and Higgs bosons where A f ij are the general 3×3 complex soft SUSY-breaking scalar trilinear coupling matrices between Higgs fields (H u ,H d ) and sfermions, in generation basis.
Several mixing effects arise in the general MSSM. In particular, the electroweak gauginos mix together with the higgsinos and give rise to the chargino and neutralino mass eigenstates. The chargino mixing matrix in the weak eigenstate (W + ,H + u ,W − ,H − d ) basis is given by where µ is the Higgs quadratic coupling and M 2 the soft SUSY-breaking wino mass. The β parameter is related to the vacuum expectation values of the two Higgs doublets present in the MSSM by The 2 × 2 unitary matrices U and V which diagonalize the chargino mass matrix M χ are defined as Their explicit expressions can be found in e.g. Refs. [16,29,30].
In the other sectors, the MFV hypothesis limits the mixing of squarks to the third generation only. This approach is still widely used in the study of the MSSM. If the MFV hypothesis is relaxed for other generations, a rich mixing dynamic arises. Concentrating on the squark sector in such a NMFV model, and starting from the Lagrangian in (2.2), one can define the super-CKM (sCKM) basis so that it rotates the (s)quarks' superfields in flavour space, making the quark mass matrices m u,d diagonal. This flavour alignment between quarks and squarks does not imply diagonal squark mass matrices, and can yield substantial flavour-changing effects. In the same manner as Ref. [29], letf be a six-component vector, wheref L ,f R are spanning generation space. We can therefore write the 6 × 6 flavour mixed squared fermion mass matrices as Collecting all sfermion mass terms in (2.2) and using (2.9) This defines the relevant mass matrices for the squark sector: . Their complete expressions can be found e.g. in Ref. [31], and a thorough analysis of the various terms at play can be found in Ref. [29]. Following Ref. [31], we define where the M 2 U , M 2 D and M 2 Q are the soft breaking squark masses defined in (2.2), and m u,d are the diagonal up-and down-type quark masses. The various D terms are given by which are obviously flavour diagonal. Finally, the T u,d terms are related to the trilinear quark-squark-Higgs couplings in (2.4) by (2.14) The final mass-ordered squark mass eigenstates are obtained by introducing the unitary transformation to the matrices in (2.11): = RqM 2 q R † q , for q = u, d, and m 2 q 1 < · · · < m 2 q 6 , (2.15) with the matrices Rũ ,d containing the flavour decomposition information of the massordered squark mass eigenstates: (2.16) Transformations between mass and flavour eigenstates are needed to perform phenomenological analyses on the model, as its parameters cannot be accessed directly from the mixed final eigenstates. The complexity of such analyses grows rapidly with the allowed mixings and free parameters. The computational challenge is such that a complete analysis of the most general MSSM with its 105 free parameters is not feasible. We propose two approaches: one within the so-called Mass Insertion Approximation with 28 free parameters, and then within a subset of the MSSM including NMFV with 42 parameters.

The Mass Insertion approach to the NMFV-MSSM
The usual approach when studying NMFV effects in the MSSM is to use the MIA approach, introduced as early as 1989 in Ref. [32]. The MIA originates as a diagrammatic technique [32,33], allowing us to choose a basis where the quark-squark-neutral gaugino couplings are flavour diagonal. The flavourchanging effects are provided by non-diagonal contributions in the sfermion propagators, as shown in Figs. 1 and 2. The new SUSY contributions (to e.g. Wilson coefficients) are then proportional to the various off-diagonal elements.
The MIA is also defined algebraically through the Flavour Expansion Theorem (FET) [34] the main features of which we will summarize in the following. All sfermion squared mass matrices M can be decomposed as a sum of a diagonal diag(M ii ) ≡ M d i and a non-diagonalM ij matrix. Calculating loop amplitudes requires evaluating Hermitian matrix functions f (M ) of the involved mass matrices, which can be expanded following the FET's conditions: where the divided difference f [k] functions are defined in Ref. [34]. This expansion expresses the loop quantities such as Wilson coefficients in terms of the flavour-violating off-diagonal entries in the squark squared mass matrices. The following dimensionless ratio is usually introduced to define the mass insertions: where M 2 f is one of the fermion soft-breaking matrices in (2.2). As the full sfermion mass matrix is actually a 6 × 6 matrix spanning both generation and chirality indices (2.9), the actual mass insertion parameter is of the form (δf ij ) AB where i, j are generation indices, and (AB) ∈ {LL, LR, RL, RR}.
In this framework, we define the relevant Mass Insertions (MI) for our study. To be consistent with the constraints from Kaon observables [29,33,35], every off-diagonal  33 , 22 , 33 , 22 .
For the δ u LL insertion, following the definition of M 2 u in (2.11), we express it in terms of the soft-breaking squark mass matrix M 2 Q as .

(3.4)
All the relevant NMFV contributions to the C 7 , C 9 and C 10 Wilson coefficients are given in Appendix A.2.

Numerical Setup
In what follows, we present a study of NMFV contributions to the b → sll processes in terms of Wilson coefficients (given in A.2) and mass insertions. The model used is an extension of the phenomenological MSSM (pMSSM), where the new contributions arise from additional flavour violation sources in the form of mass insertions. No new sources of CP violation in L SOFT with respect to the pMSSM are included, and the degeneracy between the first and second generations of squarks is kept. The third-generation trilinear interactions A t , A b and A τ are allowed to vary, while the others are set to zero. As we will see in Section 6.2.2, the slepton sector contribution to (g − 2) µ can be completely decoupled from the squark sector analysis. Therefore, no flavour-violating effect is turned on in the slepton sector, as they do not contribute to the b → sll observables. The Standard Model sector parameters are given in Table 1. The 28 input parameters for our model and their ranges (pMSSM+MIA) are given in Tables 2 and 3. They are randomly sampled from a uniform distribution. The spectrum is calculated at the electroweak scale using the code SOFTSUSY [36]. Following Ref. [30], we have introduced an average squark mass (over the first two generations), obtained from SM parameter Value

Constraints
In the following, we give the constraints considered in our study, both during and after the sampling of the parameter space. First, no tachyonic spectra are kept. This is a built-in condition in many spectrum calculators such as SOFTSUSY which is enforced during execution. We discard any spectra with a charged Lightest Supersymmetric Particle (LSP) to ensure the possibility for the LSP (often the lightest neutralino) to be a viable dark matter candidate. We impose further the latest available mass limits from supersymmetric searches given by Ref. [41]. No additional ab initio constraints are imposed, in order to keep the study as general as it can be.
As the SLHA1 [42] file format does not implement flavour mixing, the spectrum yielded by SOFTSUSY is obtained without considering the flavour-violating MIA parameters. Therefore, the spectrum considered here is pMSSM-like. The δs are considered as additional free parameters, that do not intervene in the computation of the spectrum. This can be justified a posteriori by considering constraints on the MI, as the approximation should be valid if they are small enough with respect to the diagonal mass parameters.
The following limits on the MIA parameters are also considered a posteriori: • To avoid tachyonic sparticles, all the MI parameters' ranges are reduced to • From vacuum stability arguments [30,43] The flavour-violating parameters that contribute the most to C 9 are (δ u 23 ) LL and (δ u 23 ) LR , in the chargino penguin diagrams such as the ones shown in Fig. 1, which are mainly constrained by (4.2) and (4.3). On the other hand, in thed sector, the gluino loops contribute mostly to C 7 which is already strongly limited by experimental data. Therefore, considering all double mass insertions as negligible, no constraints on the other MI parameters are imposed. A comprehensive discussion of the allowed ranges for these parameters can be found in Ref. [44]. Finally, all spectra should be considered with particular care, as flavour mixing can significantly affect the squark masses and their expected signal topologies at colliders. Also, the recast of LHC limits for general MSSM models is a non-trivial task [45,46], which goes beyond the scope of this study. Therefore, no particular limits on the sparticle masses are considered, apart from the model-independent ones present in Ref. [41].

Results and Discussion
The mass insertions allow new sources of FCNC, which give sizeable contributions to flavour observables by significantly shifting the relevant Wilson coefficients. In Fig. 3, we present the scan results with about 2 million model points. We can see an oystershaped spread of the pMSSM distribution upon turning on the NMFV contributions in the (C 9 , C 7 ) plane. In the (C 9 , C 10 ) case, we can see an isotropic spread of the pMSSM distribution in all quadrants, indicating a homogeneous behaviour of the two Wilson coefficients under flavour violation in the squark sector. On the other hand, in the (C 9 , C 7 ) case, the largest contribution to C 9 can be obtained by shifting C 7 significantly from its SM value, which is strongly constrained by the b → sγ data. However, it is clear from the impressive spread that the flavour anomalies can be given a satisfying answer using this framework, while still having reasonable values for C 7 . In Fig. 4a, a zoom in the region of interest in the (δC 9 , δC 7 ) plane is presented, together with the global best-fit patches from Ref. [47]. δC i is defined as The pMSSM distribution is shown in red, and the corresponding NMFV points are shown in blue. Imposing the constraints discussed in Section 4.1 yields Fig. 4b with 1 721 remaining points. We can see that even if the highest density of model points can be found away from the C 7 best-fit region, the presented NMFV model succeeds in proposing valid scenarios. In particular, several points seem to completely account for the flavour anomalies in the B sector, but further exploration of the full model spectrum is necessary. Also, it is clear that the pMSSM alone cannot give sufficient contributions to C 9 and C 7 : The orange bands represent the 1σ best-fit regions from Ref. [47].
the red distribution can at most account for half of the required shift in C 9 to explain the anomalies, with no other constraints imposed on the pMSSM parameters. Indeed, Fig. 5 clearly shows this feature, where we can see the spread of C 9 for both pMSSM and NMFV models, with no particular constraints on the sampled points. The pMSSM, while being able to provide compelling shifts, fails to fully account for the anomalies (best fit given in Refs. [47][48][49]) as was shown already in Ref. [22], whereas the NMFV is capable of providing hundreds of compatible scenarios if no other constraints are considered.
To examine these best-fit points, one can look at the associated mass spectra for some well-studied collider SUSY signals like electroweakinos and coloured sparticle states. In particular, in Fig. 6, we show both MI parameters and the LSP's mass distribution for our candidate models, without imposing constraints (left) and after imposing constraints (right) for some of the best points with respect to the expected C 9 shift. We see that a 10% fraction survives the tachyonic and vacuum stability constraints while still offering valid candidates for the flavour anomalies.
Similarly to what was shown in Refs. [30,50], it is mostly the top row diagrams in Fig. 1 corresponding to chargino interactions that contribute the most to C 9 , which corresponds to (δ u 23 ) LR , (δ u 23 ) LL terms in the MIA. For C 7 , the major contributions come from (δ d 23 ) LR in the gluino diagrams. However, the effect on the Wilson coefficients shown here is the result of a global effect coming from all the contributions. The correlation between the free parameters (pMSSM+MI) and the best (C 9 , C 7 ) values was not found to point towards a specific direction.
The effect of the constraints on the most important MI parameters is also shown in Figs. 6d and 6c. From left to right, the available parameter space in the (δ LL , δ u LR ) plane is reduced from This shows that vacuum stability constraints onc L −t R mixing are the most stringent ones, as expected from large average squark masses, i.e. M sq 2m t . The other MI parameters' contribute very little to C 7,9 and can be neglected and/or kept close to zero.
The results clearly show the interest of NMFV scenarios, and the need of their further exploration. Indeed, the main advantage of the MIA in our case was to easily explore the pMSSM extended with flavour violation, with direct access to the flavour-violating parameters instead of the final mass eigenstates. Also, it has the advantage of reducing the model's free parameters, if their contribution is not significant in the subject at hand, which we did by keeping fewer than 30 parameters, instead of O (50) or O(100). However, due to the obviously expected effect on the sparticle spectrum, a more general and complete approach without approximation is necessary to completely confirm the model's shown interesting features. Moreover, a complete approach should also evaluate the contribution of such models to the muon (g − 2) µ . This is precisely what is addressed in the next section.

Analytical calculations in NMFV-MSSM scenarios without approximation
In the pMSSM, analytical calculations have been performed for several one-loop quantities such as C 7 , C 9 [51] and (g − 2) µ [52].
In NMFV scenarios, some calculations have been performed at the one-loop level (see e.g. Refs. [53][54][55]), but the general contributions to C 7 , C 9 and (g − 2) µ are not known. In the following sections, we present the methods that we used to derive analytically these
We used MARTY [26,27] to calculate automatically all the involved Feynman diagrams and extract the coefficients (g − 2) µ , C 7 and C µ 9 . The number of diagrams for each contribution is presented in Table 4. As MARTY counts left and right Dirac projectors P L and P R as independent vertices, the number of diagrams is larger than what a standard counting method would imply.χ 240 96 * 24 * 24 0 C 9 /γ−penguins 240 96 * 24 * 24 0 C 9 /Z−penguins 624 1344 * 240 * 78 0 C 9 /boxes 864 13824 * 0 12 0 Table 4: Number of diagrams for each contribution calculated by MARTY. NMFV-specific contributions are the starred numbers. By definition, C 7 and (g − 2) µ only receive contributions from γ-penguin diagrams. There are in total 17949 Feynman diagrams.

Numerical evaluation
The mathematical expressions resulting from the sum of thousands of one-loop diagrams are too large for any analytical purpose. In order to obtain predictions, MARTY generates a numerical C++ library containing functions evaluating the results given a general MSSM scenario. From a set of values for the SUSY-breaking parameters presented in Eqs. (2.1) to (2.4), we are therefore able to evaluate the exact values of C 7 , C µ 9 and (g − 2) µ at the one-loop level in the library generated by MARTY.
While MARTY also generates a tree-level spectrum generator to calculate masses and mixings from the initial model parameters, loop corrections are known to be large and we therefore use SPheno [56,57] to produce a more precise spectrum including looplevel corrections and phenomenological constraints. Finally, the values of the Wilson coefficients are given to SuperIso to apply renormalization group equations and evolve the coefficients down to the b mass scale, and calculate flavour observables.

Random scan
To sample the MSSM parameter space, we used a uniform random scan in 42 dimensions with NMFV only in the squark sector to reduce the number of free parameters. Input parameter ranges are presented in  33 [−100, 100] GeV (A u/d ) 11 [−0.1, 0.1] GeV (A u/d ) 22 [−100, 100] GeV (A u/d ) 33 [−10 4 , The scan efficiency is of about 0.05%, corresponding to physical scenarios for which SPheno can calculate a spectrum. For such a low efficiency there is a large bias in the selected scenarios. Consequently, we also present some posterior distributions of the , slepton (middle) and squark (right) masses. For particle families, the distribution corresponds to the lightest particle of the family. Chargino and gluino mass distributions extend up to 3 TeV and 7 TeV respectively. spectrum in Fig. 8. The scan could be refined with better constraints on the input parameters to improve the efficiency. The following analysis is therefore more a proof of principle rather than a complete phenomenological study of the MSSM parameter space. There are two visible biases in the posterior distributions of spectrum parameters: • Charged sleptons are lighter than sneutrinos because the range for M 2 E is smaller than that of M 2 L . • The lightest neutralino is always lighter than 400 GeV contrary to the lightest chargino. This is because we impose the condition of having a neutral LSP in order to be a dark matter candidate.
To improve the scan efficiency, we considered machine learning techniques to sample the parameter space. The purpose of these techniques is to create a sampling bias toward scenarios that generate valid model points, that therefore improves the scan efficiency. However, while these techniques can be implemented without much difficulty for the pMSSM with 19 parameters, the 43-dimensional space of the NMFV scenarios we present in this paper is too large for the machine-learning-based sampling to be established. Indeed, in the absence of prior knowledge on the distribution of valid parameters, and because of the high number of dimensions, no efficient sampler could be constructed with the considered techniques such as Normalizing Flows or Hamiltonian Monte Carlo samplers. Further work is required to build efficient samplers in highly dimensional unknown and little-constrained parameter spaces with very few acceptable points, which is beyond the scope of this work.
Finally, let us stress that the current LHC limits on SUSY particle masses are not directly applicable to our study. The NMFV-MSSM being a more general model than the so-called simplified or constrained MSSM scenarios, the recasting of collider constraints on the sparticle spectrum is a non-trivial task (see e.g. Refs. [45,46] and yields weaker bounds. We nevertheless checked for points leading to significant negative contributions to C 9 that they escape the direct limits, in particular due to the degeneracy between the lightest neutralino and charigno ∆m(χ 0 1 , χ ± 1 ) ≤ 1 GeV, which makes them extremely complicated to probe experimentally.

Results
Using as input the NMFV-MSSM spectra obtained with SPheno, the numerical functions generated by MARTY evaluate the full one-loop contributions to the Wilson coefficients and (g − 2) µ . As the scan is random, we show distributions for the different quantities that we calculated for the 70282 valid model points. In the following, we study the impact on the Wilson coefficients and (g − 2) µ separately. Then, the relation between the two will be discussed.

Wilson coefficients
The distributions for the NMFV-MSSM contributions to the Wilson coefficients C 7 and C µ 9 are presented in Fig. 9. Figure 9: Distribution of the Wilson coefficients δC 7 and δC µ 9 . The 1σ best-fit regions from Ref. [47] are shown in orange.
Both distributions are centered around zero, as expected. While the majority of δC 7 points are close to zero and the best-fit region, many scenarios are already excluded because of a large shift to this coefficient. For δC µ 9 , the best fit-region is shifted by −1 from the SM value. While it is possible to obtain substantial C 9 shifts in our scenarios, only a handful of them predict δC µ 9 < 0.2. It is important to note that the best-fit region for C µ 9 should not be considered as a discriminant criterion, any scenario between the SM and the best fit can still fit better flavour observables and should be carefully considered.
A 2D distribution of (δC 7 , δC µ 9 ) is presented in Fig. 10. It is clear that the constraint on δC 7 excludes several scenarios with δC µ 9 < −0.15. It seems nevertheless possible to address both coefficients, but a larger dataset is required to explore the region with large negative δC 9 .

(g − 2) µ and combined analysis
For just over fifteen years, the anomalous magnetic moment of the muon has proven to be a persistent tension between the SM  and experimental measurements. The most recent results obtained at Fermilab [79] have not only confirmed the Brookhaven 3σ − 4σ [80] discrepancy, but raised it to the 4.2σ level with a combined experimental average of a EXP µ = 116592061(41) × 10 −11 . However, multiple questions remain as lattice QCD calculations may reduce the discrepancy to only 1.6σ [81]. In the following, we investigate whether NMFV-MSSM models can account for the observed tensions in both (g − 2) µ and the b quark flavour sector.
Our present analysis does not strictly consider NMFV parameters in the lepton sector 1 as shown in Table 5. We present the numerical results for (g − 2) µ in the following. The mass distribution for charged sleptons is around the electroweak scale, i.e. a few hundred GeV (see Fig. 8). This implies significant contributions to (g − 2) µ that are shown in Fig. 11. As the experimental deviation is very small [79], it is not hard to address (g − 2) µ alone. Figure 11: Distribution of δ(g − 2) µ . Only scenarios with a positive shift are considered and the experimental measurement with its 1σ uncertainty [79] is shown in orange.
As shown in table 4, the lepton and quark sectors are sensitive to the neutralino and 1 There is no limitation for NMFV in the lepton sector, this choice has been made to reduce the number of free parameters and concentrate on flavour observables that are more difficult to address because of the C µ 9 shift. chargino mass scales. However, while there are slepton contributions in box diagrams for C µ 9 , these contributions are small and the latter coefficient is almost independent of the slepton masses. Figure 12 shows the dependence of (g − 2) µ and C µ 9 with respect to the relative slepton mass scale. 2 Figure 12: The variation of the relative mean absolute value of C 9 and (g − 2) µ for the entire dataset is plotted as a function of the relative slepton mass scale. The initial, non-modified dataset corresponds to the point at (1, 1). This analysis shows that by rescaling the slepton masses (charged sleptons and sneutrinos), one can shift the value of (g − 2) µ and let Wilson coefficients C 7 and C µ 9 remain stable. It is therefore possible to search for a scenario that fits the flavour observables well and adjust the slepton mass scale to address (g − 2) µ .

Conclusion
We presented a first study of the pMSSM extended with non-minimal flavour-violating couplings in the context of the tensions observed in b → s + − transitions with the SM predictions, and considered the SUSY contributions to the relevant Wilson coefficients. We carry on our study first by assuming the mass insertion approximation and show that the NMFV contributions allow us to shift the Wilson coefficients sufficiently enough to fully address the anomalies. After imposing theoretical constraints on the flavour-violating parameters we still find scenarios in agreement with the experimental measurements.
While the MIA provides a direct access to the flavour-violating parameters and eases the phenomenological studies by reducing the number of free parameters, a more general approach is necessary to fully assess the impact of NMFV contributions. Hence, in a second part, we calculated for the first time the full one-loop analytical contributions in the general MSSM to the relevant Wilson coefficients, as well as to (g − 2) µ using MARTY. By scanning the MSSM parameter space randomly and setting non-zero values for some of the flavour-violating parameters, we obtained 70282 valid scenarios with their individual spectra. In these scenarios, we showed that C µ 9 can be shifted towards the best-fit region given in Ref. [47] but that we have only a few points that shift C µ 9 in the favoured direction and let C 7 come close to the SM prediction. We then discussed the scaling of (g − 2) µ with the slepton mass scale that allows us to address (g − 2) µ without modifying the predictions for flavour observables. The present analysis is limited by the small sample of scenarios. As a perspective, the scan should be optimized by searching a parameter set that is more likely to produce physical scenarios. In particular, by looking at the posterior distributions of the input parameters it is possible to refine the scan, improve the efficiency and generate more scenarios to analyse. Finally, experimental constraints could be studied more in depth to compare the obtained spectra with direct searches of SUSY particles, in particular from LHC measurements. The obtained results are nevertheless very promising and show for the first time the impact of NMFV parameters in addressing the recent anomalies.

A.1 Effective Hamiltonian at the electroweak scale
The effective Hamiltonian for the decay B → X s + − in the SM and in the MSSM is given by (neglecting the small contribution proportional to V * us V ub ), in the basis of Ref. [30] where the O i operators read with V being the CKM matrix and q L(R) = (1 ∓ γ 5 ) 2 q. This Hamiltonian is known to next-to-leading order both in the SM [82,83] and in the MSSM [84][85][86]. For the B system, the operators and coefficients of interest for the anomalies are C 7 , C 9 , and C 10 .

A.2 Wilson coefficients
We consider the contributions to the Wilson coefficients as given in Ref. [30], including the correction as suggested in Ref. [50], which we reproduce here for completeness. The In the weak eigenstates basis the chargino mass matrix is given by [29]: where the index 1 of rows and columns refers to the wino, and the index 2 to the higgsino. µ is the Higgs quadratic coupling and M 2 the soft SUSY breaking wino mass. The 3 × 3 complex matrices U and V which diagonalize M χ are introduced: Their explicit expressions can be found in e.g. Ref. [87]. All the contributions to the Wilson coefficients are evaluated at the renormalization scale µ 0 = m W .

A.2.1 Chargino contributions
In the following, we give the contributions from chargino loops, such as the two top diagrams in Fig. 1.
Z-penguin with higgsino/wino loops: This diagram is proportional to (δ u 23 ) LR , which is one of the most interesting mass insertions, and is yet to be more constrained.
Z-penguin with two wino vertices: This is the same diagram as above, but with the exchange of two winos. They differ only by the specific mass insertion and the factor λ t /g 2 . Both diagrams are null in the limit of a diagonal chargino mass matrix, and so they are negligible for large M 2 .
Gamma penguin with two wino vertices: Gamma penguin with higgsino-wino vertex: The primed operators are obtained by switching the chirality of external states. All of these contributions are used, but we are most interested in the C 9 contribution coming from the γHW vertex.
Z-penguin with two wino vertices and a double mass insertion: Even though a double mass insertion corresponds to a higher order in the perturbative expansion, it has been pointed out [30] that this particular diagram could provide enhancement in the K system. For completeness, in the B decay, the contribution is: (A.12)

A.2.2 Gluino contribution
In this part, we collect all the contributions arising from gluino loops, from the bottom diagrams in Fig. 1, with γ-penguin: (A.14) The terms proportional to Mg can be dominant over the others. However, the mass insertion which enters the diagram is strongly constrained from b → sγ [32].
Double Mass Insertion -Z −g: For completeness, we give the gluino penguin contribution with a double mass insertion:

A.2.3 Box diagrams
The following contributions come from the box diagrams in Fig 2. Box diagram with wino exchange: where: Box diagram with higgsino-bottom-stop vertex: Replacing the wino with a higgsino yields

B.1 Hypergeometric functions and integral representations
To calculate the new Wilson coefficients in the NMFV/MIA framework, the following integrals need to be evaluated : These integrals can be shown to be linear combinations of hypergeometric functions p F q (for an extensive review see Refs. [88,89]). In most cases, the integrals can be rewritten using solely 2 F 1 , which is sometimes referred to as the Gaussian or ordinary hypergeometric function. This leads to nearly arbitrary precision in the numerical implementation as the series usually converge quickly. However, in some cases, particular analytical continuation formulae have to be used, which are well known in the literature, and can be found in e.g. Ref. [88]. The integral representation of 2 F 1 is defined by Euler's formula: which is a one-valued analytic function of z, provided Re(c) > Re(b) > 0, |arg(1 − z)| < π [88].
In what follows, we will simply refer to 2 F 1 as F . Equation (B.2) provides an analytic continuation of F , which is usually defined by its series representation: where : (a) n = Γ(a + n) Γ(a) , is the (rising) Pochhammer symbol.
The generalized hypergeometric functions p F q can be defined by and a general Euler integral transform relates hypergeometric functions of higher and lower orders [89]: For completeness, let us show that in this case too, we can re-express P ij1 (a, b) using hypergeometric functions: We are interested in the set of integrals defined by P ij1 (a, b) = This holds provided Re{1 − ax − b(1 − x)} < 1 . For x ∈ [0, 1], a, b ∈ R + , this is always true.
Using the general Euler transform (B.5) we can compute the integral over 2 F 1 : The two previous results for a = b and a = b can be checked against the explicit analytical expressions given before.