Rare nonleptonic $\bar{B}_s^0$ decays as probes of new physics behind $b\to s\mu\bar\mu$ anomalies

The anomalous results of recent measurements on various $b\to s\mu^+\mu^-$ processes could be initial evidence of physics beyond the standard model (SM). Assuming this to be the case, we entertain the possibility that the underlying new physics also affects the rare nonleptonic decays of the $\bar B_s^0$ meson. We consider in particular new physics arising from the interactions of a heavy $Z'$ boson and investigate their influence on the decay modes $\bar B_s^0\to(\eta,\eta',\phi)\omega$, which receive sizable QCD- and electroweak-penguin contributions. These decays are not yet observed, and their rates are estimated to be relatively small in the SM. Taking into account the pertinent constraints, we find that the $Z'$ effects can greatly increase the rates of $\bar B_s^0\to(\eta,\phi)\omega$, by as much as two orders of magnitude, with respect to the SM expectations. We have previously shown that $\bar B_s^0\to(\eta,\phi)\pi^0$, with similarly suppressed SM rates, could also undergo substantial $Z'$-induced enhancement. These rare modes can therefore serve as complementary probes of the potential new physics which may be responsible for the $b\to s\mu^+\mu^-$ anomalies.


I. INTRODUCTION
The current data on a number of b → sµ + µ − transitions have manifested several tantalizing deviations from the expectations of the standard model (SM). Specifically, the LHCb Collaboration [1] found moderate tensions with the SM in an angular analysis of the decay B 0 → K * 0 µ + µ − , which were later corroborated in the Belle experiment [2]. Moreover, LHCb reported [3] that the ratio R K of the branching fractions of B + → K + µ + µ − and B + → K + e + e − decays and the corresponding ratio R K * for B 0 → K * 0 µ + µ − and B 0 → K * 0 e + e − decays are a couple of sigmas below their SM predictions [4]. Also, the existing measurements [5,6] on the branching fractions of B → K ( * ) µ + µ − and B s → φµ + µ − favor values below their SM estimates.
These anomalies may be harbingers of physics beyond the SM, although their statistical significance is still insufficient for drawing a definite conclusion. In fact, model-independent theoretical studies have pointed out that new physics (NP) could explain them [7,8]. This would suggest that they might be experimentally confirmed in the near future to have originated from beyond the SM. Thus, it seems timely to explore what might happen if the same underlying NP could have an appreciable influence on some other b → s processes.
Previously we have entertained such a possibility in a scenario where a new, electrically neutral and uncolored, spin-one particle, the Z boson, is behind the b → sµ + µ − anomalies [9]. In particular, we investigated the potential implications for the nonleptonic decays of theB 0 s meson which are purely isospin-violating, namelyB 0 s → (η, η , φ)(π 0 , ρ 0 ), most of which are not yet observed [6]. In the SM limit, they are not affected by QCD-penguin operators, which conserve isospin, while the effects of tree operators are suppressed by a factor |V us V ub |/|V ts V tb | ∼ 0.02 involving Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, and thus the amplitudes for these decays tend to be dominated by electroweak-penguin contributions [10]. Accordingly, their rates in the SM are comparatively small [10][11][12][13][14][15][16][17][18][19], which motivated earlier works suggesting that one or more of these decay modes could be sensitive to NP signals [19][20][21][22][23][24]. Incorporating the relevant constraints, we demonstrated in Ref. [9] that the Z influence could cause the rates of two of the modes,B 0 s → (η, φ)π 0 , to rise by up to an order of magnitude above their SM expectations. It follows that these modes could offer valuable complementary information about the NP which may be responsible for the b → sµ + µ − anomalies.
Extending our preceding analysis, the present paper coversB s → (η, η , φ)ω, which are also not yet observed [6]. Here, as in theB 0 s → (η, η , φ)(π 0 , ρ 0 ) case, the tree operators suffer from the CKM suppression, again allowing the penguin operators to become important. However, unlike the latter modes,B s → (η, η , φ)ω preserve isospin and therefore receive both electroweak-and QCD-penguin contributions. In the SM, the rates of these decays turn out to be relatively small as well [12][13][14][15][16][17], and so they could be expected to serve as additional probes of the potential NP behind the anomalies. We will show that this can indeed be realized in the aforementioned Z model, especially for the two modesB s → (η, φ)ω.
The remainder of the paper is organized as follows. In Sec. II we describe the Z interactions which impact the various processes of concern. In Sec. III we address how the Z -induced effects on the considered rare nonleptonic decays could raise some of their rates significantly. We will impose appropriate restraints on the Z couplings, including from other b → s transitions, such asB s → φρ 0 and B → πK decays. 1 Our numerical work will also involveB s → (η, φ)π 0 , which we investigated before, to see if there might be any correlation between their rate enlargement and that ofB s → (η, φ)ω. In Sec. IV we give our conclusions. An Appendix contains extra formulas.

