Muon $g-2$ and semileptonic $B$ decays in B\'elanger-Delaunay-Westhoff model with gauge kinetic mixing

In the model proposed by B\'elanger, Delaunay and Westhoff (BDW), a new sector consisted of vectorlike fermions and two complex scalars is charged under an extra Abelian symmetry $U(1)_X$. In this paper, we generalize the BDW model by introducing the kinetic mixing between the $U(1)_X$ and the standard model $U(1)_Y$ gauge fields. The new physics contributions to the muon anomalous magnetic moment and the Wilson coefficients $C_{9,10}^{(')}$ are obtained analytically. We have explored the free parameter space of the model, taking into account various constraints on the muon $g-2$ using recent data from the E989 experiment at Fermilab, the lepton universality violation in terms of $R_K$ and $R_{K^*}$, and the branching ratios of the semileptonic decays, $B^+ \rightarrow K^+ \mu^+ \mu^-$ and $B^0 \rightarrow K^{*0} \mu^+ \mu^-$, the LEP and LHC searches for sleptons and $Z'$ boson, as well as the perturbative requirement. The viable parameter regions of the model are identified. In the presence of the gauge kinetic mixing term, those regions are enlarged and significantly deformed in comparison to the case with vanishing kinetic mixing. In the near future, the E989 experiment with the projected sensitivity will be able to test significant parts of the currently allowed parameter regions.

( ) 9,10 are obtained analytically. We have explored the free parameter space of the model, taking into account various constraints on the muon g − 2 using recent data from the E989 experiment at Fermilab, the lepton universality violation in terms of R K and R K * , and the branching ratios of the semileptonic decays, B + → K + µ + µ − and B 0 → K * 0 µ + µ − , the LEP and LHC searches for sleptons and Z boson, as well as the perturbative requirement. The viable parameter regions of the model are identified. In the presence of the gauge kinetic mixing term, those regions are enlarged and significantly deformed in comparison to the case with vanishing kinetic mixing. In the near future, the E989 experiment with the projected sensitivity will be able to test significant parts of the currently allowed parameter regions.

