$B_{s(d)}-\bar{B}_{s(d)}$ Mixing and $B_s\to\mu^+\mu^-$ Decay in the NMSSM with the Flavour Expansion Theorem

In this paper, motivated by the observation that the Standard Model predictions are now above the experimental data for the mass difference $\Delta M_{s(d)}$, we perform a detailed study of $B_{s(d)}-\bar{B}_{s(d)}$ mixing and $B_s\to\mu^+\mu^-$ decay in the $\mathbb{Z}_3$-invariant NMSSM with non-minimal flavour violation, using the recently developed procedure based on the Flavour Expansion Theorem, with which one can perform a purely algebraic mass-insertion expansion of an amplitude written in the mass eigenstate basis without performing any diagrammatic calculations in the interaction/flavour basis. Specifically, we consider the finite orders of mass insertions for neutralinos but the general orders for squarks and charginos, under two sets of assumptions for the squark flavour structures (\textit{i.e.}, while the flavour-conserving off-diagonal element $\delta_{33}^\text{LR}$ is kept in both of these two sectors, only the flavour-violating off-diagonal elements $\delta_{23}^\text{LL}$ and $\delta_{i3}^\text{RR}$ ($i=1,2$) are kept in the \text{LL} and \text{RR} sectors, respectively). Our analytic results are then expressed directly in terms of the initial Lagrangian parameters in the interaction/flavour basis, making it easy to impose the experimental bounds on them. It is found numerically that the NMSSM effects with the above two assumptions for the squark flavour structures can accommodate the observed deviation for $\Delta M_{s(d)}$, while complying with the experimental constraints from the branching ratios of $B_s\to \mu^+ \mu^-$ and $B\to X_s\gamma$ decays.


