Addressing $R_{D^{(*)}}$, $R_{K^{(*)}}$, muon $g-2$ and ANITA anomalies in a minimal $R$-parity violating supersymmetric framework

We analyze the recent hints of lepton flavor universality violation in both charged-current and neutral-current rare decays of $B$-mesons in an $R$-parity violating supersymmetric scenario. Motivated by simplicity and minimality, we had earlier postulated the third-generation superpartners to be the lightest (calling the scenario"RPV3") and explicitly showed that it preserves gauge coupling unification and of course has the usual attribute of naturally addressing the Higgs radiative stability. Here we show that both $R_{D^{(*)}}$ and $R_{K^{(*)}}$ flavor anomalies can be addressed in this RPV3 framework. Interestingly, this scenario may also be able to accommodate two other seemingly disparate anomalies, namely, the longstanding discrepancy in the muon $(g-2)$, as well as the recent anomalous upgoing ultra-high energy ANITA events. Based on symmetry arguments, we consider three different benchmark points for the relevant RPV3 couplings and carve out the regions of parameter space where all (or some) of these anomalies can be simultaneously explained. We find it remarkable that such overlap regions exist, given the plethora of precision low-energy and high-energy experimental constraints on the minimal model parameter space. The third-generation superpartners needed in this theoretical construct are all in the 1-10 TeV range, accessible at the LHC and/or next-generation hadron collider. We also discuss some testable predictions for the lepton-flavor-violating decays of the tau-lepton and $B$-mesons for the current and future $B$-physics experiments, such as LHCb and Belle II. Complementary tests of the flavor anomalies in the high-$p_T$ regime in collider experiments such as the LHC are also discussed.


I. INTRODUCTION
A very likely way new physics beyond the Standard Model (SM) could show up in experiments is through anomalous features in the data that cannot be explained by any known SM physics. While some of these anomalies may well be due to statistical fluctuations and/or systematic or theoretical issues that need further understanding, it is also possible that one (or more) such deviation(s) from the SM may well be a genuine signal of some new beyond the SM (BSM) physics. Moreover, given their possible impact and global ramifications, it is worthwhile to carefully scrutinize them in light of possible underlying BSM scenarios.
Of the existing statistically persistent anomalies, particularly prominent ones are the hints of Lepton Flavor Universality Violation (LFUV) in both charged current tree-level and neutral current one-loop rare decays of Bmesons, based on the b → c −ν (with = e, µ, τ ) and b → s + − (with = e, µ) transitions, respectively. In particular, hints for LFUV are seen in the following ratios of branching ratios (BRs): where D * and K * denote excited states of D and K mesons, respectively. These ratios of BRs are interesting observables due to several reasons: (i) Different experiments with completely independent data sets, namely, BaBar [1], Belle [2][3][4] and LHCb [5,6] for R D ( * ) , as well as LHCb [7,8] and Belle [9,10] for R K ( * ) , have reported results for these observables, thus reducing the chances of a statistical fluctuation.
(ii) Such ratios are theoretically clean observables, i.e., with strongly suppressed hadronic uncertainties, thus making them less vulnerable to higher-order quantum corrections [11,12].
In this paper, we take the B-anomalies at face value and address them in a minimal R-parity violating (RPV) supersymmetric (SUSY) framework, with the third generation superpartners lighter than the other two generations (hence dubbed as 'RPV3'). The RPV3 framework was earlier proposed [51] to explain the R D ( * ) anomaly and its possible interconnection with the radiative stability of the SM Higgs boson. The basic idea behind this suggestion comes from the simple observation that the arXiv:2002.12910v1 [hep-ph] 28 Feb 2020 R D ( * ) anomaly involves b and τ , both members of the third fermion family. On the other hand, it is again another third-family fermion, namely, the top quark, that is primarily responsible for the Higgs naturalness problem in the SM. The best known candidate theory for addressing the naturalness problem is (still) low-scale SUSY. However, given the null results in direct SUSY searches at the LHC so far [49,52,53], SUSY solutions to naturalness have become less appealing. As argued in Ref. [51] (see also Ref. [54]), the RPV3 framework which assumes only the third-family fermion to be effectively supersymmetric at the low-scale, while the sfermions belonging to the first two families are decoupled from the low-energy spectrum, provides a simple and minimal solution to the naturalness issue, while being consistent with the LHC constraints so far [55,56], as well as preserving the attractive features of SUSY, such as gauge-coupling unification and existence of dark matter candidate(s).
Here we extend our previous analysis to address both R D ( * ) and R K ( * ) anomalies simultaneously in the RPV3 framework. In addition, we also examine the possibility of addressing two other intriguing and seemingly unrelated anomalies within the same RPV3 framework, namely, (i) the longstanding discrepancy between the SM prediction and experimental measurement of the muon anomalous magnetic moment [49,57], and (ii) the recent anomalous upgoing ultra-high energy (UHE) air showers seen by the ANITA balloon experiment [58,59]. Our goal in this paper is to see if we can carve out a common allowed parameter space within our RPV3 framework where the regions favored by the B-anomalies can overlap with the muon (g − 2) and ANITA-favored regions, while being consistent with all relevant theoretical and experimental constraints. For simplicity, we consider three different versions of our scenario enumerated later, based on certain symmetry arguments, and in each case, we investigate whether there is any available parameter space where all these anomalies can coexist. In one of the three scenarios we find a common overlap region at 3σ confidence level (CL) satisfying all the anomalies, while in the other two simpler scenarios not all the four anomalies could be accounted for, but a combination of either two or three of them could coexist. To the best of our knowledge, this is the first analysis of its kind to unify the B-anomalies with the muon g − 2 and ANITA anomalies in a single testable framework. In passing, let us also mention that while in the past few of years many papers  jointly discuss both R D ( * ) and R K ( * ) anomalies, only a few [15,18,22,37,39] also simultaneously address the muon (g − 2) or ANITA [32], but not both together.
In this work while we are using our RPV3 scenario to understand several of the anomalies because we think it has considerable theoretical appeal for such issues, we will in the following section also voice our concerns regarding experiments and theory pertaining to the results of interest.
The rest of the paper is organized as follows: In Sec-tion II, we give a brief description of the anomalies under consideration. In Section III, we discuss how these anomalies can be explained in our RPV3 setup. Our main numerical results are presented in Section IV. The low-energy experimental constraints used in our analysis are discussed in Section V. In Section VI, we make predictions for the LFV decays of τ and B-mesons. Complementary high-p T tests of the B-anomalies at the LHC are discussed in Section VII. We conclude in Section VIII. In Appendix A, we calculate the extra contribution to R D ( * ) from a light bino in the final state.

II. THE ANOMALIES
In this section, we critically assess the status of each of the experimental anomalies to be subsequently addressed in our RPV3 framework. Although we indulge in a BSM explanation of the anomalies using our RPV3 scenario because we think it has considerable theoretical appeal for such issues, we also voice our reservations and concerns regarding experiment and theory pertaining to the results of interest. Thus, even though the global pull of the B anomalies against the SM may appear to be over 5σ [60] (see Table I), its interpretation as robust evidence of LFUV does not seem compelling to us at this point. It is quite plausible that the resolution of some of these anomalies may well lie in fluctuation of one or more of these experimental results by a few σ. We discuss how the remaining experimental and theoretical issues may be addressed.
In Table I we summarize the anomalies and their pulls. When combining the pulls of several observables we treat all observables as independent degrees of freedom.

