Predictions for $B_s \to \bar{K}^* \ell \,\ell$ in non-universal $Z'$ models

The lepton flavor universality violating (LFUV) measurements $R_K$ and $R_{K^*}$ in $B$ meson decays can be accounted for in non-universal $Z'$ models. We constrain the couplings of these $Z'$ models by performing a global fit to correlated $b \to s \ell \ell$ and $b \to d \ell \ell $ processes, and calculate their possible implications for $B_s \to \bar{K}^*\ell \ell$ observables. For real new physics (NP) couplings, the 1-$\sigma$ favored parameters allow the corresponding LFUV ratio $R_{K^*}^{(s)}$ in $B_s \to \bar{K}^*\ell \ell$ to range between 0.8 -- 1.2 at low $q^2$. Complex NP couplings improve the best fit only marginally, however they allow a significant enhancement of the branching ratio, while increasing the range of $R_{K^*}^{(s)}$ at low $q^2$ to 0.8 -- 1.8. We find that NP could cause zero-crossing in the forward-backward asymmetry $A_{FB}$ to shift towards lower $q^2$ values, and enhancement in the magnitude of integrated $A_{FB}$. The $CP$ asymmetry $A_{CP}$ may be suppressed and even change sign. The simultaneous measurements of integrated $R_{K^*}^{(s)}$ and $A_{CP}$ values to 0.1 and 1% respectively, would help in constraining the effective NP Wilson coefficient $C_9$ in $ b \to d \mu \mu$ interactions.

