Unified explanation of flavor anomalies, radiative neutrino mass and ANITA anomalous events in a vector leptoquark model

Driven by the recent experimental hints of lepton-flavor-universality violation in the bottom-quark sector, we consider a simple extension of the Standard Model (SM) with an additional vector leptoquark $V_{\rm LQ}({\bf 3},{\bf 1},2/3)$ and a scalar diquark $S_{\rm DQ}({\bf 6},{\bf 1},4/3)$ under the SM gauge group $SU(3)_c\times SU(2)_L\times U(1)_Y$, in order to simultaneously explain the $b \to s \ell^+ \ell^-$ (with $\ell=e,\mu$) and $b \to c l^- \bar \nu_l$ (with $l=e,\mu,\tau$) flavor anomalies, as well as to generate small neutrino masses through a two-loop radiative mechanism. We perform a global fit to all the relevant and up-to-date $b \to s \ell^+ \ell^-$ and $b \to c l^- \bar \nu_l$ data under the assumption that the leptoquark couples predominantly to second and third-generation SM fermions. We then look over the implications of the allowed parameter space on lepton-flavor-violating $B$ and $\tau$ decay modes, such as $B_s \to \ell^\pm_i \ell^\mp_j$, $B \to K^{(*)} \ell^\pm_i \ell^\mp_j$, $B_s \to \phi \ell^\pm_i \ell^\mp_j$ and $\tau \to \mu \gamma (\phi)$, respectively. Minimally extending this model by adding a fermion singlet $\chi({\bf 1},{\bf 1},0)$ also explains the ANITA anomalous upgoing events. Furthermore, we provide complementary constraints on leptoquark and diquark couplings from high-energy collider and other low-energy experiments to test this model.


I. INTRODUCTION
Over the last few years, several B-physics experiments, such as the LHCb [1-10], as well as the B factories BaBar [11,12] and Belle [13][14][15][16][17] , have reported a number of deviations from the Standard Model (SM) expectations at the level of (2-4)σ [18] in the rare flavor-changing neutral-current (NC) and charged-current (CC) semileptonic B-meson decays involving the quark-level transitions b → s + − (with = e, µ) and b → cl −ν l (with l = e, µ, τ ) respectively, which provide intriguing hints of new physics (NP) beyond the SM (BSM). The non-observation of any new heavy BSM particle through direct detection at LHC experiments makes these indirect hints a powerful tool in the NP exploration. A more careful analysis of these tantalizing hints for lepton-flavor-universality violation (LFUV), taking into account the possibility of statistical fluctuations and yet unknown systematic and/or theoretical issues, is absolutely essential to confirm or rule out the possible role of NP in the B-sector. However, given their possible impact on NP searches, it is worthwhile to scrutinize these experimental results at their face values in light of possible NP scenarios.
Since the above-mentioned anomalies associated with b → s + − and b → clν l transitions probe different NP scales [19], most of the theoretical studies in the literature have attempted to address either the NC or the CC sector, but not both on the same footing. Only a few specific models, mainly those involving the color-triplet leptoquark (LQ) boson  which allows tree-level couplings between quarks and leptons, have been successful in explaining both kinds of flavor anomalies simultaneously (see also Refs. [55][56][57][58][59][60][61][62][63][64][65] for other plausible simultaneous explanations of the B-anomalies).
As discussed in Refs. [22,44,48,66], models with a single scalar LQ cannot address both these anomalies simultaneously. With the aim of understanding the experimental observations linked with both types of processes in a common framework, here we consider a simple extension of the SM by adding a single vector leptoquark (VLQ) V LQ which transforms as (3, 1, 2/3) under the SM gauge group SU (3) c × SU (2) L × U (1) Y . The existence of VLQ at low energy can be theoretically motivated from many ultraviolet (UV)-complete frameworks [67], such as grand unified theories [68][69][70][71], Pati-Salam model [72][73][74], technicolor model [75][76][77][78], etc. In the literature, the flavor anomalies have already been investigated in the VLQ scenario [20, 22-24, 27-30, 35, 38-40, 42, 45, 46, 57, 60, 79-88]. Here we update this discussion with the latest experimental data and also minimally extend the VLQ model by introducing a scalar diquark (SDQ) S DQ (6, 1, 4/3), to explain the light neutrino mass generation through a two-loop radiative mechanism. Moreover, the observation of LFUV generically implies the existence of lepton-flavor-violating (LFV) decay modes [89]. Even though some theoretical works [20,90] contradict this precept, the link between LFV and LFUV persists in several models. In this connection, we will also investigate the LFV decays of neutral and charged mesons, as well as of the tau lepton, in conjunction with the LFUV parameters for the VLQ case.
Moreover, as it turns out, minimally extending this VLQ-SDQ model with an additional fermion singlet χ(1, 1, 0), we can also accommodate the recent ANITA anomaly [91,92]. Finally, we provide complementary constraints on leptoquark and diquark couplings from collider and other low energy experiments to test this model. The organization of this paper is as follows. In Section II , we present the effective Hamiltonian in terms of dimension-six operators, describing b → s + − and b → cτν l quark-level transitions. In Section III , we discuss our model framework and the NP contributions arising due to the exchange of VLQ. The set of relevant observables that have been used to constrain the NP parameters are listed in Section IV . The numerical fit to the new Wilson coefficients from the existing experimental data on b → s + − and b → cτν τ processes is presented in Section V . Section VI contains the implication of VLQ on the LFV B and τ decay modes. In section VII , we discuss a two-loop radiative neutrino mass generation with the VLQ and SDQ particles. The SDQ signal at LHC is illustrated in Section VIII . Section IX presents an explanation of the ANITA anomaly in our model with an additional fermion singlet. Our conclusion is given in Section X . In Appendix A , we list the experimental data used in our numerical fits. Appendix B (C) contains the expressions required for B → K ( * ) i j LFV decays. The loop functions for τ → µγ are provided in Appendix D .