A. Hints for LFUV in B Decays
As alluded to above, multiple experimental results from BaBar, Belle and LHCb are pointing to nonstandard sources of LFUV in charged-current and in flavor-changing neutral current (FCNC) decays of Bmesons, based on the b → c −ν and b → s + − transitions, respectively, as measured by the ratios of the BRs, R D ( * ) and R K ( * ) [cf. Eqs. (1) and (2)]. We briefly review the current experimental results on these observables and the significance of the discrepancies with respect to the SM predictions. 1. RD, RD * and R J/ψ Measurements of R D ( * ) exist from BaBar [1], Belle [2][3][4], and LHCb [5,6]. Combining all these, we find Observable R D ( * ) , R J/ψ R K ( * ) (g − 2)µ All but (g − 2)µ All with an error correlation between R D and R D * of ρ = −38%. This is in very good agreement with the average from the Heavy Flavor Averaging Group (HFLAV) [50]. Our R D ( * ) combination is shown in the left plot of Fig. 1.
For the SM predictions we use in our analysis Note that the above uncertainties are somewhat larger than those quoted in e.g. Ref. [12], but we prefer to be conservative for reasons described below. LFUV in the same quark-level transition can also be probed by the decay B c → J/ψ ν. The corresponding experimental result from LHCb reads [63] whereas the SM prediction is [64][65][66][67][68][69] R SM J/ψ = 0.26 ± 0.02 .
The SM predictions of the individual observables disagree with the experimental results by 1.4σ (R D ), 2.5σ (R D * ), and 1.7σ (R J/ψ ). The combined discrepancy between the quoted SM predictions and our experimental average is 3.3σ, as shown in Table I. A few remarks are in order on the theoretical and experimental errors.
(i) Lattice calculations for B → D semileptonic decay are fairly mature by now with stated errors up to around 4% [61,70,71]. However, these quoted errors so far do not include corrections due to soft photons with energy below the experimental threshold; these corrections could be around a few % [72]. These calculations may also need to be corrected for electromagnetic and isospin effects, e.g difference between charged and neutral B decays etc.
(ii) For B → D * semileptonic case there appear to be more serious issues with the theory calculations. An important point that needs to be considered seriously is that since D * carries spin, its production and decay cannot be rigorously factorized. In fact in a construction of the quantum amplitude the production from B → D * ν must be correlated with the final decay, say D * → Dπ, with an appropriate spin-1 D * propagator with its width. It is quite likely that unless this effect is correctly taken into account both the extraction of V cb and R D * suffer from some inaccuracies.
(iii) Moreover, for B → D * transition, a complete lattice calculation with the full q 2 -dependent formfactors does not exist yet and from the lattice perspective given that for a vector final state there are four and not two form-factors (unlike the case of a pseudoscalar final state), it is difficult to see why the theory errors for the case of R D ( * ) should not be appreciably bigger than for R D [cf. Eqs. (5) and (6)].
(iv) There may also be a rather serious concern, at present, on the experimental side, namely, most of the experimental results so far using the leptonic τ → µνν decays involving two neutrinos seem to indicate somewhat larger deviations from theoretical expectations based on the SM compared to two recent measurements, one from LHCb [6] and the other from Belle [3] which use τ → hadron(s) + ν; see Table II. Although the error in each class of measurements is rather large so that the difference in the central values is not different by a significant amount, this difference needs to be better understood as it may originate from some important experimental systematics. Superficially, for example, τ decays involving two neutrinos in the final state appear more vulnerable to backgrounds from higher D * resonances. Theoretical estimates on such contaminations are quite unreliable and they should be subtracted by using experimental measurements, which can be quite challenging.
(v) Another issue on the experimental side that is somewhat disconcerting is that the very first experimental results on the charged-current anomaly came from BaBar [1] and they seem to indicate the most significant deviations from the SM; in contrast, all the Belle results seem to show only mild deviations [cf. Table II]. That is why excluding the BaBar results leads to a smaller pull of only 2.2σ for R D ( * ) , as shown in Table I.
The concerns regarding theory errors voiced above in (i)-(iii) on the charged-current anomaly not withstanding, we also want to stress that at this point the theory errors are subdominant and unlikely to be the sole cause of the discrepancy. Moreover, there is also an intriguing aspect of data from all three experimental groups on these semileptonic decays that is quite interesting and deserves attention. Table II shows all available results to date indicating whether the other B in the event was tagged hadronically or semileptonically and whether the τ decayed leptonically or hadronically. Table II also includes the R J/ψ ratio from similar semileptonic decays of B c to J/ψ τ ( )ν. Altogether there are 11 entries and it is quite remarkable that the experimental central value of the R-ratio for each of these is always without exception above the central value predicted by the SM. Note that the 11 experimental results in Table II are not all completely independent. In fact in some cases, these are just updates of ongoing analyses with more data. Nevertheless, many among these are independent and so the fact that so many experimental measurements are above the SM predictions is quite noteworthy. 1

RK and RK *
The most precise measurement of the LFUV observable R K comes from LHCb [8]: with the dilepton invariant mass-squared in the range 1.1 GeV 2 < q 2 < 6 GeV 2 . The SM predicts R SM K 1 with %-level accuracy [11], corresponding to a ∼ 2.5σ tension with the experimental result.
The most precise measurement of R K * is from a Run-1 LHCb analysis [7] that finds where the first and second values correspond to a q 2 range of 0.045 GeV 2 < q 2 < 1.1 GeV 2 and 1.1 GeV 2 < q 2 < 6 GeV 2 , respectively. The result for both q 2 bins are in tension with the SM prediction, R SM K * 1 [11], by ∼ 2.5σ each. Since the systematic errors here are subdominant, it is reasonable to add the deviations in these two bins in quadrature. Treating the two bins as independent observables we thus find that the deviations from the SM in R K * amounts to about 2.9σ.
Recent results for R K * and R K by Belle have sizable uncertainties and are compatible with both the SM predictions and the LHCb results. For the 1.1 GeV 2 < q 2 < 6 GeV 2 bin, Belle finds [9,10] In the right plot of Fig. 1 we show the combination of the LHCb and Belle results for R K ( * ) in the 1.1 GeV 2 < q 2 < 6 GeV 2 bin compared to the SM prediction. Combining the Belle and LHCb results, we get a net pull of 3.4σ in R K ( * ) as shown in Table I.
Unlike the charged-current semileptonic decays, in the case of FCNC decays B → K ( * ) + − , there are hardly any nagging theoretical issues. So long as the lepton pair invariant mass is larger than about 500 MeV, the SM prediction for the ratio is rather clean and unambiguous. The reservation one may have is only about light lepton invariant mass, say below 500 MeV. Then there is a concern that the electron pair may receive appreciably different radiative corrections from the muon pair [11].
The primary concerns about µ − e universality violation in FCNC is experimental. Of course the effects are only a few σ. Moreover, it is only one experiment, i.e. LHCb, and an independent confirmation by Belle II would be highly desirable. Also, if it is genuine LFUV it ought to show up irrespective of hadronic final states in B-decays. Thus one should see the corresponding b → s FCNC decays materializing into baryonic and other final states, such as Λ b → Λ + − . It also should not depend on the spectator quark. Thus charged and neutral B and also B s decays ought to exhibit similar signs of LFUV. In particular, LHCb already seems to have indications that the observed rate for B s → φµ + µ − is seemingly below "SM" expectations [74] but the absolute rate calculations may suffer from some long-distance (non-local) contaminations, so a direct test of µ − e universality via a measurement of B s → φe + e − would be very valuable.
Let us briefly add that we are primarily focusing on the LFUV anomalies as they are theoretically cleaner and for now we are choosing not to include some other possible indications of deviations from the SM, such as angular observables or absolute rate for B → K ( * ) µ + µ − [60,[75][76][77][78][79]) and also rate for B s → φµ + µ − [74] as in these cases there can be non-perturbative contributions from nonlocal effects especially in the region of low q 2 that are not under full theoretical control.
Before closing this subsection, it is worth pointing out here that the hints of LFUV are only seen in the semileptonic B-decays. Analogous semileptonic decays of charmed mesons do not show any such deviations from the SM. For instance, BESIII has recently reported a measurement of the ratio of BRs in the D + -decay [80], viz.
(14) The (g − 2) experiment at Fermilab [83] is expected to improve the experimental accuracy by a factor of about four in the next few years.
The SM prediction for a µ can be decomposed in contributions from QED, from the electro-weak interactions, from hadronic vacuum polarization and from hadronic light-by-light scattering: The QED and electro-weak contributions are known with high accuracy [84,85] a QED µ = (11, 658, 471.897 ± 0.007) × 10 −10 , (16) The hadronic vacuum polarization contribution can be determined using e + e − → hadrons data and dispersion relations. A recent such analysis gives [86] (see also Ref. [87]) where the first, second and third terms correspond to the LO, NLO, and NNLO contributions, respectively. The value is in good agreement with the findings of a hybrid approach that uses the best part of lattice results along with the best part of the experimental data and continuum dispersion relation data [88], and tends to favor the BSM interpretation of the data. This is particularly significant since in the traditional R-ratio dispersion analysis there is appreciable concern due to the discrepancy of ≈ 2σ between the BaBar data and the KLOE data [89]. Indeed the lattice hybrid approach does not use the somewhat conflicting input data from BaBar or KLOE.
Combining the results collected above leads to a discrepancy between experiment and SM prediction at 3.3σ CL [86]: For this anomaly the next year is likely to be pivotal. The new Muon (g − 2) experiment at Fermilab [83] already seems to have collected about two times the data used by the BNL experiment; the analysis of that accumulated data is expected in the next few months. How this new result compares with the previous BNL result would be crucial for the BSM interpretation.
On the lattice front, about a factor of 3 reduction in the error is anticipated in the next few months by the RBC-UKQCD Collaboration [94] and this could also have a critical bearing on the BSM interpretation. Also phenomenological approaches are pursued both for the hadronic vacuum polarization and the light-by-light scattering contribution [95][96][97][98][99]. At the moment, the so-called "hybrid" method of RBC-UKQCD [88] which uses part of the continuum dispersive calculation and in part the lattice calculation in regions which complement each other seems to tentatively favor the BSM interpretation. But it would be much better if pure lattice techniques can further reduce their error by factor of 2 to 3 so it does not use any input from experiment especially since two of the best experimental results from KLOE and BaBar have ≈ 2σ disagreement between them. Therefore pure lattice calculations with reduced errors would be very welcome in providing input for the fate of the BSM interpretation in muon g − 2. It appears we will need to wait for another year or so for this to happen.
The theory uncertainty on the hadronic vacuum polarization contribution can also be reduced by about a factor of 2 at the proposed MUonE experiment [100,101] which will make a very high-precision measurement of elastic µ − e scattering at a QED-dominated momentum exchange of q 2 = O(100 MeV) 2 . This measurement will be quite robust and invulnerable to any BSM physics that could be responsible for the muon (g − 2) anomaly [102,103].

