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

TL;DR

This paper investigates B meson mixing and decay processes within the NMSSM framework using a novel algebraic approach based on the Flavour Expansion Theorem, providing insights into experimental deviations and constraints.

Contribution

It introduces a purely algebraic mass-insertion expansion method for analyzing flavor processes in the NMSSM, avoiding traditional diagrammatic calculations.

Findings

NMSSM can explain deviations in $ riangle M_{s(d)}$

Results comply with experimental bounds on $B_s o \mu^+ u$ and $B o X_s\gamma$

Analytic expressions relate model parameters to observables.

Abstract

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.

Figures

Tables

Table 1. Table 1: Summary of part of the input parameters used throughout this paper.

QCD and EW parameters [54]
$G_F\left[{10}^{-5}{\text{GeV}}^{-2}\right]$	${\alpha }_s\left(M_Z\right)$	$M_W\left[\text{GeV}\right]$	${\sin }^2{\theta }_W$
$1.1663787$	$0.1181\left(11\right)$	$80.379$	$0.2312$
Quark masses [GeV] [54]
${\overline{m}}_b\left({\overline{m}}_b\right)$	${\overline{m}}_c\left({\overline{m}}_c\right)$	$m_t$
${4.18}_{-0.03}^{+0.04}$	${1.275}_{-0.035}^{+0.025}$	$173.0\left(4\right)$
$B$-meson parameters [55]
$M_{B_d}\left[\text{GeV}\right]$	$M_{B_s}\left[\text{GeV}\right]$	$1/{\Gamma }_H^s\left[\text{ps}\right]$
$5.280$	$5.367$	$1.609\left(10\right)$
CKM parameters [111]
${\lambda }_{\text{CKM}}$	$A$	$\overline{\rho }$	$\overline{\eta }$
$0.2251\left(4\right)$	${0.831}_{-0.031}^{+0.021}$	$0.155\left(8\right)$	$0.340\left(10\right)$

Table 2. Table 2: Two set of fixed parameters, all being defined at the scale 1 ​ TeV 1 TeV 1~{}\rm{TeV} , for the ℤ 3 subscript ℤ 3 \mathbb{Z}_{3} -invariant NMSSM with NMFV. They are given in units of “GeV” except for tan ⁡ β 𝛽 \tan\beta , λ 𝜆 \lambda , and κ 𝜅 \kappa .

	$M_1$	$M_2$	$M_{\tilde{g}}$	$\sqrt{M_S^{}{}_1{}^{}}$	${\mu }_{\text{eff}}$	$\tan \beta$	$A_{\lambda }$	$A_{\kappa }$	$\lambda$	$\kappa$
Scenario A	$500$	$1000$	$2100$	$1600$	$200$	$3$	$650$	$-10$	$0.67$	$0.1$
Scenario B	$500$	$1000$	$2100$	$1600$	$200$	$10$	$2000$	$-100$	$0.3$	$0.2$

Equations

$$\Delta M_s^{\text{SM}}=20.01\pm 1.25{\text{ps}}^{-1},\Delta M_d^{\text{SM}}=0.630\pm 0.069{\text{ps}}^{-1},$$

$$\Delta M_s^{\text{exp}}=17.757\pm 0.021{\text{ps}}^{-1},\Delta M_d^{\text{exp}}=0.5064\pm 0.0019{\text{ps}}^{-1}.$$

$$\displaystyle W_{\text{NMSSM}}=W_{\text{MSSM}}\big{|}_{\mu=0}+\lambda\,\hat{S}\,\hat{H}_{u}\cdot\hat{H}_{d}+\frac{1}{3}\,\kappa\,\hat{S}^{3},$$

$$\displaystyle-\mathcal{L}^{\text{NMSSM}}_{\text{soft}}=-\mathcal{L}^{\text{MSSM}}_{\text{soft}}\big{|}_{\mu=0}+m^{2}_{S}\,|S|^{2}+\left(\lambda\,A_{\lambda}\,S\,H_{u}\cdot H_{d}+\frac{1}{3}\,\kappa\,A_{\kappa}\,S^{3}+\text{h.c.}\right),$$

$$M_{\tilde{q}}^2=\left(\begin{array}{cc}{M}_{\tilde{q},\text{LL}}^2& M_{\tilde{q},\text{LR}}^2\\ M_{\tilde{q},\text{RL}}^2& M_{\tilde{q},\text{RR}}^2\end{array}\right),\tilde{q}=\tilde{U},\tilde{D},$$

$$M_{\tilde{U}}^2$$

$$M_{\tilde{D}}^2$$

$$M_{\tilde{U}}^2$$

$$M_{\tilde{D}}^2$$

$$M_{\chi }=\left(\begin{array}{cc}{M}_2& \sqrt{2}{m}_W\sin \beta \\ \sqrt{2}{m}_W\cos \beta & {\mu }_{\text{eff}}\end{array}\right),$$

$${M_C}_i$$

$${\delta }_{ij}^C$$

$$M_{\chi_{0}}=\left(\begin{array}[]{ccccc}M_{1}&0&-\frac{ev_{d}}{2\cos\theta_{W}}&\frac{ev_{u}}{{2}\cos\theta_{W}}&0\\ 0&M_{2}&\frac{ev_{d}}{2\sin\theta_{W}}&-\frac{ev_{u}}{2\sin\theta_{W}}&0\\ -\frac{ev_{d}}{2\cos\theta_{W}}&\frac{ev_{d}}{2\sin\theta_{W}}&0&-\mu_{\text{eff}}&-\frac{\lambda v_{u}}{\sqrt{2}}\\ \frac{ev_{u}}{{2}\cos\theta_{W}}&-\frac{ev_{u}}{2\sin\theta_{W}}&-\mu_{\text{eff}}&0&-\frac{\lambda v_{d}}{\sqrt{2}}\\ 0&0&-\frac{\lambda v_{u}}{\sqrt{2}}&-\frac{\lambda v_{d}}{\sqrt{2}}&{\sqrt{2}}\kappa v_{s}\end{array}\right),$$

$${M_N}_i=(M_{{\chi }_0}^{†}{M}_{{\chi }_0}{)}_{ii},{\delta }_{ij}^N=\frac{(M_{{\chi }_0}^{†}{M}_{{\chi }_0}{)}_{ij}}{\sqrt{{M_N}_i{M_N}_j}},$$

$$\sum _A(Z_U{)}_{iA}{l}_p(m_A^2)(Z_U{)}_{jA}^{\ast }$$

$$\sum _A(Z_D{)}_{iA}{l}_p(m_A^2)(Z_D{)}_{jA}^{\ast }$$

$$\sum _A(Z_{\chi }^{+}{)}_{iA}{l}_p(m_A^2)(Z_{\chi }^{+}{)}_{jA}^{\ast }$$

$$\sum _A(Z_{\chi }^{-}{)}_{iA}{l}_p(m_A^2)(Z_{\chi }^{-}{)}_{jA}^{\ast }$$

