Consequences of $b\, \to\, s\, \mu^+\, \mu^-$ anomalies on $B \, \to \, K^{(*)} \,\nu \, \bar{\nu}$, $B_s \to\, (\eta,\eta') \, \nu\,\bar{\nu}$ and $B_s \, \to \, \phi \, \nu \, \bar{\nu}$ decay observables

The long persistent discrepancies in $b\,\to\,s\, \ell\, \ell$ quark level transitions continue to be the ideal platform for an indirect search of new physics that lies beyond the SM. The measurements of $R_K$, $R_{K^*}$, $P_5^{\prime}$, $\mathcal{B}(B_s\, \to\, \phi\, \mu^+\,\mu^-)$ and $\mathcal{B}(B_s\, \to\, \mu^+\,\mu^-)$ persistently deviate from the standard model expectations. Similarly, the new tests of lepton flavor universality performed using the isospin partners such as $B^0 \, \to \, K_S^{0} \, \ell\,\ell$ and $B^+ \, \to \, K^{+*} \, \ell\,\ell$ exhibit the same pattern of deviation in $R_{K_S^0}$ and $R_{K^{*+}}$ with the existing results. Motivated by these anomalies we search for the patterns of new physics in the family of flavor changing neutral decays with neutral leptons in the final state undergoing $b \, \to \,s\, \nu \, \bar{\nu}$ quark level transitions. There are close relations between $b\,\to\,s\, \ell\, \ell$ and $b \, \to \,s\, \nu \, \bar{\nu}$ transitions not only in standard model but also in various beyond the standard model scenarios. For beyond the standard model physics under the $\rm SU_L(2)$ gauge symmetry one can relate the left handed charged leptons to the neutrinos. Moreover, there are several advantages of studying $b \, \to \,s\, \nu \, \bar{\nu}$ transitions over $b\,\to\,s\, \ell\, \ell$ as they are free from various hadronic uncertainties beyond the form factors such as the non-factorizable corrections and photonic penguin contributions. Hence, we explore the consequences of $b\, \to\, s\, \mu^+\, \mu^-$ anomalies on $B \, \to \, K^{(*)} \,\nu \, \bar{\nu}$, $B_s \to\, (\eta,\eta') \, \nu\,\bar{\nu}$ and $B_s \, \to \, \phi \, \nu \, \bar{\nu}$ decay observables in SMEFT platform within various 1D and 2D new physics scenarios.


I. INTRODUCTION
In the standard model (SM), the three families of leptons are identical except for their masses. More specifically, the photon, Z and W ± bosons couple to them with equal strengths. Hence the SM is lepton flavor universal. However, there exists several hints of lepton flavor universality (LFU) violation in B meson decays mediating via b → c ℓ ν charged current and b → s ℓ + ℓ − neutral current transitions reported by BABAR, Belle, and more recently by LHCb Collaboration. The recent updated measurements of R K , B(B s → ϕ µ + µ − ), B(B s → µ + µ − ) and new measurements of R K 0 S and R K * + from LHCb continue to exhibit the same pattern of deviations with respect to the SM expectations. For completeness, we report the current status of several b → s ℓ + ℓ − decay observables in Table I.
The rare semileptonic b → s ℓ + ℓ − transition processes are very interesting probes of new physics (NP) because of their sensitivity to various NP contributions that can, in principle, appear in the penguin loop diagrams or in the box diagrams. It is, however, worth mentioning that a precise SM prediction of the observables is crucial to disentangle the genuine NP contribution from the SM uncertainties that may come from meson to meson transition form factors and CKM matrix elements. In recent years, the QCD motivated approaches based on the lattice quantum chromodynamics (LQCD) and light cone sum rule (LCSR) have provided very precise value of the form factors for various b → s ℓ + ℓ − processes. For several preferred decay modes such as B → K ( * ) ℓ + ℓ − and B s → ϕ ℓ + ℓ − , a very precise value of the form factors are obtained within LCSR and LQCD [1,2] approach. Apart from these hadronic uncertainties, there exists several other challenges such as short distance contributions, non-local effects below the charmonium, the hadronic non-local effects, the non-factorizable effects arising due to spectator scattering and finite width effects [3][4][5][6][7][8] that can cause b → s ℓ + ℓ − neutral processes more difficult to access theoretically. Although these corrections tend to increase the discrepancy in the branching fractions, the normalized angular observables such as P ′ 5 and other LFU sensitive ratios are mostly insensitive to these corrections. Moreover, a global fit including all these corrections are still awaited [8,9].
Similar to the family of neutral decays with charged leptons in the final state undergoing b → s ℓ + ℓ − quark level transition, there also exist another family of flavor changing neutral transitions with two neutral leptons in the final state. Study of rare processes mediating via b → s νν quark level transitions are important for several reasons. First, these processes are theoretically cleaner than the corresponding neutral current decays with two charged leptons in q 2 bins Theoretical predictions Experimental measurements Deviation R K [1.1, 6.0] 1 ± 0.01 [10,11] 0.846 +0.044 −0.041 [20,21] ∼ 3.1σ  [10,11] 0.660 +0.110 −0.070 (stat) ±0.024 (syst) [23] ∼ 2.2 − 2.5σ 1 ± 0.01 [10,11] 0.52 +0.36 −0.26 (stat) ±0.05 (syst) [24] [1.1, 6.0] 1 ± 0.01 [10,11] 0.685 +0.113 −0.069 (stat) ±0.047 (syst) [23] 1 ± 0.01 [10,11] 0.96 +0.45 −0.29 (stat) ±0.11 (syst) [ the final state as they do not suffer from hadronic uncertainties beyond the form factors such as the non-factorizable corrections and photon penguin contributions. Second, the b → s ℓ + ℓ − and b → s νν transition decays are very closely related not only in the SM but also in beyond the standard model physics. Hence, study of these decay modes theoretically as well as experimentally will be crucial to look for potential new physics proposed to explain the anomalies present in b → s ℓ + ℓ − transition decays. Study of b → s νν transition decays are experimentally challenging because of the presence of the neutral leptons which leave no information in the detectors. There exist a few experiments that predict the upper bound of the branching ratio (BR) of B → K ( * ) νν decays. The initial experimental study on b → s νν channels was done by BaBar [38] in 2004, where the upper limit of B(B + → K + νν) and B(B + → π + νν) at 90% CL were reported using the hadronic reconstruction method. Later the results were updated in 2008 [39] and in 2013 [34], respectively. Similarly, the first result by Belle [40] was published in the year 2007. Subsequently the results were updated in 2013 and 2017. So far all the measurements used tagged approaches where the second B meson that is produced in e + e − → Υ (4S) → BB was explicitly reconstructed either in a hadronic decay or in semileptonic decay [34,36,37]. This approach of tagging, in principle, suppresses the background events and results in a low signal reconstruction efficiency which is typically below 1%. Very recently, Belle II used a novel technique based on the inclusive tagging methods and exploited the topological features of B + → K + νν decays. This inclusive tagging method has helped to identify B + → K + νν from seven dominant background processes of the generic B mesons decays. It improves the signal efficiency by 4% at the cost of higher background levels in comparison to earlier methods. An upper bound of B(B + → K + νν) < 4.1 × 10 −5 [35] at 90% CL is reported very recently. Combining this with earlier measurements from Belle and BaBar, the estimated world average for B(B + → K + νν) is reported to be (1.1 ± 0.4) × 10 −5 [35]. We summarize all the results in Table I. Our main aim is to explore the consequences of b → s ℓ + ℓ − anomalies on several b → s νν transition decays in a model independent effective theory formalism. Theoretical study on b → s νν transition decays are limited as compared to b → s ℓ + ℓ − transition decays [16,[41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58] and to the b → c ℓ ν decays [59][60][61][62][63][64][65][66][67][68][69][70][71]. It is well known that b → s ℓ + ℓ − and b → s νν transition decays are related not only in SM but also in beyond the standard model physics. In beyond the standard model physics, they are related via SU(2) L gauge symmetry and can be best exploited using SM effective field theory (SMEFT) formalism. The concept of SU(2) L gauge symmetry was established earlier in few literatures [72][73][74][75][76][77][78][79][80] to provide a simultaneous explanation of b → c ℓ ν and b → s ℓ + ℓ − anomalies. The authors in Ref. [72,73] point out that by assuming the new physics scale much larger than the weak scale, the operators can be made invariant under SU (3) C × SU (2) L × U (1) Y gauge group. There arises two consequences. First, the left-handed fermion fields must be replaced by SU(2) L doublet and second, there will be two new physics operators that are invariant under SU(2) L . As a result, these new physics operators lead to different type of contributions to the neutral current and the charged current interactions, which, in turn, can be used to explain R K and R D ( * ) anomalies simultaneously. We list out few more relevant literatures on b → s νν transition decays [74, where such connections have been addressed. More specifically, in Ref. [81], the authors did study B → K ( * ) νν decays in the SM as well as in several NP models such as MSSM, modified Z/Z ′ penguins, single scalar extension. The authors also pointed out the correlations of K → πνν, K L → πνν, B → X s ℓℓ and B s → µµ in case of right handed NP. Similarly, in Ref. [82], the authors study B → K ( * ) νν decays in the SM and in several beyond the SM models such as Z ′ model, MSSM, leptoquark model. They also use the model independent SMEFT framework and explored several NP scenarios. Very recently, in Ref. [83], the authors study the implication of b → s ℓ ℓ anomalies on several b → s νν and s → d νν decays. They also discussed the correlation between B → K ( * ) νν and K → πνν decays in the case of minimal flavor violation. In Ref. [84], the authors used the SMEFT framework and estimated the new limit on the branching ratios of B → (K, X s ) νν, B s → ϕ νν and B → (π, ρ) νν decays. In Ref. [85], the authors have explored the possibility of enhancement in the branching ratio of B → K νν using NP within scalar and vector leptoquarks and generic vector gauge boson Z ′ model assuming minimal new particle content. In the Ref. [102] the authors investigate B → K ( * ) νν decays in the context of non-standard neutrino interactions. Moreover, in Ref. [103], the authors use the SMEFT framework and perform a global fit to the R D ( * ) and R K ( * ) data. They, indeed, find a strong correlation between the C L operator of b → s νν and C V L of b → c ℓ ν decays.
In the present article, we study the implication of b → s ℓ + ℓ − anomalies on B → K ( * ) νν, B s → (η, η ′ ) νν and B s → ϕ νν decay observables within the SMEFT framework. We give predictions of the branching fractions and longitudinal polarization fraction in the SM as well as in the presence of several 1D and 2D new physics scenarios constructed from SMEFT operators. We perform a global fit to the b → s ℓ + ℓ − data to obtain the allowed new physics parameter space. Our fit analysis include the experimental measurements of R K , R K * , P ′ 5 , B(B s → ϕ µ + µ − ) and B(B s → µ + µ − ) and, in particular, we make use of the latest updated measurements of R K , B(B s → ϕ µ + µ − ) and B(B s → µ + µ − ). In addition, we also check the compatibility of the constrained new physics parameter space with the latest b → s νν experimental data.
So far we don't have many experimental results on b → s νν transition decays. The experimental techniques used for B → K ( * ) νν can be used for B s → (η (′) , ϕ) νν decays as well. Currently, Belle II can be the ideal platform to perform such analysis and predict the upper bound of the branching fractions of these decays. In contrast to the Υ (4S) resonance at Belle where it goes to BB pair, Belle II runs at Υ (5S) as well. The Υ (5S) goes into pairs of B or B s . Belle II has collected samples at the Υ (5S) resonance at an integrated luminosity of 121.4 fb −1 . By taking the cross-section for e − e + → bb and B(Υ (5S)) → B ( * )0 sB ( * )0 s = 0.172 ± 0.030 [104], one can estimate a total (7.11 ± 1.30) × 10 6 B ( * )0 sB ( * )0 s pairs at the KEKB collider. Since only a fraction of these B s decays will survive the kinematics, we expect the statistical uncertainty for B s → (η (′) , ϕ) νν to be more with respect to B → K ( * ) νν decay channel. Moreover, the missing momentum in the final state due to undetected neutrinos can cause difficulties in reconstructing these channels.
The paper is organized as follows. In Sec. II, we start with a brief overview of the standard model effective field theory and write down the effective Hamiltonian governing b → s νν and b → s ℓ + ℓ − decays. In Sec. III, we give predictions of all the observables in the SM and in several 1D and 2D NP scenarios. We conclude with a brief summary of our results in Sec. IV.

A. Standard model effective field theory
So far LHC searches do not provide any direct evidence of new particles close to the electroweak scale. It indirectly suggests the existence of NP at a scale that must lie beyond the electroweak scale. A better way to look for indirect signature of NP in a model independent basis can be attained by considering SM effective field theory (SMEFT) framework. The SMEFT Lagrangian contains all possible set of higher dimensional operators that are built out of the SM fields and are consistent with the SU (3) c × SU (2) L × U (1) Y gauge group. In SMEFT, the higher dimensional operators are suppressed by appropriate power of the NP scale. For a complete set of dimension six and dimension eight operators, we refer to Refs. [105][106][107][108][109]. It is well known that the left handed charged leptons are related to the neutral leptons via SU (2) L symmetry. In this context, the SMEFT framework can be a powerful tool to study the correlation between b → s ℓ + ℓ − and b → sνν transition decays by considering higher dimensional operators. We will consider only dimension six operators in our analysis. Moreover, it is believed that SMEFT analysis may be of great importance if no new particles are observed in LHC [105,107].
The SMEFT Lagrangian corresponding to dimension six operators is expressed as [107] where the relevant operators contributing to both b → s νν and b → s ℓ + ℓ − decays are and the operators contributing only to b → s ℓ + ℓ − decays are At low energy, we can write down the most general ∆F = 1 effective Hamiltonian governing b → s νν and b → s ℓ + ℓ − decays as [81,82] where G F is the Fermi coupling constant, V tb and V * ts are the corresponding Cabibbo Kobayashi Maskawa (CKM) matrix elements. The operators corresponding to b → s νν transition decays are represented by O L and O R with WCs C L and C R , respectively. The operators are where, P L,R = (1 ∓ γ 5 )/2 are the projection operators. In the SM, where the operators O ′ 9 and O ′ 10 exist purely in beyond the SM scenarios. After electroweak symmetry breaking, the low energy SM WCs will get contribution from the dimension six operators of SMEFT. We write C 9,10,L and C 9 ′ ,10 ′ ,R in terms of SMEFT WCs as [82] where, Hq ), c ′ Z = 1 2 ( c Hd ) and ζ ≈ 0.08 is the small vector coupling of Z to charged leptons. We refer to Refs. [82] for all the omitted details. Since we have two undetected neutrinos in the final state, we can only measure differential branching ratio as a function of q 2 for B → P νν decays, where P stands for pseudoscalar meson. Whereas, we can measure differential branching ratio and longitudinal polarization fraction F L in case of B → V νν decays, where V stands for vector meson. All the expressions pertinent for our discussion are reported in Appendix A.

A. Input parameters
For our numerical computation, we use several input parameters such as mass of mesons, quarks and leptons, CKM matrix element |V tb V * ts |, fine structure constant α, Fermi coupling constant G F and the lifetime of parent B (s) meson. For completeness, we report all the relevant input parameters taken from Ref. [110] in Table II. Similarly, for B → K form factor inputs, we use the values obtained in LQCD [2]. Again, for B → K * and B s → ϕ form factors, we use the combined LCSR and LQCD results as reported in Ref. [1]. Moreover, we use the B s → η (′) form factor input parameters from Ref. [111] that are obtained in the LCSR method.  Our main aim is to explore the consequences of b → sℓℓ anomalies on several b → s νν transition decays in a model independent SMEFT formalism. The SMEFT coefficients such as c ql and c Z corresponding to the left chiral currents appear in C 9,10 of b → sℓℓ and C L of b → s νν transitions. Similarly, the SMEFT coefficients corresponding to the right chiral currents such as c dl and c ′ Z appear in C ′ 9,10 of b → sℓℓ and C R of b → s νν transitions, respectively. We consider several NP scenarios based on NP contributions from single operators as well as from two different operators and try to find the scenario that best explains the anomalies present in b → sℓℓ transition decays. To find the best fit values of these NP WCs, we perform a naive χ 2 test with all the b → sℓℓ experimental data. The relevant χ 2 is defined as where O th i represents the theoretical value of each observable and O exp i represents measured central value of the observables. ∆O th i and ∆O exp i represent the errors associated with the theory and experimental values, respectively. We perform two different fit analysis: Fit A and Fit B. In Fit A, we include a total of five measurements for the evaluation of χ 2 , namely, In Fit B, we include only a subset of these five measurement for the evaluation of χ 2 , namely, R K , R K * and B(B s → µ + µ − ). In Table III, we report the best fit values of each SMEFT coefficients in several 1D and 2D scenarios for Fit A and Fit B. We also report the allowed 1σ range of each 1D coefficients. In addition, we report the χ 2 min /d.o.f and the Pull SM = χ 2 SM − χ 2 NP for each scenarios.
• in Fit A, we have used five measured parameters for the the evaluation of χ 2 . Accordingly, the number of degrees of freedom (d.o.f) will be 5 − 1 = 4 for 1D NP scenarios and 5 − 2 = 3 for each 2D NP scenarios. To measure the disagreement of SM with the data, we first obtain χ 2 min /d.o.f in the SM and it is found to be 10.264. The best fit value for each scenarios corresponds to the minimum χ 2 min value. The allowed range of each 1D coefficients at 95% confidence level (CL) is obtained by imposing χ 2 ≤ 37.96 constraint.
• In case of Fit B, we include only three measurement for the the evaluation of χ 2 . Accordingly, the number of d.o.f will be 3 − 1 = 2 for each 1D NP scenarios and 3 − 2 = 1 for each 2D NP scenarios. In the SM, we have found χ 2 min /d.o.f to be 7.189. The allowed range of each 1D coefficients at 95% CL is obtained by imposing χ 2 ≤ 11.98 constraint.
From Table III, it is clear that the coefficients c dl , c ′ Z and ( c dl , c ′ Z ) can not explain the anomalies present in b → sℓℓ data as the minimum χ 2 values obtained for these scenarios are as large as or in some cases larger than that of the SM χ 2 value. Hence we exclude them in the rest of our analysis. There, however, exists few 2D scenarios, namely, ( c ql , c Z ) and ( c ql , c ′ Z ) for which the Pull SM is considerably larger than the rest of the NP scenarios. Moreover, these scenarios have better compatibility with R K , The compatibility of fit results with all b → sℓℓ observables are reported in Appendix B. Again, we do not find any special features in Fit B. For some 2D scenarios, we observe that Fit A serves as a better fit to the data than Fit B. Hence, in all our future discussions, we will mainly focus on the Fit A results. We now proceed to discuss the goodness of Fit A results with the measured values of B(B → K ( * ) νν).
C. Additional constraints from B → K ( * ) νν decays We wish to determine the effect of the SMEFT coefficients on several B → K ( * ) νν decay observables, namely, B(B → K ( * ) νν), F L , R K ( * ) and R K * F L . In Table IV In Fit B, we include only a subset of these five measurement for the evaluation of χ 2 , namely, R K , R K * and B(B s → µ + µ − ). associated with each observable pertaining to B → K ( * ) νν decays in the SM and in the presence of several NP scenarios. To estimate the NP effects, we use the best fit values of the SMEFT coefficients obtained in Fit A of Table. III. In the SM, we obtain the branching fractions for both B → K ( * ) νν decays to be of O(10 −6 ). Similarly, the ratios R K , R K * and R K * F L are found to be equal to 1 in the SM. Hence any deviation from unity in these parameters could be a clear signal of beyond the SM physics. There exist a few experiments that provide the upper bound of the branching ratio of B → K ( * ) νν decays. At present, the upper bounds are found to be B(B → K νν) < 11 × 10 −6 and B(B → K * νν) < 27 × 10 −6 , respectively. Neglecting the theoretical uncertainty, we estimate the upper bound on R K ( * ) to be R K ≤ 2.75 and R K * ≤ 2.89. Our observations are as follows.
• Values of B(B → K ( * ) νν) and R K ( * ) obtained in each 1D NP scenarios with c • In case of 2D scenarios, we observe that the values of B(B → K ( * ) νν) and R K ( * ) obtained with ( c (1) ql , c ′ Z ) SMEFT coefficients are larger than the experimental upper bound. Although it can explain the anomalies present in b → s ℓ + ℓ − data, it, however, can not explain the b → sνν data simultaneously.
ql , c Z ) SMEFT coefficients, the value of B(B → K ( * ) νν) and R K ( * ) are obtained to be quite large. More precise data on B(B → K ( * ) νν) in future will put a severe constraint on these NP scenarios.
• In the SM, R K * F L = 1. Any deviation from unity is a clear signal of the presence of right handed currents. It is evident from the Table IV that the value of R K * F L remains SM like for all the scenarios with left handed currents. However, with the inclusion of right handed currents, its value seem to differ from unity. We see that the value of F L obtained in the presence of ( c Z , c ′ Z ) and ( c ql , c ′ Z ) coefficients are clearly distinguishable from SM prediction at more that 5σ level of significance.
In Fig. 1, we show the q 2 dependence of differential branching ratios and K * polarization fraction for the B → K ( * ) νν decays in the SM and for the best fit values of few selected new physics scenarios. The SM central line and the  and the three ratios R K , R K * and R K * F L in the SM and with the best fit value of each SMEFT coefficients from Fit A analysis of Table. III. FIG. 1: q 2 dependence of differential branching ratios of B → K νν (left) and B → K * νν (middle) decays and longitudinal polarization fraction of K * meson F K * L (q 2 ) (right) in the SM and in few selected NP scenarios. The green, red, orange and black lines correspond to the best fit values of ( c ql , c ′ Z ), respectively. The corresponding SM central value and the uncertainty band is shown with blue. corresponding 1σ uncertainty band is shown with blue color. The green, red, orange and black lines correspond to the best fit values of ( c ql , c ′ Z ), respectively. We observe that the new physics contributions coming from these SMEFT coefficients are quite distinct. In case of F L (q 2 ), the contribution coming from ( c (3) ql , c ′ Z ) and ( c Z , c ′ Z ) are more pronounced and they are clearly distinguishable from SM contribution. In case differential branching ratio, the deviation from the SM prediction is maximum with ( c ql , c Z ) NP scenario. The K * polarization fraction F L value, however, remains SM like as there is no right handed currents.
We wish to quantify our results in terms of the independent parameters R K , R K * and R K * F L . In the presence of ( c (3) ql , c ′ Z ), value of R K is increased by almost ∼ 30% from the SM value, whereas, value of R K * and R K * F L are decreased by almost ∼ 70% and ∼ 60% from the SM prediction, respectively. In case of ( c Z , c ′ Z ), we notice that the values of R K , R K * and R K * F L decreased by almost ∼ 20%, ∼ 80% and ∼ 50% from the SM predictions, respectively. Similarly, ql , c Z ), there is almost ∼ 80% increment in R K and R K * , whereas, the value of R K * F L remains SM like. This is because of the absence of right handed currents in this scenario. Finally, in case of ( c (1) ql + c (3) ql , c ′ Z ), value of R K decreases by almost ∼ 40%, whereas, R K * and R K * F L increase by ∼ 40% and ∼ 10%, respectively.

D. Prediction of Bs → η νν, Bs → η ′ νν and Bs → ϕ νν decay observables in SM and beyond
Study of rare B decays mediating via b → sνν quark level transition is very well motivated as they can provide complimentary information regarding NP in b → s ℓ + ℓ − quark level transition decays. To this end, we study several rare B s meson decays such as B s → η νν, B s → η ′ νν and B s → ϕ νν proceeding via b → s νν quark level transitions in a model independent SMEFT formalism. We give predictions of the branching fractions and ϕ polarization fraction in the SM and in the presence of several NP couplings. For our NP analysis, we choose four NP scenarios, namely, ( c ql , c ′ Z ), that provides the best solutions to the b → s ℓ + ℓ − anomalies. Interestingly, except ( c ql , c Z ) the rest of the scenarios include the effects from right handed currents. In Table V, we report the central values and the corresponding 1σ uncertainty associated with B(B s → (η , η ′ , ϕ)νν) and F L (ϕ) in the SM and in the presence of NP. We obtain the 1σ uncertainty associated with each of these observables by varying the input parameters such as the meson to meson form factors and the CKM matrix elements within 1σ from their central values. In addition, we also quantify the results in terms of R η , R η ′ , R ϕ and R ϕ F L . In the SM, we find the branching ratios of B s → (η , η ′ )νν decays to be of O(10 −6 ), whereas, for B s → ϕ νν decays, it is found to be of O(10 −5 ). The value of ϕ polarization fraction is obtained to be F L = 0.537 ± 0.030. The NP effects can be easily quantified in terms of R η ,η ′ ,ϕ and R ϕ F L . We observe that B(B s → (η , η ′ )νν) increases by almost ∼ 30% in the presence of ( c (3) ql , c ′ Z ), whereas, it decreases by almost ∼ 20% due to the presence of ( c Z , c ′ Z ) NP couplings. Moreover, we observe a ∼ 80% increment in the branching fraction due to ( c ql , c ′ Z ) NP couplings, it decreases by almost ∼ 40% with respect to the SM prediction. In case of B s → ϕνν channel, we notice that B(B s → ϕνν) increases by almost ∼ 80% with ( c Similarly, we observe that the branching fraction increases by almost ∼ 40% in the presence of ( c ql , c ′ Z ) NP couplings, respectively. For F L , we observe maximum deviation from the SM prediction with ( c Z , c ′ Z ) and ( c ql , c ′ Z ) NP couplings. Although, there is slight deviation observed due to ( c ql , c ′ Z ) NP couplings, the deviation from the SM prediction, however, is quite small and it is not distinguishable from the SM.  The branching ratios of B(B s → (η , η ′ , ϕ)νν) and the longitudinal polarization fraction of the ϕ meson F ϕ L in the SM and with the best fit value of few selected 2D SMEFT scenarios of Fit A. The results are also quantified in terms of R η , R η ′ ,R ϕ and R ϕ F L .

SMEFT couplings
In Fig. 2, we display the q 2 dependence of differential branching ratios and ϕ polarization fraction F L (q 2 ) for B s → (η , η ′ , ϕ)νν decays in the SM and in the presence of NP couplings. The SM central line and the corresponding uncertainty band obtained at 95% CL are represented with blue color. The green, red, orange and black lines correspond to the best fit values of ( c Table. III, respectively. Our observations are as follows.
• The differential branching ratio for B s → (η , η ′ ) νν decays is enhanced at all q 2 for ( c the SM prediction at more than 3σ and they are quite distinct from each other. The deviation from the SM prediction is more pronounced in case of ( c ql , c Z ) NP scenario.
• The differential branching ratio for B s → ϕνν decays is enhanced at all q 2 for ( c ql , c ′ Z ), whereas, it is reduced at all q 2 in case of ( c (3) ql , c ′ Z ) and ( c Z , c ′ Z ) NP scenarios. All the NP scenarios are distinguishable from the SM prediction at more than 3σ. The deviation observed is more pronounced in case of ( c ql , c Z ) NP scenarios.
• The ϕ polarization fraction F L (q 2 ) for the B s → ϕ νν decay is distinct from SM only in the presence of ( c Z , c ′ Z ), ( c ql , c ′ Z ) that includes the contribution from right handed currents. In case of ( c ql , c Z ), it is SM like. The deviation from the SM prediction observed with ( c Z , c ′ Z ) and ( c ql , c ′ Z ) is quite significant and they are distinguishable from the SM prediction at more than 5σ. A slight deviation is observed with ( c (1) ql + c (3) ql , c ′ Z ) and it is not distinguishable from the SM prediction.

IV. CONCLUSION
Motivated by the long standing anomalies in B decays with charged leptons in the final state undergoing b → s µ + µ − quark level transition, we study several B meson decays, namely, B → K ( * ) νν, B s → (η, η ′ ) νν and B s → ϕ νν mediating via b → s νν quark level transition. Our primarily goal of this study is intended to analyze the consequences of latest b → s µ + µ − anomalies on b → s νν decays in a model independent approach. We use the standard model effective field theory formalism constructed out of new operators of dimension six corresponding to the arbitrary Wilson coefficients. We study several decay observables pertaining to these decay modes in the SM and in the presence of various SMEFT coefficients in several 1D and 2D scenarios. We perform a naive χ 2 fit to the b → s µ + µ − data, namely, R K , R K * , P ′ 5 , B(B s → ϕ µ + µ − ) and B(B s → µ + µ − ), to find the best fit values of all the SMEFT coefficients in several 1D and 2D scenarios. We observe that the pull corresponding to 2D scenarios are comparatively better than the 1D scenarios. In particular, the fit results of ( c ql , c Z ) and ( c (1) ql + c (3) ql , c ′ Z ) of 2D SMEFT scenarios show better compatibility with all the five b → s ℓℓ measured data. We also check the goodness of the fit results with the additional constraints coming from the experimental upper bounds of B(B → K ( * ) νν). We observe that, although, ( c (1) ql , c ′ Z ) provides a better solution to the b → s ℓℓ data, it, however, cannot explain the existing b → sνν data. The estimated value of B(B → K ( * ) νν) with the best fit value of ( c (1) ql , c ′ Z ) exceeds the experimental upper bounds of B(B → K ( * ) νν).
In case of branching ratio, we observe a significant deviation from the SM prediction in all the four NP scenarios. All the NP scenarios are distinguishable from the SM prediction at more than 3σ significance. The deviation observed with ( c ql , c Z ) are more pronounced. Similarly, for F L , The deviation from the SM prediction observed with ( c Z , c ′ Z ) and ( c ql , c ′ Z ) is quite significant and they are distinguishable from the SM prediction at more than 5σ. Study of B s → (η, η ′ ) νν and B s → ϕ νν decay modes are very well motivated theoretically as well as experimentally as they can, in principle, provide complementary information regarding NP in b → s ℓ + ℓ − decays. Experimental investigations of decay observables in b → s νν in the future will definitely help us in identifying the possible new physics Lorentz structures in b → s ℓ + ℓ − decays. In particular, measurement of F L will be very crucial to not only examine the effects of right handed currents but also to distinguish between various new physics models.

Acknowledgments
NR would like to thank CSIR for the financial support in this work.
Appendix A: Differential decay distribution for B → (P, V ) νν decays The differential decay distribution for B → P νν decays, where P denote pseudoscalar meson, can be written as [81,82] dΓ(B → P νν) Similarly, for B → V νν, it can be written as where, |A ⊥ |, |A ∥ |, |A 0 | are the B → V transversity amplitudes which can be expressed in terms of form factors and Wilson coefficients as Here N is the normalization factor and q 2 is the invariant mass of the neutrino-antineutrino pair. The factor λ is defined as λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ca). The B → P and B → V form factors are defined in terms of f + (q 2 ), V (q 2 ), A 1 (q 2 ), A 2 (q 2 ), respectively. Similarly, longitudinal polarization fraction of the final vector meson can be written as In addition to the differential branching ratio and polarization fraction, we define R P , R V and R V F L where, P and V represent pseudoscalar and vector mesons, respectively. They are expressed in terms of the three real parameters ϵ, η and κ η as [82]. where, where, ρ i is rescaled form factors. It is important to note that the value of R P is independent of decay mode as it only depends on the WCs C L,R . However, R V and R V F L do depend on the decay mode through the factor κ η . The contribution from κ η is observed to be very tiny for B → K * νν and B s → ϕ νν decays.
Appendix B: Best estimates of RK , RK * , P ′ 5 , B(Bs → ϕ µ + µ − ) and B(Bs → µ + µ − ) in the presence of several 1D and 2D SMEFT coefficients from Fit A and Fit B analysis.   Table. III. In the first row we report the experimental central value and the corresponding 1σ and 2σ range of each of these observables.