B s → φ µ + µ − [5] and angular observable P 5 in B → K * µ + µ − decay [6][7][8]. Hence it is natural to try accounting for the discrepancies in all the above measurements by assuming new physics only in the muon sector.
These anomalies may be addressed in a model-independent way using the framework of effective field theory, where the effects of beyond-SM physics are incorporated by adding new operators to the SM effective Hamiltonian. Various groups have performed global fits to all available data in the b → s + − sector in order to identify the Lorentz structure of possible new physics operators [9,[27][28][29][30][31][32][33]. Some of these new physics operators can be generated in Z [13][14][15][16] or leptoquark models [21][22][23][24][25][26]. It has been shown that several models with Z , either light or heavy, can help account for the anomalies in b → s µ + µ − sector [17][18][19][20].
Since the Z boson would in general couple to all generations, the imprints of such a Z would be seen in other flavor sectors as well. Therefore, it is worth extending this model to include other related decays. This will provide further insights into the NP flavor structure. In this work, we consider possible observable effects of Z models on decays induced by the quark-level transition b → d µ + µ − .
A large number of b → dµ + µ − decays, at the level of thousands or tens of thousands, would be observed after the LHC upgrade. For example, about 17000 B + → π + µ + µ − events are expected to be observed after collection of the full 300 fb −1 dataset. For B s →K * µ + µ − decays, the full angular analysis is expected to be possible after the LHCb Upgrade-II dataset, where around 4300 events could be observed [39]. This would enable the measurements of angular observables in B s →K * µ + µ − decays with a precision even better than the existing measurements of angular distributions in B d → K * µ + µ − decay.
Currently, as there are not many measurements in the b → d sector, a model-independent analysis would not be very useful in constraining new physics. However, in the context of specific models (like Z ), some of the couplings can be constrained from the b → s sector and neutrino trident production. Therefore we choose this approach to constrain the effective couplings in the b → d µ + µ − sector, and identify potential observables in the B s →K * µ + µ − decay where large new physics effects are possible.
The paper is organized as follows. In section II, we introduce the Z model considered in this work and indicate how it can be constrained by available measurements. We then describe our fit methodology in Section III. The fit results along with predictions of various B s →K * µ + µ − observables are presented in Section IV. We summarize in Section V.

II. THE Z MODEL AND SOURCES OF CONSTRAINTS
In the non-universal Z model that we consider, the Z boson is associated with a new U (1) symmetry. It couples to both left-handed and right-handed muons but not to leptons of other generations. It couples to both left-handed and right-handed quarks, however we assume its couplings to right-handed quarks to be flavor-diagonal, thereby avoiding contribution of new chirality flipped operators to flavor changing neutral current (FCNC) decays [66,67]. The change in the Lagrangian density due to the addition of such a heavy Z boson is where The right-hand side in eq. (2) includes only the terms contributing to FCNC processes. Here P L(R) = (1 ∓ γ 5 )/2, Q i is the i th generation of quark doublet, and L = (ν µ , µ) T is the second generation doublet. Further, g µµ L(R) are the left-handed (right-handed) couplings of the Z boson to muons, and g bq L to quarks. One can integrate out the heavy Z and get the relevant terms in the effective four-fermion Hamiltonian as, In eq. (3), the first (third) term corresponds to b → s(d)µ + µ − transitions, the second (fourth) terms give rise to B s -B s (B d -B d ) mixing, whereas the fifth term contributes to the neutrino trident production ν µ N → ν µ N µ + µ − (N = nucleus). Consequently, the products g bs L g µµ L,R (g bd L g µµ L,R ) are constrained by the b → s(d)µ + µ − data, and individual magnitudes |g bs | (|g bd |) from the B s -B s (B d -B d ) mixing. The neutrino trident production puts limits on the individual muon couplings g µµ L,R . We now discuss constraints on the Z couplings arising from each of the above measurements.
The effective Hamiltonian for b → qµ + µ − transition in the SM is  6,8) contribute to these processes through the modifications C 7,9 (µ) → C eff 7,9 (µ, q 2 ), where q 2 is the invariant mass-squared of the final state muon pair. We drop the superscript "eff" from here on for the sake of brevity. Addition of the new Z boson to the SM particle spectrum modifies the WCs as C bq 9,10 → C bq,SM 9,10 + C bq,NP 9,10 , where In the Z models, (C bs,NP 9 and C bs,NP 10 are in general independent. Two of the one-parameter scenarios, C bs,NP 10 = 0 (popularly known as C bs,NP 9 < 0) and C bs,NP 9 = −C bs,NP

10
, can be realized by substituting g µµ L = g µµ R and g µµ R = 0, respectively.
The dominant contribution to B q −B q mixing within the SM comes from the virtual top quark in the box diagram. The Z boson contributes to B q −B q mixing at the tree-level. The combined contribution to M q 12 , the dispersive part of the box diagrams responsible for the mixing, is where Here η B =0.84 is the short-distance QCD correction calcualated at NNLO [40], f Bq is the decay constant, and B Bq is the bag factor. The mass difference ∆M q = 2|M q 12 | is while the relevant weak phase φ q is

Neutrino trident production
Within the Z models, the modification of the cross section σ for neutrino trident production, where v = 246 GeV and s W = sin θ W .

III. FIT METHODOLOGY
We now determine favored values of the new physics couplings g bs L , g bd L , g µµ L and g µµ R . We nominaly take the mass of the Z boson to be M Z = 1 TeV. Note that since M Z only appears through the combination g 2 /M 2 Z , the constraints on couplings can be appropriately scaled with the actual value of M Z .
All the observables in the b → sµ + µ − sector put constraints on the combinations g bs L g µµ L and g bs L g µµ R . For the fit related to b → sµ + µ − , we closely follow the methodology of Ref. [27]. The χ 2 function for all the b → sµ + µ − observables listed above is calculated as where C i = C bs,NP 9,10 . Here O th (C i ) are the theoretical predictions of b → sµ + µ − observables calculated using flavio [54], and O exp are the corresponding experimental measurements. The total covariance matrix C is obtained by adding the individual theoretical and experimental covariance matrices.
We now turn to B q −B q mixing. Here we consider constraints from ∆M d , ∆M s , and the two CP-violating phases. Using f B d B B d = (225 ± 9) MeV [58], along with other input parameters from ref. [59], eq.  [57], the contribution of ∆M d to χ 2 is where we denote the experimental mean value of an observable X by X exp,m , and the uncertainty in the observable by σ X . In order to obtain σ X , we add the experimental and theoretical uncertainties in quadrature. Here, σ ∆M d is dominated by the theoretical uncertainty.
In order to minimize the impact of theoretical uncertainties, we use ∆M s constraints through the ratio M R = ∆M d /∆M s . In the SM, where [56] and |V td /V ts | = 0.2088 +0.0016 −0.0030 [55], we obtain M SM R = 0.0297 ± 0.0009, where we have added the errors in quadrature. Wherever there are asymmetric errors, we take a conservative approach and use the larger of the errors on two sides.
The value of M exp R = 0.0285 ± 0.0001 [57], so the contribution to χ 2 due to this ratio is The observables ∆M d and M R constrain |g bs L | and |g bd L |.
The CP -violating constraints from J/ψφ and J/ψK S decays contribute to the χ 2 as where Here we have taken the measurements to be S exp J/ψφ = 0.02 ± 0.03 and S exp J/ψK S = 0.69 ± 0.02 [59]. For the constraints from neutrino trident production, we use the quantity R ν = σ/σ SM , whose theoretical expression is given in eq. (10). We have taken R exp ν ≡ 0.82 ± 0.28 [60,75]. The contribution to the total χ 2 is This observable constraints g µµ L and g µµ R . The b → dµ + µ − decays are CKM-suppressed as compared to b → sµ + µ − . In our analysis, we include constraints from the branching ratios of B + → π + µ + µ − and B d → µ + µ − decays. We do not include the measurements of observables in B s →K * µ + µ − decay in our fit, since we are interested in obtaining predictions for these.
The branching ratio of B d → µ + µ − in our model is given by and the contribution to χ 2 is We [56], and other inputs from [59].
Finally, combining all the above constraints, we obtain In the next section, we present our fit results, along with predictions of several observables in The blue curve is the boundary of the 1σ-favored region due to constraints from measurments in b → d sector and neutrino trident production. The pink shaded region represents the 1σ-favored parameter space after including additional constraints from b → sµ + µ − data and B s −B s mixing.

IV. FIT RESULTS AND PREDICTIONS
In a model-independent analysis, it is hard to obtain useful constraints on the new physics couplings for b → d decays, since there are only a few measurements in this sector. However, within the context of a Z model, one can obtain useful constraints using correlated b → s and b → d processes. This can be seen from Fig. 1, which depicts the allowed (g bd L , g µµ L = g µµ R ) parameter space corresponding to a Z model which generates the 1D scenario C bd,NP 10 = 0. The elliptical region represents the 1σ-favored parameter space with constraints only from b → d sector, i.e., branching ratios of B + → π + µ + µ − and B d → µ + µ − decays, B d −B d mixing, and neutrino trident production. The two shaded regions represent the 1σ-favored parameter space obtained by including additional constraints from all relevant measurements related to b → sµ + µ − decays and B s −B s mixing. It can be seen that the allowed range of NP couplings, in particular g µµ L,R , reduces considerably after including constraints from the b → s sector. Therefore, it is worth studying implications of several measurements in the b → s sector on the observables in b → dµ + µ − decays.
Performing a fit to the relevant observables in b → s and b → d sectors, we determine the 1σ-favored parameter space of the couplings g bd L , g bs L , g µµ L and g µµ R , considering g bs L and g bd L to be (i) real, (ii) complex. These can be used to find constraints on the NP Wilson coefficients (C bd,NP 9 , C bd,NP 10 ), and to put limits on the allowed NP in the following observables in B s →K * µ + µ − decay: differential branching ratio, the LFUV ratio R (s) K * , muon forward-backward asymmetry A F B , longitudinal polarization fraction F L , and direct CP asymmetry A CP .
The matrix element for the decay amplitude of B s →K * µ + µ − can be written as where C bd,SM 9 and C bd,SM 10 are taken from Ref. [70] and Ref. [69] respectively. The matrix elements appearing in eq. (21) have been calculated using form factors obtained by a combined fit to lattice calculations and QCD sum rules on the light cone (LCSR) [68]. We also include the non-factorizable corrections due to soft gluon emission and charmonium resonance, which have been computed for 71,72], and parameterized as corrections to C SM 9 . These effects are assumed to be roughly the same for B s →K * due to flavor symmetry [65].
The decay B s →K * µ + µ − may be described in terms of the four-fold distribution as [65] d where q 2 is the lepton invariant mass, θ V and θ l are the polar angles, and φ is the angle between the dimuon plane and K * decay plane. The relevant observables can be obtained from the four-fold distribution as where the functions I i can be expressed in terms of the transversity amplitudes [64]. HereB corresponds to the decay modeB s → K * µ + µ − .
We present our results for the above observables at four benchmark NP scenarios as given in   A. Real couplings In order to quantify how well the Z model is able to account for all data in the b → s and b → d sectors, we define ∆χ 2 = χ 2 SM − χ 2 NP , where the minimum χ 2 in the SM, and in the presence of NP Z couplings, is denoted by χ 2 SM and χ 2 NP , respectively. For the case of real couplings, we find the best fit values to be g bd L = ± 0.3 × 10 −3 , g µµ L = ∓ 0.4, and g µµ R = ∓ 0.2. The value of χ 2 SM ≈ 221 and ∆χ 2 ≈ 41, so that the SM point corresponding to g µµ L = g µµ R = g bd L = 0 is highly disfavoured. The 1σ-favored parameter space of the couplings (g bd L , g µµ L g µµ R ) is shown in Fig. 2. It can be seen from (g bd L , g µµ L ) and (g bd L , g µµ R ) planes that, while g µµ R = 0 is barely disfavored within 1σ, a rather wide strip |g µµ L | ≤ 0.25 lies beyond the 1σ-favored region. This is because the anomalies in  Table I . b → sµµ decays need a non-zero value of C bs,NP 9 , which in turn require a non-zero value of g µµ L or g µµ R . Furthermore, the scenario C bs,NP 9 = −C bs,NP 10 [27], which provides a good fit, favors g µµ R = 0, thus requiring g µµ L to be away from zero. Note that the results in the (g µµ L , g µµ R ) plane indicate the class of favored solutions that lie along g µµ L = g µµ R , corresponding to C bs,NP