Introduction
It is known that Supersymmetric (SUSY) extensions of the Standard Model (SM) are well motivated by being able to provide a unification of the SM gauge couplings at high scales, a solution of the hierarchy problem, and a viable dark matter candidate [1,2]. As one of the low-energy realizations of SUSY, the Minimal Supersymmetric Standard Model (MSSM) [3], carrying the minimal field content consistent with observations, has received over years of continuous attentions; see e.g. Refs. [4,5] for reviews. Despite having many advantages [4,5], the MSSM needs to be extended because of the following two main motivations. The first one is due to the existence of the "µ-problem" in the MSSM [6], where µ is a dimensionful parameter set by hand at the electroweak (EW) scale before spontaneous symmetry breaking. The second one is driven by the recent discovered 125 GeV SM-like Higgs boson [7,8], which has imposed strong constrains on the parameter space of MSSM [9][10][11]. Among the various non-minimal SUSY models, the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [12][13][14][15], the simplest extension of the MSSM with a gauge singlet superfield, not only has the capability to fix these shortcomings of the MSSM, but also alleviates the tension implied by the lack of any evidence for superpartners below the EW scale [16]. Specifically, this model can solve the "µ-problem" of the MSSM elegantly by generating an effective µ-term at the SUSY breaking scale [17,18]. Restrictions on the Higgs sector can also be relaxed in the NMSSM, because the Higgs field can acquire a larger tree-level mass with a low SUSY scale, making accordingly the quantum corrections to the observed 125 GeV Higgs boson small [14,15].
The strength of B s(d) −B s(d) mixing is described by the mass difference ∆M s(d) between the two mass eigenstates of neutral B s(d) mesons, and the uncertainties of the bag parameters and decay constants, which quantify the hadronic matrix elements of local four-quark operators between |B s(d) and |B s(d) states, make up the largest uncertainty by far in the prediction of ∆M s(d) [35,36]. These parameters can be determined either by lattice simulations [37][38][39][40][41][42][43] or with Heavy Quark Effective Theory (HQET) sum rules [44][45][46][47][48], with their results being now compatible in precision with each other [48,49]. Depending on the values of these input parameters as well as the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements [50,51], the SM prediction of ∆M s(d) varies from being consistent with to being larger than the experimental data; see e.g. Refs. [33,52,53] for detailed discussions. In particular, taking the latest lattice averages for these non-perturtive parameters from the Flavour Lattice Averaging Group (FLAG) [49], which are dominated by the values of Fermilab Lattice and MILC (FNAL/MILC) collaboration in 2016 [42], one gets [42,52]  This observation has profound implications for New Physics (NP) models [33,52,53]. As detailed in Refs. [27,33,52,53], this implies particularly that the constrained minimal flavour violating (CMFV) models, in which all flavour violations arise only from the CKM matrix, have difficulties in describing the current data on ∆M s(d) . Thus, in order to reconcile such a tension, one has to resort to the scenarios with non-minimal flavour violation (NMFV), which involve extra sources of flavour-and/or CP-violation and can, therefore, provide potential negative contributions to ∆M s(d) [56][57][58]. This motivates us to investigate whether the Z 3 -invariant NMSSM with NMFV, in which the extra flavour violations arise from the non-diagonal parts of the squark mass matrices related to the soft SUSY breaking terms, can accommodate the observed deviation for ∆M s(d) .
In SUSY models, a physical transition amplitude is more conveniently calculated in the interaction/flavour basis in which all gauge interactions are flavour diagonal and the flavourchanging interactions originate from the off-diagonal entries of the mass matrices in the initial Lagrangian before diagonalization and identification of the physical states, than in the mass eigenstate (ME) basis in which the amplitude is expressed in terms of the physical masses and mixing matrices. This can be achieved with the following two different methods. The first one is based on the well-known diagrammatic technique called the Mass Insertion Approximation (MIA) [59][60][61][62]. Here the diagonal elements of the mass matrices are absorbed into the definition of the (un-physical) massive propagators and the amplitude is, at each loop order, expanded into an infinite series of the off-diagonal elements of the mass matrices, commonly referred to as the mass-insertion (MI) parameters. The second one is based on the Flavour Expansion Theorem (FET) [63], according to which an analytic function about zero of a Hermitian matrix can be expanded polynomially in terms of its off-diagonal elements with coefficients being the divided differences of the analytic function and arguments the diagonal elements of the Hermitian matrix. Being a purely algebraic method, it offers an alternative derivation of the MIA result directly from the amplitude calculated in the ME basis, without performing the tedious and error-prone diagrammatic calculations with MIs in the interaction/flavour basis [63]. Even in the case where there is no clear diagrammatic picture, the FET expansion can still give a consistent MIA result. This method has also been automatized in the package MassToMI [64], facilitating the expansion of an ME amplitude to any user-defined MI order.
In the Z 3 -invariant NMSSM with NMFV, we further assume that the direct mixing between the first two generations of squarks is absent in the mass squared matrices 2 , which is strongly suppressed by the data on K 0 −K 0 mixing and rare kaon decays [30,[68][69][70]. For example, under the constraint from K 0 −K 0 mixing and even with a rough estimate of the hadronic matrix elements, the MI parameters of the first two generations are bounded to be smaller than 0.04 in both the LL and RR sectors [61,62,70]. As we focus only on the flavour mixing effects observed in B s(d) −B s(d) mixing as well as B s → µ + µ − and B → X s γ decays [22][23][24][25][26][27][28][29][30][31][32][33][34], it is natural to assume that the third-generation squarks can mix with the other two generations simultaneously. With such specific squark flavour structures, we shall then adopt the FET procedure [63] where W MSSM µ=0 is the MSSM superpotential but without the µ term [78,79], andŜ denotes the Higgs singlet superfield, The dimensionless parameters λ and κ can be complex in general, but are real in the CP-conserving case. After the scalar component ofŜ gets a non-zero vacuum expectation value (VEV), S = v s / √ 2, the second term in Eq. (2.1) generates an effective µ term, with µ eff = λv s / √ 2, which then solves the "µ-problem" of the MSSM [14,15].
With the scalar components of the Higgs doublet and singlet superfields being denoted by H u , H d , and S, respectively, the soft SUSY breaking Lagrangian of the Z 3 -invariant NMSSM is then given by [76,77] where L MSSM soft µ=0 corresponds to the MSSM part but with the µ-related term removed [78,79]. While the mass parameter m 2 S is real, the trilinear couplings A λ and A κ are generally complex, but are also real in the CP-conserving case, as is assumed throughout this paper.

Flavour structures of the Z 3 -invariant NMSSM
Firstly, we focus on the up-and down-squark mass squared matrices M 2 U and M 2 D , which can be written in their most general 2 × 2-block form as [77] in the so-called super-CKM basis [62]. In the NMFV paradigm [29,30], these two mass matrices are not yet diagonal and can introduce general squark flavour mixings that are usually described by a set of dimensionless parameters δ AB ij , with A, B=L, R referring to the left-and right-handed superpartners of the corresponding quarks and i, j = 1, 2, 3 the generation indices [61].
Throughout this paper, we assume that the third-generation squarks can mix with the other two generations simultaneously, while the direct mixing between the latter two, which is strongly suppressed by the data on Kaon physics [30,[68][69][70], is absent in the squark mass squared matrices. It is also noted that the off-diagonal elements in the LR and RL sectors are severely restricted by the dangerous charge and colour breaking (CCB) minima and unbounded from below directions (UFB) in the effective potential [80][81][82][83]. Here we set all the flavour-violating off-diagonal elements in the LR and RL sectors to be zero, because the CCB and UFB bounds on them are always stronger than the ones from flavour-changing neutral-current processes, and even do not decrease when the SUSY scale increases [80][81][82]. These upper bounds on the flavour-conserving off-diagonal elements are, however, weaker than on the flavour-violating counterparts, which are both proportional to the related running quark masses [83,84]. So we only keep δ LR 33 non-zero and set all the other flavour-conserving off-diagonal elements to be zero, as the upper bounds on the latter are much stronger than on δ LR 33 . In addition, a non-zero δ LR 33 is needed to reproduce the observed 125 GeV SM-like Higgs boson [30,68,69]. All the above observations promote us to consider the following two sets of squark flavour structures: while the flavour-conserving off-diagonal element δ LR 33 is kept in both of these two sectors, only the flavour-violating off-diagonal elements δ LL 23 and δ RR i3 (i = 1, 2) are kept in the LL and RR sectors, respectively. Then, the two mass squared matrices M 2 U and M 2 D in case I, in which the flavour violation resides only in the LL sector, are given, respectively, by have also assumed that the first two generations of squarks are nearly degenerate in mass [4].
In case II with the flavour violation arising only from the RR sector, on the other hand, the two mass squared matrices are given, respectively, by where δ 46 ≡ δ RR 13 and δ 56 ≡ δ RR 23 . Here, for simplicity, we have assumed that all the δ AB ij parameters are real and hence δ AB ij = δ BA ji , due to the hermiticity of the squark mass matrices. The mass matrix for charginos in the interaction basis reads [79] where M 2 is the wino mass, and β = tan −1 (v u /v d ) is the mixing angle of the two Higgs doublets, M P i and the MI parameters δ C ij , δ P ij are defined, respectively, by where i = j and the summation is not applied for the same indices here.
The neutralino mass matrix is given in the basis (B,W 3 ,H 0 d ,H 0 u ,S) T by [14,57] where M 1 is the bino mass, and θ W is the weak mixing angle. Such a mass matrix indicates that the singlinoS couples only to the HiggsinosH 0 d andH 0 u , but not to the gauginosB and W 3 [14]. The squared masses M N i and the MI parameters δ N ij are defined, respectively, by where i = j and the summation is also not applied for the same indices. Diagonalization of Eq. (2.11) is rather involved and has to be in practice performed numerically [14].

FET expansion with different MI order
Before performing the FET expansion, one has to write down the transition amplitude in the ME basis [63]. All the relevant Feynman rules are taken from Refs. [14,15,79,86,87]. Then, the procedure of FET expansion with general/finite MI order includes the following three steps.
Step 1. One transforms the amplitude written in the ME basis into the intermediate result expressed in terms of the blocks L X (i, j), which are defined, respectively, as [64] A for the up and down squarks behaving as scalar fields, for the charginos that behave as Dirac fermions, and for the neutralinos behaving as Majorana fermions. Here l p (m 2 A ) represents symbolically part of the transition amplitude that depends on the mass m A of an internal physical particle in a Feynman diagram, at tree or loop level; for example, l p (m 2 can be a propagator, with q µ being the momentum of the particle A. The unitary transformation matrices Z U , Z D , We use A, B, · · · and i, j, k, · · · to represent the flavour indices of field multiplets in the ME and the interaction/flavour basis, respectively. After applying the transformation rules specified by Eqs. (2.13)-(2.21), one can see that the blocks L X (i, j) depend only on the matrix elements of some functions with arguments being the Hermitian mass squared matrices, and can be expanded as [63] where L X (n; i, j) represents the n-th term in the MI order of the blocks L X (i, j).
Step 2. During our calculation, we also encounter the case in which a Feynman diagram contains two lines involving the same particle. In such a case, the product of two blocks with MI orders specified by m and n, L X (m; i, j)L X (n; i, j), should be firstly combined into a single term with fixed MI order, such as L UU (2n; 3, i; 3, j) = n n 1 =1 L U (2n 1 − 1; 3, i)L U (2n − 2n 1 + 1; 3, j), where L XX (n; i, j; i , j ) denotes the n-th term in the MI order of the block L XX (i, j; i , j ), with (2.23) All the non-zero block terms L X (n; i, j) and L XX (n; i, j; i , j ) for squarks and charginos, with the flavour structures specified in Sec. 2.2, can be easily derived and are listed in the appendix.
As the neutralino mass matrix M χ 0 , given by Eq. (2.11), has many non-zero elements, one can use the following recursive formulas to represent the corresponding blocks , and can be re-expressed in terms of l r N (i x ) (x = 0, 1, · · · , n) using the "divided difference" method [64].
Step 3. One now needs to perform the loop-momentum integration over the products of blocks L X (n; i, j), L XX (n; i, j; i , j ), and l r N (i x ) introduced in the last step. As only one-loop amplitudes are involved throughout this paper, this can be done using iteratively the operation ∂ M SQ Loop X , where M SQ denotes symbolically the squared mass, equaling to M 2 for scalars and to M † M or M M † for fermions, and Loop X is the n-point one-loop integrals in the Passarino-Veltman (PV) basis [88], such as D 2n , C 2n and B 0 introduced in Ref. [22].
Following these three steps, one can then transform an amplitude written initially in the ME basis into an expansion in terms of the MI parameters, up to any user-defined MI order [63,64].

MI-order estimation
Although the FET procedure can provide with us the result expanded to any MI order, an optimal cutting-off should be applied due to the time costing of the programme running. Here we illustrate an efficient MI-order estimation method to get the appropriate order which was usually set by hand (order 2 or 4) in most of recent works [57,[65][66][67].
Taking the block term L UU (2n; 3, 2, 3, 2) in case I, q 2 −M S 2 , as an example, and using the inequality [63] PV (n+1) 0 satisfied by the PV-integrals with vanishing external momenta that are defined by , (2.27) with the assumption that n 3 to avoid divergent integrals, we can obtain where L oth represent the blocks related to the other types of particles, such as L CC . Then, for a given small constant 0 < c 0 < 1, only when indices i 0 and j 0 , the summation over the MI index can be terminated to the n-th order.
In this section, we shall apply the FET procedure with general/finite MI order to B s(d) −B s(d) mixing and B s → µ + µ − decay, within the Z 3 -invariant NMSSM with NMFV.

B s(d) −B s(d) mixing
The strength of B s(d) −B s(d) mixing is described by the mass difference ∆M s(d) , defined by [89] where M q 12 denotes the off-diagonal element in the neutral B q -meson mass matrix, and the effective weak Hamiltonian can be written in a general form as [21] H ∆B=2 Within the Z 3 -invariant NMSSM with NMFV, the following eight operators, as defined in Ref. [22], are all found to be relevant: where α and β are the colour indices, σ µν = 1 2 [γ µ , γ ν ], and P L,R = (1 ∓ γ 5 )/2. To the lowest order in the EW theory, the corresponding Wilson coefficients C i , at the matching scale, of the operators Q i are obtained by evaluating the various one-loop box diagrams mediated by heavy particles appearing in the SM and beyond 3 . Within the SM, only C VLL 1 gets a non-negligible contribution from the one-loop box diagrams with up-type quarks and W bosons circulating in the loops [89], and the next-to-leading-order perturbative QCD corrections to C VLL 1 are also known [92]. In the context of Z 3 -invariant NMSSM with NMFV, on the other hand, all the eight Wilson coefficients C i can get non-zero contributions from the additional one-loop box diagrams mediated by: 1) charged Higgs, up-quarks; 2) chargino, upsquarks; 3) gluinos, down-squarks; 4) neutralinos, down-squarks; 5) mixed gluino, neutralino, down-squarks [22,27,93]. With the aid of the packages FeynArts [94] and FeynCalc [95], all these Feynman diagrams can be calculated and the resulting Wilson coefficients are expressed in terms of the rotation matrices Z U , Z D , Z + χ , Z − χ , and Z χ 0 , as well as the blocks L X (i, j) and L XX (i, j; i , j ). Our results for the Wilson coefficients agree with the ones given in Refs. [27,29,57]. Then, following the procedure detailed in Sec. 2.3, we can transform these Wilson coefficients given in the ME basis into the FET results. Here we have made full use of the hierarchies among the CKM parameters to simplify the final results. For example, when calculating ∆M s in case I, we encounter a term As L UU (i, j), with i = j, does not vanish only when (i, j) = (2, 3) or (3,2), and because of Once the initial conditions for the Wilson coefficients are obtained, we need to include the renormalization group (RG) running effects from the matching scale down to the low-energy scale, at which the hadronic matrix elements are evaluated by the lattice groups [42,49]. All the relevant ingredients for this RG running can be found in Ref. [21]. In this way, we can obtain the final result of the off-diagonal element M q 12 and hence the mass difference ∆M s(d) . For convenience of later discussions, the total contributions to M q 12 are split into the following different parts:

