$B$ meson anomalies and large $B^{+}\to K^{+}\nu\bar{\nu}$ in non-universal $U(1)^\prime$ models

In view of both the latest LHCb measurement of $R_{K^{(*)}}$ and the new $2.8\sigma$ deviation reported by Belle II on $B^{+}\to K^{+}\nu\bar{\nu}$ decays, we present a fit to the $B$ meson anomalies for various one and two dimensional hypothesis including complex Wilson coefficients. We show in a model-independent way that the generic non-universal $U(1)^{\prime}$ extensions of the SM, without flavour violation, fail to simultaneously fit those observables and corroborate that they can modify $\mathrm{BR}(B^{+}\to K^{+}\nu\bar{\nu})$ up to only a $10\%$. In view of this deficit, we propose a new way in which those models can accommodate the data at tree level by introducing lepton flavour violating couplings and non-diagonal elements of the charged lepton mixing matrix, with implications in future charged lepton flavour violation searches.


Introduction
A set of different discrepancies (or so called anomalies) between the theoretical predictions of the Standard Model (SM) and their corresponding experimental measurements for different flavour observables, has attracted the attention of the particle physics community for a long time now.The largest deviations are from semi-leptonic quark level decays involving b → s transitions and the anomalous magnetic moment of the muon a µ .Here we will focus on b to s transitions 1 , where the long standing anomalies in b → sµ + µ − have been recently augmented by the first evidence of B + → K + ν ν decays from the Belle II collaboration.Preliminary results were presented at the EPS 2023 conference2 , give a measurement of (2.4 ± 0.5 +0.5  −0.4 ) × 10 −5 for the branching fraction of B + → K + ν ν decays, which deviates by 2.8σ from the SM prediction [2].Combining this with previous data they find a new world average of (1.4 ± 0.4) × 10 −5 .
Within the context of new physics (NP) effects, the most recent global fits to b → sµ + µ − data prefer a vector-like lepton universal (LU) coupling with a pull of around 6σ with respect to the SM [3], though this was before the new b → sν ν measurement.For a recent discussion on the hadronic effects interpretation see [4].Other flavour anomalies are the B meson decays related to b → cℓν ℓ transitions, the Cabibbo angle anomaly (CAA), leptonic τ decays of the form τ → µνν and non-resonant di-electrons [5] which we will not address here.
In general, it is difficult to find suitable NP models able to simultaneously accommodate all of the flavour observables, so using the most precise measurements as discriminators could help at excluding many of those NP scenarios.The simplest NP scenarios for explaining all the flavour anomalies are leptoquark models [6,7], two Higgs doublet models (2HDMs) [8][9][10] and Z ′ models [11,12].In this paper, we will focus on the latter class of models, namely gauged non-universal U (1) ′ extensions of the SM, which are well motivated and have been extensively studied in the literature (for a review see e.g., [13]).One common feature to all non-universal U (1) ′ models is the emergence of new exotic fermions as a result of requiring a theory free of chiral anomalies [14,15], adding then extra degrees of freedom that can alleviate most of the tensions with the experiments 3 .In particular, those models are able to solve the B meson anomalies [18], dark matter [19,20], the fermion mass hierarchy problem [21] and the muon g − 2 in the gauged L µ − L τ model [22].However previous analyses of Z ′ models found only small contributions to B + → K + ν ν decays [23][24][25][26].Here we corroborate that in general at most such models can only explain 10% of the new world average for B + → K + ν ν.In view of this deficit, we propose a new way in which a lepton flavour violating (LFV) non-universal U (1) ′ model could simultaneously fit the B meson neutral anomalies and the latest Belle II measurement.This paper is structured as follows: in section 2 we provide an update for the modelindependent fits to the neutral anomalies with both real and complex Wilson coefficients.Then in section 3 we review the basic ingredients for a generic non-universal U (1) ′ model and in section 4 we present our results and predictions.Finally, we summarise our conclusions in section 5.