C. Anomalous ANITA Events
A recent anomalous observation in ultra-high energy cosmic ray (UHECR) air showers made by the ANITA collaboration has also hinted at some BSM physics. Two anomalous upward-going events with deposited shower energies of 0.6 ± 0.4 EeV [58] and 0.56 +0. 3 −0.2 EeV [59] (1 EeV = 10 9 GeV) have been reported. Both these events originate from well below the horizon, with large negative elevation angles of (−27.4 ± 0.3) • and (−35.0 ± 0.3) • , respectively. They do not exhibit phase inversion (opposite polarity) due to Earth's geomagnetic effects, and hence, are unlikely to be downgoing UHECR air showers reflected off the Antarctic ice surface, although there is some uncertainty in modeling the roughness of the surface ice [104][105][106]. Potential background events from anthropogenic radio signals that might mimic the UHECR characteristics, or unknown processes that might lead to a non-inverted polarity on reflection from the ice cap are estimated to be ≤ 0.015, resulting in a 3σ evidence for direct upward-moving Earth-emergent UHECR-like air showers above the ice surface [59]. This poses considerable difficulty for interpretation of such events within the SM framework due to the low survival rate ( 10 −6 ) of EeV-energy neutrinos over long chord lengths in Earth∼ 7000 km, even after accounting for the probability increase due to ν τ regeneration [58,107,108]. Moreover, as pointed out in earlier studies [106,109], the strength of isotropic cosmogenic neutrino flux needed to account for the two events is in severe tension with the upper limits set by Pierre Auger [110,111] and IceCube [112,113]. Therefore, a BSM explanation with an anisotropic astrophysical source with some exotic generation and propagation mechanism of upgoing events is desirable to solve the ANITA anomaly, provided it stands further scrutiny after more data release from future ANITA flights.