II. Z INTERACTIONS
In our Z scenario of interest, the mass eigenstates of the u, d, s, and b quarks have nonstandard interactions described by [9] where the constants ∆ sb L,R are generally complex, while ∆ µµ V and ∆ uu,dd L,R are real due to the Hermiticity of L Z , and P L,R = (1 ∓ γ 5 )/2. As in Ref. [9], we suppose that any other possible couplings of the Z to SM fermions are negligible and that it does not mix with SM gauge bosons but is not necessarily a gauge boson. 2 Moreover, for simplicity we concentrate on the special case in which where ρ L,R are real numbers.
For the Z being heavy, the couplings to bs and µμ in Eq. (1) contribute to the effective Lagrangian where α e and G F are the usual fine-structure and Fermi constants, and C 9µ = C sm 9 + C np 9µ and C 9 µ = C np 9 µ are the Wilson coefficients, with C sm 9 being the flavor-universal SM part ( = e, µ, τ ) and [9] According to model-independent analyses [7], one of the best fits to the anomalous b → sµ + µ − data corresponds to C np 9µ ∼ −1.1 and C np 9 µ ∼ 0.4, with no NP in b → se + e − . The bsZ couplings in L Z above also affect B s -B s mixing at tree level and hence need to satisfy the restrictions inferred from its data. As elaborated in Ref. [9], the requirements from b → sµ + µ − processes and B s -B s mixing together imply that the left-handed bsZ coupling must be roughly ten times stronger than the right-handed one, and so ρ L ∼ 10ρ R . This will be taken into account later on.
In the SM, their amplitudes proceed from b → s four-quark operators O u 1,2 , O 3,4,5,6 , and O 7,8,9,10 derived from charmless tree, QCD-penguin, and electroweak-penguin diagrams, respectively. 3 In models beyond the SM, new ingredients may alter the Wilson coefficients C j of O j and/or generate extra operatorsÕ j which are the chirality-flipped counterparts of O j . In our Z case, at the W -mass (m W ) scale only C 3,5,7,9 andC 3,5,7,9 get Z contributions given by [20][21][22]27] where 3,5,7,9 are the Wilson coefficients with [9] C Z 3,5 = we have assumed that renormalization group evolution (RGE) between the m Z and m W scales can be neglected, and e u = −2e d = 2/3. At the b-quark mass (m b ) scale, all the penguin coefficients acquire Z terms via RGE, which we treat in the next section.