II. GENERAL EFFECTIVE HAMILTONIAN
The effective Hamiltonian responsible for the CC b → cτν l quark level transitions is given by [93] H CC eff = where G F is the Fermi constant, V cb is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element and C l X are the Wilson coefficients, with X = V 1,2 , S 1,2 , T , which are zero in the SM and can arise only in the presence of NP. The corresponding dimension-six effective operators are given as where f L(R) = P L(R) f are the chiral fermion (f ) fields with P L(R) = (1 ∓ γ 5 )/2 being the projection operators.
The effective Hamiltonian mediating the NC leptonic/semileptonic b → s + − processes can be written as [94,95] Here V tb V * ts is the product of CKM matrix elements, C i 's are the Wilson coefficients [96] and O i 's are the dimension-six operators, expressed as where α em is the electromagnetic fine structure constant. The SM has vanishing contribution from primed as well as (pseudo)scalar operators, which can be generated only in the BSM theories.

III. MODEL FRAMEWORK
We build a simple model by extending the SM by a color-triplet, SU (2) L -singlet vector leptoquark V LQ (3, 1, 2/3) for explaining the flavor anomalies (see Section V). We also add a color-sextet, SU (2) Lsinglet SDQ S DQ (6, 1, 4/3) to explain the neutrino masses by radiative mechanism (see Section VII), with some interesting collider signatures as well. Finally, we add a fermion singlet χ(1, 1, 0) to account for the ANITA anomaly (see Section IX). The relevant interaction Lagrangian is given by where Q L (L L ) is the left-handed quark (lepton) doublet, u R (d R ) is the right-handed up (down) quark singlet and R is the charged lepton singlet. Here λ L(R) αβ are the VLQ coupling strengths to left (right) handed quarks and leptons, with α, β being the generation indices, λ S is the coupling constant of SDQ S DQ with up type quarks and µ S represents the strength of V LQ V LQ S DQ three-point interaction. We also include a coupling λ χ between the VLQ, singlet fermion and right-handed uptype quarks for the ANITA phenomenology. Note that the coupling µ S in the Lagrangian (5) softly breaks lepton number by two units while the baryon number is conserved, so there is no proton decay in this model, while a nonzero Majorana neutrino mass can be induced (see Section VII). The inclusion of these new fields can be realized in gauged B − L extensions of SM or in UV-complete models. For illustration, one such UV-completed scenario is the asymmetric left-right extension of the SM with gauge group SU (3) C × SU (2) L × U (1) R × U (1) B−L in which the electric charge relation is defined as Q = T 3L + T 3R + (B − L)/2 ≡ T 3L + Y where T 3L and T 3R are the third components of isospin generators corresponding to the gauge groups SU (2) L and U (1) R respectively, and B − L is the difference between baryon and lepton numbers. Apart from these usual quarks and leptons, these extra fields like the VLQ, SDQ as well as the singlet fermion are transforming under this asymmetric left-right gauge symmetry as χ(1 C , 1 L , 1/2 R , −1 B−L ). However, in this work we will not focus on any specific model details.
Instead, we work with the effective Lagrangian (5) and discuss its phenomenology in subsequent sections.
After expanding the SU (2) indices in Eq. (5) and performing the Fierz transformation, we obtain the new Wilson coefficients for the process b → cτν l [cf. Eq. (1)] as [79].
where V k3 denotes the CKM matrix element. There are also additional contributions from C It should be noted here that the SU (2) L -singlet VLQ does not provide any additional tensor-type contribution to either b → cτν l or to b → s channels. The tree level Feynman diagram for b → cτν τ (left panel) and b → s (right panel) processes mediated via VLQ are shown in Fig. 1 .
After having the idea about the NP contributions to the Wilson coefficients for both b → s and b → cτν l , we now move forward to constrain these new parameters. For this purpose, we classify the new parameters into the following four scenarios: • Scenario-I (S-I): Includes C LQ V 1 for b → cτν τ and C LQ 9 = −C LQ 10 for b → s (contains only LL couplings).
• Scenario-III (S-III): These new couplings for various scenarios are constrained by performing a global fit (as discussed in Section V), to the relevant experimental observables as listed in the next section.