Model independent fits
We make use of the flavio package [27] with the same observables as in [28] 4 including the latest LHCb measurements for R K ( * ) [29,30], now in agreement with the SM prediction well within 1σ.In addition, we also include the new Belle II measurement and the related upper limits reported by Belle in [31], for a total of 180 observables.In particular, 96 of those observables are individual bins reported by LHCb [32], CMS [33] and ATLAS [34] related to the coefficients in the angular distributions of B 0 → K * 0 µ + µ − and B + → K * + µ + µ − decays, described by the di-lepton invariant mass squared q 2 and three angles θ K , θ ℓ and ϕ [35], + S 4 sin 2θ K sin 2θ ℓ cos ϕ + S 5 sin 2θ K sin θ ℓ cos ϕ where A FB is the forward-backward asymmetry of the di-muon system and F L is the fraction of longitudinal polarisation of the K * 0 meson.S i are form factor dependent observables which can be related to the so called optimised P ′ basis as [36], We choose to work in the S basis for all angular observables when possible.This is motivated by remarks in [37] pointing out that the optimised observables are clean from hadronic uncertainties only at leading order, besides being non-linearly correlated with each other and therefore, not offering major advantages for precision in phenomenological analysis.
To simplify the analysis, none of the low bins q 2 < 1.0 GeV 2 sensitive to both electromagnetic and chromomagnetic operators in the b → sℓ + ℓ − angular observables and branching ratios were included in the fit.Besides that, and from the latest CMS-ATLAS-LHCb combination of BR(B s → µ + µ − ) presented in [28], scalar operators are strongly constrained and are taken here to be zero, allowing NP effects to come only from vector and axial-like operators.Furthermore, couplings to right-handed quarks are assumed to be negligibly small compared to the left-handed ones and therefore, the prime operators (P L ↔ P R ) will not be taken into account.
In this way, the effective Hamiltonian responsible for b → sℓ + ℓ − transitions can be written as where ℓ = e, µ, G F is the Fermi constant, V tb(s) are elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix and the C ℓℓ i are the effective Wilson coefficients (WCs) at the scale of the bottom quark mass, and are the flavour-changing neutral currents (FCNC) local operators encoding the low-energy description of the high energy physics that has been integrated out.The WCs can be written as where C SM i is the SM contribution to the ith WC and Table 1: Left: Common one dimensional hypothesis for real NP Wilson coefficients.Right: Two dimensional hypothesis for real Wilson coefficients with correlation ρ.The notation for the last row follows from [38] where the short-distance Wilson coefficients can contain two types of NP contribution, where the first contribution is lepton universal and the second one is lepton flavour violating.
C ℓℓ,NP i is the NP contribution.For simplicity of notation, from now on we will drop the NP index unless specified otherwise.
Here, we first review and corroborate the most common 1D and 2D hypothesis presented in [3,28,39] as well as some new 2D scenarios for the WCs from b → sℓ + ℓ − transitions.In Table 1 we provide the best fit values with 1σ uncertainties and correlations for both the 1D and 2D cases, with the Pull SM metric calculated as in [40,41], where F is the χ 2 cumulative distribution function and n dof is the number of degrees of freedom.
We see in Table 1 that one of the most common hypothesis, namely the {C µµ 9 , C µµ 10 } 2D hypothesis, is now disfavoured by 1.6σ with respect to the scenarios fitted by C µµ 9 and C ee 9 , in accordance to the updated LHCb measurements for R K ( * ) [29,30], now SM-like 5 .Particularly interesting is the new 2D scenario {C µµ 9 , C ee 9 = −C ee 10 } (Fig. 1) with a pull of 6.1σ with respect to the SM.
In Fig. 2 we also present an update for the common 2D case {Re(C µµ 9 ), Im(C µµ 9 )} of complex WCs and compare it to the new scenario {Re(C µµ   in [39].We corroborate the findings of [39,43,44] and observe no improvement over the real WCs case, although we can see a slight preference from the Belle+BaBar data for a large negative imaginary part of C µµ 9 .Regarding b → sν ν, if no right-handed neutrinos are included, or equivalently, if their mixing with active neutrinos is suppressed, the effective Hamiltonian at the bottom quark scale is defined as where with x t = m 2 t /m 2 W and the function X(x t ) is defined at the next-to-leading order in QCD in [45,46], giving C SM LL = −6.38 ± 0.06 [47].In this way, we will have at first 9 WCs coming from the purely left-handed 10 operators contributing to the branching ratio [25], with the SM prediction given by BR(B + → K + ν ν) SM ≃ 4.4 × 10 −6 from the flavio package.
In the most general case for real WCs, we will have four independent flavour conserving coefficients associated to the O ℓℓ 9 and O ℓℓ 10 operators, summing up to a total of 13 WCs for the fit.In Table 2 we can see again the well-known preference for the muonic WC C µµ 9 ≃ −1.0 accompanied in this case by a small positive C µµ 10 and our data shows a new interesting correlation C ee 9 ≃ −C ee 10 = −0.5 for the electronic WCs.As for the Belle II measurement, it gets fitted by a LU violating coupling with di-electron neutrinos equal to C νeνe LL ≃ −7 with order 7% correlations with the muonic WCs.Taking these last values and correlations into account, we proceed now to map the WCs into the parameter space for a generic non-universal Z ′ model.

Non-universal U (1) ′ models
These kinds of models assume that the dominant NP contribution to b → sℓ + ℓ − and b → sνν transitions comes from tree level exchange of a SU (2) L singlet Z ′ gauge boson associated to an extra U (1) ′ symmetry [23].This symmetry is broken spontaneously at the TeV scale by the vacuum expectation value (VEV) of a new neutral scalar χ, making the Z ′ boson to acquire a mass term M Z ′ = g X υ χ with g X the new gauge coupling.The Z ′ will also mix with the SM Z boson in the kinetic terms of the Lagrangian, however we will consider such mixing negligible 6 .The purely new physics interaction Lagrangian with fermions in the flavour basis is defined then by, where, under the same assumptions as before, the relevant couplings will be g bs L , g ℓℓ L,R and g νν L , and the WCs generated at tree level by the effective Hamiltonian in Eq.(2.4) are given by, where the lepton couplings are defined as, and can be in general non-universal, determined by both the U (1) ′ charges and by the n×n rotation matrix V l L which transforms flavour eigenstates to mass eigenstates with n = 3 + E, where E is the number of exotic charged leptons.For concreteness, we define the SM part of V l L as where α L , δ and ϵ are free parameters of the rotation matrix constrained by both unitarity and the masses of the charged leptons.This texture is justified by recent searches of lepton flavour violation in Higgs boson decays from both ATLAS [49] and CMS [50,51], reporting upper limits on the branching ratios of H → τ µ, H → τ e and H → µe decays, consistent with decoupling of 1-2 and 1-3 rotations.On the other hand, we will take the matrix V l R to be diagonal, i.e., with lepton universal right-handed couplings, Regarding the b → sνν process, the WC in question according to Table 2 is given by, V PMNS with V PMNS the Pontekorvo-Maki-Nakagawa-Sakata (PMNS) matrix, meaning that depending on the ordering (inverted or normal) for the neutrino mixing parameters, the lepton couplings entering in Eq.(3.7) may differ.

Results
Equipped with the required WCs both from the model-independent fit and for the U (1) ′ models, we define the quadratic approximation to the likelihood function where C = {C µµ 9 , C µµ 10 , C ee 9 , C ee 10 , C νeνe LL } is defined in terms of the WCs in Eqs.(3.2-3.7),C bf is given by the central values in the second column of Table 2, and Cov is the covariance matrix or Hessian associated to the correlation matrix ρ also in Table 2 obtained by using the MIGRAD minimization algorithm.With this function at hand, a random generator in Mathematica is requested to find points inside the ellipsoids defined by ∆χ 2 ≤ σ, 2σ for 2 degrees of freedom and boundaries defined for two scenarios: 1. Scenario with lepton flavour universality violation (LFUV).

Scenario with lepton flavour violation (LFV) in addition to LFUV.
For the first scenario we scan over the following couplings where the ranges are chosen based on the results from [48], where a global analysis was performed constraining the lepton couplings from various leptonic processes mediated at tree level by a generic Z ′ and the muonic (electronic) couplings were found to be at most of order 2 (0.2). Regarding the values for the Z ′ mass, we follow [52] and [53] and set the lower and upper mass values for the scan as 1 and 10 TeV respectively.We also evade more stringent constraints from colliders given that, as in [23], we do not specify the flavourconserving couplings to first and second generation quarks and therefore the high-mass Drell-Yan bounds will be suppressed assuming that the Z ′ boson couples dominantly to bottom quarks [28] and to leptons [48].
For the LFV scenario, we additionally scan over the following parameters, for i = e, µ, τ , X R is the single flavour universal charge for the right-handed leptons and the flavour non-universal charges for the left-handed lepton couplings are related to the lepton couplings by Eq.(3.4).As a consequence of lepton flavour violation in this scenario, the couplings of Eq.(3.4) will be non-diagonal allowing for processes as µ → eγ, which we take explicitly into account as a constraint from the latest measurement made by the MEG collaboration [54], with the theory prediction given by [48], where Γ µ is the decay width of the muon and with c 12 R obtained from c 12 L by interchanging L and R. We can see in Fig. 3 that the LFUV scenario (red) can generate at most C νeνe LL ≃ −0.7 at the 1σ level while simultaneously fitting the rest of the B meson anomalies, meaning a relative increment of ∼ 10% for BR(B + → K + ν ν) with respect to the SM prediction 7 .This deficit with respect to the Belle-II measurement is alleviated in our LFV scenario (blue), obtaining the required value for the Wilson coefficient C νeνe LL ≃ −7, interfering constructively with the SM contribution.Assuming bi-large mixing in the PMNS matrix (s 13 → 0) and normal ordering for the neutrino parameters [56], the best fit point for the LFV scenario is given by  L , X τ L } plane using the same colours as the left panel and in the {X R , X e L } using green and light green.The purple regions show how these regions shrink when we additionally include constrain our samples with µ → eγ by post-processing the data.We also include neutrino trident constraints but these do not have any impact on the contours in these planes.
In Fig. 4 we consider additional constraints on the model.In both panels we additionally constrain our samples using the ∆M s mass difference from B s − B s mixing [42,53,57], and show the remaining 1σ and 2σ regions in blue and light blue respectively, while the constraint itself is indicated by the red shading in the left panel, where we present the results in the plane of the gauge boson mass and the g bs L coupling.In the right panel of Fig. 4 we show these results (including the ∆M s constraint) in planes featuring the lepton charges.We also consider constraints from from µ → eγ using Eq.(4.4) and from the neutrino trident production cross section given by [58], which we required to be within σ exp /σ SM = 0.83 ± 0.18 [57].However we find that the neutrino trident constraint has no impact on our samples, which is in agreement with [58], where the neutrino trident production was found to be less restrictive for TeV masses.On the other hand, the µ → eγ constraint can eliminate many scenarios in the full parameter space, with a clear impact (purple regions on the right panel of Fig. 4) on the projected 1σ contours in the {X µ L , X τ L } and {X R , X e L } planes.Note that this is from requiring µ → eγ is fulfilled at the 1σ level, whereas the effect would be very mild if we relaxed this to 2σ since our parameter ranges were chosen to avoid lepton flavour constraints at the 2σ level.
In both plots in Fig. 4 we observe an allowed region for the Z ′ mass in the 2-8 TeV range at the 1σ level, although requiring large values for the X µ,τ L couplings at the best fit point, near the perturbativity limit of √ 4π and within reach of the future Muon Collider (MuC) [59,60].In contrast, rather small values for both the left-handed electron coupling (compatible with zero at the 1σ level) and for the right-handed universal coupling X R are preferred by the fit.The values of these lepton couplings will in general depend in the NP scenario chosen from the model-independent fits.However our fit results could provide a useful tool for constraining the quantum numbers for the SM charged leptons, and subsequently, for finding the values for the charges of the exotic fermions needed in the cancellation of the chiral anomalies.Such an algorithm for constructing anomaly-free tree-level Z ′ models directly from experimental data is left for a follow-up work.One of the main consequences from the LFV scenario is the prediction of the charged lepton flavour violating (CLFV) decays B s → ℓ ± ℓ ′∓ and B + → K + ℓ ± ℓ ′∓ .In particular, using the theoretical expressions from [62], in Fig. 5 we show the generated values from the allowed parameter space for two of those CLFV modes, obtaining a good discovery potential associated to the Z ′ lepton flavour violating scenario in the B s → e ± µ ∓ channel, within reach of future searches at both Belle II and the LHCb Upgrade II.

Summary and conclusions
We presented a model-independent update to some of the most common one and two dimensional scenarios fitting the neutral B meson anomalies and corroborated that the {C µµ 9 , C µµ 10 } NP scenario is now disfavoured with respect to the scenarios fitted by C µµ 9 and C ee 9 .After this, and in view of the new Belle II measurement, we performed a 5 dimensional fit to both the b → sℓℓ and b → sν ν observables and found a pull of Pull SM = 5.9σ.Based on this last result and using a quadratic approximation to the χ 2 function, we then mapped out the resultant Wilson coefficients into the parameter space of two U (1) ′ scenarios, one with LFUV and the second with LFV.We found that the LFUV scenario can generate an increment to the B + → K + ν ν branching ratio of 10%, insufficient to explain the Belle II data.In the LFV scenario we showed that such branching ratio, including constraints from B s − B s mixing and neutrino trident production, can be large enough once O(0.4) off-diagonal elements of the charged lepton rotation matrix were included.Finally, we presented predictions for CLFV decays, the most promising being the B s → e ± µ ∓ decay, within reach of future searches at both Belle II and the LHCb Upgrade II.

9 )Figure 1 :
Figure 1: Global fits showing 1σ (lighter colours) and 2σ (darker colours) regions of the b → sℓ + ℓ − data available from LHCb (bright blue), CMS (green), ATLAS (orange), Belle and BaBar (both combined and displayed in blue).Given that the vast majority of observables come from LHCb, this experiment drives the total fit (red), almost indistinguishable from the LHCb data alone.Left: {C µµ 9 , C µµ 10 } scenario now disfavoured with respect to the scenarios fitted by C µµ 9 and C ee 9 .Right: New physics 2D hypothesis involving electronic Wilson coefficients.

Figure 4 :
Figure 4: Constraints to the parameter space of the LFV scenario fitting the b → sℓ + ℓ − data simultaneously with BR(B + → K + ν ν).Left: Projected 1 and 2σ regions (blue and light blue respectively) in the {M Z ′ , g bs L } plane constrained at 1σ level by ∆M s with the region allowed by ∆M s shown in red.Right: The projected 1 and 2σ regions constrained at 1σ level by ∆M s is shown in the {X µL , X τ L } plane using the same colours as the left panel and in the {X R , X e L } using green and light green.The purple regions show how these regions shrink when we additionally include constrain our samples with µ → eγ by post-processing the data.We also include neutrino trident constraints but these do not have any impact on the contours in these planes.