III. Z EFFECTS ON RARE NONLEPTONICB s DECAYS
To estimate the Z impact onB s → (η, η , φ) ω, following Ref. [9] we employ the soft-collinear effective theory (SCET) [16][17][18]28]. For any one of them, we can write the SCET amplitude at leading order in the strong coupling α s (m b ) as [17] AB where f M denotes the decay constant of meson M , the ζs are nonperturbative hadronic parameters extractable from experiment, the T s represent hard kernels containing the Wilson coefficients C j andC j at the m b scale, and φ M (ν) is the light-cone distribution amplitude of M normalized as 1 0 dν φ M (ν) = 1. We collect the hard kernels, from Refs. [16][17][18], in Table I, where the flavor states η q ∼ uū + dd / √ 2 and η s ∼ ss are linked to the physical meson states η and η by η = η q cos θ − η s sin θ and η = η q sin θ + η s cos θ with mixing angle θ = 39.3 • [17,18,29]. We note that the so-called charming-penguin contribution is absent from AB s→M1M2 , which is one of the reasons why these decays have low rates [16,17].
In the presence of NP which also gives rise toÕ j , the quantities c k and b k in Table I depend not only on C j andC j but also on the final mesons M 1,2 besides the CKM factors λ u,t . The dependence on M 1,2 is due to the fact that, in view of the nonzero kernels in this table, for each 4-quark operator the contraction of theB s → M 1 and vacuum → M 2 matrix elements in the amplitude can lead to an overall negative or positive sign for the contribution of the operator, the sign being fixed by the chirality combination of the operator and by whether M 1,2 are pseudoscalars (P P ), vectors (V V ), or P V . Thus, forB s → P P andB s → φω we have where N c = 3 is the color number, C − j = C j −C j , and b 2,3,5,6 , which enter T 2J,2Jg (ν), are also functions of ν because ω 2 = νm Bs and ω 3 = (ν − 1)m Bs [17]. However, forB s → (η q , η s )ω we need to make the sign change C − j → C + j = C j +C j in c 2,3,5,6 and b 2,3,5,6 .
The formulas in Eq. (9) generalize the SM ones from Refs. [16][17][18], which also provide the C sm j values at the m b scale, C sm 1,2 = (1.11, −0.253) and C sm 7,8,9,10 = (0.09, 0.24, −10.3, 2.2) × 10 −3 , calculated at leading-logarithm order in the naive dimensional regularization scheme [30] with the prescription of Ref. [31]. We will incorporate these numbers into c k and b k . The Z -generated coefficients in Eq. (6) contribute to Eq. (9) 3,5,7,9 due to RGE from the m W scale to the m b scale and δC i are analogously related toC Z 3,5,7,9 . To evaluate AB s→M1M2 , in light of Table I, we employ the decay constant f ω = 192 MeV [16] and treat the integral in Eq. (8) with the aid of [17,18]. Furthermore, for the ζ's in AB s→(ηq ,ηs)ω , we adopt the two solutions from the fit to data performed in Ref. [17]: It is worth remarking here that, since they were the outcome of fitting to B → P P , P V data with SM Wilson coefficients [17], using these solutions in an investigation of NP is justifiable provided that the impact of the NP on the channels which dominate the fit is small compared the SM contribution. This requisite can be met by our Z scenario, as we set the mass of the Z to be O(1 TeV) and ensure that its couplings comply with the various constraints described in this study. From Eq. (10), we then have ζ Bηq,s (J) = ζ P (J) and ζ Bηq,s (J)g = ζ (J)g , assuming flavor-SU(3) symmetry [17]. Other input parameters are CKM matrix elements from Ref. [32] as well as the meson masses m η = 547.862, m η = 957.78, m ω = 782.65, m φ = 1019.461, and m Bs = 5366.89, all in units of MeV, and the B s mean lifetime τ Bs = 1.509 × 10 −12 s, which are their central values from Ref. [6]. For the third (φω) mode, we choose the CKM and SCET parameters supplied recently in Ref. [16].
Before dealing with the Z influence onB s → (η, η , φ)ω numerically, with the above SCET prescription we arrive at their SM branching fractions, listed in Table II. ForB s → (η, η )ω, the entries in the last two columns correspond to the two solutions of SCET parameters in Eq. (10). The central values of the SCET predictions agree with those in Refs. [16,17] (10). The second and third columns exhibit numbers computed with QCDF [12] and PQCD [14,15].
have added the errors shown. For comparison, in the second and third columns we quote results from the QCD factorization (QCDF) [12] and perturbative QCD (PQCD) [14,15] approaches. Evidently, the SCET numbers can be compatible with the QCDF and PQCD ones within the sizable errors. One concludes that for NP to be unambiguously noticeable in the rates it would have to amplify them relative to their SM ranges by significantly more than a factor of two.
In the presence of the Z contributions, we find the changes of the Wilson coefficients C j at the m b scale to be where in each coefficient we have kept only the Z term with the biggest numerical factor, upon applying the RGE at leading-logarithm order [30] with α e = 1/128, α s (m Z ) = 0.119, m b = 4.8 GeV, and m t = 174.3 GeV [18]. Furthermore, at the m b scaleC j = δC j are related tõ C Z 3,5,7,9 in an analogous manner. Combining the SM and Z portions, for m Z = 1 TeV and the central values of the input parameters we derive the amplitudes forB s → (η, η , φ) ω to be, in units of 10 −9 GeV, where the superscripts (1) and (2) refer to SCET Solutions 1 and 2, respectively, Z terms with numerical factors below 0.005 in size are not displayed, and Given that δ ± and ρ ± participate in the amplitudes forB 0 s → (η, η , φ)(π 0 , ρ 0 ) as well [9], as Eqs. (A1)-(A3) in the Appendix show, it is germane to include them in this analysis. What is more, as discussed in Ref. [9], the LHCb finding B B s → φρ 0 exp = (0.27 ± 0.08) × 10 −6 [6,33], which is in line with some of its SM estimates albeit within large errors [12-14, 16, 19, 21], translates into an important constraint on the Z couplings. Additionally, treating all these rare decays at the same time would allow us to see if there might be correlations among their rate increases/decreases compared to the SM expectations. Such correlations would constitute Z predictions potentially testable in upcoming experiments.
The couplings in Eq. (13) also affect other nonleptonic b → s processes which have been observed and hence need to respect the restrictions implied by their data. Here we focus on the well-measured decays B − → π 0 K − , π −K 0 andB 0 → π 0K 0 , π + K − plus their antiparticle counterparts. Their rates in the SM, with ∼ 40% errors [18], agree with their measurements [6]. Incorporating the Z terms, we have calculated the B → πK amplitudes and collected their expressions in Eqs. (A6) and (A7).
To illustrate how the Z interactions contribute to the decays of interest, we put together 5,000 randomly generated benchmarks fulfilling the following conditions. We impose which is the 2σ range of B B s → φρ 0 exp . For the B → πK requirement, since the SM predictions have uncertainties of around 40% and are compatible with their data, we demand that the Z effects alter the B → πK rates relative to their SM values by no more than 20%. 4 For the Z parameters, we select ρ R = 0.1 ρ L as in Ref. [9] and let the products of ρ L and the other nonleptonic couplings vary within the intervals having already set m Z = 1 TeV in Eq. (12) and in the other amplitudes written down in the Appendix. We present the results in the two figures below which depict two-dimensional projections of the benchmarks for a number of quantities.
In Fig. 1 we display the distributions of the enhancement factor of theB s → M M rate with respect to its SM prediction for a few pairs of final states M M . As the top-left plot reveals, R(ηω) and R(φω) go up or down simultaneously and can reach roughly 50 and 150 (270 and 170), respectively, for Solution 1 (2). It follows that, in light of the SCET central values in Table II, the Z influence can boost the branching fractions ofB s → ηω and B s → φω to ∼ 2 × 10 −6 and ∼ 6 × 10 −6 , respectively, for both solutions. Accordingly, these decay channels are potentially sensitive to NP signals, and moreover the correlation between R(ηω) and R(φω) may be experimentally checked.
As regardsB s → η ω, for which we do not provide any graphs, with Solution 1 (2) we get at most R(η ω) ∼ 80 (only 2.5), which translates into B(B s → η ω) 0.08 (0.5) × 10 −6 based 4 The CP asymmetry difference ∆A exp [6] might be another restraint. It excludes the SM central values ∆A SM CP −0.01 for the two SCET solutions, but the theoretical error is large, ∼ 0.15 [18], implying that the predictions are still consistent with ∆A exp CP . Similarly, although our Z benchmarks yield −0.025 ∆A SM+Z CP 0.002, they are not in conflict with ∆A exp CP , considering the substantial theoretical uncertainty.  Table II. However, the corresponding upper error in this table for Solution 1 suggests that B(B s → η ω) might still undergo a Z -mediated boost to the 10 −6 level, thereby offering an additional window to the Z interactions.
The top-right plot in Fig. 1 indicates that R(ηπ 0 ) and R(φπ 0 ), like R(ηω) and R(φω), increase/decrease at the same time, although the former two can rise to only about 8.0 and 4.5 (10 and 7.3), respectively, for Solution 1 (2). Nevertheless, as elaborated in Ref. [9], such enhancement factors are sufficiently sizable to makeB s → (η, φ)π 0 promising as extra tools in the quest for the potential NP behind the b → sµ + µ − anomalies. The correlation between R(ηπ 0 ) and R(φπ 0 ) is obviously a testable prediction as well. We notice that the preceding Solution-2 numbers are roughly 20% less than their counterparts (12 and 9.1) in Ref. [9], mostly due to the aforementioned B → πK requisite.
From the bottom plots in Fig. 1, unlike the top ones, it is not evident if there is a connection between R(ηω) or R(φω) and R(ηπ 0 ). The former two also do not seem to have clear correlations with R(φπ 0 ), although this is not illustrated here. We will ignore possibly related consequences forB 0 s → η π 0 , (η, η )ρ 0 because the Z impact on their rates is only modest [9]. Information about relationships between R(M M ) and the Z couplings is highly valuable for examining the latter if one or more of these decays are observed. For our decay channels of greatest interest, it turns out that there are a few relationships that are more or less plain, which we exhibit in Fig. 2. As might be expected, the curves in the fourth plot resemble the corresponding ones in Ref. [9].
It is worth noting that the restriction we imposed above from the B → πK sector is of significance to some extent, although how stringent the condition should be is unclear in view of the 40% uncertainties of the SM rate predictions [18]. For illustration, making it stricter so that the Z -induced changes to the B → πK rates in the SM not exceed 10%, we arrive at the graphs in Fig. 3. In this instance, R(ηπ 0 ) and R(φπ 0 ), especially the latter, become somewhat less remarkable than before, but R(ηω) and R(φω) are still considerable, and so all these decays remain useful for probing the Z effects. We further find, however, that if the B → πK rates are allowed to deviate from the SM expectations by up to 25% or higher, the impact on the maxima of these Rs will start to diminish and they can have the wider spreads depicted in Fig. 4. Lastly, we comment that the Z couplings in our benchmarks are consistent with collider constraints, as discussed in Ref. [9], including those from LHC searches for new high-mass phenomena in the dilepton final states [34]. This is partly because the products of the muon-Z coupling ∆ µµ V and the quark-Z flavor-diagonal couplings (δ ± , ∆ ± ) can be rendered small enough to evade the bounds by sufficiently increasing the size of ρ L while maintaining ∆ µµ V ρ L in Eq. (4) and (δ ± , ∆ ± )ρ L to stay within their desired respective ranges.