IV. OBSERVABLES USED FOR NUMERICAL FIT
In our analysis, we consider the following most relevant flavor observables to constrain the new parameters.
In the b → sµµ sector, we include the following observables and their corresponding experimental data.
• R K and R K * : The lepton flavor universality violation ratios R K and R K * are defined as In 2014, the measurement on the LFUV parameter R K , in the low q 2 ∈ [1, 6] GeV 2 region by the LHCb experiment [4]: (where the first uncertainty is statistical and the second one is systematic) has attracted a lot of attention, as it amounted to a deviation of 2.6σ from its SM prediction [97] (see also [98]) The updated LHCb measurement of R K in the q 2 ∈ [1.
also shows a discrepancy at the level of 2.5σ.
In addition to these LHCb results, Belle experiment has recently announced new measurements on R K [17] and R K * [15] in several other bins: R Belle One can notice that the Belle results have comparatively larger uncertainties than the LHCb measurements on R K * ; therefore, we do not include the Belle results for R K ( * ) in our fit for constraining the new parameters.
• B + c → τ + ν τ : This channel has not been measured yet, but indirect constraints on BR(B + c → τ + ν τ ) 30% have been imposed using the lifetime of B c [121]. A stronger constraint of 10% was obtained from LEP data at the Z peak [122]. However, it assumes that the B c hadronization fraction measured in proton-proton collisions can be simply translated to e + e − collisions and it uses this method to predict the number of B c mesons produced at LEP. However, B c production has not been observed at LEP, so there is a large uncertainty in this number, which was not considered in Ref. [122]. Therefore, we will use the more conservative bound of 30% on the B + c → τ + ν τ branching ratio.
To estimate the SM values of the above-discussed observables, we use all the particles masses and lifetime of B q mesons from PDG [125]. The SM results of B s → µ + µ − (τ + τ − ) processes are taken from Ref. [101]. The B → K form factors evaluated in the light cone sum rule (LCSR) approach [126] are considered to estimate B → K processes in the SM. For B (s) → K * (φ) decay modes, we use the form factors from Refs. [127,128]. The decay constant of B c meson is considered as f Bc = 489 MeV [129] to compute branching ratio of B c → τ ν τ .
Since the singlet (3, 1, 2/3) VLQ does not provide additional contributions to b → sν ν type decay modes at tree level due to charge conservation violation, the branching ratio of B → K ( * ) ν ν remains SM-like. Though the charge current D meson decays mediated by c → s ν transitions such as D + s → + ν , D + → K 0 + ν , D 0 → K ( * )− + ν can also play a pivotal role in constraining VLQ couplings, however, they provide very weak bounds on these couplings. Thus, we do not consider these decay modes in our analysis. We further assume that the NP couplings associated with first-generation down-type quark and leptons are negligible. However, the coupling to up-type first generation quark can be non-vanishing via CKM matrix. Since we are mainly interested in the new couplings associated with second and third generation fermions, we do not consider the constraints coming from leptonic/semileptonic K(D) meson decay modes and the We also do not consider the decays like B u → τ ν τ which require new couplings to the first-generation fermions to have the b → uτ ν τ transition, which can be chosen to be small without affecting the b → cτ ν τ transitions, we are interested in.
The VLQ also contributes to loop-level flavor-changing processes, such as the B s − B s mixing, radiative b → sγ and b → sνν decays, as well as Z → l ilj decays. However, the simple VLQ model considered here is, by itself, non-renormalizable, which undermines the predictivity of these looplevel processes, unless some UV-complete framework generating the VLQ mass is explicitly specified; see e.g. Refs. [36-39, 59, 130, 131]. Therefore, in the numerical analysis discussed below, we have considered only those processes which occur at tree-level through the exchange of a VLQ to derive constraints on our simplified model parameter space. However, we will consider a few loop-level processes for tau LFV prediction (see Section VI E) and neutrino mass generation (see Section VII), which should be used with caution due to this caveat.

V. NUMERICAL FITS TO MODEL PARAMETERS
In this section, we consider the NP contributions to both b → s and b → cτν τ processes, and fit the NP parameters by confronting the SM predictions with the observed data. The expression for χ 2 used in our analysis is given by where  [29,35,60]. For concreteness, we fix the VLQ mass at M V LQ = 1.2 TeV in the following analysis, which is consistent with the current LHC constraints [132].
We consider various possible sets of data to fit different scenarios of new Wilson coefficients.
These different cases are further classified as follows. C-III : Includes measurement on B decay modes, which decay either to third generation or second In Fig. 2 , we show the constraints on new leptoquark couplings by using different data sets of above discussed observables for Scenario-I (see Section III), which includes only LL type operators, i.e. C LQ V 1 contribution from b → cτν τ and C LQ 9,10 from b → s + − . Here, the constraint plots for the new couplings for C-Ia (left), C-Ib (middle) and C-II (right) cases are presented in the top panel.
The bottom panel of Fig. 2 represents the constraint plots for C-III in the λ L 33 − λ L 32 (left) and λ L 32 − λ L 22 (right) panels. In each plot of Fig. 2 , different colors represent the 1σ, 2σ, and 3σ contours and the black dot stands for the best-fit value. The corresponding best-fit values obtained for various cases are presented in Table I . In this Table, we have also provided the χ 2 min,VLQ+SM /d.o.f as well as the pull= χ 2 SM − χ 2 best−fit values. For C-Ia case, we have 4 observables with two parameters for fit, thus the number of degrees of freedom (d.o.f.) is 2. Here we find χ 2 min /d.o.f = 2.3/2 = 1.15, which implies the fit is acceptable. The χ 2 min /d.o.f for C-Ib case is found to be 0.575 i.e., the singlet VLQ can explain both b → cτν τ and b → sτ + τ − data simultaneously. We find χ 2 min,VLQ+SM /d.o.f < 1 for both C-II and C-III cases, which implies the VLQ can accommodate b → sτ τ (µµ) anomalies as well as the issues in both b → sτ τ (µµ) and b → cτν τ very well. This analysis implies that the presence of only LL type VLQ couplings can illustrate the B anomalies associated with both b → cτν τ and b → s kind of processes on equal footing. In Scenario-II with the new leptoquark couplings of RR operator type, the constraint on the new couplings associated with right-handed quark and lepton singlets is depicted in Fig. 3 . Since the VLQ has no λ R ij type coupling contribution to b → cτν τ , so we fit the new λ R ij parameters from only b → sµµ data (C-II case of our analysis). In Table I Table I . The χ 2 min /d.o.f is found to be 4.235 (2.117) for C-Ia (C-Ib), which means the fit is rather poor. The χ 2 min /d.o.f for both C-II and C-III cases are slightly greater than 1. Thus, the presence of only pseudo(scalar) type couplings arising due to the exchange of VLQ is not good enough to explain the anomalies in both b → cτν τ and b → s processes.  Table I . We notice that χ 2 min /d.o.f is slightly larger than one. Though the fit is not as good as Scenario-I, but still it is acceptable. The total effective Hamiltonian for b → sl − i l + j processes in the VLQ model can be written as Different colors represent the 1σ, 2σ, and 3σ contours and the black dot stands for the best-fit value.

Scenarios Cases Couplings
Best-fit Values where the vector-axial vector (VA) and scalar-pseudoscalar (SP) parts are given by τ → µγ (bottom panel) processes mediated via the VLQ.
This leads to the following LFV processes: The branching ratio of the LFV B s → l − i l + j decay process in the presence of VLQ is given as where f Bs is the B s decay constant and is the so-called triangle function.
The differential branching ratio of B → Kl − i l + j process is given as [81] dBR where The coefficients J i in Eq. (30) are given in Appendix B.
Including the VLQ contribution, the differential branching ratio of B → K * l − i l + j decay process is given by [133] dBR where the expressions for the angular coefficients I i (q 2 ) are provided in Appendix C . The same expression can be applied to B s → φl i l j process with appropriate changes in the particles masses and the lifetime of B s meson.
Assuming that there is no NP contribution to the first-generation fermions, here we compute the branching fractions for the LFV decay modes of B meson to second and third generation leptons only. One can also notice that the leptoquark couplings required to investigate the above-defined and B s → τ ± µ ∓ [136] for which we find our predictions for the branching ratios are well below the current 90% CL experimental upper limits. Our predictions on branching ratio of B s → τ ± µ ∓ process is = 1.82 × 10 −9 for S − III , which is much lower than the current experimental limit at 90% C.L. [136]: Our estimated branching ratios of the LFV processes B (s) → (K, K * , φ)µ − τ + (µ + τ − ) even for Scenario-III are reasonable and within the reach of future B-physics experiments, such as LHCb upgrade [137] and Belle-II [138].

Decay modes
Predicted values for S-I Predicted values for S-III Experimental Limit (90% CL) In the following two subsections, we study the LFV decay modes of the τ lepton.

D. τ → µφ
The Feynman diagram for τ → µφ LFV decay process is presented in the middle panel of Fig.   6 . The branching ratio of τ → µφ channel is given by [66] BR where f φ is the decay constant of φ meson. Using the values of various masses and lifetime of τ from PDG [125], f φ = (238 ± 3) MeV from Ref. [141] and best-fit values of required new parameters for C-III case of Scenario-I and Scenario-III from Table I , we estimate the branching fraction of τ → µφ as shown in Table II . We find that the τ − → µ − φ branching ratio is substantially enhanced in S-I scenario; it is just below the current experimental upper limit [139] and within the reach of Belle-II experiment.
The right panel of Fig. 6 represents the one loop Feynman diagram for radiative τ → µγ channel.
The effective Hamiltonian for radiative τ − → µ − γ decay mode can be expressed as [142] H Here σ µν is the photon field strength tensor and the Wilson coefficients C L(R) generated due to VLQ exchange are given as where x b = m 2 b /M 2 V LQ , N c = 3 is the color factor and the expression for the loop functions f 1,2,3,4 (x b ) andf 1,2,3,4 (x b ) are given in Appendix D . The branching ratio of this process is [142] where τ τ is the lifetime of τ lepton. In the presence of VLQ, the predicted branching ratio of τ → µγ for C-III of Scenario-I is given in Table II which is roughly an order of magnitude below the current experimental limit [140] . It should be noted that, except C-III of S-I, none of the scenarios can provide the required new parameters to study the τ → µγ process.
Though the muon anomalous magnetic moment gets an additional contribution through one-loop diagram with internal VLQ and down-type quark in the loop, the observed discrepancy cannot be accommodated by using our predicted allowed parameter space.

VII. RADIATIVE NEUTRINO MASS GENERATION
With the particle content of the model discussed in Section III , there are no tree level contributions to light neutrino masses as well as no one-loop level contributions. However, there is a two-loop contribution to light neutrino masses, similar to a variant of the Zee-Babu model [143,144] proposed in Ref. [145] where the lepton doublet is replaced by up-type quark while the singly and doubly charged scalars are replaced by VLQ and SDQ, respectively. The corresponding Feynman diagram for two-loop neutrino mass generation is displayed in Fig. 8 .
The two-loop contribution to light neutrino masses in the flavor basis is given by where I jk is the finite part of the two-loop integral, given by The evaluation of this loop-integral can be done following Ref. [145]. Assuming that the VLQ and SDQ are much heavier than the SM quarks in the loop, the loop function can be reduced to [146]  where the function I(x) has closed-form analytic expression in the following limits: Note that here we have assumed the VLQ loops in Fig. 8  Ref. [24] considered two VLQs (instead of a VLQ and a SDQ as in our case) to cancel the remaining infinities contained in the Passarino-Veltman function by summing over both VLQs. In their case, the neutrino mass can be generated at one-loop level by Higgs-VLQ mixing.
Since for the flavor anomalies, we have only considered couplings to third and second generation fermions, and therefore, do not have full information on all the λ L αj couplings needed to fit the 3-neutrino oscillation data using Eq. (41), we will only derive an order-of-magnitude estimate for the neutrino mass constraint on the model parameters. For illustration, let us take the Scenario-I Case-III which provides the best-fit to both b → cτν τ and b → s anomalies, as discussed in Section V . In this case, the best-fit values of the relevant λ L couplings, (λ L 33 , λ L 23 ) = (0.435, 0.45) can be read off from Table I . Also recall that we have fixed the VLQ mass at M V LQ = 1.2 TeV for the flavor anomalies. We still have three unknowns in Eq. (41) , namely, the trilinear mass term µ S , Yukawa coupling λ S and the SDQ mass M S DQ . As we will see in Section VIII , the λ S coupling cannot be arbitrarily large due to collider constraints from diquark searches. Similarly, the

VIII. SCALAR DIQUARK AT THE LHC
The new TeV-scale colored particles in our model predominantly couple to third and secondgeneration fermions, and offer rich phenomenology at current and future hadron colliders, such as the LHC and its high-luminosity phase, as well as future hadron colliders. The collider phenomenology of color-triplet VLQs coupling to third-generation fermions has been extensively studied; see e.g.
On the other hand, the color-sextet SDQ introduced here to explain the neutrino mass is not constrained by the B-physics sector. In this section, we explore the collider constraints on the SDQ and its future detection prospects. At the LHC, the SDQ can be either singly produced by the annihilation of two up-type quarks, or pair-produced via the gluon-gluon annihilation. The single production has the advantage for relatively heavy SDQ due to the s-channel resonance [158,159], so we focus on this channel only. The single production cross section is governed by the Yukawa coupling λ S in Eq. (5), which in general has a flavor structure. For simplicity, we assume λ S to be proportional to the identity matrix, so that it couples with equal strength to uu, cc and tt. Note that for the neutrino mass generation, it might suffice to have a nonzero coupling to tt and cc only; however, for its production at LHC, an SDQ coupling to uu is desirable due to the large u-quark PDF inside a proton.
Once produced on-shell, the SDQ will decay back to the diquark final states. For M S DQ 2m t , the branching ratios to all quark flavors is roughly the same; for M S DQ close to the 2m t threshold, one has to include the phase space suppression factor of (1 − 4m 2 t /M 2 S DQ ) 3/2 in the S → tt partial decay rate. In our model, for M S DQ > 2M V LQ , the SDQ can also decay into a pair of VLQs; however, we will choose the corresponding coupling strength µ S to be small, so that the diquark decay modes are the dominant ones. A small µ S is also favored by the neutrino mass constraint, if we allow for relatively large λ S values.
In Fig. 10 , we show the SDQ single production cross section (normalized to |λ S | 2 = 1) times branching ratio into dijet (uu + cc) and ditop (tt) final states at √ s = 13 TeV LHC. These numbers were obtained at parton level using MadGraph5 [160] at the leading order. We have used NNPDF3.1 PDF sets [161] with default dynamical renormalization and factorization scales. Also shown is the current 95% CL upper limit on the dijet cross section times branching ratio times acceptance from a recent CMS analysis [151]. Comparing the dijet cross section with the corresponding CMS upper limit, one can derive an upper limit on the coupling λ S as a function of the SDQ mass, as shown by the blue shaded region in Fig. 9. We find that the dijet constraint requires λ S 0.01 − 0.1 for a multi-TeV SDQ.
The same-sign top pair (tt) final state offers a promising test of the SDQ in this model, since the SM background is very small [162][163][164][165][166]. The current experimental limit at 95% CL from a recent ATLAS analysis [167] is shown in Fig. 10 , which only goes till 3 TeV resonance mass. The corresponding constraint on λ S turns out to be weaker than the dijet constraint derived above. However, we expect the ditop sensitivity to improve in the high-luminosity phase and/or in the future hadron colliders.

IX. ANITA ANOMALY
Recently, the ANITA collaboration has reported two anomalous upward-going ultra-high energy cosmic ray (UHECR) air shower events with deposited shower energies of 0.6 ± 0.4 EeV [91] and 0.56 +0.3 −0.2 EeV [92]. This is difficult to explain 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 [91,168,169]. Moreover, as pointed out in Refs. [170,171], 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 [172] and IceCube [173]. Therefore, a NP explanation with an anisotropic astrophysical source with some exotic generation and propagation mechanism of upgoing events is desirable to solve the ANITA anomaly.
As already pointed out in Ref. [45], the observed ANITA events can be explained in our VLQ model framework by including a fermion singlet field χ (1, 1, 0), which couples to the VLQ as given by the last term in Eq. (5). Note that this is one of the handful models that admit LQ coupling to a singlet fermion (aka sterile neutrino) [44,67,174]. This new coupling leads to the production 11. VLQ mediated Feynman diagram for neutrino-quark interaction resulting into production of a long-lived particle χ (left-panel) and the decay of this long lived particle χ into D + s − (right panel).
of χ in the neutrino interactions with up-type quarks (u, c) in Earth matter mediated by the VLQ, which can be resonantly enhanced for TeV-scale VLQ. This occurs for the incoming neutrino energy GeV is the nucleon mass and x is the Bjorken scaling variable, which has an average value of ∼ 10 −3 for these deep inelastic scattering processes. The χ being a SM-singlet can in principle be long-lived and traverses the required chord length before decaying via the same interactions into a D s meson and a charged lepton; see Fig. 11 . We will assume that the charged lepton produced from the χ decay is a tau lepton, which comes from the λ L 23 coupling of the VLQ that is also relevant for the B-anomalies discussed above.
After integrating out the VLQ and performing Fierz transformation, the relevant interaction Lagrangian obtained from Eq. (5) becomes where u = u, c, t is any up-type quark, d = b, s down-type quark and = µ, τ charged leptons. Here i, j are the indices of down-type quarks and charged leptons respectively. Using Eq. (45), the rate of χ → τ D s decay mode is given by The masses of D + s meson and τ lepton are taken from PDG [125] and the decay constant f D + s = 257.86 MeV. We simulate the production of χ by implementing our model file into MadGraph5 [160] at the leading order and using the NNPDF3.1 PDF sets [161]. This is followed by the decay of χ given by Eq. (46) to estimate the event rate at ANITA [170]: FIG. 12. The 2σ and 3σ preferred region in the (M χ , λ χ ) parameter space to explain the ANITA anomalous events in our VLQ model. The green-shaded region is allowed by the D 0 − D 0 mixing constraint. In the vertical grey-shaded region, the χ decay shown in the right panel of Fig. 11 is not kinematically allowed.
where T = 53 days for the total effective exposure time, Φ ν = 2 × 10 −20 (GeV · cm 2 · s · sr) −1 for the cosmic neutrino flux, 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 χ and the cross section for the χ production; see Ref. [170] for the explicit expression for an analogous bino production in supersymmetry. From Eq. (47), we know that the overall event number N is a function of m χ and λ χ for a given VLQ mass. 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. This is shown in Fig. 12 , where the dark and light blue-shaded regions can explain the ANITA events at 2σ and 3σ CL respectively for M V LQ = 1.2 TeV. The vertical grey-shaded region is kinematically forbidden for the χ decay shown in the right panel of Fig. 11. Note that our results for the ANITApreferred region in Fig. 12 are slightly different from those given in Ref. [45]. We also include the D 0 − D 0 mixing constraint, as discussed below.
In presence of the singlet χ, there will be a new contribution to the D 0 − D 0 mass difference from the box diagrams with the VLQ and χ flowing in the loop, as shown in Fig. 13. The effective Hamiltonian for D 0 − D 0 mixing in the presence of VLQ is where the NP Wilson coefficient is given as with x χ = M 2 χ /M 2 V LQ and the loop function The SM contribution to the mass difference is negligible and the corresponding measured value is given by [125] ∆M D = 0.0095 +0.0041 −0.0044 ps −1 .
The green-shaded region in Fig. 12 shows the allowed parameter space from this constraint.
The presence of χ also leads to an additional contribution to the B + c → τ + ν process, via a diagram similar to the right-panel of Fig. 11 (with s replaced by b), since the χ practically behaves like a neutrino for the energies involved in the B-decays. The corresponding branching ratio is given by For M Bc = 6.25 GeV and m τ = 1.77 (m µ = 0.106) GeV [125], the maximum mass value of χ for which the B + c → + χ process can occur kinematically is M χ = 4.47 (6.144) GeV. However, from Fig. 12 , we see that the overlap between the ANITA and ∆m D preferred regions occurs only between M χ = [6,30] GeV. Hence, the B + c → τ + χ decay is not relevant here.

X. CONCLUSION
The recently observed various flavor anomalies in the CC and NC mediated semileptonic B meson decays may be considered as one of the most compelling hints of NP at the TeV scale. To explain these intriguing set of discrepancies in a coherent manner using a single framework is a challenging task,  Table II) can be used to test this scenario in the future B-physics experiments, such as LHCb upgrade and Belle-II.
In addition, augmenting the VLQ model with a color-sextet SDQ can explain the neutrino mass at two-loop level (see Section VII). We discussed the LHC constraints on the SDQ mass and Yukawa coupling with up-type quarks, and identified the same-sign top-pair production as an excellent probe of this scenario for a multi-TeV SDQ in the future high-energy collider experiments, such as highluminosity LHC (see Section VIII). Further, adding a GeV-scale SM-singlet fermion to the VLQ model can also explain the ANITA anomalous upgoing events. It was shown to be consistent with the D 0 − D 0 mixing constraint (see Fig. 12). In summary, we have proposed a unified explanation of the flavor anomalies, radiative neutrino mass and ANITA events. Different aspects of the model can be tested in future collider and B-physics experiments.
Appendix D: τ → µγ The loop functions required to compute τ → µγ decay mode in the presence of VLQ are [142] f (D1)