Predictions for
The top left panel of Fig. 3 shows predictions for the differential branching ratio corresponding to real Z couplings, for the SM as well as two benchmark scenarios NP1 and NP2 from Table. I. These scenarios roughly correspond to the maximum deviation on either side from the SM predictions in the 1σ favored NP parameter space. The maximum enhancement (suppression) in the differential branching ratio corresponds roughly to a maximum positive (negative) value of C bd,NP

9
. It can be seen from the figure that only a marginal enhancement or suppression over the SM value is possible in the differential branching ratio. A clean distinction among the predictions of different scenarios is difficult owing to the large uncertainties (about 20%) arising from the form-factors.
A measurement of the LFUV ratio R (s) K * (q 2 ) in a few q 2 bins would be possible with the LHCb upgrade-II data set [39]. The predictions for this quantity in the benchmark scenarios NP1 and NP2 are shown in the top right panel of Fig. 3. In the SM, R (s) K * (q 2 ) is unity in the entire low-q 2 region, while an enhancement up to 1.3 and a suppression up to 0.8 is allowed. The maximum enhancement (suppression) roughly corresponds to the maximum positive (negative) value of C bd,NP
Within the SM, the forward-backward asymmetry A F B (q 2 ) is predicted to vanish around q 2 ≈ 3.5 GeV 2 , and the zero-crossing is from negative to positive, as can be seen from the bottom left panel in Fig. 3. The maximum value of A F B (q 2 ) in the SM is ≈ 10 %. The positive (negative) value of C bd,NP 9 also shifts the zero-crossing towards lower (higher) q 2 value. The integrated value of A F B over q 2 = (1 − 6) GeV 2 bin is (−0.6 ± 1) % within the SM . The predictions for integrated A F B for the benchmark scenarios NP1 and NP2 are (3.1 ± 1.4)% and (−5 ± 1.7)%, respectively.
The predictions for longitudinal polarization fraction F L (q 2 ) are shown in the bottom right panel of Fig. 3. Within the SM, the peak value of F L (q 2 ) is ≈ 0.9 around q 2 ≈ 1.8 GeV 2 . The shape of F L (q 2 ) does not change with NP and only a marginal deviation from SM is allowed for the benchmark NP scenarios considered here.
Thus, in the case of real couplings, R (s) K * (q 2 ) is useful to distinguish the predictions of the two benchmark NP scenarios from the SM expectation, while the predictions for the differential branching ratio, A F B (q 2 ), and F L (q 2 ) may not have distinct NP signatures, owing to the large form factor uncertainties. The results obtained for integrated dB/dq 2 and R (s) K * (q 2 ) over the q 2 = (1 − 6) GeV 2 bin are presented in Fig. 4. These results are depicted in the (C bd,NP

B. Complex couplings
We would now like to see how the predictions for the above observables in the B s →K * µ + µ − decay would change if the couplings g bd L and g bs L are allowed to be complex. Note that since the leptonic current in eq. (3) is self-conjugate, g µµ L and g µµ R must be real. We also study the impact of these complex couplings on the direct CP asymmetry in this decay. , g µµ L , and g µµ R . The minimum χ 2 in the presence of the complex NP couplings is χ 2 NP ≈ 178, so that ∆χ 2 ≈ 43, thereby providing a slightly better fit as compared to the case of real couplings (∆χ 2 = 41). The corresponding best fit values are Re[g bd L ] = ± 2.7 × 10 −3 , Im[g bd L ] = ∓ 3.8 × 10 −3 , g µµ L = ∓ 1 and g µµ R = ∓ 0.255. As χ 2 NP for complex couplings is lower compared to that for real couplings, the 1σ-favored parameter space shifts further away from the SM point, which has Re[g bd L ] = 0. A larger parameter space is allowed for the muon couplings compared to the real case, |g µµ L | ≤ 2.5 and |g µµ R | ≤ 1.4. The 1σ favored region encompasses g µµ R = 0, whereas a rather large region around g µµ L = 0 (i.e. |g µµ L | ≤ 0.4) is disfavoured within 1σ. The allowed range of Im[g bd L ] is qualitatively similar to that of Re[g bd L ]. Note that the complex nature of g bd L is constrained only from B d −B d mixing measurements, since no CP -violating measurements are currently available in the b → dµµ sector.
The predictions for differential branching ratio and R (s) K * (q 2 ) for the Z model with complex couplings are shown in the top panel of Fig. 6 , for SM as well as the benchmark scenarios NP3 and NP4 in Table I Table I . in NP3, which could be useful in identifying deviations from the SM. A large enhancement is also possible in the LFUV ratio R (s) K * (q 2 ) in the NP3 scenario, with the maximum value of R (s) K * = 1.8 at q 2 = 6 GeV 2 . While the scenario NP4 cannot be distinguished from the SM using only the branching ratio, the value of R (s) K * (q 2 ) in this scenario can be as low as 0.85. Therefore, R (s) K * (q 2 ) would be useful to identify deviations from the SM.
A marginal enhancement in A F B (q 2 ) is possible for the scenario NP3, which would also display zero-crossing at much lower q 2 values (q 2 ≈ 2.5 GeV 2 ) compared to that in the SM (q 2 ≈ 3.5 GeV 2 ).
A marginal suppression in F L (q 2 ) is also possible in NP3. The scenario NP4, on the other hand, does not show significant deviations from the SM for these two observables.

Direct CP asymmetry
The direct CP asymmetry in the b → d µ + µ − sector is expected to be about an order of magnitude larger than b → s µ + µ − . As direct CP violation in b → s µ + µ − sector is expected to be ∼ 0.1%, its experimental observation would be possible only if some new physics provides an order of magnitude enhancement to bring it up to the level of a few percent. In b → d µ + µ − decays, the A CP in SM itself is at the level of a few per cent, and can be within experimental reach. Fig. 7 shows A CP (q 2 ) in the low-q 2 region for the decay B s →K * µ + µ − , considering the benchmark scenarios NP1 and NP2 (real couplings), as well as NP3 and NP4 (complex couplings). It can be seen from the left panel of the figure that for real couplings, A CP (q 2 ) is either marginally below the SM prediction or almost consistent with it.
For complex couplings, however the suppression in A CP (q 2 ) can be quite large. It can even lead to A CP (q 2 ) falling below a per cent level, hence making its measurement extremely difficult.
In some scenarios (e.g. NP3), it is even possible for A CP (q 2 ) to be negative for very low q 2 values. After scanning over the 1σ-favored parameter space, we find no significant enhancement in A CP (q 2 ). So an NP signal can be established if the measurements put an upper bound which is firmly below the SM prediction of A CP (q 2 ).   ]. An enhancement in integrated R ]. We find that the measurements of integrated R