IV. CONCLUSIONS
We have explored the possibility that the anomalies detected in the current b → sµ + µ − data arise from physics beyond the SM and that the same new physics also affects the rare nonleptonic decays of theB s meson, most of which are not yet observed. Since the rates of these modes in the SM are comparatively low, one or more of them may be sensitive to NP signals. Adopting a scenario in which the NP is due to the interactions of a heavy Z boson, we investigate the implications for the rare decaysB 0 s → (η, η , φ)(π 0 , ρ 0 , ω). Taking into account the pertinent restraints, we demonstrate that the Z effects can hugely amplify the rates ofB 0 s → (η, φ)ω above their SM predictions, by as much as two orders of magnitude. The corresponding enhancement factors forB 0 s → (η, φ)π 0 could be substantial as well, up to an order of magnitude, in line with our previous work. Thus, these four rare modes are potentially very consequential should future experiments establish that the b → sµ + µ − anomalies are really manifestations of NP.
We have derived the main formulas for the Z contributions toB s → (η, η , φ)(π 0 , ρ 0 ) in Ref. [9] under the SCET framework. Therein we did not include Solution 1 in the evaluation of the Z effects and neglected renormalization-group evolution for simplicity. In the present paper, we include the latter and give results for both Solutions 1 and 2. Thus, summing the SM and Z terms, with the central values of the input parameters and m Z = 1 TeV, we calculate the amplitudes to be, in units of 10 −9 GeV, Bs→ηρ 0 2.56 + 0.77i + (6.35 − 0.12i)δ + ρ + , where the superscripts (1) and (2) refer to Solutions 1 and 2, respectively, Z terms with numerical factors below 0.005 are not displayed, and δ ± and ρ ± were already defined in Eq. (13). We note that in SCET at leading order the rates ofB 0 s → (η, η , φ)(π 0 , ρ 0 , ω) are equal to their antiparticle counterparts [16][17][18]. We have checked that under the same requirements as in Ref. [9] the resulting maximal enhancement factor of theB 0 s → ηπ 0 φπ 0 rate is around 12 (8.6), which is almost identical to (5% below) what we determined earlier [9] ignoring RGE. Imposing also the B → πK condition as discussed in Sec. III may lead to smaller enhancement factors, depending on how stringent it is.