III. RPV EXPLANATION OF THE ANOMALIES
As we suggested before [51], RPV SUSY is a particularly interesting theoretical framework to address the flavor anomalies. For one thing, for the charged-current tree level indication of BSM physics, RPV is a natural candidate and if LFUV is involved then this is especially so. Moreover, since members of the third family, namely, b and τ are involved in B → D ( * ) τ ν, it may well be that this anomaly is a hint that it is related to the issue of the radiative stability of the Higgs mass which is an important persistent problem of the SM. Motivated by the naturalness arguments and to keep the RPV SUSY scenario minimal, for reasons of simplicity, we have suggested that it may well be that the third generation superpartners are the lightest. In that scenario proton stability issues are less relevant and for that reason too R-parity breaking is a viable option [54]. Lastly, we have shown that even with such an economical setup involving effectively only one generation of superpartners a very attractive feature of SUSY, namely unification, is retained. Finally we also want to remark that our objective is to use the latest experimental data with the current set of indications to constrain as best we can the parameters of this interesting theoretical construct.
We start from the LQD part of the RPV SUSY Lagrangian that contains the λ couplings which are relevant for an explanation of R D ( * ) [51,[114][115][116][117][118][119][120] and R K ( * ) [116,117,119,[121][122][123][124]: As we will see below, for explanations of the R K ( * ) anomaly and the (g − 2) µ anomaly it is useful to also include the effect of the LLE part of the RPV SUSY Lagrangian which contains the λ couplings [117]: One thing to keep in mind is that the λ couplings are antisymmetric in the first two indices: λ ijk = −λ jik . Also note that the simultaneous presence of λ and λ couplings is consistent with proton decay constraints, as long as we do not switch on the λ (U DD-type) couplings. Following Ref. [51], throughout this work we will assume, for minimality, that the third-generation squarks, sleptons and sneutrinos are considerably lighter than the first and second generation ones. Integrating out the heavier SUSY particles we therefore can neglect the first and second generation sfermions, as their effect is suppressed by a higher mass scale in the RPV3 scenario. Outλ of the 27 independent RPV couplings λ ijk in Eq. (21) and the 9 independent λ ijk in Eq. (22), there are 19 λ -type and 7 λ-type couplings that involve light third generation sfermions, namely, the right-handed sbottom b R , lefthanded stop t L , left-handed tau-sneutrino ν τ and both left-and right-handed staus τ L,R . We will treat these five masses as free parameters in our numerical analysis in Section IV. In addition, we require a light long-lived bino ( χ 0 1 ) for the ANITA anomaly. As for the choice of couplings, we first analyze each of the experimental anomalies discussed above in the RPV-SUSY context and show the dependence of the observables on the relevant couplings. Then in the following Section IV, we present three different scenarios for our parameter set-up and the corresponding fit results.
A. Explanation of RD and RD * In Ref. [51] we had identified BSM contributions to b → cτ ν transitions in the RPV setup, which can arise at the tree level from sbottom exchange [cf. Fig. 2(a)]. The sbottom exchange leads to contributions to the decay amplitude that have the same chirality structure as the SM contribution and thus modify R D and R D * in a universal way. Here we note that in the presence of the LLE couplings, also diagrams with light sleptons, in particular a light left-handed stau, can contribute to the decays [cf. Fig. 2(b)]. However, in the scenarios we will consider below, the left-handed stau will be fairly heavy and the corresponding contributions will be negligible. We will therefore focus on the sbottom contribution from the diagram in Fig. 2

(a).
It is important to note that R D and R D * measured by BaBar and Belle correspond to ratios of the tauonic decay modes to an average of the muonic and electronic modes, while the LHCb measurements are ratios of tauonic to muonic modes. Using the notation from Ref. [117], we find in our setup v = 246 GeV is the Higgs VEV and V ij are the CKM matrix elements.
In the case of the B-factories, we instead have where ξ e parameterizes the relative weight of the electronic and muonic decay modes in the R D ( * ) measurements at the B-factories. We note that ξ e can in principle be different for each experimental analysis but we expect ξ e ∼ 50% (see e.g. [125]). We explicitly checked that varying ξ e has no significant impact on our results. This is due to the fact that µ − e universality in b → c ν decays is observed with high accuracy. Translating the results from Ref. [126] into our RPV scenario, we have Therefore, it is an excellent approximation to combine the LHCb and B-factory results as done in Section II A 1.
In that case we find 2 both for the LHCb and the B-factory expressions [cf. Eqs. (23) and (25)].
1. Implications of the observed q 2 distribution and of the D * polarization Recently, Ref. [127] in an interesting study have included q 2 (where q is the 4-momentum carried by the leptonic pair) and also the longitudinal polarization of the D * in addition to the integrated rates in order to discriminate against models. To analyze the data in a model independent manner they allow all possible current structures in the weak Hamiltonian subject only to the constraint that only left-handed neutrinos are involved in the interaction; thus, with the operators The parameter space explaining the R D ( * ) data automatically explains the R J/ψ data, because the underlying transition is the same b → c ν. Therefore, we do not discuss the R J/ψ fits separately.
and weighted by the corresponding Wilson coefficients C i . In this representation, the operator O V L is of special significance as it encapsulates the SM interaction. In their study of the existing experimental data, Ref. [127] find that the simplest solution to the chargecurrent anomaly is with a small non-vanishing value of C V L ≈ 0.08, with all other C's equal to zero. This has the important consequence that the polarization of the D * or for that matter of the τ will not be different from the SM. Recently Belle collaboration reported, for the longitudinal polarization of the D * [128] which is in mild tension of about 1.6 σ with the SM which predicts [129][130][131] F L (D * ) SM = 0.46 ± 0.03 .
In the past ≈ 2 years, Belle collaboration has also attempted to measure the polarization of the τ and found [3] At this point this result on tau polarization within its large errors is consistent with the SM expectations of [132] The fact that the experimentally observed q 2 distribution in the semileptonic B → D ( * ) decays supports a small non-vanishing value, C V L ≈ 0.08 is also very significant for our RPV3 BSM scenario. One can see from Eq. (21) that as long as only the LQD interactions are relevant, in RPV3 the dimension-6 effective interaction for the semileptonic decays is essentially identical to the (V − A) × (V − A) structure of the SM effective Hamiltonian with the difference being just in the overall coefficient. Whereas in the SM the overall coefficient is G F × V cb / √ 2, RPV3 has the overall coefficient λ ×λ /m 2 b . Thus the coefficient, C V L ≈ 0.08 is consistent with m b ≈ 2 TeV for λ 0.5, as we will explicitly see below in the numerical fits.

Bino Contribution
There is an additional contribution to the B → D ( * ) ν decays that can arise in our RPV scenario. If the bino, χ 0 1 , is extremely light and has a very long lifetime (as motivated by an explanation of the ANITA anomaly, see Section III D below), then the decays B → D ( * ) χ can be open and mimic the B → D ( * ) ν decays. In this case, we could have the B → D ( * ) χ processes via either left-handed stau or right-handed sbottom exchange which effectively give contributions to operators of the form (cP R b)(τ P R χ) and (cσ µν P R b)(τ σ µν P R χ). Details are given in Appendix A. Evaluating these contributions, we find that the effect is rather small: BR(B → Dτ χ)/BR(B → Dτ ν) SM 1% and thus this extra channel cannot significantly affect R D ( * ) .
We also did a similar analysis with regard to the possible contribution from the extra bino-channel to the longitudinal polarization fraction F L (D * ) of B → D * ν. 3 We do expect a non-zero correction to F L (D * ) coming from the extra Bino channel because of the different operators that are involved. However, we find that the effect is tiny ∆F L (D * ) 8 × 10 −5 which is not significant given the large uncertainties in the current experimental value [cf. Eq. (30)] and the SM value [cf. Eq. (31)].

B. Explanation of RK and RK *
The BSM contributions to the rare decays B → Kµ + µ − and B → K * µ + µ − are conveniently described by shifts in the Wilson coefficients of semileptonic fourfermion operators in the effective Hamiltonian [60] and Q 9,10 are obtained from Q 9,10 by replacing P L → P R . Recall that in the SM, the Wilson coefficients are universally for all = e, µ, τ . Fits of R K and R K * show that the observed pattern can be accommodated with BSM in the coefficients (C 9 ) e , (C 10 ) e , (C 9 ) µ , (C 10 ) µ , as long as BSM in the primed coefficients is subdominant, otherwise it leads to an anti-correlated effect in R K and R K * , contradicting the current data.
Global fits of all relevant data on rare B decays find a particular consistent BSM picture which is characterized by non-standard effects in muonic coefficients in the combination of Wilson coefficients (C 9 ) µ = −(C 10 ) µ [60] (see also [75][76][77][78][79]). As we will see below, our RPV SUSY scenario will generate contributions to both (C 9 ) µ = −(C 10 ) µ and (C 9 ) µ = −(C 10 ) µ . Such a scenario provides an excellent fit to the data for the following values [60] Note that the combination (C 9 ) µ −(C 10 ) µ corresponds to BSM that mainly affects left-handed muons. All other coefficients are compatible with zero at the 2σ level. The correction to the SM values of the Wilson coefficients C SM 9 −C SM 10 4 is at the level of −15% for the muon flavor, while for the electron flavor the corrections vanish. The above BSM values for the coefficients explain not only the observed values for R K and R K * , but also other (theoretically less clean) anomalies in rare B decays, like the angular observable P 5 or the branching ratio of B s → φµµ (see Refs. [60,[75][76][77][78][79]).
Note that in our RPV setup the simultaneous presence of muon and electron couplings would likely lead to extremely stringent constraints from searches for µ → e transitions, like the µ → eγ decay, or µ → e conversion in nuclei [133]. We therefore focus on muonic couplings only.
In the considered RPV scenario, contributions to b → s transitions arise both at the tree level and the loop level. Tree-level exchange of stops (see Fig. 3a) gives contributions to the wrong chirality Wilson coefficients. In agreement with Ref. [122] we find where α em is the fine structure constant. The abovediscussed preferred ranges for these coefficients in Eq. (40) translate into the approximate bound In addition, there are various classes of 1-loop contributions to the b → sµµ decays that we consider (see Fig. 3b-f). There are loops with right-handed sbottoms and W bosons (Fig. 3b), with two right-handed sbottoms (Fig. 3c), as well as with stops and sneutrinos (Fig. 3d). 4 These contributions are all governed by the λ RPV couplings. In the presence of the λ RPV couplings there are additional 1-loop effects (as first pointed out by Ref. [117]). We take into account loops with right-handed sbottoms and staus (Fig. 3e), as well as with left-handed sneutrinos (Fig. 3f). All those diagrams give contributions to the left-handed Wilson coefficients and therefore can in principle explain the anomalies in R K and R K * . Summing up all these contributions we get [117,122,123] ( where the X and X factors are the following combinations of RPV couplings: cussed in Ref. [123], assuming that winos are sufficiently heavy in our RPV3 scenario. Note that this does not necessarily spoil the gauge coupling unification in RPV3 [51], as the renormalization group (RG) running is logarithmic, and O(10 TeV) winos (and similar mass for the gluino to satisfy the stringent LHC constraints), along with light bino (and Higgsinos), are acceptable. X bs = λ 331 λ 321 + λ 332 λ 322 + λ 333 λ 323 , X µµ = |λ 213 | 2 + |λ 223 | 2 + |λ 233 | 2 , X µµ = |λ 231 | 2 + |λ 232 | 2 + |λ 233 | 2 , X bµ = λ 331 λ 231 + λ 332 λ 232 + λ 333 λ 233 , X sµ = λ 321 λ 231 + λ 322 λ 232 + λ 323 λ 233 , It is intriguing that the RPV setup produces BSM contributions that follow the (C 9 ) µ = −(C 10 ) µ pattern that is preferred by the data. Note that the first term in (43) arises from the sbottom-W boxes and has the wrong sign, i.e. it always worsens the agreement with data. The coupling combinations that enter in the other terms are constrained for example by B s mixing and B → Kνν. The last two terms in (43) involve both the λ and λ couplings (the last one was not included in Ref. [117]). These additional terms provide more freedom to explain the R K ( * ) anomalies in the context of RPV SUSY. An explanation of the anomalies requires negative C 9 . Given that V ts −0.04, this in turn requires some of the λ or λ couplings to be negative.
C. Explanation of (g − 2)µ The contributions to (g − 2) µ can arise in RPV SUSY both from the λ and λ couplings. The diagrams involving λ are shown in Figs. 4a-d and those involving λ (with sleptons and leptons in the loop switched to squarks and quarks) are shown in Figs. 4e-h. In our RPV3 setup, these contributions can be summarized as [134] We find that the net contribution from the λ-dependent terms is typically dominant, as the relevant λ couplings tend to be less constrained than the λ couplings (cf. Table IV). It is worth noting here that the electron g − 2 also has a ∼ 2.4σ discrepancy between the experimental measurement [135] and SM prediction [136], due to a new measurement of the fine structure constant [137]: It is difficult to explain the opposite sign with respect to ∆a µ using RPV couplings only. However, within the minimal supersymmetric SM (MSSM), it is possible to explain ∆a e by either introducing explicit lepton flavor violation [138] or using threshold corrections to the lepton Yukawa couplings [139] or arranging the bino-slepton and chargino-sneutrino contributions differently between the electron and muon sectors [140]. Since this is independent of the RPV sector, we do not include the electron (g − 2) in our subsequent discussion.

D. Explanation of ANITA Upgoing Events
We interpret the ANITA upgoing anomalous events [58,59] as signals from the decay of long-lived bino in RPV SUSY, produced by interactions between UHE neutrinos and nucleons/electrons inside Earth matter via exchange of a TeV-scale sparticle mediator. As first discussed in Ref. [109], the whole process could be divided into four sub-processes, namely, the generation of the bino on the far-side of Earth, its propagation through Earth matter, followed by its decay in the atmosphere and signal detection at ANITA. The generation and the decay of bino could both be described by Fig. 5 with one of the vertices coming from either λ or λ sector, while the other being U (1) Y gauge coupling g . The contribution from the λ sector involving the ν − e interactions turns out to be sub-dominant in our case due to the choice of small λ i13 and the lower probability to have an s-channel resonance for ν − e interactions as compared to ν − q interactions, since all three down-type quark PDFs are sizable at EeV energies [109]. After the bino is generated, it is required to have a long lifetime to travel through a chord length of ∼ 7000 km, as inferred from the ANITA events. The decay width of bino is parameterized by its mass m χ 0 1 , the mediator sbottom or stau mass and the λ, λ , g couplings as: As mentioned above, the λ-contribution is subdominant and we will only keep the λ terms in Eq. (46). The longevity of the bino in our model comes from two things: First, it is electrically neutral and interacts with the SM particles very weakly. Secondly, it is produced with a very high Lorentz boost factor of γ ∼ 10 6 . So as long as Eq. (46) gives a mean lifetime of ∼ 10 ns in the bino rest-frame, which translates to a lifetime ∼ 0.01 s in the lab frame, it can safely propagate through a chord length of ∼ 7000 km without losing much energy. This happens for a relatively light bino with m χ 0 1 ∼ a few GeV. After propagating the chord length of a few thousand km, as it reaches near the surface of Earth, it undergoes a 3body decay back to quarks (or leptons) and neutrinos, followed by hadronization of the quarks, producing an extensive air shower due to the Askaryan effect [141]. The radio signal from the air showers is then detected by the ANITA balloon detector.
The expected number of events can be estimated as follows [109]: where we have taken T = 53 days for the total effective exposure time, Φ ν (E ν ) = 2 × 10 −20 (GeV · cm 2 · s · sr) −1 for the cosmic neutrino flux, 5 and A eff · ∆Ω is the effective area integrated over the relevant solid angle, averaged over the probability for interaction and decay to happen over the specified geometry. The effective area contains all the information of the geometry, decay width of bino and the cross section for the bino generation process; see Ref. [109] for the explicit expression. From Eq. (47), we know that the overall event number N is a function of m χ 0 1 , m b R and λ ij3 for our RPV3 scenario. Therefore, comparing the simulated event numbers with the ANITA observation of two anomalous events gives us the best-fit parameter region at a given CL.

IV. NUMERICAL RESULTS
After examining Eqs. (24), (43), (44) and (46), all the relevant parameters contributing to the anomalies discussed above in our RPV3 scenario are summarized in Table III. For the purpose of convenience, we also collect the dominant terms in the expressions for anomalies in Table III. The same is done for the relevant experimental constraints in Table IV which we discuss in detail in the following Section V).
As mentioned before, in our RPV3 setup, there are six free mass parameters relevant for the anomalies, namely, As for the choice of RPV couplings shown in Table III, we apply certain symmetry rules to reduce the number of parameters. We consider the following three different cases and present our numerical fit results in each case.

A. Case 1: CKM-like Structure
This symmetry is inspired by the observed hierarchy in the CKM mixing matrix in the quark sector. This is brought out most clearly in the Wolfenstein parameterization of the CKM-matrix [142], where the first generation plays the central role. The coupling of first-to-first generation quarks are of order one, whereas the coupling of the first to the second carries a suppression factor of λ sin θ C ≈ 0.23. Similarly, the coupling of second generation to the third carries a suppression of λ 2 , and the coupling of first generation to the third carries a suppression factor of λ 3 . Inspired by this structure, in our RPV 5 This is consistent with the recent upper bound for transient sources, based on a joint analysis of ANITA detection and Ice-Cube non-detection results [113]. To be more specific, our transient anisotropic flux value Φν integrated over the small solid angle ∆Ω corresponding to the uncertainty in the observed elevation angles for the ANITA events is Φ int = 4.9 × 10 −24 (GeV · cm 2 · s) −1 at 0.5 EeV, to be compared with the upper bound on Φ int ≤ 8 × 10 −24 (GeV · cm 2 · s) −1 for the steady analysis [113].
scenario which is third-generation-centric, we postulate the λ -couplings to be of the form with λ 333 ∼ O(1) and each time any of the three indices {i, j, k} differs from 3, we pay an appropriate factor of , which is a tunable small parameter in the model. A similar rule is applied to the λ sector, where we choose for the nonzero λ's: 6 where i < j and λ 233 ∼ O(1). This setup reduces the number of couplings from 27 (λ ijk )+9 (λ ijk )=36 to only 3, namely, In Fig. 6, we show a benchmark scenario for Case 1 in the (m b R , ) plane, while fixing the other free parameters as follows: The two coupling values are mainly chosen to simultaneously maximize the overlap region where the anomalies can be explained, as well as to evade the current existing bounds. A particularly stringent constraint comes from τ → νν (see Section V G) which involves both λ 333 and λ 233 couplings, and the masses of right-handed stau m τ R and right-handed bottom, m b R . Thus we need to change λ 333 and λ 233 together so that their overall effect mostly cancels to give a narrow allowed window from τ → νν. These two couplings are set as large as possible so that the cancellation takes place, and meanwhile gives a maximized overlap region as long as the other constraints do not become too strong. The masses chosen here are consistent with the 13 TeV LHC constraints [49]. The stau mass is chosen close to the experimental limit of 900 GeV to obtain the maximally allowed parameter space, while satisfying the bound from τ → νν, i.e. choosing a larger stau mass will shrink the available parameter space shown in Fig. 6, while a smaller stau mass will shrink the window of the allowed region from τ → νν. As for the choice of the sneutrino mass, from Table IV we could see that the term involving m ντ contributes dominantly to the B s − B s bound and thus to alleviate this bound, we set m ντ at a relatively larger value of 15 TeV. We choose m τ L to be 10 TeV to suppress the possible contribution to R D ( * ) from LLE couplings. Also, m t L is set at Observable Parameter dependence Relevant terms 4 TeV to suppress the tree-level contribution to b → s as mentioned in Eq (42). The favored regions for explaining the R D ( * ) , R K ( * ) and ANITA anomalies are shown in Fig. 6 by cyan, red and orange-shaded regions, with the 2σ and 3σ regions depicted by thin and thick solid contours respectively. The parameter is required to take negative values in order to find overlap between R D ( * ) and R K ( * ) regions. This is due to the fact that we need C 9 < 0 to fit the data [cf. Eq. (39)] and since C 9 is composed of odd powers of with positive definite factors [cf. Eq. (43)], this inevitably sets negative. On the other hand, the R D ( * )favored regions are divided into two different branches due to the polynomial dependence of λ ijk and λ ijk upon [cf. Eq. (24)]. As for the ANITA-favored region, it is mostly governed by the bino mass which is set at 2.0 GeV, apart from the sbottom mass and λ couplings.
Other shaded regions in Fig. 6 with dashed/dotted boundaries are the relevant experimental constraints; see Section V and Table IV for details. The main constraints come from B s − B s mixing [50] and B → Kνν [25,143,144] measurements. Note that the B s -meson mixing bound has a branch-cut feature which is due to the cancellation between the terms in Eq. (65). Somewhat less constraining bounds come from B → τ ν [50], D − D mixing [145], τ → νν [49], and Z → ¯ data [49]. Finally, the vertical shaded region below m b R < 1.0 TeV is excluded from direct sbottom searches at the LHC [49].
The overlap region between R D ( * ) , R K ( * ) and ANITA is highlighted by the green shaded region in Fig. 6 around (m b R , ) = (2.2 TeV, −0.015). This is remarkable, given how simple the coupling choice is, even though it occurs only at the 3σ CL. However, a major drawback of this scenario is that the (g − 2) µ -favored region lies around − ∼ O(10), which is far away from our CKM-like assumption that | | 1; therefore, it is not shown in Fig. 6.

B. Case 2: Flavor Symmetry
The second benchmark point we study is inspired by a U (2) q × U (2) flavor symmetry proposed in Ref. [117]. In this case, the values of λ ijk and λ ijk couplings are decided by the specific flavon VEVs in the model. They have the generic structure λ ijk ∼ c ijk and λ ijk ∼ c ijk , where the and values may differ for each coupling, while c ijk and c ijk are O(1) free parameters. Here we choose (pink), and τ → νν (blue). The overlap region simultaneously explaining the R D ( * ) and R K ( * ) anomaly is shown by the green shaded region, and the region also explaining the ANITA anomaly along with R D ( * ) and R K ( * ) is shown by the green shaded region with thick boundaries. a simplified version of this model and assume that c ijk and c ijk are strictly equal to the overall scales of λ and λ respectively, i.e. λ ijk ∼ λ and λ ijk ∼ λ with and fixed by the flavor structure parameters as indicated in Ref. [117]. Moreover, to accommodate R K ( * ) , we choose λ 333 to be negative and set it as a free parameter to be fit numerically. All other λ ijk values are fixed by the overall scale λ , i.e. λ 1jk = λ 211 = λ 231 = λ 213 = λ 311 = λ 331 = λ 313 0 , where q ≈ m s /m b 0.025, q ≈ q m d /m s 0.005 and 1 [117]. Similarly, all λ ijk values are fixed by the overall scale λ, i.e.
where 0.004 and S 0.06 [117]. Therefore, this choice is equivalent to taking 3 free parameters for the couplings, i.e.
{λ 333 , λ , λ} , The choice of the combination of λ, λ 333 , m τ R and m τ L is mainly due to the consideration of enlarging the overlapping region and avoiding current constraints. Larger magnitude of λ and λ 333 will push R K ( * ) region downwards and R D ( * ) upwards giving a larger overlap. However, both B s − B s mixing, and the B → Kνν and τ → νν decays are sensitive to the choice of these four parameters (see Table IV) and most of them become stronger as we increase the couplings. The more complicated relation comes from τ → νν which involve λ, λ 333 , m τ R and m b R . As described in Eq. (75), the two dominant terms of τ → νν, involving λ, m τ R and λ 333 , m b R respectively cancel each other. Thus we choose m τ R = 3.0 TeV to maintain a window in the right range of m b R ∼ 2.5 TeV where R K ( * ) , R D ( * ) and ANITA overlap. A smaller m τ R will shrink the window and move it to the left, but choosing m τ R to be larger will cause the R K ( * ) region to shrink, due to nonlinear dependence on m τ R . Meanwhile, we increase λ, λ 333 simultaneously so that their effects on τ → νν window mostly cancel. To avoid RG running problems (i.e. hitting the Landau pole too close to the TeV-scale), λ 333 is set at its largest possible magnitude of −3.5. This large coupling results in severe B s − B s mixing bound and to alleviate this, we choose m ντ to be 9 TeV. m τ L is chosen, different from m τ R , at 10 TeV, as mentioned in the previous case to suppress the possible contribution to R D ( * ) from LLE couplings. The color scheme for the shaded regions is the same as in Fig. 6. Now we also show the 2σ (3σ) preferred region for (g − 2) µ at the upper left corner of Fig. 7 by the thin (thick) blue line with the arrow pointing into the allowed region. The horizontal hatched region is theoretically disfavored from perturbativity constraint on λ ≤ √ 4π. The location and shape of the favored regions for R D ( * ) and R K ( * ) anomalies are different from Case 1 mainly due to the fact that the parameter planes are different. In Fig. 6, the y-axis shows the parameter which plays the role as the relative scale between two λ or two λ couplings, while in Fig. 7 the y-axis shows the overall scale for the λ -couplings. Generally speaking, the overall scale could be larger but the relative scale should be heavily suppressed due to the polynomial dependence. Therefore, the overlap region in Fig. 7 has λ ∼ 0.8, as compared to that in Fig. 6 which has ∼ −0.01.
Also note that in Case 2, the 3σ allowed region for ANITA shrinks dramatically, in both m b R and λ 233 directions, which is mainly due to the structure of the λ couplings in Eq. (53). The favored region shrinks in the m b R direction because there are larger λ couplings and thus the simulated number of events for ANITA gets more sensitive to change of m b R . Shrinking in the λ direction is a combined effect of the structural change of the λ s and the change of y-axis from relative scale ( in Case 1) to overall scale (λ 233 in Case 2).
The overlap region of R D ( * ) , R K ( * ) and ANITA anomalies is marked by the green block around (m b R , λ ) ∼ (2.5 TeV, 0.8). No overlap could be achieved with (g−2) µ region in this parameter setup. We find that (g − 2) µ is most sensitive to m ντ and we have tried an extreme case of setting m ντ at the current LHC lower bound of 900 GeV [49], which does expand the (g − 2) µ region downward but not enough to have an overlap while in the meantime B s meson mixing bound becomes much severe and rules out the whole parameter region. Thus in this case (g − 2) µ cannot be accounted for.
The bounds also appear differently in Case 2 than in Case 1 due to the change of y-axis. The most stringent bounds in this case are τ → νν [49] and B s meson mixing processes [50]. Similar to Fig. 6, the branch-cut feature in the B s -meson mixing bound is due to the cancellation between the terms in Eq. (65).

C. Case 3: No Symmetry
In this final benchmark scenario, we do not invoke any symmetries. Instead, we adopt a pragmatic approach to choose our parameters so that we maintain the necessary freedom to explain all the anomalies while satisfying all experimental constraints. At the same time, we want to keep the total number of free parameters the same as in the other two cases, i.e. six mass parameters and three couplings. Thus, we try to equalize the non-zero parameters as much as possible. We end up with the following 3 free coupling parameters, with all the other λ and λ couplings are set to be very small (essentially zero in practice). As shown in Fig. 8, our benchmark point in this sce- while we vary the remaining two parameters λ and m b R to find the common overlap region for R D ( * ) , R K ( * ) , (g − 2) µ and ANITA. We are able to do so around (m b R , λ ) = (3.0 TeV, 0.3). The overlap region is highlighted as the green block in Fig. 8. In this parameters setup, R D ( * ) and R K ( * ) are brought together mainly by setting a large negative λ 223 = −1.5. When combined with setting λ 333 = 0, this setup results in R D ( * ) being dominated by −X µ c ∼ −λ 223 λ /m b R , which gives a positive contribution as we want. Meanwhile, for R K ( * ) , the dominant term is the second term from Eq. (43) , which gives a negative contribution as required. The (g − 2) µ -favored region in this case is vastly expanded compared to Case 2, and covers pretty much the entire parameter space shown in Fig. 8. This is mainly due to the choice of small m ντ and the multiple O(1) λs, where we choose λ to be −3.5, which give larger overlap compared to the positive value due to the dominant λ term contribute to the denominator of R D ( * ) . This setting guarantees the dominant contribution to be the λ terms in Eq. (44) and thus the subdominant λ terms could have a much larger range. In this case, the effect of m τ R on g − 2 and R K ( * ) is gone due to the van-  Fig. 7. In addition, the D 0 → µ + µ − constraint is shown by the blue shaded region (marked by the dashed blue boundary). The 2σ (g − 2)µ region covers almost the entire shown parameter space, so the 3σ region is not shown. Also, as in Fig. 7, the horizontal hatched region is theoretically disfavored from perturbativity constraint on λ ≤ √ 4π.
ishing couplings λ k23 . So the only influence of m τ R is on D − D mixing bound, which inversely depends on m 2 τ R (see Section V D). Therefore, we simply set m τ R =2 TeV, same as in Case 1. On the other hand, from the same consideration of reducing the effect of LLE coupling on R D ( * ) like previous two cases, we set m τ L =10 TeV.
The relevant bounds, including B → τ ν [50], D − D mixing [145], B s − B s mixing [50], B → Kνν [25,143,144] and D 0 → µ + µ − [146], are also shown in Fig. 8 by dark blue, magenta, gray, yellow and blue shaded regions respectively, while the LHC bound on sbottom mass is shown by the vertical brown-shaded region. In this case, the most stringent constraints come from B s −B s mixing and D → µµ which shrink the overlap region substantially. The B s − B s mixing, as mentioned earlier in Case 1 and Case 2, is a typical bound for our RPV3 model trying to explain the B-anomalies since the relevant couplings λ i33 , λ i23 and λ i32 all contribute to B-meson mixing. The branch-cut feature of the B s −B s mixing bound seen in Figs. 6 and 7 is absent in Fig. 8 because in this case there is no cancellation in Eq. (65), as the third term dominates due to the choice of small sneutrino mass. On the other hand, the D → µµ bound is crucial mainly due to the important role of λ 223 in this particular Case 3. Note that in this case the τ → νν bound is not relevant due to vanishing couplings λ 333 = λ 233 = 0; see Section V G for more details.

V. CONSTRAINTS
For the record let us briefly mention that just before the advent of the two asymmetric B-factories, the general perception was that RPV had so many parameters and that it was so completely unconstrained that it can accommodate just about anything; see e.g. p.921, Table  13.6 in Ref. [147]. On the contrary, what we will show here is that the situation now has dramatically improved, thanks to the enormous experimental and theoretical progress in the past two decades. In fact, despite the many parameters our RPV3 scenario is remarkably wellconstrained as we discuss below so much so that more accurate measurements of say R D ( * ) preserving the central value could have appreciable adverse consequences at least for the version of RPV that we are now finding to be favorable.
In this section, we discuss all relevant constraints on our RPV3 scenario shown in Figs. 6, 7 and 8, with the parameter dependence and dominant terms in the corresponding expressions summarized in Table IV. A. B → τ ν and Bc → τ ν In the notation of Ref. [117], for B ± → τ ν, we have where the sum over l is for all flavors of neutrinos in final state, and which includes processes involving both LLE and LQD vertices. Notice that the extra factor in front of the second term is due to the difference between vector and pseudoscalar current. The B → τ ν channel has been experimentally measured and the most updated results is reported in Ref. [50]: with a SM prediction of [51]: Comparing these numbers for the experimental measurement and SM calculation, a constraint could be imposed on the combination of RPV couplings and masses of sparticles in Eq. (59). In Figs. 6 and 8, this constraint has been shown by the blue shaded region with dashed dark blue boundary. The constraint turns out to be not relevant for the parameter choice in Fig. 7. Similarly, the decay B c → τ ν also gets a contribution from Eq. (60) with V uj /V ub replaced by V cj /V cb . This channel has not been measured and may not be measured in the near future. Previously, constraints have been imposed using the life time of B c , τ Bc = 0.51 × 10 −12 s [49]) and a 10%-40% estimate on the maximal allowed BR(B c → τ ν) [148][149][150]. We do not use this channel as a constraint, since we find that in our scenarios B → τ ν gives always stronger bounds. For completeness, we provide the predictions for BR(B c → τ ν) for our benchmark points: 25.6% (Case 1), 0.9% (Case 2), and 2.0% (Case 3). The corresponding ratio of the BR(B c → ντ ) between the RPV3 scenario and SM is found to be BR(Bc→τ ν)RPV3 BR(Bc→τ ν)SM = 34.2 (Case 1), 1.2 (Case 2), and 2.7 (Case 3).

B. B → K ( * ) νν and B → πνν
Tree-level exchange of sbottoms contributes to the decays B → Kνν and B → K * νν. Taking into account decay modes into different neutrino flavor combinations we get for the branching ratios: with the top loop function X t = 1.469 ± 0.017 [151] and s w being the weak mixing angle. Note that we consider both b L and b R exchanges, a feature only valid for final state with two neutrinos. Depending on the chosen benchmark, this equation simplifies into different forms and we use m b L = m b R for numerical purposes. A bound for this ratio has been given by [25,144] R B→K ( * ) νν < 5.2 at 95% CL, which is adopted for our parameter setting and indicated in Figs. 6, 7 and 8 as the yellow-shaded regions with dashed yellow boundary. An analogous expression holds for the decays B → πνν and B → ρνν: However, the experimental bounds on those decays are much weaker than the B → K ( * ) νν bounds and are always satisfied for the parameter choice we have, and hence, are not shown in Figs. 6, 7 and 8.

Constraint
Parameter dependence Relevant terms

C. Bs − Bs Mixing
Here, RPV contributions can arise at the tree level from sneutrino exchange, as well at the one-loop level from box diagrams with sbottoms, sneutrinos, or stops. Based on the derivation from Ref. [117], we have: where we update the hadronic P factors from Ref. [152] with the latest lattice input from Ref. [153] giving: The mass difference ∆M Bs in neutral B s meson mixing is measured with excellent precision, ∆M Bs = (17.757 ± 0.021) ps −1 [50]. The SM prediction, on the other hand, has sizable uncertainties stemming mainly from the hadronic matrix elements and the CKM matrix element V cb [154]. For the SM prediction, we use the latest lattice average of hadronic matrix elements from Ref. [155] (see also Refs. [153,156,157] where f Bs is the B s decay constant, andB Bs a so-called bag parameter. For the CKM matrix element we use |V cb | = (41.0 ± 1.4) × 10 −3 , which is the conservative PDG average of recent inclusive and exclusive determinations [49]. We find ∆M SM Bs = (19.3 ± 1.7) ps −1 . This is in good agreement with the experimental value and a recent SM prediction based on light cone sum rule calculations [158]. Combining our SM prediction with the experimental result we obtain the following bound at 95% C.L.
The combination of the experimental result and the SM prediction puts a bound on the parameter space by confining the possible contribution to ∆M Bs from RPV. This bound is indicated as the grey-shaded region in Figs. 6, 7 and 8.

D. D − D Mixing
Here contributions from stau or sbottom loops could arise from λ ijk ≡ λ ilk V jl couplings. The effective Hamiltonian for the loop-level contributions to D − D mixing from RPV is described as [159] Using this we can derive a bound on the RPV parameters and relate them to x D ≡ ∆M D /Γ D (where Γ D is the mean decay width of D meson): Combining this with the experiment result: x expt D = (5.8 ± 1.9) × 10 −3 [145], we get the bound for (m b R , λ ), which is denoted by the pink-shaded region in Fig 8. E. D 0 → µ + µ − Tree-level contributions from sbottom exchange to this rare D 0 decay width can be expressed as [116]: In case 3, this bound become most important and the expression reduces to a function of λ , λ 223 and m b R . An updated upper bound on this branching ratio [146] is set at 7.6 × 10 −9 at 95% CL and the corresponding bound is shown as the light purple area in Fig 8. In other two cases, this bound is subdominant and is not shown in Fig  6 and Fig 7. F. Z → ¯ This process gets modified by top-sbottom loops. More specifically, Z will decay to an off-shell tt pair both of which through LQD give a lepton and a sbottom. The sbottom line closes to form a loop, while the product of the process becomes a ll final state. Due to different i index in λ i33 , we may have different flavor final states (Z → τ τ or Z → µµ) and even flavor violating final states such as Z → τ µ. A change in the Z decay process from the SM prediction will affect the ratio of the vector and axial-vector couplings of the Z boson with different lepton flavors. Experimental measurements on these couplings are given in Ref. [49] as: The contributions to these ratios from RPV model are given by [117] where δg i j is a simplification of Eq. (30) in Ref. [160] where we only keep the top Yukawa-related terms. It is denoted as follows: Taking i, j both equal to 3 and using Eqs. (71) and (72), we derive a bound on the parameter space, as shown by the vertical pink-shaded region in Figs. 6 and 7. This bound is not shown in Fig. 8 because the choice λ 333 ∼ 0 makes it irrelevant for Case 3.
In principle, bounds could also be put on λ 333 λ 233 by evaluating the experimental bounds on the LFV branching ratio of Z → τ µ process. However, the current experimental bound for BR(Z → τ µ) is of order of 10 −5 while the contribution to this branching from RPV is typically < 10 −7 [160]. Therefore, no substantial bound can be put from the flavor violating Z coupling. Also worth noting is that the W couplings could also be altered by RPV loop processes. However, such bounds from the W coupling variations are not shown here since they are not as strong as the bound from τ → νν process [160], which is described in the next subsection.

G. τ → νν
The LLE coupling will result in the change of the decay rate of τ → eνν and of τ → µνν via the exchange of τ R . This effect could be tested by the ratio: Based on the derivation from Ref. [117], in the SM+RPV case, we have: This can be used to put constraints on the parameter space when combined with the experimental values [161]: (R τ /e τ ) exp = 1.0060 ± 0.0030 .
The corresponding bound is displayed in Figs. 6 and 7 as dark blue region, while it is not shown in Fig. 8 because in Case 3 this bound becomes irrelevant due to λ 333 ∼ λ 323 ∼ 0.
However, as pointed out by Refs. [163][164][165][166], BSM effects from both R-parity conserving and violating terms could contribute to this channel either directly via oneloop diagrams or indirectly via RG running. Considering the direct RPV contribution only, we take the bound in Ref. [163] adopting it to the updated measurement [50], which gives: Substituting the benchmark mass values for our three cases we find that the constraints are λ ij2 λ ij3 1.64, 1.15, 1.00 and λ i2j λ i3j 7.14, 25.94, 1.01 respectively for Case 1, 2 and 3, while the actual values of these coupling products we have for all these cases are O(0.01). Thus the b → sγ constraint is always satisfied for all our benchmark points. The weakness of this bound could be understood both from the partial cancellation between the two terms in Eqs. (80) and (81), and from the dependence of the upper bounds on the sparticle masses.

I. Neutrino Mass
The trilinear RPV couplings in Eqs. (21) and (22) contribute to neutrino masses at one-loop level through the quark-squark and lepton-slepton loops, respectively. Using the general expression [167][168][169] and dropping the terms involving the first two generation sfermions, we obtain: where ( m d LR ) 2 and ( m e LR ) 2 are the left-right squark and slepton mixing matrices respectively, given by (and similarly for ( m e LR ) 2 in terms of A e and y e ), where A d,e are the soft trilinear terms, y d,e are the Yukawa couplings, and tan β = v u /v d is the ratio of the VEVs of the two Higgs doublets in the MSSM.
In the basis in which the charged lepton masses and the down quark masses are diagonal, it is customary to assume that the A-terms are proportional to the Yukawa couplings, i.e. A d 33 = A b y b and A e 33 = A τ y τ . With this substitution, Eq. (82) simplifies to where m b and m τ are the average sbottom and stau masses. We must ensure that the trace of the M ν matrix in Eq. (84) (i.e. the sum of its eigenvalues m νi ) should satisfy the cosmological bound on the sum i m νi 0.1 eV [170]. For the three cases discussed earlier, we find that this requires (A b,τ −µ tan β) O(0.5 MeV) for Cases 1 and 2, while for Case 3, the upper bound is relaxed to about a GeV. With this choice, the neutrino mass constraint can be readily satisfied, and therefore, we do not include it in Figs. 6, 7 and 8.

VI. LFV PREDICTIONS
In this section, we make predictions for LFV decay modes of the τ -lepton and rare decays of the B-mesons for our three benchmark cases, anticipating that future experiments like Belle II [171] might be able to test some of these predictions.
For illustration purpose, let us first consider the simple Case 1 with CKM-like coupling structure. Concretely, we plan to implement the third-generation centric rotations due to RPV interactions in complete analogy with the SM. We just have to bear in mind that in RPV3 we interchange the role of the first and third generations compared to the SM. Moreover, as in the SM, the order parameter, λ 0.23 in the Wolfenstein representation [183] can be used for flavor rotations in our RPV3 set up. In particular, when RPV interactions τ → u and µ → u are involved, in a similar fashion, these can be accompanied by suppression factors, say, 31 21 , where 31 ≈ λ 2 and 21 ≈ λ 3 . In line with our thinking that superpartners of third generation quarks are the lightest, these rotations may be analogous to V ub and V cb respectively with the product causing a suppression in the rate of order λ 10 ≈ 4 × 10 −7 . Thus, with a mediator mass of M 1.6 TeV (20 times heavier than W ), this can re-sult in a branching ratio of O(10 −12 ) and be completely consistent with the current bounds.
So clearly there is significant model dependence involved at this stage and we will just need to dig the appropriate effects of these rotations in flavor space from the experimental data. In this third-generation centric RPV3 model of ours, it would seem that τ → µss final states may be less suppressed than those with uu, dd and sd. The τ → µss process gives rise to distinctive final states such as τ → µφ [K + K − ]. Making the ad hoc assumption that these couplings go as 32 ≈ λ ≈ 0.23, a mediator mass of 1.6 TeV can lead to Yet another simple mode where we can make a statement about the branching ratio is τ → µK 0 . This can be normalized conveniently to the SM mode τ → νK + which has a branching ratio of about 7 × 10 −3 . Note that as above the τ → s RPV vertex will carry a suppression of λ. The µ → d vertex couples second generation to first; thus this is analogous to V cb in the SM and the rate goes as (λ 3 /|V cb |) 2 ≈ λ 2 . Putting all the factors together, one finds BR(τ → µK 0 ) ≈ 5 × 10 −10 .
Another interesting example is τ → µµµ. This arises at tree level via use of LLE couplings of RPV [cf. Eq (22)]. We again assume a suppression of 32 λ ≈ 0.23. Then again for a mediator mass of 1.6 TeV, we can get where the index k in λ 32k is either 2 or 3. In this calculation we have assumed that when the third-generation sneutrino couples to two muons which are from second generation, there is a suppression of O(λ 2 ) in the vertex.
In Table V, we summarize the above-mentioned treelevel LFV decay modes of τ , with the dominant coupling dependence in our RPV3 setup and the model predictions in each of the three cases discussed above, corresponding to the parameters in the overlap regions shown in Figs. 6, 7 and 8. Also shown are the current experimental constraints on each channel. As can be seen, all the three benchmarks are consistent with the current bounds, while some of the predictions might be accessible at future B-factories. Note that the tree-level BRs in Case 1 turn out to be much smaller than our naive estimate discussed above, because we have used the value of | | = 0.02 for the overlap region in this case (cf. Fig. 6), which is a factor of 10 smaller than the simple choice of | | λ ≈ 0.23.

B. LFV via Loop Decays of τ
There are interesting LFV loop decays of τ that we can estimate quite easily by using existing calculations of b → sγ [184] and of b → s + − [185]. These calculations are relevant as the virtual top quark dominates in b decay as well as in τ decays because of the simple picture of mixing angles that we have adopted. For τ → µγ, we find where I ∼ 0.075 is a loop function [186] that depends on the ratio of top mass to sbottom mass. Thus, we estimate that BR(τ → µγ) ∼ 10 −9 .
In an analogous fashion, in the loop decays τ → µ + − (for = µ, e), the virtual top-quark dominates as in the case of b → s + − . This leads one to the estimate, Thus, we conclude that the loop contribution to BR(τ → 3µ) is about a hundred times smaller compared to the tree contribution estimated above. Our RPV3 predictions for the loop-level τ LFV decays are also summarized in Table V for all three cases, along with the corresponding experimental bounds.

C. LFV Decays of B-mesons
We briefly discuss here some illustrative examples of distinctive LFV decays of B-mesons that proceed via tree processes in our RPV3 scenario and can be estimated readily.
First example we want to discuss is b → sτ µ. In our third-generation centric setup with flavor rotations as explained above and with the reminder that when it comes to flavor rotations accompanying RPV interactions, quarks and leptons are treated on the same footing. With that in mind b → sτ µ results in the exclusive modes, B → K( * )µτ [187] or B s → φµτ . In this case the b → τ vertex has no suppression but s → µ, both being second generation fermions, carry a suppression of λ 2 . Thus once again taking the mediator mass of 1.6 TeV and taking the normalizing weak decay B → νX c with BR ≈ 11% which involves a suppression factor of V 2 cb results in BR(b → sµτ ) ≈ 7 × 10 −7 . The BRs of the corresponding exclusive manifestations are likely a factor of 5 smaller as indicated in Table V. Another related extremely interesting example is B s → τ µ. Let us normalize this to the SM mode of B → τ ν. In this case though the LFV BR(B s → τ µ) carries a suppression of λ 4 , that is more than compensated by V 2 ub factor in the normalizing mode. Thus, again for a mediator mass of 1.6 TeV, we get BR(B s → τ µ) ≈ 8.4 × 10 −8 .
It is to be stressed that these BRs of flavor violations involving τ µ final states of B and B s are rather large and future experiments like Belle II and LHCb should be able to constrain them quite well.
For completeness, we also list in Table V our RPV3 predictions in the lepton flavor conserving FCNC decay modes, such as B s → + − . Although in all our benchmark cases, the model predictions are quite small for these channels, it is conceivable that in a less restrictive setup, the RPV contributions could be within reach of upcoming experiments.

VII. HIGH-pT PREDICTIONS AT THE LHC
As shown in Ref. [51], simple crossing symmetry arguments can be used to establish a high-p T modelindependent test of the R D ( * ) anomaly in CMS and AT-LAS experiments; see also Refs. [188][189][190]. The basic idea is that the underlying quark-level process for R D ( * ) is b → cτ ν, which by crossing symmetry also implies the existence of the processes like gc → bτ ν and gb → cτ ν, which can be searched for in the high-p T LHC experiments. We do not wish to repeat the same analysis here, but would like to stress the point that similar modelindependent tests can be done for the R K ( * ) anomaly as well; see also Refs. [191,192].
Specifically, the underlying parton-level process for R K ( * ) is b → s + − (with = e, µ), and by crossing symmetry, the following processes must also occur in the pp collisions at the LHC: (i) bs → + − , (ii) gb → s + − and (iii) gs → b + − (we generically denote both quarks and anti-quarks by q). So if the R K ( * ) anomaly were true, we must also have an anomaly in these channels, which might be observable depending on the signal to background ratio.
The signal in each case can be analyzed in the four-fermion setup with the vector operators defined in Eq. (34) and an effective mass scale of O(TeV). Scalar and tensor operators do not work here, unlike in the R D ( * ) case, because of the B s → µ + µ − constraint, which can be helicity-suppressed only in the vector case. Moreover, as the LHC center-of-mass energy is comparable to the mass scale of the effective operator being studied here, it is more accurate to use an explicit mediator. To be concrete, we use the RPV3 model as our benchmark, where one of the squarks serves as the mediator for the processes listed above and couples via the λ -type LQD interaction. Note that the λ-type LLE interactions do not enter here at the leading order, since we must have b and s quarks in the external legs to relate to the R K ( * ) anomaly.
For case (i), we find that the SM background is overwhelmingly large, mostly coming from Z → + − . Imposing an invariant mass cut on M to exclude the Zmass window helps, but still we find it difficult to achieve a signal significance of more than 3σ. Similarly, for case (ii), the signal cross section is suppressed due to the bottom-quark PDF in the proton. So the best case scenario is case (iii), where the additional b jet in the final state provides a better handle on the signal over background. Some simple kinematic distributions for the corresponding signal and background are shown in Figure 9.
Here we have used the minimal trigger cuts: p j,b T > 20 GeV, |η j,b, | < 2.5 and ∆R > 0.4 and an average btagging efficiency of 70%. From the distributions, we find that although the M distribution can distinguish the RPV signal from the SM background to some extent, the striking signature of RPV comes in the distribution of the invariant mass M b with the correct lepton combination. This is because in RPV3, we have the pro-cess gs → t R µ − → bµ + µ − through the λ -couplings [cf. the penultimate term in Eq. (21)]. Thus, if kinematically allowed, the stop can be produced on-shell in the s-channel, followed by its decay into bµ + , thereby giving a resonance peak in the M bµ + distribution, as shown in the bottom right panel of Fig. 9 for a representative stop mass of 1 TeV. Using this resonance feature, it is possible to achieve more than 3σ signal significance for the overlap region shown in Fig. 8 at the 14 TeV LHC with an integrated luminosity of 3 ab −1 .
Similarly, the RPV3 explanation of the ANITA anomaly can be independently tested at colliders. The key thing to note here is that we require a light, longlived bino with a rest lifetime of about 10 ns to explain the ANITA anomaly [109] (see Section III D). The bino can be produced at the LHC from either gluino or squark decay through gauge interactions, followed by the 3-body decay of bino into two quarks and a lepton (through LQD coupling) or into three leptons (through LLE coupling). For TeV-scale gluinos and squarks, a GeV-scale bino will have a boost factor γ ∼ 10 3 at the LHC and will have a decay length of ∼ 100 m in the lab frame. This leads to distinct displaced vertex signatures [193], which should be accessible to dedicated long-lived particle searches at the LHC [194,195].

VIII. CONCLUSION
Taking the reported B-physics anomalies, as well as the muon anomalous magnetic moment and ANITA anomalous upgoing air shower events at face value, we examine the exciting possibility that these anomalies in vastly different systems could actually be connected by a single underlying BSM framework. In particular, we point out that the origin of these anomalies might be related to a third-generation-centric BSM scenario, which could also address the SM Higgs naturalness issue, while preserving all the good features of a generic supersymmetric framework. In a promising minimalist approach, we consider the so-called 'RPV3' scenario, wherein only the superpartners of the third-generation SM fermions are relatively light, at (sub)TeV scale, whereas all other sparticles (except the lightest neutralino) are much heavier and do not play a significant role in explaining the anomalies.
We have considered three benchmark cases for this RPV3 setup and analyzed the reduced parameter space to carve out the regions favored by each of the above mentioned anomalies, while making sure that all relevant experimental constraints are satisfied. We find that some combination of these indication(s) of deviations from the Standard Model can be explained in all three cases, but finding an allowed overlap region between all of them may only be possible in one of the three cases studied here. Nevertheless, it seems remarkable to us that such an overlap region exists at all (see Fig. 8), given the stringent experimental constraints from a large number of low and high energy processes on the masses and couplings. As shown in Fig 10, there are two possible diagrams for the b → cτ χ 0 1 decay process. In theory, the vertex with bτ c may also give rise to a third diagram but we only consider the third generation of sparticles to be light, so that diagram will not be considered here.
The effective Hamiltonian for the two diagrams could be written respectively as: where λ ijk = λ ilk V jl , A τ χ is a linear combination of the SU (2) L and U (1) Y gauge couplings in the sparticle sector, and B b χ is related to the U (1) Y gauge coupling g ; for the exact definitions of A and B, see Eqs. (8.88a) and (8.88c) in Ref. [197]. Via Feirz transformation and using the property of charged conjugate operator, we could rewrite H 2 in a similar fashion as H 1 : [4(cRb)(τ Rχ) − (cRσ µν b)(τ σ µν χ)] .

(A3)
To compare with the SM contributions, we can rewrite both H 1 and H 2 into the form of standard Wilson coefficients multiplying the corresponding operators [cf. Eq. (28)]: with the following Wilson coefficients which are the dominant ones in our RPV3 scenario: and the operators defined as: With these expressions, we could easily compare the contribution to the decay width from the extra channel and from the SM. We define the following ratio: where H s V,0 , H s V,t , H s S and H s T are the helicity amplitudes defined in the same way as in Appendix B of Ref. [127]. Taking the Wilson coefficients from Eqs. (A5) and (A6), masses of the final state particles as in our benchmark cases, helicity amplitudes and form factors for B and D mesons from Refs. [127,198], we find that for m χ 0 1 = 2 GeV, R extra int = 0.6%, which is insignificant compared to the SM and typical RPV contributions discussed in Section III A.