Introduction
Although the standard model (SM) predictions have been verified in many experiments showing an excellent agreement with data, new physics seems to be around the corner due to unanswered questions. On the one hand, one of the most important precision tests of the SM is the muon anomalous magnetic moment (g − 2) whose value was determined accurately as [1,2] a exp µ = (11659206.1 ± 4.1) × 10 −10 . ( This is the average value taking into account the new result from the E989 experiment at Fermilab [2] which confirmed the previous E821 measurement at Brookhaven National Laboratory. However, the SM prediction for the muon g − 2 is presently [3] a SM µ = (11659181.0 ± 4.3) × 10 −10 , corresponding to a 4.2σ deviation from the above world-average experimental value. According to the projected sensitivity of the E989 experiment [4] as well as the future experiment (E34) at J-PARC [5], the precision will be improved by a factor of four that will shed light on this deviation. If the experimental center value of the muon g − 2 remains unchanged, the above deviation will raise up to about more than 5σ [6], evidently indicating the existence of new physics coupled to the lepton sector [7].
On the other hand, with the improvement of experimental accuracy, rare decays of B-mesons are useful as the precision tests of the SM. The small branching ratios of these processes make them good probes to search for new physics beyond the SM. In fact, anomalies have been observed in the rare semileptonic B decays related to the quark transition process b → s + − .
In particular, the measured relative branching ratios for 1.1 GeV < q 2 < 6.0 GeV [8,9] R K = BR(B + → K + µ + µ − ) BR(B + → K + e + e − ) = 0.846 +0.044 −0.041 , deviate from the corresponding SM predictions being close to unity [10,11,12] at the levels of more than 3σ and 2.4σ, respectively [13,14]. This may be a signature of the lepton universality violation implying the existence of new physics beyond the SM. Recently, the updated LHCb results on the angular analysis of the decay process B 0 → K * 0 (892)µ + µ − show a 3.3σ deviation from the SM prediction [15,16] that is slightly increased in comparison to the previous observation [17]. This is due to the higher statistics regarding to the inclusion of 2016 data for 13 TeV collisions.
proach with the Hamiltonian [18] where α and G F are the fine structure constant and the Fermi constant, have been performed to determine the new physics contributions to the relevant Wilson coefficients [19] C corresponding to the operators The Wilson coefficients C ( ) S,P and C ( ) 7 are strictly constrained by the leptonic decay process B s → µ + µ − [20] and the radiative B decays [21], respectively. In the case of semileptonic decays of B mesons, the fitting results showed that scenarios with new physics contributions to the Wilson coefficients is much more favored than the pure SM, especially for the coefficient C µ 9 . While the new physics contributions to muonic Wilson coefficients (C µ( ) 9 , C µ( ) 10 ) play an important role, the electronic coefficients (C e( ) 9 , C e( ) the gauge kinetic mixing, we calculate analytically the new physics contributions to the muon g − 2, and the Wilson coefficients C ( ) 9,10 which are used to compute the B-meson semileptonic decays. The most updated data from the Heavy Flavor Averaging Group [32] and the LHCb Collaboration [8,9] are used in our consideration to constrain the model's relevant parameters.
Using the projected sensitivity of the muon g − 2 experiment E989, we study its ability to test the model in the near future.
The structure of the paper is as follows. In Section 2, we briefly review the structure of the BDW model, then generalize it by introducing the gauge kinetic mixing. In Section 3, the analytic results of the new physics contributions to the muon g − 2 and the Wilson coefficients calculations are presented. In Section 4, the phenomenological constraints are used to specify the allowed parameter space. Finally, Section 5 is devoted to conclusions.

The model
The new particles introduced in the BDW model beside the SM particles are the vectorlike lepton and quark doublets of the gauge group SU (2) L , and two complex scalars, χ and φ, that are singlets under the SM gauge groups. The symmetry of this model is an extension of the SM symmetry by adding an extra Abelian gauge group, transformation, while the new particles transform nontrivially with the U (1) X charges given in Table 1 together with other properties.  [31].
The Lagrangian include the SM part and the part involving new physics: where Here, the SM left-handed lepton and quark doublets are denoted as The explicit form of the scalar potential V 0 (χ, φ) is given by Among the new scalars, we assume that only φ can develop a vacuum expectation value where with H = 174 GeV being the VEV of the SM Higgs field. Due to the nonzero VEV, φ , the gauge group U (1) X is spontaneously broken, leading to a massive Z boson with a mass where g X is the U (1) X gauge coupling.
By decomposing the complex scalar field φ into their real and imaginary components, the masses of these fields are found to be respectively. Note that ϕ i is a massless Nambu-Goldstone boson that can be absorbed by Z in the unitary gauge. For the case of χ, after the decomposition the mass matrix for these real scalar fields is obtained: where In the simple case where the coupling r is real, the matrix M 2 χ is diagonal, and the masses of the particles χ r and χ i are respectively Since the field χ does not develop a nonzero VEV, there is no mass mixing between the SM leptons and the vectorlike ones. However, the situation for quarks is more involved because the VEV of φ generates mass mixing terms via the new Yukawa interactions with the couplings w = (w 1 , w 2 , w 3 ) in Eq. (14). To diagonalize the quark mass matrices, M u and M d , we need to use four 4 × 4 unitary matrices to transform the quark gauge eigenstates, (u 1 , u 2 , u 3 , U ) and (d 1 , d 2 , d 3 , D), into the mass eigenstates, (u, c, t, U) and (d, s, b, D): The diagonal mass matrices of up-type and down-type quarks then read

The gauge kinetic mixing
In a model with two Abelian gauge symmetries, the gauge kinetic mixing term of the form − k 2 F 1 µν F 2µν can be introduced in the Lagrangian without violating any symmetry [33]. Here, F 1 µν , F 2µν and k are respectively the field strength tensors of the U (1) Y and U (1) X gauge fields and the kinetic mixing coefficient. In fact, the gauge kinetic mixing term is always generated radiatively [34], even though it is set to zero at high energy scales [35]. In the presence of such term, the gauge kinetic part of the Lagrangian relating to the Abelian groups can be written as L gauge kinetic ⊃ − By an appropriate transformation in the space of the Abelian gauge fields, the kinetic Lagrangian can be made canonical. In the new basis, the covariant derivative is then expressed where the new charge X and the new gauge coupling g X are determined by the original quantities and the kinetic mixing coefficient k as Here, Y (X) and g Y (g X ) are the charge and the gauge coupling corresponding to the Abelian group U (1) Y (U (1) X ), respectively. The nonzero kinetic mixing coefficient k induces a shift in the U (1) X charge and modifies the relevant gauge coupling.
After the electroweak symmetry breaking by the VEV of the SM Higgs field, H , the kinetic coefficient generates a mass mixing between the Z and Z bosons, leading to a non-diagonal mass matrix for these particles: To diagonalize the above mass matrix, we use the following orthogonal rotation: where Z and Z are the mass eigenstates, and the mixing angle α Z is determined as It is worth noticing that, in the limit of no kinetic mixing (k = 0), the pure BDW model is recovered, namely X = X, g X = g X , and α Z = 0.
From the covariant derivative of muon, we find the interaction terms between the muon and the Z and Z bosons as where In the limit of no kinetic mixing, muons interact only with the Z boson just as in the SM. Similarly, the interaction terms between the new charged lepton E R and these two massive neutral gauge bosons are obtained as where For the new scalar sector, the derivation of Z and Z couplings with χ r,i from their covariant derivatives leads to where X χ = X χ = −1.
The flavor changing neutral currents (FCNCs) in the quark sector are induced at the tree level due to the mixing between the SM quarks and the vectorlike quarks. We parameterize the couplings between the b-quark, the s-quark and the massive neutral gauge bosons as follows where and the parameters A bs , B bs and C bs characterize the FCNCs in the b → s transition. At the tree level, these parameters are determined as Assuming that the mixing between the vectorlike quarks and the first generation quarks is negligible for simplicity, the parameter C bs can be approximated as At the loop-level, these parameters are considered to be the effective couplings encoding the new physics relevant to the quark sector.
3 New physics contributions to muon g − 2 and Wilson coefficients 3.1 Muon g − 2 In this model, new physics contributes to the muon anomalous magnetic moment via the gauge interaction associated with the massive boson Z , and the new Yukawa interaction between the scalar χ, the right-handed lepton E R and the muon. The Feynman diagrams corresponding to the leading new physics contributions to the muon g − 2 are shown in Figure 1. The diagrams (a) and (b) in this figure are due to the Yukawa coupling y µ in Eq. (14). The contribution related to the diagram (c) is generated from the gauge kinetic mixing effect.
From the matrix elements of these one-loop diagrams, after performing some algebraic calculation, we obtain the following expression for the new physics contributions to the muon where The first term in Eq. (54) with the squared brackets corresponds to the new physics contributions in Figures 1a and 1b. Here, the loop function F g (x) is defined as This term is in agreement with the result in Ref. [31] for the case of no kinetic mixing The second term in Eq. (54) with the integral results from the effect of the gauge kinetic mixing via the diagram in Figure 1c. It vanishes in the limit k = 0. Since the second term is suppressed by the factor β as well as the kinetic mixing coefficient, the sign of ∆a N P µ is determined by the sign of the first term in Eq. (54) that is always positive. Since the scalar and pseudoscalar Wilson coefficients (C S,P ) are severely restricted by the leptonic decay B s → µ + µ − [20], in this subsection we consider the Wilson coefficients C ( ) 9,10 to be used for the calculation of the semileptonic branching ratios of B mesons. The leading new physics contributions to the Wilson coefficients C ( ) 9 and C ( ) 10 are depicted by the Feynman diagrams in Figure 2. Here, we consider the diagrams whose the matrix elements are proportional to the second order of the gauge couplings g 2 or g X .

Wilson coefficients
In Figure 2a, the new physics enters this tree-level diagram only via the coupling of the bsZ vertex resulting from the mixing between the SM quarks and the vectorlike quarks. The contribution according to the diagram 2b is due to the gauge kinetic mixing effect. The contributions due to the diagrams (c)-(f) stem from the gauge interactions of the vectorlike charged lepton E R with both Z and Z bosons. The cases for the diagrams (g)-(j) relevant to the gauge interactions of the scalars χ r,i are more involved. While the contributions from the diagrams (h) and (j) always exist, those from the diagrams (g) and (i) are due to the ZZ mixing that only emerges when k = 0.
The new physics contributions to the Wilson coefficients C ( ) 9,10 according to the diagrams in Figure 2 are expressed as In these above formulas, the loop functions f A , g A , f B , and g B are given by as the results of the Feynman parameterization. In these formulas, is one of the SM charge leptons {e, µ, τ }.

Numerical analysis
In the numerical analysis, we assume, for simplicity, that the new vectorlike leptons have sizable coupling with muons (y µ ), while the corresponding coupling with electrons (y e ) is negligible.
Therefore, the set of relevant input parameters are m χr , m Z , k, g X , y µ , τ, δ, A bs , B bs .
In this section, we consider the phenomenological constraints including the muon anomalous magnetic moments, and the rare semileptonic decays of B mesons. The current deviation between the experimental value and the SM prediction of muon g − 2 is [1,2,3] The on-going E989 experiment will be able to reach a precision of 140 parts-per-billion [4].
Assuming the same center value of ∆a µ as the current result (72), the projected difference between the SM prediction and the experimental value reads corresponding to a 5.5σ deviation.
For the B-meson semileptonic decays, we take into account the branching ratios of the and the observables, R K and R K * , characterizing the violation of lepton flavor universality. For the muon invariant mass in the region q 2 = [1.1, 6.0] GeV 2 , the 2σ allowed ranges for the following B-meson observables are: [39,9] according to the updated results from the LHCb experiment. The branching ratios of the Bmeson decay processes are calculated using the methods described in Refs. [40,41,42]. The relevant form factors from Ref. [43] are employed in our calculation.
The searches for the Z boson have been carried out at the LHC run II in various channels where the resonance decays into dilepton [44,45,46], diquark [47,48,49,50,51], and Zh [52] for a wide range of m Z up to 6 TeV. These analyses set the upper limits on the production cross section times branching ratios of the Z boson that in turn impose constraints on the Z mass and its couplings to the SM particles. In the BDW model, to address the muon g − 2 anomaly, the loop-induced effective coupling between Z and muons needs to be large enough.
As a consequence, the Z boson decays dominantly to µμ and ν µνµ [31]. Therefore, among these constraints, those from the dimuon searches is the most severe for the BDW model.
The constraints on the kinetic mixing coefficient k and m Z from electroweak precision tests and other various experimental data from channels like the h → ZZ and h → Z Z decays, the Drell-Yan Z production were studied in Ref. [53]. According to that, for a wide range of the Z -boson mass below O(1) TeV, the current limit for the gauge kinetic mixing coefficient is . For small mass region of the U (1) X gauge boson below 10 GeV, the analyses by the BarBar Collaboration [54] and the KLOE-2 Collaboration [55] indicate the most stringent upper bound on the kinetic mixing coefficient. Recently, the CMS Collaboration investigated the muon pair production channel at the LHC with the center of mass energy √ s = 13 TeV in the search for a narrow resonance [56]. Similar analysis was also studied by the LHCb Collaboration [57]. These results show a severe constraint on the kinetic mixing coefficient with a Z boson lighter than 200 GeV. The approximated upper bound on k for various range of m Z is summarized as follows At the LEP experiment, the search for sleptons used the channels with the same final states as those coming from the vectorlike leptons. Therefore, they provide a lower limit for the mass of charged vectorlike leptons [58]: Similarly, the LHC constraint on the vectorlike lepton masses can be derived from the data of the slepton searches at the ATLAS and CMS experiments at 13 TeV [59]. According to that, the vectorlike leptons must satisfy either m L O(1) TeV, or m L − m χr 60 GeV.
To explain the muon g − 2 while keeping the coupling y µ in the perturbative regime, the vectorlike leptons must be light enough. Therefore, the scenario with a small gap between m L and m χr is preferable. For m χr ∼ O(100) GeV, the condition (80) implies τ ∼ O(1). Assuming that the particle χ r is stable, the condition (80) becomes 0 < m L − m χr 60 (GeV).
The parameter region with such compressed mass spectra is subjected to a constraint from the recent analysis by the ATLAS Collaboration [60].

Constraints on m χr :
Since the mass of the scalar field χ r only appears as an independent parameter in Eq. (54), the bounds (72) induce a constraint on the parameter m χr to explain the observed muon g − 2.
In Figure 3, we show the dependence of ∆a NP µ on m χr for fixed values of other inputs, τ = 1.78, δ = 1, y µ = 3, and k = 0. In this case, the current bounds on the muon g − 2 yield the 2σ allowed range for the χ r mass to be 94 GeV m χr 157 GeV. 2 In the near future, when the measurement at the E989 experiment is completed, we expect that this range will be improved.
The projected bounds in Eq. (73) imply more severe 2σ limits for this parameter that are 98 GeV m χr 143 GeV. In the subsequent analysis regarding the muon g − 2, we choose m χr = 120 GeV as a representative value.
The strips corresponding to each of these constraints appear approximately in the elliptical For the right viable region, the corresponding limits are 208.5 × 10 −5 A bs 213.7 × 10 −5 , and −12.55 × 10 −5 B bs −4.46 × 10 −5 . From the magnitudes of these two parameters, it is clearly that all the relevant FCNC processes are very much suppressed. In Figure 7, the left regions satisfying all the considered constraints from B-meson decays are depicted for different scenarios with the kinetic mixing coefficient to be k = −0.001, 0, 0.001, and 0.002. They are shown as the yellow, red slashed, purple back-slashed, and cyan regions in the plot. We observe that the viable region shifts to the left toward smaller values of A bs , and increases its area when the kinetic mixing coefficient becomes larger. As the consequence, the windows for the parameters A bs and B bs become more relaxed for larger values of k. In Figure   8, the situation for the right regions is shown when changing the value of the kinetic mixing coefficient, namely k = −0.001, 0, 0.001, and 0.002. We observe that the lower bound of B bs in the right region (∼ 12.55 × 10 −5 ) remains almost unchanged for various values of k, while its upper bound slightly increases for larger k. Regarding to the parameter A bs , a larger value of k shifts its allowed range toward the left. Meanwhile, the width of this range becomes slightly larger when increasing k. Since A bs in the right region is about O(10) times larger than that in the left region, it enhances the cross section σ(pp → Z → µμ) via the bs/sb annihilations roughly by a factor of O(10 2 ). Taking into account the constraint in Ref. [45], we find that the right region is excluded while the left one is allowed. Constraints on (τ, δ) plane: The parameters τ and δ involve in all the considered observables. Since δ is defined by the Eq. (56), it must satisfy the theoretical condition δ −1. In Figure 9, we show how the constraints on ∆a NP τ , R K , R K * , BR(B + → K + µ + µ − ), and BR(B 0 → K * 0 µ + µ − ) at the level of 2σ affect the (τ, δ) plane in the case with a vanishing kinetic mixing coefficient and fixed values of other inputs: m χr = 120 GeV, m Z = 300 GeV, y µ = 3, g X = 3, A bs = 24.2 × 10 −5 , and B bs = −11.5 × 10 −5 . We observe that most of the interested parameter region satisfies the constraint on BR(B 0 → K * 0 µ + µ − ). The region with −0.5 δ 0.5 results in a more frequent decay of B 0 mesons into K * 0 and a pair of muons that is above the upper limit. Therefore, this region is excluded. The bounds on new physics contributions to the muon g − 2 require 0.3 τ 3.2 when δ 3.7. It is because, for such large values of δ, the second term in the squared bracket in Eq. (54) is negligible, and ∆a NP µ depends mostly on the first term. As a consequence, the bounds on ∆a NP µ specify the allowed range for τ that is almost independent of δ. However, this region with large δ is roughly excluded by the experimental data on R K * 0 .
For δ smaller than about 3.7, the dependence of the allowed range of τ on δ becomes clearer.
From this figure, beside the constraint on the muon g − 2, we see that those on the semileptonic branching ratios of the B 0 and B + mesons play an important role in determining the 2σ allowed parameter region.
In Figure 10, the regions satisfying all these current bounds is shown in the cyan color.
When taking into account the projected result from the E989 experiment, the viable parameter regions significantly reduce. These regions are depicted by the hatched areas in Figure 10. For m χr = 120 GeV, the LEP constraint in Eq. (79) is automatically satisfied when we assume the LHC constraint in Eq. (81). The latter leads to the constraint on τ that is shown as the parameter region between the two vertical red dashed lines in Figure 10. These bounds for τ from Eq. (81) are actually more severe than those expected at the E989 experiment in the near future. We see that there are two separated viable regions corresponding to positive and negative values of δ. For the chosen set of other inputs, the allowed ranges for the two parameters in the former region are 1.00 < τ 2.25 and 0.73 δ 1.72, while those ranges in the latter one are 1.82 τ 2.25 and −0.62 δ −0.53. The search for the electroweak production of supersymmetric particles with compressed mass spectra at the LHC 13 TeV [60] excludes the following range for the τ parameter given that m χr = 120 GeV: Therefore, while the above negative-δ region remains intact, Eq. (82) leads to two distinct positive-δ regions for the τ parameter: (i ) 1.00 < τ 1.02, and (ii ) 1.43 τ 2.25.
For nonzero kinetic mixing coefficients, the allowed parameter space gradually changes. In Figure 11, the allowed regions on the (τ, δ) plane are plotted for various values of the kinetic mixing coefficient, namely k = −0.001, 0, 0.001, and 0.002. Similar to Figure 10, the region between the two vertical red dashed lines satisfies the constraint in Eq. (81). Here, we see that when increasing k, the viable range for |δ| moves toward smaller values, leading to a narrower gap between the positive-δ and the negative-δ regions. It is noteworthy that the area of the positive-δ region is much larger than the negative-δ one implying that the scenario with lighter χ r is more favorable than the one with lighter χ i . Relevant to the parameter τ , when the value of k increases, the range of τ in the positive-δ region remains unchanged due to the constraint (81), while the range of τ in the negative-δ region expands to the left.
Constraints on the (g X , y µ ) plane: In Figure 12, the phenomenological constraints are plotted on the (g X , y µ ) plane in the case of no kinetic mixing and fixed values of other inputs, namely m χr = 120 GeV, m Z = 300 GeV, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, and δ = 1. Since the muon g − 2 does not depend on the coupling g X , the corresponding allowed region has the form of a horizontal band with 2.27 |y µ | 3.84. The parameter regions satisfying the constraints for other observables (R K , R K * , and the branching ratios of B + and B 0 decays) have the hyperbolic forms. This is due to the Wilson coefficients C ( )NP 9,10 that define the correlation between g X and y µ for given values of the branching ratios. The viable region determined by all of the five phenomenological constraints is shown separately as the cyan region in Figure 13. We see that the overlapping region is a thin strip in the hyperbolic form. Here, the viable range for the parameter g X is 2.35 g X 3.97. In the near future, the E989 experiment will impose a more severe constraint on the (g X , y µ ) plane which is shown in the hatched region. To explain the anomalies on the B-meson decays, the new Yukawa coupling, y µ , and the U (1) X gauge coupling, g X , are required to have large values of O(1). Taking into account the perturbation limits for these two parameters shown as the horizontal dot-dashed and the vertical dashed lines in this figure, we observe that this theoretical condition excludes a large portion of the allowed parameter region. Moreover, the perturbation condition on g X indirectly determines the lower bound for y µ , and the perturbation condition on y µ indirectly determines the lower bound for g X . When combining with the B-meson decay constraints, these two theoretical conditions lead to even more severe bounds for g X and y µ than those expected at the E989 experiment. As the result, we have 2.53 g X √ 4π, and 2.53 |y µ | √ 4π.
For the case of nonzero kinetic mixing, since the effect of the kinetic mixing on the muon g − 2 is small, the allowed range determined from the ∆a NP µ constraint remains almost intact. Although the constrained regions by the B-meson decay processes (the green, the yellow, the blue hatched and the red hatched regions) still have the hyperbolic shape as before, they slightly shift downward to the region with smaller values of |y µ | (while the values of g X are fixed) in comparison to the case of vanishing kinetic mixing. This is due to the additional contributions of the diagrams with the Z -boson exchange at the tree level to the Wilson coefficients when k = 0. As a result, the above thin strip of allowed parameter region in Figure 13 slightly changes with respect to the change of the kinetic mixing coefficient. To magnify this behavior, in Figure 14, we plot the viable parameter regions in the (g X , y µ g X ) plane with four benchmark values of the kinetic mixing coefficient k = −0.001, 0, 0.001, and 0.002. It is observed that, when increasing k, the viable parameter region shift downward indicating that smaller values of the product |y µ |g X are preferable for larger values of k. In other words, for a given U (1) X gauge coupling, g X , smaller values of |y µ | are more favored for larger values of k. In this g X |y μ | × g X Figure 14: Viable parameter regions on the (g X , y µ g X ) plane for the case of m χr = 120 GeV, m Z = 300 GeV, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, δ = 1, and various values of the kinetic mixing coefficient k = −0.001, 0, 0.001, 0.002. The dot-dashed and dashed lines correspond to the perturbation limits on y µ and g X , respectively.
figure, the perturbation conditions on y µ and g X are represented by the dot-dashed and dashed lines, respectively. The parameter region satisfying these conditions stays in between these two straight lines. Constraints on the (g X , m Z ) plane: The phenomenological constraints on the (g X , m Z ) plane are presented in Figure 15 for the case of vanishing kinetic mixing and fixed values of other parameters as y µ = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, and δ = 1. For given values of the branching ratios, the parameters g X and m Z are linearly dependent. This is because these two parameters always come in the term of . However, there is only one overlapping region satisfying all of these four constraints. Furthermore, this allowed region is much more severe than the overlapping region determined by only two constraints from R K and R K * (the blue hatched and the green regions). Therefore, additional consideration of the branching ratios of the decay processes B + → K + µ + µ − and B 0 → K * 0 µ + µ − (the red hatched and the yellow regions) is crucial. The viable parameter region from all of these four constraints is plotted separately in Figure 16. From this, we can determine the allowed range for the ratio m Z g X to be 98.2−100.8 Figure 15: Phenomenological constraints on the (g X , m Z ) plane for the case of y µ = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, δ = 1, and k = 0. GeV. Taking into account the perturbation limit for g X , we find the upper bound for m Z to be approximately 354 GeV in this case. When the gauge kinetic mixing is switched on, the linear correlation between g X and m Z is deformed due to nonzero values of k in the Wilson coefficients. In Figure 17, we demonstrate the phenomenological constraints on the (g X , m Z ) plane for k = 0.002. Here, the deformation appears when m Z ≈ 85 GeV is due to the sign flipping of the ZZ mixing angle, α Z , in Eq. (37). In the case with nonzero gauge kinetic mixing, the constraint from the muon g −2 measurement need to be considered since ∆a NP µ depends on m Z via the β term in Eq. (54). For k = 0.002, we find the lower bound for the Z -boson mass from this constraint to be m Z 0.68 GeV.
In Figure 18, we plot the viable parameter regions for various values of the kinetic mixing coefficients (k = −0.001, 0, 0.001, and 0.002) in the (g X , m Z g X ) plane. Comparing these areas, we see that, for larger values of k, the allowed region in this plane shifts upward implying that larger values of m Z are more favored for a given value of g X . Similar to Figure 16, the perturbation limit for g X is important in ruling out the region with large Z -boson mass. It sets the upper limit for m Z for fixed values of other parameters. It is seen that larger values of the kinetic mixing coefficient k lead to slightly larger upper bounds on the Z -boson mass.
Constraints on the (y µ , m Z ) plane: The phenomenological constraints on the (y µ , m Z ) plane are plotted in Figure 19 where the values of other fixed inputs are m χr = 120 GeV, g X = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 ,  Figure 15, the B-meson decay width leads to an approximately linear dependence between these two parameters (y µ and m Z ) on the plane due to the factor |yµ| m Z in the Wilson coefficients when there is no kinetic mixing. In particular, each of the constraints on R K , R K * , BR(B + → K + µ + µ − ), and BR(B 0 → K * 0 µ + µ − ) determines two allowed ranges for the values of the ratio |yµ| m Z . However, there is only one overlapping region that satisfies all of the considered 2σ bounds. This region is shown separately in Figure 20 with the cyan color. Since the allowed region is a thin strip, we find that the correlation among the two considered parameters is m Z g X ∼ 100 GeV for the given set of other inputs. From this, we can extract the limits for the Z -boson mass at 95% C.L. as 226 GeV m Z 381 GeV. In the near future, when the E989 experiment get the full data, assuming the center value of the muon g − 2 remains the same, we can expect the new limits for the Z -boson mass to narrow down to 248 GeV m Z 370 GeV as illustrated by the hatched region. Taking into account the perturbation condition for y µ , the upper bound for m Z is even more severe than the one imposed by the projected E989 result, namely, m Z 355 GeV. When the gauge kinetic mixing is introduced, the constrained regions are strongly distorted Figure 19: Phenomenological constraints on the (y µ , m Z ) plane for the case of m χr = 120 GeV, g X = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, δ = 1, and k = 0. as depicted in Figure 21 for the case of k = 0.002. On the one hand, regarding the constraint on ∆a NP µ (the pink region) we see that the new Yukawa coupling, y µ , can be small in the presence of the kinetic mixing term as the Z boson is light enough with its mass of less than O(1) GeV. On the other hand, the constraint on the branching ratio of the decay process B + → K + µ + µ − imply that m Z must be larger than about 8 GeV. Therefore, small values of |y µ | are forbidden in the case k = 0.002. The blue hatched regions (the green regions, the yellow regions) determined by the bounds on R K (R K * , BR(B 0 → K * 0 µ + µ − )) which contain two separated ranges of the ratio |yµ| m Z as in Figure 19, become connected in this case.
The red hatched regions determined by the bounds on BR(B + → K + µ + µ − ) also experience strong deformations particularly when 8 GeV m Z 90 GeV. Especially, all the constrained regions are deformed when the Z -boson mass is close to about 85 GeV where the mixing angle α Z changes its sign. The viable parameter regions on the plane (y µ , m Z yµ ) satisfying all the considered constraints are presented in Figure 22 for different values of the gauge kinetic mixing coefficient k = −0.001, 0, 0.001, and 0.002. We observe that the viable region shifts up when increasing the coefficient k. This indicates that larger values of m Z are more favored for larger values of k as y µ is fixed. Meanwhile, the allowed range for y µ is almost independent on the gauge kinetic mixing coefficient k.
The phenomenological constraints on the gauge kinetic mixing coefficient k are plotted in Figure 23 when other parameters are fixed as m χr = 120 GeV, m Z = 300 GeV g X = 3, y µ = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, and δ = 1. We observe that the constraint on ∆a NP µ only excludes large values of k being close to ±1. The R K and R K * constraints are quite severe since they rule out large portions of the possible range of k. The most stringent limits on k are given by the constraints on the branching ratios of the decay processes B + → K + µ + µ − and B 0 → K * 0 µ + µ − . There are two narrow ranges of k satisfying each of these two constraints.
This observation once again reveals how important these branching ratios are beside R K and R K * . By overlapping all the ranges of k allowed by the considered constraints, we find the viable range of the kinetic mixing coefficient to be −0.002 k 0.003 when other parameters are fixed. This range is magnified in the figure by its inset. This result shows that the upper bound for the kinetic mixing coefficient in the case m Z = 300 GeV is comparable with that derived from other experimental data in the case m Z 200 GeV as in Eq. (78). Therefore, the B-meson decays are important channels to constraint the parameter space. At the oneloop level, the kinetic mixing coefficient is estimated [61] to be at the magnitude of about O(10 −1 ) − O(10 −2 ) that is already larger than the current upper bound. Therefore, in this case, k must be nonzero at the tree level with a similar order of magnitude and an opposite sign so that the total effective value of this parameter is consistent with the experimental bounds.
In Table 2, as a demonstration, the relevant observables are calculated for four benchmark values of the gauge kinetic mixing k = −0.002, 0, 0.002, and 0.003 while other parameters are chosen the same as those in Figure 23. All of them satisfy the corresponding experimental bounds at the 2σ level. We notice that, in this table, while the values of ∆a NP µ only show the differences from the eighth significant digits when changing k, differences for the values of other observables can be seen from the third or the forth significant digits. It implies that the gauge kinetic mixing coefficient has a much smaller effect on the muon g − 2 than those on the semileptonic B decays and the violation of lepton universality. The FCNC parameters in the model, A bs and B bs , induce the tree-level contributions to the B 0 s −B 0 s mixing. Due to the constraints on the rare B-meson decays, these parameters are required to be very small of about O(10 −4 ). Therefore, such contributions to the mixing observables, ∆m s and ∆Γ s , are negligible in comparison to the SM contributions 3 that are consistent with the experimental values [63]. Beside the mixing of vectorlike and SM quarks, the gauge kinetic mixing induces additional tree-level contributions to the couplings between the Z boson and all SM flavors that are approximately proportional to k O(10 −3 ). Therefore, the Z production cross sections at hadron colliders, as well as its decay widths are slightly modified by small amounts roughly proportional to k 2 O(10 −6 ) in comparison to the case of vanishing kinetic mixing. Assuming the existence of the kinetic mixing, we have calculated the Z production cross section times branching ratio to dimuon of pp collisions at √ s = 13 TeV. It is found to be about O(1) fb that is of the same order as the experimental limit in Ref. [45]. Thus, the benchmark point is marginally acceptable. In the near future, more precise analyses at the LHC will be able to test the model. 1.05028 × 10 −7 0.76737 1.95104 × 10 −7 0.55434 Table 2: The considered observables for the case of m χr = 120 GeV, m Z = 300 GeV g X = 3, y µ = 3, A bs = 24.2 × 10 −5 , B bs = −11.5 × 10 −5 , τ = 1.78, δ = 1, and four benchmark values of the gauge kinetic mixing coefficient k.
The particle χ r in this model is stable and neutral under the SM gauge groups. It can be a candidate for dark matter. In the original BDW model, the leading contribution to the spin-independent cross section between χ r and nucleon was estimated to be σ p SI ∼ O(10 −50 ) cm 2 [31]. In our analysis, the nonzero kinetic mixing slightly enhances the chance of the spin-independent scattering by allowing additional one-loop contributions. However, since the kinetic mixing is limited to be k ∼ O(10 −3 ), such new one-loop contributions to the spinindependent cross section is suppressed. As a result, the total cross section is still smaller than the coherent neutrino-nucleus scattering background [64]. With the chosen parameter sets, the pair annihilation process of χ r into a pair of leptons (µ or ν µ ) is effective due to the large y µ coupling. In addition, the coannihilation between the vectorlike lepton and χ r also reduces its relic density since their masses are relatively close. Therefore, the relic density of χ r is smaller than the observed dark matter abundance as obtained in Ref. [31]. Once the kinetic mixing is switched on, the above pair annihilation and coannihilation processes happen more frequently in the early universe leading to a smaller value of Ω χr h 2 than that in the case of vanishing kinetic mixing. In order to account for the dark matter relic density measured by the Planck Collaboration [65], an additional dark matter candidate is necessary.

Conclusion
The BDW model with additional vectorlike lepton and quark doublets and two complex scalars charged under the U (1) X gauge group is well motivated due to its ability in explaining various anomalies at the same time. In this paper, we have generalized this model by introducing the gauge kinetic mixing term. The new physics contributions to the muon anomalous magnetic moment and the Wilson coefficients (C ( ) 9,10 ) have been calculated analytically. We have investigated the parameter space of the model taking into account the phenomenological constraint on the muon g − 2, the updated LHCb data on lepton universality violation (R K and R K * ), the branching ratios of the semileptonic rare decays (B + → K + µ + µ − and B 0 → K 0 µ + µ − ), the LEP data on slepton searhces, and the LHC 13 TeV data on both slepton and Z -boson searches, as well as requirements from the perturbative theory. The viable parameter regions satisfying all the considered constraints at the level of 2σ have been identified. The results indicate that the FCNC parameters, A bs and B bs , are small enough to be consistent with experiment data. In the presence of the gauge kinetic mixing term, the allowed regions are shifted and deformed in comparison to the case with k = 0. Especially, the kinetic mixing coefficient plays an important role in extending the viable parameter regions. The analysis also shows that the rare B-meson decays are important channels that provide important constraint on the gauge kinetic mixing beside the Z resonance searches. In the near future, with the projected sensitivity, the E989 experiment will be able to test certain parts of the free parameter space and put a more severe constraint on the acceptable parameter regions of the model. Sine the Z production cross section times branching ratio to dimuon at the LHC is of the same order as the current limit for certain parameter regions, the model can be tested in more precise analyses of this channel.