$$\sum _A(Z_{\chi }^{-}{)}_{iA}{m}_A{l}_p(m_A^2)(Z_{\chi }^{+}{)}_{jA}^{\ast }$$

$$\sum _A(Z_{\chi }^{+}{)}_{iA}{m}_A{l}_p(m_A^2)(Z_{\chi }^{-}{)}_{jA}^{\ast }$$

$$\sum _A(Z_{{\chi }_0}{)}_{iA}{l}_p(m_A^2)(Z_{{\chi }_0}{)}_{jA}^{\ast }$$

$$\sum _A(Z_{{\chi }_0}{)}_{iA}{m}_A{l}_p(m_A^2)(Z_{{\chi }_0}{)}_{jA}$$

$$\sum _A(Z_{{\chi }_0}{)}_{iA}^{\ast }{m}_A{l}_p(m_A^2)(Z_{{\chi }_0}{)}_{jA}^{\ast }$$

$$L_{\text{X}}(i,j)=\sum _{n=0}^{\infty }{L}_{\text{X}}(n;i,j),$$

$$L_{\text{XX}}(i,j;i^{\prime },j^{\prime })=\sum _{n=0}^{\infty }{L}_{\text{XX}}(n;i,j;i^{\prime },j^{\prime }).$$

$$L_{\text{N}}(n;i_0,i_n)=$$

$$\times \left({\delta }_{i_1{i}_2}^N\sqrt{{M_N}_{i_1}{M_N}_{i_2}}\right)\cdots \left({\delta }_{i_{n-1}{i}_n}^N\sqrt{{M_N}_{i_{n-1}}{M_N}_n}\right),$$

$$L_{\text{UU}}\left(2n;3,2;3,2\right)=\frac{n{M_S}_1{M_S}_2{\delta }_{23}^2}{(q^2-{M_S}_1{)}^2(q^2-{M_S}_2{)}^{n+1}}{\Delta }_1^{n-1},$$

$$\left|{\text{PV}}_0^{(n+1)}(m_1^2,m_2^2,\cdots ,m_n^2,m_{n+1}^2)\right|⩽\frac{1}{m_{n+1}^2}\left|{\text{PV}}_0^{(n)}(m_1^2,m_2^2,\cdots ,m_n^2)\right|,$$

$${\text{PV}}_0^{(n)}(m_1^2,m_2^2,\cdots ,m_n^2)=-i(4\pi {)}^2\int \frac{d^{4}q}{(2\pi {)}^4}\frac{1}{{\Pi }_{j=1}^n(q^2-m_j^2)},$$

$$\left|\int d^{4}q{L}_{\text{UU}}(2n;3,2;3,2)L_{\text{oth}}\right|⩽n{\left({\delta }_{23}^2+{\delta }_{36}^2\right)}^{n-1}\left|\int d^{4}q{L}_{\text{UU}}(2;3,2;3,2)L_{\text{oth}}\right|,$$

$$(n+1){\left({\delta }_{23}^2+{\delta }_{36}^2\right)}^n<c_0,$$

$$\left|\int d^{4}q{L}_{\text{N}}(n;i_0,j_0)L_{\text{oth}}\right|⩽\sum _{i_1,i_2,\cdots ,i_{n-1}}\left|\int d^{4}q{l}_{\text{N}}^r(i_0)L_{\text{oth}}\right|\left|{\delta }_{i_0{i}_1}^N{\delta }_{i_1{i}_2}^N\cdots {\delta }_{i_{n-1}{j}_0}^N\right|\sqrt{\frac{{M_N}_{i_0}}{{M_N}_{j_0}}},$$

$$\left|\int d^{4}q{L}_{\text{N}}(n;i_0,i_n){(M_{{\chi }_0})}_{i_n{j}_0}{L}_{\text{oth}}\right|⩽\sum _{i_1,i_2,\cdots ,i_n}\left|\int d^{4}q{l}_N^r(i_0)L_{\text{oth}}\right|$$

$$\times \left|{\delta }_{i_0{i}_1}^N{\delta }_{i_1{i}_2}^N\cdots {\delta }_{i_{n-1}{i}_n}^N{(M_{{\chi }_0})}_{i_n{j}_0}\right|\sqrt{\frac{{M_N}_{i_0}}{{M_N}_{i_n}}}.$$

$$\Delta M_q=2|M_{12}^q|=2|⟨B_q|{\mathcal{H}}_{\mathrm{eff}}^{\Delta B=2}|{\bar{B}}_q⟩|,q=s,d,$$

Full text

1 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 “$\mu$-problem” in the MSSM [6], where $\mu$ 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 “$\mu$-problem” of the MSSM elegantly by generating an effective $\mu$-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].

Along with the dedicated direct searches at the Large Hadron Collider (LHC) for SUSY particles [16], it is also interesting and even complementary to investigate the virtual effects of these hypothesized particles on the low-energy processes, such as the neutral-meson mixings as well as the CP violation and rare decays of various hadrons [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. In this paper, we shall focus on the $\mathbb{Z}_{3}$-invariant NMSSM [14, 15], a simpler scenario of NMSSM with a scale-invariant superpotential, and study its effects on $B_{s(d)}-\bar{B}_{s(d)}$ mixing, $B_{s}\to\mu^{+}\mu^{-}$ and $B\to X_{s}\gamma$ decays.

The strength of $B_{s(d)}-\bar{B}_{s(d)}$ mixing is described by the mass difference $\Delta 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)}\rangle$ and $|\bar{B}_{s(d)}\rangle$ states, make up the largest uncertainty by far in the prediction of $\Delta 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 [49, 48]. Depending on the values of these input parameters as well as the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements [50, 51], the SM prediction of $\Delta 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 [52, 42]555While the 2016 FNAL/MILC derived values of the bag parameters still need an independent cross-check by other lattice groups or by other methods, we shall adopt in this paper its direct results for the matrix elements, i.e. $f_{B_{s(d)}}\sqrt{B_{s(d)}}$, which profit from the full set of correlations among them and are, therefore, more suitable for phenomenological applications [42].

[TABLE]

which are about $1.8~{}\sigma$ above their respective experimental values [54, 55]

