Connecting $b\to s\ell\bar\ell$ anomalies to enhanced rare nonleptonic $\bar{B}_s^0$ decays in $Z'$ model

The present data on a number of observables in $b\to s\mu^+\mu^-$ processes manifest some tensions with the standard model (SM). Assuming that these anomalies have a new physics origin, we consider the possibility that a $Z'$ boson is responsible for them. We further assume that its interactions with quarks also affect rare nonleptonic decays of the $\bar B_s^0$ meson which are purely isospin-violating and tend to be dominated by electroweak-penguin contributions, namely $\bar B_s^0\to(\eta,\eta',\phi)(\pi^0,\rho^0)$. Most of these decays are not yet observed, and their rates are expected to be relatively small in the SM. Taking into account constraints from various measurements, including the evidence for $\bar B_s^0\to\phi\rho^0$ recently seen by LHCb, we find that the $Z'$ effects on $\bar B_s^0\to(\eta,\phi)\pi^0$ can make their rates bigger than the SM predictions by up to an order of magnitude. For $\bar B_s^0\to\eta'\pi^0,(\eta,\eta')\rho^0$, the enhancement factors are at most a few. Since the $Z'$ contributions to the different channels depend on different combinations of its couplings, observations of more of these decays in future experiments, along with improved $b\to s\mu^+\mu^-$ data, will probe this $Z'$ scenario more thoroughly.


I. INTRODUCTION
The latest measurements of various b → sµ + µ − processes have turned up some intriguing discrepancies from the expectations of the standard model (SM) of particle physics. Specifically, the LHCb Collaboration in its angular analysis of the decay B 0 → K * 0 µ + µ − found tensions with the SM at the 3.4σ level [1]. This was later confirmed in the Belle experiment on the same process, but with lower statistical confidence [2]. Furthermore, LHCb findings [3,4] on the ratio R K of the branching fractions of B + → K + µ + µ − and B + → K + e + e − decays and on the corresponding ratio R K * for B 0 → K * 0 µ + µ − and B 0 → K * 0 e + e − decays are all below their SM predictions [5][6][7] by 2.1σ to 2.6σ. In addition, the current data [8][9][10] on the branching fractions of B → K ( * ) µ + µ − and B s → φµ + µ − favor values less than their SM estimates.
Although the statistical significance of the aforesaid deviations from SM expectations is still too low for a definite conclusion, they may be early clues about interactions beyond the SM in b → s transitions. Recent model-independent theoretical analyses have in fact demonstrated that new physics (NP) could account for these anomalies [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. In view of the possibility that these tentative hints of NP will be confirmed by upcoming experiments, it is of interest to explore the potential implications for other b → s processes.
Among them are the nonleptonic decaysB s → ηπ 0 ,B s → η π 0 ,B s → φπ 0 ,B s → ηρ 0 , B s → η ρ 0 , andB s → φρ 0 . These transitions fully break isospin symmetry, implying that their amplitudes receive no contributions from QCD-penguin operators and arise instead from charmless tree and electroweak-penguin operators [26,27]. The product of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements in the tree contributions is suppressed compared to that in the electroweak-penguin ones, and the suppression factor is |V us V ub |/|V ts V tb | ∼ 0.02. Consequently, although the Wilson coefficients in the former are much bigger, the latter turn out to dominate these channels and lead to relatively low rates [26][27][28][29]. Most of them are not yet observed, the exception beingB s → φρ 0 . Evidence for the latter decay was detected by LHCb last year [30] with a branching fraction of (0.27 ± 0.08)×10 −6 [10], which agrees with some of its estimates in the SM within sizable errors [31][32][33][34][35].
The smallness of the rates ofB s → (η, η , φ)(π 0 , ρ 0 ) in the SM implies that they may serve as probes of physics beyond it. This has been considered to varying extents in the contexts of different models [34][35][36][37][38]. In this paper, we treat these rare nonleptonicB s decays along similar lines and suppose that the NP influencing them also causes the aforementioned b → sµ + µ − anomalies. We adopt in particular a scenario where an electrically neutral and uncolored spinone particle, the Z boson, is responsible for the new interactions in these two sets of b → s transitions. We assume that it couples nonuniversally to SM fermions and does not mix with SM gauge bosons, but it is not necessarily a gauge boson and could even be composite.
The rest of the paper is organized as follows. In Sec. II, we address the Z contributions to b → sµ + µ − and take into account constraints from the relevant empirical information, including that on B s -B s mixing. In Sec. III, we examine the impact of the Z interactions with SM quarks onB s → (η, η , φ)(π 0 , ρ 0 ). To evaluate their amplitudes, we employ the soft collinear effective theory [39][40][41][42][43][44][45]. We show that in the Z presence the rates ofB s → η π 0 , (η, η )ρ 0 can increase by as much as factors of a few with respect to the predictions in the SM, while the rates of B s → (η, φ)π 0 can exceed their SM values by up to an order of magnitude. We make our conclusions in Sec. IV.
Global analyses [18,19] have found that some of the best fits to the most recent anomalous b → sµ + µ − measurements result from effective interactions given by where C 9µ = C sm 9 + C np 9µ and C 9 µ = C np 9 µ are the Wilson coefficients, α e = 1/133 is the fine structure constant at the b-quark mass (m b ) scale, and G F is the Fermi constant. The same SM part C sm 9 occurs in the electron and tau channels b → s(e + e − , τ + τ − ), but they are not affected the NP. In Fig. 1, for later use, we display the 2σ (cyan) region of C np 9 µ versus C np 9µ permitted by the data, from the global fit carried out in Ref. [18].
In the literature, many models possessing some kind of Z particle with different sets of fermionic couplings have been studied in relation to the b → sµ + µ − anomalies . In the Z scenario considered here, the interactions responsible for C np 9µ,9 µ in Eq. (1) are described by where the constants ∆ sb L,R are generally complex and ∆ µµ V is real due to the Hermiticity of L Z . Any other possible Z couplings to leptons are taken to be negligible. To simplify the analysis, hereafter we focus on the special case in which where ρ L,R are real numbers, and so they do not supply any new CP -violation phase. Accordingly, for a heavy Z with mass m Z we obtain In Fig. 2, we illustrate the ranges of ρ L and ρ R corresponding to the allowed C np 9 µ -C np 9µ (cyan) region in Fig. 1 for m Z = 1 TeV and some sample choices of ∆ µµ V , namely ±0.03 (red), ±0.05 (orange), ±0.1 (yellow), and ±0.3 (green). We note that these ∆ µµ V values contribute positively to the SM muon anomalous magnetic moment, but with m Z = 1 TeV are too small to explain the disagreement with its measurement [84].
The Z couplings in Eq. (2) also induce tree-level effects on ∆M s = 2|M s 12 |, which pertain to B s -B s mixing and has been measured to be ∆M exp s = (17.757 ± 0.021)/ps [10]. We can express the sum of the SM and Z contributions as [85] M s 12 = M s,sm where [85]r = 0.985 for m Z = 1 TeV is a QCD factor, g 2 sm = 1.7814 × 10 −7 GeV −2 , the SM loop function S 0 2.35 for a top-quark mass m t 165 GeV, and To apply restrictions on ρ L,R from the B s -B s mixing data, we impose which is the 95% confidence level (CL) range from the latest UTfit global analysis [86]. Since some of the numbers quoted in the last paragraph have uncertainties up to a few percent, we let κ LR vary by up to 10% from its central value when scanning the parameter space for ρ L,R values which conform to Eq. (7). For m Z = 1 TeV, we incorporate the scan result into Fig. 2, represented by the blue area. Thus, in this figure each overlap of the blue area with one of the other colored ones of a particular ∆ µµ V value corresponds to the parameter space that can explain the b → sµ + µ − anomalies and simultaneously satisfies Eq. (7). With smaller choices of |∆ µµ V |, such overlaps could be found at larger |ρ L,R | values. This graph also reveals that in the absence of the right-handed coupling, ρ R = 0, the allowed range of ρ L would be rather narrow, indicating the importance of nonvanishing ρ R for gaining bigger viable parameter space [48].

III. Z CONTRIBUTIONS TO RARE NONLEPTONICB s DECAYS
Given thatB s → (η, η , φ)(π 0 , ρ 0 ) change both strangeness and isospin, in the SM their amplitudes proceed from b → s four-quark operators O u 1,2 and O 7,8,9,10 which are derived from charmless tree and electroweak-penguin diagrams, respectively. In contrast, the QCD-penguin operators O 3,4,5,6 , which preserve isospin symmetry, do not affect these processes. 1 In many models beyond the SM, new interactions may modify the Wilson coefficients C i of O i and/or give rise to extra operatorsÕ i which are the chirality-flipped counterparts of O i . A flavor-violating Z boson may contribute to some of them, depending on the details of its properties.
In our scenario of interest, besides its couplings in Eq. (2), the Z has flavor-conserving interactions with the u and d quarks via the constants ∆ uu,dd L,R being real, but does not couple flavor-diagonally to other quarks. From Eqs. (2) and (8), we can derive tree-level Z -mediated diagrams contributing to nonleptonic b → s reactions. For a heavy Z , these diagrams yield after applying Eq. (3). It is straightforward to realize that these bring about modifications to the coefficients of the QCD-and electroweak-penguin operators O 3,5,7,9 in the SM and also generate their chirality-flipped partnersÕ 3,5,7,9 [87]. We can express them in the effective Lagrangian for b → s transitions as where where As O 3,5 andÕ 3,5 do not break isospin, only C Z 7,9 andC Z 7,9 contribute toB s → (η, η , φ)(π 0 , ρ 0 ). To estimate the Z impact on these decays, we make use of the soft collinear effective theory (SCET) [39][40][41][42][43][44][45], similarly to what was done in Ref. [38] in the case of a leptophobic-Z model. For any one of them, the SCET amplitude at leading order in α s (m b ) can be written as [45] AB where f M is the decay constant of meson M , the ζ's are nonperturbative hadronic parameters which can be fixed from experiment, the T 's are hard kernels which are functions of the Wilson coefficients C i andC i , and φ M (ν) is the light-cone distribution amplitude of M which is normalized as 1 0 dν φ M (ν) = 1. The so-called charming-penguin term, which in this case conserves isospin, is absent from AB s→M1M2 . The hard kernels for the decays of concern are available from the literature [33,44,45] and have been listed in Table I, where the flavor states η q ∼ uū + dd / √ 2 and η s ∼ ss are related to the physical meson states η and η by η = η q cos θ − η s sin θ and η = η q sin θ + η s cos θ with mixing angle θ = 39.3 • [44,45,89].
In the presence of NP which also generates the extra operatorsÕ i , the quantities c 2,3 and b 2,3 in Table I depend not only on C i andC i , but also on the final mesons M 1 and M 2 , as well as on the CKM factors λ t and λ u = V * us V ub . The dependence on M 1 and M 2 arises from the fact that, with regard to the nonzero kernels in this table, for each 4-quark operator theB s → M 1 and vacuum → M 2 matrix elements and their contraction in the amplitude can lead to an overall negative or positive sign for the contribution of the operator, the sign being determined by the chirality combination of the operator and by whether the final mesons are pseudoscalars (P P ), vectors (V V ), P V , or V P . Thus, forB s → (η q , η s )π 0 andB s → φρ 0 we have 2 where N c = 3 is the color number and b 2,3 , which are contained in T 2J (ν) and T 2Jg (ν), are also functions of ν via [45] ω 2 = νm Bs and ω 3 = (ν − 1)m Bs . However, forB s → (η q , η s )ρ 0 and B s → φπ 0 we need to make the sign change −C i → +C i in c 2,3 and b 2,3 .
For numerical computation of AB s→M1M2 , in view of Table I, the meson decay constants which we need are only f π = 131 MeV and f ρ = 209 MeV, and the integral in Eq. (13) can be treated with the aid of the relations M for M = π, ρ, in which cases χ −1 π = 3.3 and χ −1 ρ = 3.45 [44,45]. Moreover, for the ζ's we adopt the two solutions derived from the fit to data done in Ref. [45]: From these, we can obtain ζ Before addressing the Z influence onB s → (η, η , φ)(π 0 , ρ 0 ), we provide the SM predictions for their branching fractions, which are collected in Table II. For the first five modes, the SCET numbers have been evaluated with the preceding formulas and parameter values, and the last two columns correspond to the two solutions of SCET parameters in Eq. (15). For the sixth (φρ 0 ) mode, the SCET entry has been computed with the CKM and SCET parameters supplied very recently in Ref. [33]. The two errors in each of the SCET predictions are due to flavor-SU(3)-breaking effects which we have assumed to be 20% and due to the errors in the ζs from the fits to data, respectively, the latter errors being given in Refs. [33,45]. The SCET numbers forB s → (η, η )ρ 0 , φπ 0 B s → φρ 0 are close to the corresponding ones determined in The errors of the SCET predictions are due to assumed 20% flavor-SU(3)-breaking effects and the errors in the ζs from fits to data, respectively. For comparison, the second and third columns contain results calculated in the frameworks of QCDF [31] and PQCD [32].
Ref. [45] ( [33]). 3 For comparison, in the second and third columns we quote numbers calculated with QCD factorization (QCDF) [31] and perturbative QCD (PQCD) [32]. Evidently, these two methods produce results comparable to those of SCET, especially with its Solution 2 in the case of the first five modes, considering the errors in the predictions. The entries forB s → φρ 0 are also compatible with the new measurement B B s → φρ 0 exp = (0.27 ± 0.08)×10 −6 [10] mentioned earlier. An important implication of what we see in this table is that NP would not be easily noticeable in the rates of these decays unless it could enhance them by more than a factor of 2. This possibility may be unlikely to be realized in the case ofB s → φρ 0 which has been detected having a rate consistent with SM expectations. Nevertheless, as we demonstrate below, substantial enhancement can still occur in some of the other channels. Now we include the Z contributions from Eq. (11) in order to examine their impact on these decays. As Table II indicates that the predictions of the SCET Solution 1 forB s → η π 0 , ρ 0 are comparatively quite suppressed, from this point on we employ only Solution 2 parameters in our treatment ofB s → (η, η )(π 0 , ρ 0 ), φπ 0 . Thus, summing the SM and Z terms for m Z = 1 TeV, with the central values of the input parameters, we find the amplitudes (in units of GeV) for the π 0 channels to be and for the ρ 0 channels 10 9 AB s→ηρ 0 2.56 + 0.77i + (6.32 − 0.12i)(ρ L + ρ R )(δ L + δ R ) , where δ L,R are defined in Eq. (12).
We notice that the amplitudes in Eqs. (16) and (17) do not all have the same dependence on ρ L,R and δ L,R . Therefore, although B B s → φρ 0 exp implies a restraint on the values of (ρ L − ρ R )(δ L + δ R ) in AB s→φρ 0 , the amplitudes for the other channels, which have different combinations of ρ L,R and δ L,R , may generally still be altered considerably with respect to their SM parts. However, in our particular Z case ρ L,R must satisfy ρ R ∼ 0.1 ρ L , as can be inferred from Fig. 2. Hence, based on Eq. (17), we may expect that the amplitudes forB s → (η, η )ρ 0 do not deviate hugely from their SM values. To look into this more concretely, for definiteness we take ρ R = 0.1 ρ L and impose 3 The SCET predictions in Table II differ from those obtained in [38] because some of the input parameters used in our two papers are not the same.
which is the 2σ range of B B s → φρ 0 exp . From the allowed values of the product ρ L (δ L + δ R ) we can assess how much the branching fractions ofB s → (η, η )ρ 0 are modified compared to the central values of their respective SM predictions in Table II under SCET Solution 2. We show the results in Fig. 3, which also depicts Eq. (18) relative to the SM prediction. We further find that the ranges ρ L (δ L + δ R ) ∈ [−0.99, −0.71] and [−0.23, 0.048] fulfill Eq. (18). Within these ranges, represented by the horizontal portions of the unshaded areas in this figure, we learn that B B s → ηρ 0 (red solid curve) can reach up to ∼ 2.7 times its SM value, whereas forB s → η ρ 0 (blue solid curve) the enhancement is at most about 1.9 times. 4 ForB s → (η, η , φ)π 0 , the Z -induced terms in Eq. (16) are proportional to δ L −δ R . Therefore, these channels are not subject to the condition in Eq. (18), and their amplitudes may be affected by the Z contributions more than their ρ 0 counterparts. To examine this more quantitatively, we set ρ R = 0.1 ρ L as in the previous paragraph and subsequently compute the branching fractions of these π 0 modes for −1 ≤ ρ L (δ L − δ R ) ≤ 1. In Fig. 4 we present the results divided by the central values of their respective SM predictions in Table II under SCET Solution 2. We observe that over most of the ρ L (δ L − δ R ) > 0 region covered in this plot the Z effects can cause the branching fractions ofB s → (η, φ)π 0 to exceed their SM counterparts by at least a factor of 2 and up to about an order of magnitude. Moreover, theB s → (η, φ)π 0 rates tend to be enhanced together with roughly similar enlargement factors. One also notices that the Z impact could instead bring about substantial reduction of their rates. In Table III Table II under SCET Solution 2, versus the product ρ L (δ L − δ R ) in the case where ρ R = 0.1 ρ L and m Z = 1 TeV. It is worth remarking that the Z -generated coefficients C Z 3,5,7,9 andC Z 3,5,7,9 in Eq. (11) enter the amplitudes for nonleptonic b → s decays which are not dominated by the contributions of the electroweak penguin operators, such as B s → K ( * )K ( * ) and B → πK ( * ) . Since these transitions have been observed, their data imply restrictions on the size of ρ L (δ L ± δ R ), as the requirement ρ R ∼ 0.1 ρ L in our Z scenario implies that the role ofC Z 3,5,7,9 is minor. Our choices |ρ L (δ L ± δ R )| ≤ 1 above correspond to C Z 7 ± C Z 9 ≤ 0.0202 2|C sm 9 |, where C sm 9 = −0.0103 as quoted before. We have checked that the changes to the rates of those decays due to |C Z 7,9 | |C sm 9 | are less than the uncertainties of the SCET estimates in the SM, which are typically around 20% to 40% [33,44,45]. As for the influence of C Z 3,5 , it can be minimized by adjusting the extra free parameters ∆ dd L,R in Eq. (11). 5 For comparison, earlier studies [34,36,37] concerning potential NP inB s → φ(π 0 , ρ 0 ) pointed out that rate enhancement factors of a few to an order of magnitude could still occur in the φπ 0 mode and that |C np j /C sm 9 | 2 was not yet disfavored.
Finally, we would like to mention that the Z coupling parameters of interest are separately consistent with constraints which may be pertinent from collider measurements. We illustrate this with the specific examples in Table IV for different sets of ρ L (δ L ± δ R ) and ∆ µµ V values in the aforesaid case where ρ R = 0.1 ρ L and m Z = 1 TeV. The choice ρ L , ∆ µµ V = (0.8, 0.03) is evidently within the region covered in Fig. 2, while ρ L , ∆ µµ V = (1.1, 1.2), 0.02 lie in the extension thereof. In this table, the displayed numbers for δ L,R can comfortably comply with the condition |δ L,R | ≤ 1.0 1 + (1.3 TeV) 2 /m 2 Z m Z /(3 TeV) inferred in Ref. [90] from the study on LHC bounds in Ref. [91]. Furthermore, based on the results of Refs. [75,92], the selected ∆ µµ V values are compatible with LEP data on Z-boson decays into lepton pairs [10]. Lastly, as theēeZ interaction is supposed to be vanishing, restraints from LEP II measurements on e + e − → ff can be evaded.  The quark-Z coupling parameters ρ L and δ L,R corresponding to a few sample sets of ρ L (δ L ± δ R ) and ∆ µµ V in the ρ R = 0.1 ρ L and m Z = 1 TeV case.

IV. CONCLUSIONS
We have explored the possibility that the recently observed anomalies in several b → sµ + µ − processes are attributable to the interactions of a Z boson which also contribute to rare nonleptonic decays of theB s meson, namelyB s → (η, η , φ)(π 0 , ρ 0 ). Given that these purely isospinviolating decays are controlled mainly by the electroweak-penguin operators, their decay rates are expected to be relatively small in the SM, making these modes potentially sensitive to signals beyond the SM. The Z couplings are subject to various restrictions, particularly from the data on B s -B s mixing and the new experimental finding onB 0 s → φρ 0 , besides the measurements of b → sµ + µ − transitions. We showed that, within the allowed parameter space, the Z impact onB 0 s → (η, φ)π 0 can cause their rates to grow up to an order of magnitude greater than their expectations in the SM. On the other hand, the enhancement factors forB 0 s → η π 0 , (η, η )ρ 0 are at most a few. The different enlargement factors of these different channels depend not only on the combinations of the Z couplings occurring in their amplitudes, but also on how the SM and Z terms in the amplitudes interfere with each other. Therefore, the observations of more of these decays in future experiments, together with improved upcoming data on b → sµ + µ − , will test our Z model more comprehensively.