B s → µ + µ − decay
The rare decay B s → µ + µ − proceeds dominantly via the Z-penguin and W -box diagrams within the SM and its branching ratio is highly suppressed [89]. In the context of Z 3 -invariant NMSSM with NMFV, there are in general three types of one-loop diagrams that contribute to this decay, including the various box, the Z-penguin, and the neutral-Higgs-penguin diagrams [22,[96][97][98][99][100][101][102].
where the vector (O V XY ) and scalar (O SXY ) operators are defined, respectively, by The branching ratio of B s → µ + µ − decay is then calculated to be [101,102] where τ Bs is the B s -meson lifetime, and the squared matrix element is given by [100][101][102][103] (4π with the scalar, pseudo-scalar, and axial-vector form factors defined, respectively, by [101] where f Bs is the B s -meson decay constant, and m b(s) denotes the b(s)-quark running mass.
Our main task is then to calculate the Wilson coefficients C V XY and C SXY . Within the SM, only C V LL gets a non-negligible contribution 5 , and we shall use the fitting formula in Eq. (4) of Ref. [106] to get the numerical result for it. The additional contributions to C V XY and C SXY from the Z 3 -invariant NMSSM with NMFV are calculated by ourselves, with the aid of FeynArts and FeynCalc packages. Then, the FET procedure with general/finite MI order is applied for the NMSSM contributions, in exactly the same way as for the B s(d) −B s(d) mixing. 4 Here we need not consider the tensor operators O T X = (bσ µν P X s)(μσ µν µ). While also receiving contributions from these three types of one-loop diagrams in the Z 3 -invariant NMSSM with NMFV, they do not contribute to this process due to the vanishing matrix elements, 0|bσ µν P X s|B s = 0.
5 When the small external momenta and masses are considered, the SM W -box and Z-penguin diagrams also generate non-zero contributions to C SXY , besides the ones from the Higgs-penguin diagrams [104,105].
It should be noted that the branching ratio given by Eq. (3.9) is the so-called "theoretical" branching ratio, which corresponds to the decay time t = 0, while the experimentally measured branching ratio is time-integrated, which is given by [107][108][109] where A ∆Γs is a time-dependent observable [107][108][109], and y s is related to the decay-width difference ∆Γ s between the two B s -meson mass eigenstates, defined by with Γ s L (Γ s H ) denoting the lighter (heavier) eigenstate decay width and Γ s = τ −1 Bs the average decay width of B s meson. In the absence of beyond-SM sources of CP violation, which is assumed throughout this paper, both A ∆Γs and y s will take their respective SM values [107][108][109][110], and the two branching ratios are then related to each other via a simple relation which holds to a very good approximation [106].

Numerical results and discussions
After getting the analytic FET results for the B s(d) −B s(d) mixing and B s → µ + µ − decay, we now proceed to analyze numerically the parameter space of Z 3 -invariant NMSSM with NMFV that is allowed under these experimental constraints.

Choice of input parameters
Firstly, we collect in Table 1    where Mg is the gluino mass. In this paper, we shall consider two sets of fixed parameters that are collected in Table 2. Scenario A is characterized by a large λ and a small tan β, to avoid suppressing the NMSSM-specific tree-level contributions to the SM-like Higgs mass [14,15], while scenario B is featured by negligible NMSSM-specific effects on the SM-like Higgs mass [112]. The remaining parameters M S 2 , δ LL 23 , δ LR 33 , δ RR 13 , and δ RR 23 can vary freely in the ranges considered.
The set of fixed parameters in scenario A is similar to that of the scenario TP3 in Ref. [112], with λ being close to the perturbative limit but still avoiding running into a Landau pole well below the GUT scale [15]. The one in scenario B is, however, featured by a large A λ , which is closely related to the charged-Higgs mass [15]. In scenario A, the predicted lightest and nextto-lightest neutralino masses are given, respectively, by mχ0 1 ∼ 74 GeV and mχ0 2 ∼ 211 GeV, which comply with the limits set by the most recent search for electroweak production of charginos and neutralinos, via the most promising channel pp →χ ± 1χ 0 2 with a 100% branching fraction ofχ 0 2 → H SMχ 0 1 , by the CMS collaboration [113]. In scenario B, on the other hand, mχ0 1 ∼ 179 GeV and mχ0 2 ≈ mχ± 1 ∼ 207 GeV, being compatible with the mass bounds from the same channel but with a 100% branching fraction ofχ 0 2 → Zχ 0 1 [113]. The mass splitting ∆m(χ 0 2 ,χ 0 1 ) = mχ0 2 − mχ0 1 in both of these two scenarios is also compatible with the bound set by the CMS collaboration [113,114]. In addition, the choice µ eff = 200 GeV not only complies with the lower bounds on chargino and neutralino masses, but also ensures that no tree-level fine-tuning is necessary to achieve the EW symmetry breaking [15].
The LHC direct searches have also led to stringent limits on the masses of stops, sbottoms, and gluinos [115][116][117][118][119][120][121][122]. The recent ATLAS result [122] has shown that, for pair produced stops decaying into top quarks, stop masses up to 940 GeV are already excluded in the phenomenological MSSM with a wino-like next-to-lightest supersymmetric particle, while the excluding limit can be up to 860 GeV in scenarios with a Higgsino-like lightest supersymmetric particle.
Obviously, being based on simplified models, these bounds are obtained without considering the most general squark flavour structures, and can be relaxed when taking into account these mixings [123]. Furthermore, there exist possible NMSSM-specific effects with a light singlino in the searches for squarks, which could modify these limits [124][125][126][127][128][129]. Here we simply assume that the mass of the lightest squark is above 940 GeV, and fix the gluino mass at 2100 GeV [130].
Scenario A will give a light singlet scalar with mass around 90 GeV and a 125 GeV SM-like Higgs, while the charged-Higgs mass is around 650 GeV, which may be beneficial for describing the branching ratio of B → X s γ decay [131,132]. The SM-like Higgs predicted in scenario B is, on the other hand, the lightest among the neutral scalars, and the charged Higgs with mass around 2 TeV can make its effect on the B → X s γ decay marginal. Taking together with the values of tan β, we can say that both of these two scenarios make the Higgs-penguin effects negligible for both B s(d) −B s(d) mixing [57] and B s → µ + µ − decay [101].

FET result with optimal MI order 4.2.1 Cutting-off MI-order estimation
We firstly estimate the optimal cutting-off MI order for neutralinos. Here the MI parameters δ N ij depend on the six Z 3 -invariant NMSSM parameters M 1 , M 2 , µ eff , tan β, λ, and κ. Keeping for the moment λ and κ as free variables, but with tan β = 3 and 10 corresponding respectively to the two scenarios defined in Table 2, we find that δ N 15 = δ N 25 = 0 and δ N 12 = −0.007. The dependence of the other MI parameters δ N ij on λ and κ are displayed in Figure 2. From the magnitudes of these δ N ij shown and based on the criterion specified in Sec. 2.4, one can see that the effects of δ N 12 , δ N 13 , and δ N 34 are negligible with the MI order higher than 1, and that of δ N 14 and δ N 23 can be neglected starting from the third MI order (even from the second MI order for δ N 23 in scenario B); while the MI orders higher than 3 should be kept for δ N 24 , δ N 35 , and δ N 45 . When the parameters λ and κ are fixed at the values in scenarios A and B, it is further found that the values of (M χ 0 ) in 3 / M N in , for i n = 1, 2, 3, are much smaller than for i n = 4, 5, and all the terms involving L N (n; 3, i n )(M χ 0 ) in3 can be, therefore, discarded safely for i n = 1, 2, 3.
As the MI parameters for charginos are all less than 0.6 in scenarios A and B, the optimal cutting-off MI order for chargino and squark is determined to be 8 using Eq. (2.29), when the squark MI parameters are less than 0.6 and the given small parameter c 0 is set to be 0.1. in Figure 3 shows the variation of ∆M s with respect to δ LL 23 by considering the cutting-off MI order for squarks from 2 to 8 in case IA, in which the set of fixed parameters in scenario A is used under the case I assumption for the squark flavour structures (and similar definitions apply to the cases IB, IIA, and IIB, to be mentioned below). Here the results with MI order of 2 are similar to what are usually considered in the MIA method. One can see that ∆M s varies with respect to δ LL 23 only slowly for MI order of 2 but decreases obviously for higher MI orders. The results with MI orders of 6 and 8 are nearly identical, which justifies the validity of our MI-order estimation, and hence the cutting-off MI order 8 or 6 is optimal for the considered range of δ LL 23 . In order to show the convergence of the FET results in the case II assumption for squark flavour structures, we consider the ratio which gives the difference between ∆M s by considering the squark MI orders of 6 and 8, respectively. As shown in the right panel of Figure 3, O ∆Ms is nearly 1 in the whole area in case IIA, implying that the convergence has been verified. Similar observations are also made for ∆M s in case IIB and for ∆M d in all the four cases.

MI-order comparison for
For our prediction of the time-integrated branching ratio B(B s → µ + µ − ) in the Z 3 -invariant NMSSM with NMFV, the slepton mass squared matrices are set to be diagonal, with all diagonal elements being given by M S 1 . As an example of convergence checking, we show in Figure 4 the variation of B(B s → µ + µ − ) with respect to δ LL 23 for δ LR 33 = −0.15 in case IA. Being affected by δ RR 13 and δ RR 23 quite weakly [57,140], the convergence of B(B s → µ + µ − ) in case II is not shown here. It can be seen from Figure 4 that the FET result obtained with squark MI order of 2 has no significant deviation from that with higher MI orders, which is obviously different from what is observed in the case of B s(d) −B s(d) mixing.
When the off-diagonal element δ LL 23 is zero, there exists no squark MI contribution and the NMSSM contributions come only from the charged Higgs and the diagonal part of M 2 U . Even in this case, B(B s → µ + µ − ) can be obviously enhanced, putting it to be larger than the LHCb measurement, (3.0 ± 0.6 +0. 3 −0.2 ) × 10 −9 [141]. As a reference, our prediction within the SM is (3.6 ± 0.3) × 10 −9 , obtained by using the fitting formula in Eq. (4) of Ref. [106] with the updated input parameters listed in Table 1 and the decay constant f Bs from Ref. [54].

FET vs. mass diagonalization
From the previous analyses, we can see that, in some regions of the NMSSM parameter space, the usually adopted MIA is not adequate by considering only the first one or two MI orders, and higher orders in the MI expansion must be considered. It is also shown that the FET results with optimal cutting-off MI orders do demonstrate good convergence. In this subsection, we show that these results also agree well with the ones calculated in the ME basis with exact diagonalization of the mass matrices that is achievable only numerically [63,67].
As an example, we show in Figure 5 the NMSSM contributions to Re(2M   [141], is in reasonable agreement with the SM prediction, (3.6 ± 0.3) × 10 −9 . This will impose much more stringent constraints on NP [110].

(4.4)
Here the photon energy cutoff E γ > 1.6 GeV is defined in the decaying meson rest frame.
The charged-Higgs contribution in the NMSSM (or MSSM) belongs to the case of Model-II considered in Ref. [132], and always interferes with the SM one in a constructive manner. Then, the extra one-loop SUSY contributions should involve a cancellation with that from the charged Higgs, so as to comply with the B → X s γ constraint [22,70]. As the SM and charged-Higgs contributions to B dγ , the branching ratio of B → X d γ decay, are both suppressed by the CKM factor |K 31 /K 32 | 2 with respect to B sγ , while the contributions from the squark MI parameters are not affected by this factor, it is expected that δ LL 13 in M 2 U will be strongly constrained by B dγ . As a result, we have set δ LL 13 to be zero in Sec. 2.2 from the very beginning. In this work, we use the public code SUSY FLAVOR [73][74][75] to get the NP contribution to B sγ numerically, by adapting the initial values of input parameters to our case; especially, one must change the related parameters in the file to make sure that the photon energy cutoff E γ > 1.6 GeV.
In the following, we shall exploit the 95% C.L. bounds from ∆M s , ∆M d , B(B s → µ + µ − ), and B sγ , to set the allowed regions for the NMSSM parameters M S 2 , δ LL 23 , δ LR 33 , δ RR 13 , and δ RR 23 , during which only the experimental and the SM theoretical uncertainties are considered. Explicitly, we firstly construct an allowed range for each observable by taking into account both the experimental data and the SM prediction with their respective 1.96σ uncertainties added in quadrature, and then require the NMSSM central values to lie within the range, to get the allowed regions for the NMSSM parameters.

Results in scenario A
Firstly, we show in Figure 6 the allowed regions for the squark MI parameters δ LL 23 , δ LR 33 , δ RR 13 and δ RR 23 in scenario A. Here, to match the SM-like Higgs mass, as shown in Figure 1, three choices of M S 2 from 1100 2 to 1500 2 GeV 2 are made in both cases IA and IIA, and the squark MI parameter δ LR 33 is set to be −0.4 in case IIA. The scenario A is featured by the observation that the charged-Higgs contribution to ∆M s is positive and large for small charged-Higgs mass and tan β, while its contribution to B sγ , being also positive and large for small charged-Higgs mass, is less sensitive to the variation of tan β. The bound from B(B s → µ + µ − ) is, however, not shown in Figure 6, because on the one hand this bound is satisfied in the whole parameter regions displayed in Figure 6 for case I, and on the other hand B(B s → µ + µ − ) is nearly not affected by δ RR 13 or δ RR 23 in case II [57,140]. One can see clearly that the excluded area by the lower bound of 940 GeV on squark masses based on a rough estimate of the LHC direct searches for squarks [115][116][117][118][119][120][121][122] becomes reduced with increasing M S 2 .
The allowed region for δ LL 23 from ∆M s also becomes reduced when M S 2 increases, and only the one with large magnitudes of δ LL 23 exists in case IA. It is particularly observed that, for M S 2 = 1100 2 GeV 2 , there exists no allowed region for δ LL 23 and δ LR 33 in this case. In case IIA, on the other hand, the bounds from ∆M s and ∆M d become stronger for a larger M S 2 , and the one from ∆M s is so strong that it is no longer compatible with that from the lower bound of 940 GeV on squark masses. The bound from B sγ also shows that a larger magnitude of δ LL 23 or δ LR 33 is required in case IA, but almost the whole area of δ RR 13 and δ RR 23 is allowed in  case IIA. All these features can be explained by the following observations: in scenario A, as mentioned before, there exist considerably large and positive contributions to ∆M s and B sγ from the charged Higgs; then a large magnitude of δ LL 23 or δ RR 23 (δ LR 33 ) is needed to provide large but negative contributions to ∆M s (B sγ ), so as to cancel the charged-Higgs effects [22,70].
Especially in case IIA, the chosen value of −0.4 for δ LR 33 is already large enough to cancel the charged-Higgs effect on B sγ , and the flavour-violating contribution from the RR sector is small, making the whole area of δ RR 13 and δ RR 23 being allowed by this observable. Compared to the case from ∆M s , the allowed regions for the squark MI parameters from ∆M d are relatively large, because the charged-Higgs contribution to ∆M d is suppressed by the CKM factors in both cases IA and IIA, with some areas being not allowed in case IIA due to the effect of δ RR 13 . After taking into account the 95% C.L. bounds from ∆M s , ∆M d , B sγ , B(B s → µ + µ − ), as well as the SM-like Higgs mass, we find numerically that the squark MI parameters δ LL 23 > 0.45 and |δ LR 33 | ∼ 0.15, and the allowed region is severely small in case IA. In case IIA, on the other hand, there exists no allowed region for δ RR 13 and δ RR 23 .

Results in scenario B
In Figure 7, we show the allowed regions for δ LL 23 , δ LR 33 , δ RR 13 and δ RR 23 in scenario B. To match the SM-like Higgs mass, as shown in Figure 1, the squark MI parameter δ LR 33 is set to be 0.2 in case IIB, and three choices of M S 2 from 1700 2 (1800 2 ) to 1900 2 (2000 2 ) GeV 2 are made in case IB (IIB). The scenario B is featured by the observation that the positive contributions to ∆M s and B sγ with a large charged-Higgs mass are small, and considerably negative contributions from squarks and gauginos are not needed to cancel the charged-Higgs effects. As already observed in scenario A, the excluded area by the lower bound of 940 GeV on squark masses also becomes reduced with increasing M S 2 , and the bound from B(B s → µ + µ − ) is satisfied in the whole parameter regions shown here.
In case IB, the allowed region for the squark MI parameters from ∆M s(d) is large because of the small charged-Higgs contribution, and is almost independent of M S 2 ; while the bound from B sγ indicates that a smaller δ LL 23 is favored. In case IIB, on the other hand, while the allowed region from ∆M d stays nearly the same, the one from ∆M s shrinks slowly when M S 2 increases.
Under the bound of B sγ , an allowed region with a positive δ RR 23 begins to emerge only when M S 2 is about 1900 2 GeV 2 .
After considering the 95% C.L. bounds from ∆M s , ∆M d , B sγ , as well as the SM-like Higgs  Figure 6. mass, the parameter δ LL 23 can be small, and the allowed region is relatively larger in case IB compared to the one in case IA. In case IIB, on the other hand, the allowed area of δ RR 13 and δ RR 23 exists only when M S 2 is larger than about 1900 2 GeV 2 .

Conclusion
In this paper, motivated by the observation that the SM predictions are now above the experimental data for the mass difference ∆M s(d) , and the usual CMFV models have difficulties in reconciling this discrepancy, we have investigated whether the Z 3 -invariant NMSSM with NMFV, in which the extra flavour violations arise from the non-diagonal parts of the squark mass squared matrices, can accommodate such a deviation, while complying with the experimental constraints from the branching ratios of B s → µ + µ − and B → X s γ decays.
We have calculated the NMSSM contributions to ∆M s(d) and the branching ratio B(B s → µ + µ − ), using the recently developed FET procedure, which allows to perform a purely algebraic MI expansion of a transition amplitude written in the ME basis without performing the tedious and error-prone diagrammatic calculations in the interaction/flavour basis. Specifically, we have considered the finite MI orders for neutralinos but the general MI orders for squarks and charginos, under the following two sets of assumptions for the squark flavour structures: while the flavour-conserving off-diagonal element δ LR 33 is kept in both of these two sectors, only the flavour-violating off-diagonal elements δ LL 23 and δ RR i3 (i = 1, 2) are kept in the LL and RR sectors, respectively. In this way, our analytic results are polynomials with the MI parameters as the variables and are expressed directly in terms of the initial Lagrangian parameters in the interaction/flavour basis, making it easy to impose the experimental bounds on them and, at the same time, allowing for a more transparent understanding of the qualitative behaviour of the transition amplitude. We have also presented an efficient method to estimate the optimal cutting-off MI orders for neutralinos, charginos, and squarks.
For the numerical analyses, we have considered two sets of NMSSM parameters that are denoted, respectively, by scenarios A and B in Table 2. Both of them can match the SM-like Higgs mass, and also make the Higgs-penguin effects negligible for B s(d) −B s(d) mixing and B s → µ + µ − decay. Together with the two assumptions for the squark flavour structures, there are totally four different cases, IA, IB, IIA, and IIB, to be discussed. Firstly, after getting the optimal cutting-off MI orders for neutralinos, charginos, and squarks with our estimation rules, we have verified the convergence of the FET results obtained with the corresponding MI orders, even in the case when the MI parameters are large. Then, taking Re(2M s 12 ) in cases IA and IIA as an example, we have demonstrated that the FET results with optimal cutting-off MI orders agree quite well with the ones calculated directly in the ME basis with an exact diagonalization of the mass matrices that is usually achievable only numerically. This also proves the necessity to consider the optimal cutting-off MI orders when performing an NMFV expansion, as the expansion at leading order is usually insufficient and could easily be misleading.
Finally, after considering the 95% C.L. bounds from the observables ∆M s , ∆M d , B(B s → µ + µ − ), B sγ , as well as the SM-like Higgs mass, we have discussed the allowed regions for the parameters M S 2 , δ LL 23 , δ LR 33 , δ RR 13 , and δ RR 23 . It is found that only large values of δ LL 23 , with |δ LR 33 | ∼ 0.15, are allowed in case IA, and the allowed region for δ LL 23 becomes reduced when M S 2 increases from 1100 2 to 1500 2 GeV 2 . In case IB, on the other hand, with δ LR 33 being fixed at about 0.2, relatively smaller magnitude of δ LL 23 is found to be allowed, and the allowed region for δ LL 23 becomes almost independent of M S 2 . In case IIA, with δ LR 33 being fixed at −0.4, there exists no allowed region for δ RR 13 and δ RR 23 , because of the strong bound from ∆M s . In case IIB, on the contrary, with δ LR 33 being fixed at about 0.2, the allowed region for δ RR 13 and δ RR 23 exists only when M S 2 is larger than about 1900 2 GeV 2 .
As a final remark, we should mention that, with the experimental progress in direct searches for SUSY particles as well as the more and more precise theoretical predictions for these observables, the NMSSM effects on these low-energy flavour processes can be further exploited.

Appendix: Block terms for squarks and charginos
In this appendix, we list all the non-zero block terms for squarks and charginos. For convenience, we introduce the notations ∆ 1 ≡ in case II. In the following, n = 1, 2, 3, · · · , denotes the MI-order index.
For up-type squarks, the non-zero block terms are given as in case II.