[TABLE]

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 $\Delta 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 $\Delta M_{s(d)}$ [56, 57, 58]. This motivates us to investigate whether the $\mathbb{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 $\Delta 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 flavour-changing 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. See e.g. Refs. [63, 65, 66, 57, 67] for recent applications of this method.

In the $\mathbb{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 matrices666Note that, even in such a case, the simultaneous presence of non-zero $1\text{-}3$ and $2\text{-}3$ mixings can induce the $1\text{-}2$ mixing starting at the second order in the MI expansion. This effect is, however, too small to violate the constraint from Kaon physics [30, 68, 69, 70]., which is strongly suppressed by the data on $K^{0}-\bar{K}^{0}$ mixing and rare kaon decays [30, 68, 69, 70]. For example, under the constraint from $K^{0}-\bar{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 [70, 62, 61]. As we focus only on the flavour mixing effects observed in $B_{s(d)}-\bar{B}_{s(d)}$ mixing as well as $B_{s}\to\mu^{+}\mu^{-}$ and $B\to X_{s}\gamma$ 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] to calculate the mass difference $\Delta M_{s(d)}$ and the branching ratio of $B_{s}\to\mu^{+}\mu^{-}$ decay. It should be noted that the mass matrices of squarks (see Eqs. (2.4)-(2.7)) and charginos (see Eq. (2.8)) can be diagonalized analytically, and only the neutralino mass matrix (see Eq. (2.11)), being a $5\times 5$ matrix, has to be done numerically. For the squarks and charginos sectors, both the full FET approach and the calculation directly in the ME basis can be used to get an analytical result of the observables considered. But for the neutralino sector, an analytical result can be obtained only in the FET approach. As a purely algebraic method, the FET approach allows for a systematic expansion of a transition amplitude in terms of the MI parameters. Thus, the result obtained is a polynomial with the MI parameters as the variables, and hence allows for a more transparent understanding of the qualitative behaviour of the transition amplitude. In addition, being expressed directly in terms of the initial Lagrangian parameters in the interaction/flavour basis, they can be conveniently and easily exploited to put the experimental bounds on the parameters of the initial Lagrangian in this way. These features are, however, lost in the result obtained directly in the ME basis, in which the MI parameters are usually contained in the denominators and the arguments of various logarithms. On the other hand, one should note that the FET results agree very well with the ones calculated directly in the ME basis, as will be demonstrated in Sec. 4.2.4. So, for consistency, we shall adopt the FET approach throughout this paper. Specifically, we consider the general MI orders for squarks and charginos but the finite MI orders for neutralinos. Although the general MI orders have been considered in Refs. [71, 72], only one kind of MI parameter is kept in the whole “fat propagators”. In our case, however, there exist two kinds of MI parameters in each line and the mixed arrangement of them is required. For concreteness, we call our procedure the FET expansion with different MI order and, by checking if the FET results agree with the ones calculated directly in the ME basis, test our estimation for the optimal cutting-off MI orders. For the branching ratio of $B\to X_{s}\gamma$ decay, the public code SUSY_FLAVOR [73, 74, 75] is used.

Our paper is organized as follows. In Sec. 2, after specifying the flavour structures assumed throughout this paper, we introduce the FET procedure with different MI order, which is then used to calculate the $B_{s(d)}-\bar{B}_{s(d)}$ mixing and $B_{s}\to\mu^{+}\mu^{-}$ decay in Sec. 3. Detailed numerical results and discussions are then presented in Sec. 4. Our conclusions are finally made in Sec. 5. For convenience, the block terms of squarks and charginos are listed in the appendix.

2 $\mathbb{Z}_{3}$-invariant NMSSM and the FET procedure

2.1 Lagrangian of the $\mathbb{Z}_{3}$-invariant NMSSM

At the Lagrangian level, the $\mathbb{Z}_{3}$-invariant NMSSM differs from the MSSM by the superpotential and the soft SUSY breaking part. The scale-invariant superpotential of NMSSM reads [76, 77]

[TABLE]

where W_{\text{MSSM}}\big{|}_{\mu=0} is the MSSM superpotential but without the $\mu$ term [78, 79], and $\hat{S}$ denotes the Higgs singlet superfield, while $\hat{H}_{u}=(\hat{H}_{u}^{+},\hat{H}_{u}^{0})^{T}$ and $\hat{H}_{d}=(\hat{H}_{d}^{0},\hat{H}_{d}^{-})^{T}$ are the two Higgs doublet superfields, with the convention $\hat{H}_{u}\cdot\hat{H}_{d}=\hat{H}_{u}^{+}\hat{H}_{d}^{-}-\hat{H}_{u}^{0}\hat{H}_{d}^{0}$. The dimensionless parameters $\lambda$ and $\kappa$ can be complex in general, but are real in the CP-conserving case. After the scalar component of $\hat{S}$ gets a non-zero vacuum expectation value (VEV), $\langle S\rangle=v_{s}/\sqrt{2}$, the second term in Eq. (2.1) generates an effective $\mu$ term, with $\mu_{\text{eff}}=\lambda v_{s}/\sqrt{2}$, which then solves the “$\mu$-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 $\mathbb{Z}_{3}$-invariant NMSSM is then given by [76, 77]

[TABLE]

where \mathcal{L}^{\text{MSSM}}_{\text{soft}}\big{|}_{\mu=0} corresponds to the MSSM part but with the $\mu$-related term removed [78, 79]. While the mass parameter $m^{2}_{S}$ is real, the trilinear couplings $A_{\lambda}$ and $A_{\kappa}$ are generally complex, but are also real in the CP-conserving case, as is assumed throughout this paper.

2.2 Flavour structures of the $\mathbb{Z}_{3}$-invariant NMSSM

Firstly, we focus on the up- and down-squark mass squared matrices $M_{\tilde{U}}^{2}$ and $M_{\tilde{D}}^{2}$, which can be written in their most general $2\times 2$-block form as [77]

[TABLE]

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 $\delta_{ij}^{\text{AB}}$, 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 $\delta_{33}^{\text{LR}}$ 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 $\delta_{33}^{\text{LR}}$. In addition, a non-zero $\delta_{33}^{\text{LR}}$ 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 $\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 LL and RR sectors, respectively. Then, the two mass squared matrices $M_{\tilde{U}}^{2}$ and $M_{\tilde{D}}^{2}$ in case I, in which the flavour violation resides only in the LL sector, are given, respectively, by

[TABLE]

where $\delta_{23}\equiv\delta_{23}^{\text{LL}}$ and $\delta_{36}\equiv\delta_{33}^{\text{LR}}$. In the LL sectors, which satisfy the relation $(M_{\tilde{D}}^{2})_{\text{LL}}=K^{\dagger}(M_{\tilde{U}}^{2})_{\text{LL}}K$ (with $K$ being the CKM matrix) due to the $SU(2)_{L}$ gauge invariance [62], we have neglected safely the ${\cal O}(\lambda_{\text{CKM}}^{2})$ terms in Eq. (2.5) (and also in Eq. (2.7)), where $\lambda_{\text{CKM}}=|V_{us}|$ is the expansion parameter in the Wolfenstein parameterization [85] of the CKM matrix. Here we 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

[TABLE]

where $\delta_{46}\equiv\delta_{13}^{\text{RR}}$ and $\delta_{56}\equiv\delta_{23}^{\text{RR}}$. Here, for simplicity, we have assumed that all the $\delta_{ij}^{\text{AB}}$ parameters are real and hence $\delta_{ij}^{\text{AB}}=\delta_{ji}^{\text{BA}}$, due to the hermiticity of the squark mass matrices.

The mass matrix for charginos in the interaction basis reads [79]

[TABLE]

where $M_{2}$ is the wino mass, and $\beta=\tan^{-1}(v_{u}/v_{d})$ is the mixing angle of the two Higgs doublets, defined in terms of their VEVs $v_{u}=\sqrt{2}\,\langle H_{u}\rangle$ and $v_{d}=\sqrt{2}\,\langle H_{d}\rangle$. The squared masses ${M_{C}}_{i}$, ${M_{P}}_{i}$ and the MI parameters $\delta^{C}_{ij}$, $\delta^{P}_{ij}$ are defined, respectively, by

[TABLE]

where $i\neq j$ and the summation is not applied for the same indices here.

The neutralino mass matrix is given in the basis $(\tilde{B}$, $\tilde{W}^{3}$, $\tilde{H}_{d}^{0}$, $\tilde{H}_{u}^{0}$, $\tilde{S})^{T}$ by [57, 14]

[TABLE]

where $M_{1}$ is the bino mass, and $\theta_{W}$ is the weak mixing angle. Such a mass matrix indicates that the singlino $\tilde{S}$ couples only to the Higgsinos $\tilde{H}_{d}^{0}$ and $\tilde{H}_{u}^{0}$, but not to the gauginos $\tilde{B}$ and $\tilde{W}^{3}$ [14]. The squared masses ${M_{N}}_{i}$ and the MI parameters $\delta^{N}_{ij}$ are defined, respectively, by

[TABLE]

where $i\neq 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].

2.3 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. [79, 15, 14, 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_{\text{X}}(i,j)$, which are defined, respectively, as [64]

[TABLE]

for the up and down squarks behaving as scalar fields,

[TABLE]

for the charginos that behave as Dirac fermions, and

[TABLE]

for the neutralinos behaving as Majorana fermions. Here $l_{p}(m_{A}^{2})$ 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_{A}^{2})=1/(q^{2}-m_{A}^{2})$ can be a propagator, with $q^{\mu}$ being the momentum of the particle $A$. The unitary transformation matrices $Z_{U}$, $Z_{D}$, $Z_{\chi}^{+}$, $Z_{\chi}^{-}$, and $Z_{\chi_{0}}$ are introduced to diagonalize the Hermitian mass squared matrices $M_{\tilde{U}}^{2}$, $M_{\tilde{D}}^{2}$, $M_{\chi}^{\dagger}M_{\chi}$, $M_{\chi}M_{\chi}^{\dagger}$, and $M_{\chi_{0}}^{\dagger}M_{\chi_{0}}$, respectively. We use $A,B,\cdots$ and $i,j,k,\cdots$ 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_{\text{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]

[TABLE]

where $L_{\text{X}}(n;i,j)$ represents the $n$-th term in the MI order of the blocks $L_{\text{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_{\text{X}}(m;i,j)L_{\text{X}}(n;i,j)$, should be firstly combined into a single term with fixed MI order, such as $L_{\text{UU}}(2n;3,i;3,j)=\sum_{n_{1}=1}^{n}L_{\text{U}}(2n_{1}-1;3,i)L_{\text{U}}(2n-2n_{1}+1;3,j)$, where $L_{\text{XX}}(n;i,j;i^{\prime},j^{\prime})$ denotes the $n$-th term in the MI order of the block $L_{\text{XX}}(i,j;i^{\prime},j^{\prime})$, with

[TABLE]

All the non-zero block terms $L_{\text{X}}(n;i,j)$ and $L_{\text{XX}}(n;i,j;i^{\prime},j^{\prime})$ 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_{\chi_{0}}$, given by Eq. (2.11), has many non-zero elements, one can use the following recursive formulas to represent the corresponding blocks

[TABLE]

where $l^{r}_{\text{N}}(i_{0},i_{1},\cdots,i_{n})=\frac{1}{(q^{2}-{M_{N}}_{i_{0}})(q^{2}-{M_{N}}_{i_{1}})\cdots(q^{2}-{M_{N}}_{i_{n}})}$, and can be re-expressed in terms of $l^{r}_{\text{N}}(i_{x})$ $(x=0,1,\cdots,n)$ using the “divided difference” method [64].

Step 3. One now needs to perform the loop-momentum integration over the products of blocks $L_{\text{X}}(n;i,j)$, $L_{\text{XX}}(n;i,j;i^{\prime},j^{\prime})$, and $l^{r}_{\text{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 ${\partial}_{M_{\text{SQ}}}{\text{Loop}}_{\text{X}}$, where $M_{\text{SQ}}$ denotes symbolically the squared mass, equaling to $M^{2}$ for scalars and to $M^{{\dagger}}M$ or $MM^{{\dagger}}$ for fermions, and ${\text{Loop}}_{\text{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].

2.4 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 [65, 57, 66, 67].

Taking the block term $L_{\text{UU}}(2n;3,2,3,2)$ in case I,

[TABLE]

where $\Delta_{1}\equiv\frac{{M_{S}}_{1}{M_{S}}_{2}\delta^{2}_{23}}{q^{2}-{M_{S}}_{1}}+\frac{{M^{2}_{S}}_{2}\delta^{2}_{36}}{q^{2}-{M_{S}}_{2}}$, as an example, and using the inequality [63]

[TABLE]

satisfied by the PV-integrals with vanishing external momenta that are defined by

[TABLE]

with the assumption that $n\geqslant 3$ to avoid divergent integrals, we can obtain

[TABLE]

where $L_{\text{oth}}$ represent the blocks related to the other types of particles, such as $L_{\text{CC}}$. Then, for a given small constant $0<c_{0}<1$, only when

[TABLE]

can the terms starting from the $(2n+1)$-th MI order in the series expansion of the block $L_{\text{UU}}(3,2;3,2)$ be safely neglected. Thus, the cutting-off MI order for $L_{\text{UU}}(3,2;3,2)$ should be $2n$ at least. The same method can be applied for other blocks, and the final cutting-off MI orders for squarks and charginos can be determined accordingly.

For the neutralino blocks, Eq. (2.26) still works for estimating the required MI order. In this case, we obtain

[TABLE]

So, when $\left|\delta^{N}_{{i_{0}}{i_{1}}}\,\delta^{N}_{{i_{1}}{i_{2}}}\cdots\delta^{N}_{{i_{n-1}}{j_{0}}}\right|<c_{0}$ and $\left|\delta^{N}_{{i_{0}}{i_{1}}}\,\delta^{N}_{{i_{1}}{i_{2}}}\cdots\delta^{N}_{{i_{n-1}}{i_{n}}}\,{(M_{\chi_{0}})}_{i_{n}j_{0}}\right|\frac{1}{\sqrt{{M_{N}}_{i_{n}}}}<c_{0}$ for fixed indices $i_{0}$ and $j_{0}$, the summation over the MI index can be terminated to the $n$-th order.

3 $B_{s(d)}-\bar{B}_{s(d)}$ mixing and $B_{s}\to\mu^{+}\mu^{-}$ decay

In this section, we shall apply the FET procedure with general/finite MI order to $B_{s(d)}-\bar{B}_{s(d)}$ mixing and $B_{s}\to\mu^{+}\mu^{-}$ decay, within the $\mathbb{Z}_{3}$-invariant NMSSM with NMFV.

3.1 $B_{s(d)}-\bar{B}_{s(d)}$ mixing

The strength of $B_{s(d)}-\bar{B}_{s(d)}$ mixing is described by the mass difference $\Delta M_{s(d)}$, defined by [89]

[TABLE]

where $M_{12}^{q}$ 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]

[TABLE]

Within the $\mathbb{Z}_{3}$-invariant NMSSM with NMFV, the following eight operators, as defined in Ref. [22], are all found to be relevant:

[TABLE]

where $\alpha$ and $\beta$ are the colour indices, $\sigma_{\mu\nu}=\frac{1}{2}[\gamma_{\mu},\gamma_{\nu}]$, and $P_{L,R}=(1\mp\gamma_{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 beyond777Here we do not consider the double-penguin diagrams, which involve the exchange of CP-even and CP-odd scalars, and can give significant contributions only for large values of $\tan\beta$ [22, 90, 91]. This is justified by our choices of the two sets of SUSY parameters collected in Table 2, with $\tan\beta$ being fixed at $3$ and $10$, respectively.. Within the SM, only $C_{1}^{\rm VLL}$ 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_{1}^{\rm VLL}$ are also known [92]. In the context of $\mathbb{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, up-squarks; 3) gluinos, down-squarks; 4) neutralinos, down-squarks; 5) mixed gluino, neutralino, down-squarks [27, 93, 22]. 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_{\chi}^{+}$, $Z_{\chi}^{-}$, and $Z_{\chi_{0}}$, as well as the blocks $L_{\text{X}}(i,j)$ and $L_{\text{XX}}(i,j;i^{\prime},j^{\prime})$. Our results for the Wilson coefficients agree with the ones given in Refs. [57, 29, 27]. 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 $\Delta M_{s}$ in case I, we encounter a term

[TABLE]

As $L_{\text{UU}}(i,j)$, with $i\neq j$, does not vanish only when $(i,j)=(2,3)$ or $(3,2)$, and because of $|K^{*}_{23}K_{32}|\ll|K^{*}_{33}K_{22}|$, we can neglect safely the term with $(i,j)=(2,3)$, to get

[TABLE]

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_{12}^{q}$ and hence the mass difference $\Delta M_{s(d)}$. For convenience of later discussions, the total contributions to $M_{12}^{q}$ are split into the following different parts:

[TABLE]

where $M_{12}^{(q)\,\text{SM}}$, $M_{12}^{(q)\,\text{CH}}$, $M_{12}^{(q)\,\text{C}}$, and $M_{12}^{(q)\,\text{NG}}$ represent contributions from the SM, the charged Higgs, the charginos, as well as the neutralinos and gluinos, respectively.

3.2 $\boldsymbol{B_{s}\to\mu^{+}\mu^{-}}$ decay

The rare decay $B_{s}\to\mu^{+}\mu^{-}$ 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 $\mathbb{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]. The relevant effective weak Hamiltonian reads [101, 102]888Here we need not consider the tensor operators $\mathcal{O}_{TX}=(\bar{b}\sigma^{\mu\nu}P_{X}s)(\bar{\mu}\sigma_{\mu\nu}\mu)$. While also receiving contributions from these three types of one-loop diagrams in the $\mathbb{Z}_{3}$-invariant NMSSM with NMFV, they do not contribute to this process due to the vanishing matrix elements, $\langle 0|\bar{b}\sigma^{\mu\nu}P_{X}s|B_{s}\rangle=0$.

[TABLE]

where the vector ($\mathcal{O}_{VXY}$) and scalar ($\mathcal{O}_{SXY}$) operators are defined, respectively, by

[TABLE]

The branching ratio of $B_{s}\to\mu^{+}\mu^{-}$ decay is then calculated to be [101, 102]

[TABLE]

where $\tau_{B_{s}}$ is the $B_{s}$-meson lifetime, and the squared matrix element is given by [101, 102, 103, 100]

[TABLE]

with the scalar, pseudo-scalar, and axial-vector form factors defined, respectively, by [101]

[TABLE]

where $f_{B_{s}}$ is the $B_{s}$-meson decay constant, and $\overline{m}_{b(s)}$ denotes the $b(s)$-quark running mass.

Our main task is then to calculate the Wilson coefficients $C_{VXY}$ and $C_{SXY}$. Within the SM, only $C_{VLL}$ gets a non-negligible contribution999When 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]., 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_{VXY}$ and $C_{SXY}$ from the $\mathbb{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)}-\bar{B}_{s(d)}$ mixing.

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]

[TABLE]

where $A_{\Delta{\Gamma}_{s}}$ is a time-dependent observable [107, 108, 109], and $y_{s}$ is related to the decay-width difference $\Delta\Gamma_{s}$ between the two $B_{s}$-meson mass eigenstates, defined by

[TABLE]

with $\Gamma^{s}_{L}$ ($\Gamma^{s}_{H}$) denoting the lighter (heavier) eigenstate decay width and $\Gamma_{s}=\tau_{B_{s}}^{-1}$ 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_{\Delta{\Gamma}_{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

[TABLE]

which holds to a very good approximation [106].

4 Numerical results and discussions

After getting the analytic FET results for the $B_{s(d)}-\bar{B}_{s(d)}$ mixing and $B_{s}\to\mu^{+}\mu^{-}$ decay, we now proceed to analyze numerically the parameter space of $\mathbb{Z}_{3}$-invariant NMSSM with NMFV that is allowed under these experimental constraints.

4.1 Choice of input parameters

Firstly, we collect in Table 1 part of the input parameters used throughout this paper. For the $B_{s(d)}$-mixing matrix elements and the decay constants $f_{B_{s(d)}}$, we take the values provided by the FNAL/MILC collaboration [42] and the averages by the Particle Data Group [54], respectively. The relevant model parameters of the $\mathbb{Z}_{3}$-invariant NMSSM with NMFV include

[TABLE]

where $M_{\tilde{g}}$ 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 $\lambda$ and a small $\tan\beta$, to avoid suppressing the NMSSM-specific tree-level contributions to the SM-like Higgs mass [15, 14], while scenario B is featured by negligible NMSSM-specific effects on the SM-like Higgs mass [112]. The remaining parameters ${M_{S}}_{2}$, $\delta_{23}^{\text{LL}}$, $\delta_{33}^{\text{LR}}$, $\delta_{13}^{\text{RR}}$, and $\delta_{23}^{\text{RR}}$ 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 $\lambda$ 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_{\lambda}$, which is closely related to the charged-Higgs mass [15]. In scenario A, the predicted lightest and next-to-lightest neutralino masses are given, respectively, by $m_{\tilde{\chi}_{1}^{0}}\sim 74$ GeV and $m_{\tilde{\chi}_{2}^{0}}\sim 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\to\tilde{\chi}_{1}^{\pm}\,\tilde{\chi}_{2}^{0}$ with a $100\%$ branching fraction of $\tilde{\chi}_{2}^{0}\to H_{\text{S}M}\,\tilde{\chi}_{1}^{0}$, by the CMS collaboration [113]. In scenario B, on the other hand, $m_{\tilde{\chi}_{1}^{0}}\sim 179$ GeV and $m_{\tilde{\chi}_{2}^{0}}\approx m_{\tilde{\chi}_{1}^{\pm}}\sim 207$ GeV, being compatible with the mass bounds from the same channel but with a $100\%$ branching fraction of $\tilde{\chi}_{2}^{0}\to Z\,\tilde{\chi}_{1}^{0}$ [113]. The mass splitting $\Delta m(\tilde{\chi}_{2}^{0},\tilde{\chi}_{1}^{0})=m_{\tilde{\chi}_{2}^{0}}-m_{\tilde{\chi}_{1}^{0}}$ in both of these two scenarios is also compatible with the bound set by the CMS collaboration [114, 113]. In addition, the choice $\mu_{\text{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\rightarrow X_{s}\gamma$ 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\rightarrow X_{s}\gamma$ decay marginal. Taking together with the values of $\tan\beta$, we can say that both of these two scenarios make the Higgs-penguin effects negligible for both $B_{s(d)}-\bar{B}_{s(d)}$ mixing [57] and $B_{s}\to\mu^{+}\mu^{-}$ decay [101].

The parameters $\delta_{33}^{\text{LR}}$ and ${M_{S}}_{2}$ are chosen to get the SM-like Higgs mass in the range $122~{}{\rm GeV}\leqslant m_{H_{SM}}\leqslant 128~{}{\rm GeV}$ [133]. The allowed regions for $\delta_{33}^{\text{LR}}$ and ${M_{S}}_{2}$, obtained with the aid of the package NMSSMCALC [76, 134, 135, 136, 137, 138, 139], are shown in Figure 1. It can be seen that $940^{2}~{}\text{GeV}^{2}\leqslant{M_{S}}_{2}\leqslant 2000^{2}~{}\text{GeV}^{2}$ and $-0.5\leqslant\delta_{33}^{\text{LR}}\leqslant 0.5$ in scenario A, while $1500^{2}~{}\text{GeV}^{2}\leqslant{M_{S}}_{2}\leqslant 2000^{2}~{}\text{GeV}^{2}$ and $\delta_{33}^{\text{LR}}$ is only allowed to be around $0.2$ in scenario B. These bounds will be taken into account in the following numerical analyses.

4.2 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 $\delta^{N}_{ij}$ depend on the six $\mathbb{Z}_{3}$-invariant NMSSM parameters $M_{1}$, $M_{2}$, $\mu_{\text{eff}}$, $\tan\beta$, $\lambda$, and $\kappa$. Keeping for the moment $\lambda$ and $\kappa$ as free variables, but with $\tan\beta=3$ and $10$ corresponding respectively to the two scenarios defined in Table 2, we find that $\delta^{N}_{15}=\delta^{N}_{25}=0$ and $\delta^{N}_{12}=-0.007$. The dependence of the other MI parameters $\delta^{N}_{ij}$ on $\lambda$ and $\kappa$ are displayed in Figure 2. From the magnitudes of these $\delta^{N}_{ij}$ shown and based on the criterion specified in Sec. 2.4, one can see that the effects of $\delta^{N}_{12}$, $\delta^{N}_{13}$, and $\delta^{N}_{34}$ are negligible with the MI order higher than $1$, and that of $\delta^{N}_{14}$ and $\delta^{N}_{23}$ can be neglected starting from the third MI order (even from the second MI order for $\delta^{N}_{23}$ in scenario B); while the MI orders higher than $3$ should be kept for $\delta^{N}_{24}$, $\delta^{N}_{35}$, and $\delta^{N}_{45}$. When the parameters $\lambda$ and $\kappa$ are fixed at the values in scenarios A and B, it is further found that the values of ${{(M_{\chi_{0}})}_{i_{n}\,3}}/{\sqrt{{M_{N}}_{i_{n}}}}$, for $i_{n}=1,\,2,\,3$, are much smaller than for $i_{n}=4,\,5$, and all the terms involving $L_{\text{N}}(n;3,i_{n}){(M_{\chi_{0}})}_{i_{n}3}$ 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$.

4.2.2 MI-order comparison for $\boldsymbol{B_{s(d)}-\bar{B}_{s(d)}}$ mixing

After obtaining the optimal cutting-off MI orders, we now check the convergence of the FET results with different MI orders. Let us firstly discuss the $B_{s(d)}-\bar{B}_{s(d)}$ mixing. The left panel in Figure 3 shows the variation of $\Delta M_{s}$ with respect to $\delta_{23}^{\text{LL}}$ 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 $\Delta M_{s}$ varies with respect to $\delta_{23}^{\text{LL}}$ 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 $\delta_{23}^{\text{LL}}$.

In order to show the convergence of the FET results in the case II assumption for squark flavour structures, we consider the ratio

[TABLE]

which gives the difference between $\Delta M_{s}$ by considering the squark MI orders of $6$ and $8$, respectively. As shown in the right panel of Figure 3, $O_{\Delta M_{s}}$ is nearly $1$ in the whole area in case IIA, implying that the convergence has been verified. Similar observations are also made for $\Delta M_{s}$ in case IIB and for $\Delta M_{d}$ in all the four cases.

4.2.3 MI-order comparison for $\boldsymbol{B_{s}\to\mu^{+}\mu^{-}}$ decay

For our prediction of the time-integrated branching ratio $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ in the $\mathbb{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 $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ with respect to $\delta_{23}^{\text{LL}}$ for $\delta_{33}^{\text{LR}}=-0.15$ in case IA. Being affected by $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ quite weakly [140, 57], the convergence of $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ 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)}-\bar{B}_{s(d)}$ mixing.

When the off-diagonal element $\delta_{23}^{\text{LL}}$ is zero, there exists no squark MI contribution and the NMSSM contributions come only from the charged Higgs and the diagonal part of $M_{\tilde{U}}^{2}$. Even in this case, $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ can be obviously enhanced, putting it to be larger than the LHCb measurement, $(3.0\pm 0.6^{+0.3}_{-0.2})\times 10^{-9}$ [141]. As a reference, our prediction within the SM is $(3.6\pm 0.3)\times 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_{B_{s}}$ from Ref. [54].

4.2.4 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 $\text{Re}(2M_{12}^{(s)\,i})$ from different gauginos (see Eq. (3.6) for their respective definitions) obtained with these two methods. Here only the NMSSM contribution to $\text{Re}(2M_{12}^{(s)\,\text{NG}})$ is shown in case IIA, because only the parameters $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ from $M_{\tilde{D}}^{2}$ are involved in this case, and they do not contribute to $\text{Re}(2M_{12}^{(s)\,\text{C}})$. It can be seen clearly that the FET results with squark MI order of $8$ agree quite well with the ones calculated directly in the ME basis. One can also find from the left panel that, as the MI parameter $\delta_{23}^{\text{LL}}$ increases with the other parameters fixed in case IA, the NMSSM contribution to $\text{Re}(2M_{12}^{(s)\,\text{C}})$ becomes larger, while both $\text{Re}(2M_{12}^{(s)\,\text{NG}})$ and $\text{Re}(2M_{12}^{(s)\,\text{C}+\text{NG}})$ become smaller; in addition, their dependence on the parameter ${M_{S}}_{2}$ becomes weaker when $\sqrt{{M_{S}}_{2}}$ varies from $1100$ GeV to $1500$ GeV. In case IIA, on the other hand, the variations of $\text{Re}(2M_{12}^{(s)\,\text{NG}})$ with respect to $\delta_{23}^{\text{RR}}$ and $\delta_{13}^{\text{RR}}$ depend heavily on the chosen values of $\sqrt{{M_{S}}_{2}}$: while the contours for negative $\text{Re}(2M_{12}^{(s)\,\text{NG}})$ are similar, the trends for the extremum values of $\text{Re}(2M_{12}^{(s)\,\text{NG}})$ with $\sqrt{{M_{S}}_{2}}=1100$ GeV are quite different from the ones with $\sqrt{{M_{S}}_{2}}=1300$ GeV or $1500$ GeV. All these observations can be more conveniently and easily understood in terms of the FET result for $\text{Re}(2M_{12}^{(s)\,i})$, which is essentially a polynomial with the MI parameters as the variables.

4.3 Constraints on the $\mathbb{Z}_{3}$-invariant NMSSM parameters

As mentioned in the Introduction, the SM predictions for the mass differences $\Delta M_{s}$ and $\Delta M_{d}$ are now larger than their respective experimental values. For the time-integrated branching ratio $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$, on the other hand, the 2017 LHCb measurement101010This decay has been searched for by CDF [142] and D0 [143], and was observed for the first time by LHCb [144] and CMS [145], with their combined average for the branching ratio given in Ref. [146]. Searches for this decay have also been performed by ATLAS [147]. The 2017 LHCb measurement includes the LHC Run 2 data and represents the first single-experiment observation of this decay, with a $7.8\sigma$ significance [141]., $(3.0\pm 0.6^{+0.3}_{-0.2})\times 10^{-9}$ [141], is in reasonable agreement with the SM prediction, $(3.6\pm 0.3)\times 10^{-9}$. This will impose much more stringent constraints on NP [110].

The inclusive radiative decay $B\to X_{s}\gamma$ provides also important constraints on NP scenarios with extended Higgs sectors or SUSY theories [148, 149, 150, 151]. Measurements of its CP- and isospin-averaged branching ratio by the BaBar [152, 153, 154] and Belle [155, 156] collaborations lead to the following combined result [132]

[TABLE]

which is in excellent agreement with the state-of-the-art SM prediction [157]

[TABLE]

Here the photon energy cutoff $E_{\gamma}>1.6\,\text{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\to X_{s}\gamma$ constraint [70, 22]. As the SM and charged-Higgs contributions to $\mathcal{B}_{d\gamma}$, the branching ratio of $B\to X_{d}\gamma$ decay, are both suppressed by the CKM factor $\left|K_{31}/K_{32}\right|^{2}$ with respect to $\mathcal{B}_{s\gamma}$, while the contributions from the squark MI parameters are not affected by this factor, it is expected that $\delta_{13}^{\text{LL}}$ in $M_{\tilde{U}}^{2}$ will be strongly constrained by $\mathcal{B}_{d\gamma}$. As a result, we have set $\delta_{13}^{\text{LL}}$ 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 $\mathcal{B}_{s\gamma}$ 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_{\gamma}>1.6\,\text{GeV}$.

In the following, we shall exploit the $95\%$ C.L. bounds from $\Delta M_{s}$, $\Delta M_{d}$, $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$, and $\mathcal{B}_{s\gamma}$, to set the allowed regions for the NMSSM parameters ${M_{S}}_{2}$, $\delta_{23}^{\text{LL}}$, $\delta_{33}^{\text{LR}}$, $\delta_{13}^{\text{RR}}$, and $\delta_{23}^{\text{RR}}$, 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\sigma$ 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.

4.3.1 Results in scenario A

Firstly, we show in Figure 6 the allowed regions for the squark MI parameters $\delta_{23}^{\text{LL}}$, $\delta_{33}^{\text{LR}}$, $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ 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}~{}\text{GeV}^{2}$ are made in both cases IA and IIA, and the squark MI parameter $\delta_{33}^{\text{LR}}$ is set to be $-0.4$ in case IIA. The scenario A is featured by the observation that the charged-Higgs contribution to $\Delta M_{s}$ is positive and large for small charged-Higgs mass and $\tan\beta$, while its contribution to $\mathcal{B}_{s\gamma}$, being also positive and large for small charged-Higgs mass, is less sensitive to the variation of $\tan\beta$. The bound from $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ 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 $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ is nearly not affected by $\delta_{13}^{\text{RR}}$ or $\delta_{23}^{\text{RR}}$ in case II [140, 57]. 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 $\delta_{23}^{\text{LL}}$ from $\Delta M_{s}$ also becomes reduced when ${M_{S}}_{2}$ increases, and only the one with large magnitudes of $\delta_{23}^{\text{LL}}$ exists in case IA. It is particularly observed that, for ${M_{S}}_{2}=1100^{2}~{}\text{GeV}^{2}$, there exists no allowed region for $\delta_{23}^{\text{LL}}$ and $\delta_{33}^{\text{LR}}$ in this case. In case IIA, on the other hand, the bounds from $\Delta M_{s}$ and $\Delta M_{d}$ become stronger for a larger ${M_{S}}_{2}$, and the one from $\Delta 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 $\mathcal{B}_{s\gamma}$ also shows that a larger magnitude of $\delta_{23}^{\text{LL}}$ or $\delta_{33}^{\text{LR}}$ is required in case IA, but almost the whole area of $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ 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 $\Delta M_{s}$ and $\mathcal{B}_{s\gamma}$ from the charged Higgs; then a large magnitude of $\delta_{23}^{\text{LL}}$ or $\delta_{23}^{\text{RR}}$ ($\delta_{33}^{\text{LR}}$) is needed to provide large but negative contributions to $\Delta M_{s}$ ($\mathcal{B}_{s\gamma}$), so as to cancel the charged-Higgs effects [70, 22]. Especially in case IIA, the chosen value of $-0.4$ for $\delta_{33}^{\text{LR}}$ is already large enough to cancel the charged-Higgs effect on $\mathcal{B}_{s\gamma}$, and the flavour-violating contribution from the RR sector is small, making the whole area of $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ being allowed by this observable. Compared to the case from $\Delta M_{s}$, the allowed regions for the squark MI parameters from $\Delta M_{d}$ are relatively large, because the charged-Higgs contribution to $\Delta 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 $\delta_{13}^{\text{RR}}$.

After taking into account the $95\%$ C.L. bounds from $\Delta M_{s}$, $\Delta M_{d}$, $\mathcal{B}_{s\gamma}$, $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$, as well as the SM-like Higgs mass, we find numerically that the squark MI parameters $\delta_{23}^{\text{LL}}>0.45$ and $|\delta_{33}^{\text{LR}}|\sim 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 $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$.

4.3.2 Results in scenario B

In Figure 7, we show the allowed regions for $\delta_{23}^{\text{LL}}$, $\delta_{33}^{\text{LR}}$, $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ in scenario B. To match the SM-like Higgs mass, as shown in Figure 1, the squark MI parameter $\delta_{33}^{\text{LR}}$ 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}$) $\text{GeV}^{2}$ are made in case IB (IIB). The scenario B is featured by the observation that the positive contributions to $\Delta M_{s}$ and $\mathcal{B}_{s\gamma}$ 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 $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ is satisfied in the whole parameter regions shown here.

In case IB, the allowed region for the squark MI parameters from $\Delta M_{s(d)}$ is large because of the small charged-Higgs contribution, and is almost independent of ${M_{S}}_{2}$; while the bound from $\mathcal{B}_{s\gamma}$ indicates that a smaller $\delta_{23}^{\text{LL}}$ is favored. In case IIB, on the other hand, while the allowed region from $\Delta M_{d}$ stays nearly the same, the one from $\Delta M_{s}$ shrinks slowly when ${M_{S}}_{2}$ increases. Under the bound of $\mathcal{B}_{s\gamma}$, an allowed region with a positive $\delta_{23}^{\text{RR}}$ begins to emerge only when ${M_{S}}_{2}$ is about $1900^{2}~{}\text{GeV}^{2}$.

After considering the $95\%$ C.L. bounds from $\Delta M_{s}$, $\Delta M_{d}$, $\mathcal{B}_{s\gamma}$, as well as the SM-like Higgs mass, the parameter $\delta_{23}^{\text{LL}}$ 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 $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ exists only when ${M_{S}}_{2}$ is larger than about $1900^{2}~{}\text{GeV}^{2}$.

5 Conclusion

In this paper, motivated by the observation that the SM predictions are now above the experimental data for the mass difference $\Delta M_{s(d)}$, and the usual CMFV models have difficulties in reconciling this discrepancy, we have investigated whether the $\mathbb{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}\to\mu^{+}\mu^{-}$ and $B\to X_{s}\gamma$ decays.

We have calculated the NMSSM contributions to $\Delta M_{s(d)}$ and the branching ratio $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$, 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 $\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 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)}-\bar{B}_{s(d)}$ mixing and $B_{s}\to\mu^{+}\mu^{-}$ 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 $\text{Re}(2M_{12}^{s})$ 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 $\Delta M_{s}$, $\Delta M_{d}$, $\overline{\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$, $\mathcal{B}_{s\gamma}$, as well as the SM-like Higgs mass, we have discussed the allowed regions for the parameters ${M_{S}}_{2}$, $\delta_{23}^{\text{LL}}$, $\delta_{33}^{\text{LR}}$, $\delta_{13}^{\text{RR}}$, and $\delta_{23}^{\text{RR}}$. It is found that only large values of $\delta_{23}^{\text{LL}}$, with $|\delta_{33}^{\text{LR}}|\sim 0.15$, are allowed in case IA, and the allowed region for $\delta_{23}^{\text{LL}}$ becomes reduced when ${M_{S}}_{2}$ increases from $1100^{2}$ to $1500^{2}~{}\text{GeV}^{2}$. In case IB, on the other hand, with $\delta_{33}^{\text{LR}}$ being fixed at about $0.2$, relatively smaller magnitude of $\delta_{23}^{\text{LL}}$ is found to be allowed, and the allowed region for $\delta_{23}^{\text{LL}}$ becomes almost independent of ${M_{S}}_{2}$. In case IIA, with $\delta_{33}^{\text{LR}}$ being fixed at $-0.4$, there exists no allowed region for $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$, because of the strong bound from $\Delta M_{s}$. In case IIB, on the contrary, with $\delta_{33}^{\text{LR}}$ being fixed at about $0.2$, the allowed region for $\delta_{13}^{\text{RR}}$ and $\delta_{23}^{\text{RR}}$ exists only when ${M_{S}}_{2}$ is larger than about $1900^{2}~{}\text{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.

Acknowledgements

We would like to thank Michael Paraskevas, Janusz Rosiek, and Xing-Bo Yuan for valuable discussions and be grateful to Michael Paraskevas for his guidance on SUSY_FLAVOR. This work is supported by the National Natural Science Foundation of China under Grant Nos. 11675061, 11775092, 11521064, and 11435003. X.L. is also supported in part by the self-determined research funds of CCNU from the colleges’ basic research and operation of MOE (CCNU18TS029). Q.H. is also supported by the China Postdoctoral Science Foundation (2018M632896).

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 $\Delta_{1}\equiv\frac{{M_{S}}_{1}{M_{S}}_{2}\delta^{2}_{23}}{q^{2}-{M_{S}}_{1}}+\frac{{M^{2}_{S}}_{2}\delta^{2}_{36}}{q^{2}-{M_{S}}_{2}}$ and $\Delta^{\prime}_{1}\equiv\frac{{M_{S}}_{1}{M_{S}}_{2}(1+\lambda^{2}_{\text{CKM}})\delta^{2}_{23}}{q^{2}-{M_{S}}_{1}}$ in case I, and $\Delta_{2}\equiv\frac{{M_{S}}_{1}{M_{S}}_{2}(\delta^{2}_{46}+\delta^{2}_{56})}{q^{2}-{M_{S}}_{1}}$ in case II. In the following, $n=1,2,3,\cdots$, denotes the MI-order index.

For up-type squarks, the non-zero block terms are given as

[TABLE]

in case I.

For down-type squarks, the non-zero block terms are given as,

[TABLE]

in case I, and

[TABLE]

in case II.

The non-zero block terms for charginos include

[TABLE]

where $(i,i^{\prime})=(1,2)$ or $(2,1)$,

[TABLE]

[TABLE]

where $(i,j)=(1,1),(1,2),(2,1)$, or $(2,2)$,

[TABLE]

where $i=1$ or $2$, and

[TABLE]

Here the subscripts X and Y can be C or P.