V. SUMMARY AND CONCLUSIONS
In non-universal Z models, instrumental in accounting for the flavor anomalies, the observables in b → s and b → d processes would be correlated. In this paper, we study the constraints on the couplings of a non-universal Z model from the measurements in b → qµµ (q = s, d) decays, B q −B q mixing, and neutrino trident production. These couplings give rise to new additional contributions to the Wilson coefficients C bq 9 and C bq 10 . Using the above constraints, we perform a global fit to determine 1σ-favored regions in the parameter space of the Z couplings g bd L , g bs L , g µµ L , and g µµ R . We analyze the cases when quark-Z couplings g bd L and g bs L are (i) real, and (ii) complex. We also present our predictions for some important observables in B s →K * µµ decays -the differential branching ratio dB/dq 2 , the LFUV ratio R (s) K * , the angular observables A F B and F L , and the CP asymmetry A CP -for some benchmark scenarios.
It is observed from our analyses that the Z model improves the global fit over the SM by ∆χ 2 ≈ 41 (real couplings) and ∆χ 2 ≈ 43 (complex couplings). The favored regions in the parameter space lie along g µµ L ≈ g µµ R , corresponding to C bd,NP 9 ≈ 0, while the region around g µµ L = 0 is disfavored. These are mainly dictated by the R K and R K * anomalies in b → s sector.
For the observables in B s →K * µ + µ − decays, when the couplings are real, we find that the enhancement and suppresion in dB/dq 2 cannot be cleanly identified due to the large uncertainties in the SM prediction. However, the value of R . There is no significant deviation from SM in the predictions of F L (q 2 ), and the predictions of A CP (q 2 ) also stay close to the SM expectation for all the favored values of NP Wilson coefficients. Further, we find that a measurement of integrated R (s) K * in the low-q 2 bin with a precision of ∼ 0.1 can help narrow down the ranges of (C bd,NP 9 , C bd,NP
In the case of complex couplings, a larger NP parameter space is allowed, leading to larger possible deviations in the B s →K * µµ observables. In particular, a ∼ 50% enhancement in dB/dq 2 is allowed. Moreover, the LFUV ratio R (s) K * (q 2 ) can be enhanced up to 1.8 in scenarios with large positive and negative Im[C bd,NP