Scalar leptoquark effects on $B \to \mu \bar\nu$ decay

Purely leptonic $B$ meson decays provide unique probes for physics Beyond the Standard Model. We study the impact of a scalar leptoquark, $S_1$, on $B \to \mu \bar\nu$ decay. We find that, for $m_{S_1}\sim 1$ TeV, the $S_1$ leptoquark can modify the $B \to \mu \bar\nu$ rate significantly. Such a leptoquark can in prinicple also alter the $B \to \tau \bar\nu$ rate. However, current searches from LHC and low energy physics provide some constraints on the parameter space.


I. INTRODUCTION
Purely leptonic B − meson decays provide clean probes to new physics beyond the Standard Model (SM). The prime example is B → τν [1], where the experimental observation [2] has provided one of the strongest constraints on parameters of a charged Higgs boson, especially in the so-called two Higgs doublet model (2HDM) type II [3] that automatically arises with supersymmetry. The B → µν decay is further helicity suppressed, and has not been observed so far. It was pointed out recently that, while the ratio B(B → µν)/B(B → τν) is predicted to be the same for both SM and the popular 2HDM type II [1], the value could deviate [4,5] from the SM expectation in the more general 2HDM (g2HDM) that allows extra Yukawa couplings.
The Belle experiment has recently measured B(B → µν) = (6.46 ± 2.22 ± 1.60) × 10 −7 [6] using the full dataset of 711 fb −1 , finding 2.4σ significance. This should be compared with the Standard Model (SM) expectation around 3.92 × 10 −7 [5]. In anticipation of Belle II data, the measurement was further updated to B(B → µν) = (5.3 ± 2.0 ± 0.9) × 10 −7 [7], with the aim of improving the systematics error. Despite a slight drop in central value, the significance moved up from 2.4σ to 2.8σ [7]. Assuming the SM rate, the Belle II experiment should be able to observe B → µν decay with ∼ 5 ab −1 [8] in its early running. If the rate is actually larger than SM, observation would come sooner.
It is well known that leptoquarks (LQ) can also affect semileptonic and purely leptonic meson decays [9], which has been of some interest lately. The impact of a scalar LQ (SLQ) on B → τν decays and the associated constraints have been discussed in Refs. [10,11]. In this paper we study the effect of the SLQ, S 1 , on B → µν decay and discuss the possible constraints. For sake of comparison, the impact of S 1 on B → τν decay is also considered. The experimental searches for SLQs at the LHC (see e.g. Refs. [12,13]) are based on the minimal Buchmüller-Rückl-Wyler model [14]. In our study, we allow S 1 to couple to different generations of quarks and leptons [9]. This work is therefore complementary to the previous study of Ref. [5] on H + effects in g2HDM, where the effective 4-Fermi operator approach was adopted to match the experimental presentation [7]. Our starting point would be from this New Physics (NP) Wilson coefficient language.
We find that deviations of the B → µν, τν decay rates from their SM expectations are constrained in particular by direct searches at the LHC, as well as several low energy measurements. The direct search constraints arise primarily from S 1 pair production, followed by S 1 → qµ ± decay [15,16], which cuts into the parameter space allowed by B(B → µν). Although B(B → τν) is less constrained, the case of complex S 1 Yukawa couplings is constrained by the electric dipole moment (EDM) of the neutron [11].
The paper is organized as follows. We give the formalism in Sec. II, then present our results on Wilson coefficients in Sec. III. We discuss possible constraints on SLQ Yukawa couplings in Sec. IV, and summarize with some discussions in Sec. V.

II. FORMALISM
Let us consider the SLQ, S 1 , which has quantum numbers (3, 1, 1/3) under the SM gauge group. The relevant Lagrangian of S 1 interacting with SM quarks and leptons can be written as 1 In the above, Q L (L L ) denotes the left-handed quark (lepton) doublet under SU (2) L , while u R ( R ) denotes the right-handed up-type quark (charged lepton) singlet, and i, j are generation indices. These fermion states are written in the down-type quark mass basis. After rotating to mass eigenbasis via the transformations where V and U are the CKM and PMNS matrices, respectively, the Lagrangian of Eq. (1) can be expanded in the mass eigenbasis as where L, R = (1 ∓ γ 5 )/2.
1 Concerning the stability of the proton, we turn off the coupling of SLQ to di-quarks by imposing appropriate symmetry. We also do not consider right-handed neutral leptons in this paper.

arXiv:1909.00403v2 [hep-ph] 25 Nov 2019
For purely leptonic B − decays, the effective Hamiltonian in neutrino flavor-basis is given by where The SM contributes only to the V −A interaction via Wboson exchange, where the Wilson coefficients are written as The leptoquark contributions, at the m S1 scale, are given by where we have approximated the factor V ud k y L * ki , with k summed over, as y L * 1i , i.e. we assume the other y L * 2i , y L * 3i factors do not overpower the CKM suppression of V us and V ub , respectively. The branching ratio for B → ν in SM is well known Adding the SLQ contributions, the branching ratio can be expressed as where the undetected anti-neutrino flavor in the final state is summed over. Eq. (7) is essentially the same form used by Belle [7], except that we allow for the phase(s) where m b (m b ) and m b (µ 0 ) are the MS running masses evaluated at m b and m S1 , respectively. Defining the ratio in SM one finds which is relatively precise (subject to mild QED corrections) since it involves only lepton masses. The prediction for 2HDM type II is the same [1,4,5]. The B → ν ( = µ, τ ) decay rates can deviate from SM predictions in the presence of the SLQ S 1 . We shall ignore y L(R) i1 couplings and set them to zero for simplicity [17]. This is similar to the treatment of Ref. [18] where both y L(R) i1 and y L(R) 1i were assumed to be zero. Here we set y L(R) i1 = 0, and note that these parameters enter in the electron EDM at one loop, hence receive severe constraints from the recent ACME result [19]. But y L(R) 1i are less constrained. In the following, we will first present the Belle measurements of B → µν and B → τν decays in terms of C S(V )L Wilson coefficients, then discuss in some detail the constraints on relevant y Yukawa couplings, including other possible processes.

III. CONSTRAINING WILSON COEFFICIENTS
From Eq. (7) one can see that both the Wilson coefficients C SL and C V L can alter B → ν decay rates. However, the former is more efficient, as C SL receives the m 2 B /m m b enhancement factor compared with the C V L term. This is especially true for B → µν decay, where m 2 B /m µ m b ∼ 60. For B → τν decay, the C τ SL mechanism does not get large enhancement because 1/m τ is smaller than 1/m µ , but there is still some advantage over C τ V L by m 2 B /m τ m b ∼ 4. Hence, from here on we will primarily focus on C SL , and touch only briefly on the C V L mechanism.
A. B → µν decay We first focus on B → µν decay. Depending on the type of anti-neutrino flavor in the final state, there are different C SL Wilson coefficients that can modify the B → µν rate. For muon anti-neutrino in final state, C µµ SL interferes with the SM contribution, i.e. δ µµ = 1 in Eq. (7), while C µe SL and C µτ SL effects add in quadrature for electron and tau anti-neutrino emission. But as already mentioned, C µe SL receives stringent constraints, hence we ignore this Wilson coefficient for simplicity.
We plot B(B → µν) in the |C µµ SL | vs φ S µµ plane in the left panel of Fig. 1, setting all other Wilson coefficients to zero, while in the right panel we give the dependence of B(B → µν) on |C µτ SL |, where aν τ is emitted. In generating Fig. 1, we used B(B → µν µ )| SM 3.92×10 −7 , which arises from utilizing f B = 190 MeV from FLAG [20], and the exclusive value |V ub | excl. = 3.70 × 10 −3 [2].
Taking a closer look, Fig. 1  central value from Belle is illustrated by the black dotted line, with green dark (light) shaded regions illustrating the 1σ (2σ) range [7]. Forν τ emission, which is not distinguished by experiment, the dependence of B(B → µν) on |C µτ SL | is given by blue dashed line in Fig. 1(right), which can be only constructive as it adds in quadrature. At this level of discussion, as one is using operator language with Wilson coefficients that follow the presentation by Belle, the plots bear similarity to those in Ref. [5] that treat H + effects in g2HDM.
Before turning to B → τν, let us compare the C SL and C V L mechanisms. We see from Fig. 1 that, to account for the Belle central value for B(B → µν) [7], one needs while one finds from Eq. (7) that which are much larger in value. Similarly, |C µτ SL | = 0.009 is sufficient to explain the Belle central value, but |C µτ V L | would need to be 0.568 to produce the same effect. The required value of |C µµ SL | for φ S µµ = 0 (π) is about a factor of m 2 B /m µ m b ∼ 60 smaller than that of |C µµ V L | for φ V µµ = π (0), similarly for |C µτ SL | versus |C µτ V L |. This illustrates that C SL provides a much more efficient mechanism to modify B(B → µν) compared with C V L . Fig. 2(left) we give the contours of B(B → τν) in the |C τ τ SL | vs φ S τ τ plane, whereν τ is emitted, as well as dependence on |C τ µ SL | in Fig. 2(right), forν µ emission. We have used the Belle average value of B(B → τν) = (9.1 ± 2.2) × 10 −5 from PDG, while the SM expectation is 8.73 × 10 −5 [5]. These are given by black dotted and red solid lines, respectively. The 1σ and 2σ allowed ranges are illustrated by the dark and light cyan shaded regions, while the blue dashed line depicts the dependence of B(B → τν) on |C τ µ SL | for Fig. 2(right). Unlike In the left panel, the black and green colored contours are for φµµ = 0 and π respectively. In the right panel the dotted and dashed contours are plotted in black, however, it should be understood that these contours do not depend on the overall phase, and hence black colors does not correspond to any phase labeling. The red solid lines in the respective panels illustrate SM expected value when the magnitude of the any one of the couplings is zero. The purple shaded regions are excluded by ATLAS SLQ pair production [16]. See text for further discussion. The contours types and color schemes are basically same as in Fig. 3 except for one additional blue dotted contour, which represents the Belle central value for purely imaginary phase (φττ = π/2). The orange shaded region and the thick orange line illustrate the current and future reach of constraint from neutron EDM for purely imaginary φττ , as it probes CP violation. The light gray shaded regions in both panels are excluded by ATLAS pp → τ τ search [48]. See text for further discussion.
B → µν decay, the |C τ τ,τ µ SL | mechanisms do not have large enhancement factors over |C τ τ,τ µ V L |, as can be seen from Eq. (7). As a result, larger Wilson coefficient values are needed compared with the B → µν case of Fig. 1.

IV. CONSTRAINTS ON LEPTOQUARK YUKAWA COUPLINGS
Our presentation so far is not so different from the H + study [5] in g2HDM with extra Yukawa couplings, as we follow the effective Hamiltonian approach used by Belle [7], except allowing the Wilson coefficients to be complex. The underlying physics is, however, quite different. To find the allowed parameter space for S 1 lepto-quark, we turn to study B(B → µν) and B(B → τν) in terms of the relevant S 1 Yukawa couplings y R 12 , y R 13 , y L 32 and y L 33 .
A. Constraints from B → µν and B → τν For simplicity, in the left (right) panel of Fig. 3, we assume y R 12 and y L 32 , i.e. equivalently y R uµ and y L bνµ (y R 12 and y L 33 , i.e. equivalently y L bντ ) are the only nonvanishing couplings. Similarly, for Fig. 4(left), we set all couplings to zero except y R 13 and y L 33 , i.e. equivalently y R uτ and y L bντ , whereas we take y R 13 and y L 32 , i.e. equivalently y L bνµ as the only nonzero couplings for Fig. 4(right). As we focus in this subsection on the |C SL | mechanism, we drop the S superscript from φ from here on.
Let us understand Figs. 3 and 4 better. We have chosen m S1 = 1.2 TeV for illustration. In both figures, the central value of Belle measurement, +2σ and −2σ ranges, are denoted by dotted, dashed and dotdashed contours respectively. The solid curves illustrate the SM expectation, but correspond to rather sizable S 1 Yukawa couplings, or when one of the couplings vanishes, which would be elucidated later. The left (right) panel of Fig. 3 corresponds essentially to Fig. 1 left (right), where we set all C SL = 0 except C µµ SL (C µτ SL ). However, we note that there is a negligibly small contribution (proportional to V ub ) in Fig. 3(left) from C µµ V L if y L 32 is non-zero, which we neglect. A similar procedure is followed for Fig. 4.
For the detailed respective contours, in Fig. 3(left) we have two sets of contours for Belle central values of B(B → µν): one for φ µµ = 0, the other for φ µµ = π, which are denoted as black and green colors respectively. The former requires larger values of |y R 12 | and |y L 32 | couplings than the latter. This can be understood from Eqs. (7) and (8), where the S 1 contribution for φ µµ = π interferes constructively with SM in Eq. (8), while destructively for φ µµ = 0. Similarly, there are two dashed contours for Belle +2σ range, the black one is associated with φ µµ = 0, while the green one is for φ µµ = π. For Belle −2σ range, however, destructive interference with the SM contribution would be needed, which is possible only for φ µµ = 0, but not for φ µµ = π. In this case one has two black dotdashed contours, illustrating the quadratic solutions (see Eq. (7)) for |y R 12 | and |y L 32 | couplings to give the Belle −2σ range. A similar explanation goes for B(B → τν) in Fig. 4(left), where there are two dotted and dashed contours each for the Belle central value and +2σ range, corresponding to φ τ τ = 0 and π, shown by black and green respectively. But for the Belle −2σ range there are again two black dotdashed contours for φ τ τ = 0, due to quadratic solutions for |y R 13 | and |y L 33 | couplings.
The black and red solid contours in Figs. 3(left) and 4(left) correspond to the SM expectations for B(B → µν) and B(B → τν), respectively, which require some explanation. One would obviously recover the SM value when Yukawa couplings vanish, as illustrated by the red solid straight lines in each of these figures. But for φ µµ , φ τ τ = 0 where the SLQ interferes destructively, large Yukawa couplings can overpower the SM effect to reach the SM value, which are the black solid lines displayed in Figs. 3(left) and 4(left).
In contrast, the right panels of Figs. 3 and 4 are straightforward: the SLQ contributions can only add in quadrature, therefore one has black dotted and dashed lines, respectively for Belle central values and +2σ ranges.
We have also plotted the Belle central value for φ τ τ = π/2 in Fig. 4(left), to compare with the constraint from neutron EDM discussed later. In this case, the SLQ contribution is purely imaginary, and adds in quadrature to the SM effect. Note that the Wilson coefficients are generally complex, as illustrated in Figs. 1 and 2. So, the actual interpretation in terms of SLQ Yukawa couplings are more complex than what is presented here.

B. b → cµν and b → cτ ν observables
The presence of y L 32 and y L 33 can induce b → cµν and b → cτ ν transitions via vector Wilson coefficients respectively. Most notably, such Wilson coefficients contribute e.g., to B(B → D ( * ) µν)/B(B → D ( * ) eν) and, B(B → D ( * ) τ ν)/B(B → D ( * ) ν) (with = e, µ) ratios respectively. The latter ratios are refereed to as R(D ( * ) ) ratios, where some tensions are found between the experimental measurements and the SM predictions [21]. Let us first focus on the parameter ranges for y L 33 in the context of R(D ( * ) ) anomalies. The latest SM predictions [21] are: R SM (D) = 0.299 ± 0.003 and R SM (D * ) = 0.258 ± 0.005, whereas the world averages of the experimental measurements from HFLAV [21] are R(D) = 0.340 ± 0.027 ± 0.013 and R(D * ) = 0.295 ± 0.011 ± 0.008. The combination of R(D)-R(D * ) world averages deviate 3.1σ from SM predictions; which are together as the so called "R(D ( * ) ) anomalies". To find the constrain on |y L 33 |, we utilize the 1D global fit value of Ref. [22] which includes BaBar, Belle and, LHCb data on B → Dτ ν, B → D * τ ν and B c → J/Ψτ ν decay observables. The Wilson coefficient C cb;τ τ V L is found to be [0.04, 0.11] at 2σ for the matching scale µ = 1 TeV [22], which indicates some tension with SM. Taking this range as ballpark value for S 1 LQ with m S1 = 1.2 TeV, we find 1.9 |y L 33 | 3.1 at 2σ. While finding this constraint, we assumed b → ceν and b → cµν transitions to remain SM like.
This illustrates, current 1D global fit favors rather large |y L 33 |, which is beyond the plotted ranges in Figs. 3 (right) and 4 (left). However, such large |y L 33 | would also induce C τ τ SL (C µτ SL ) if y R 13 (or y R 12 ) are non-vanishing and, could be sensitive to direct search limits [23]. We defer discussion regarding this constraint for Sec. 4.6. In such scenarios, to open up the parameter space for smaller |y L 33 |, one may require other non-zero couplings such as y R 23 and y L 23 , or possibly more leptoquark couplings as discussed in Refs. [11,24]. A detailed study with all three y R 12 (or y R 13 ), y R 23 , y L 23 non-vanishing couplings would be interesting in the light of new flavor and direct search results. We leave out such analysis for future. If future measurements of LHCb and Belle-II support the current tension in R(D ( * ) ) and the anomalies become more prominent, the global fit value of Wilson coefficient C cb;τ τ V L would deviate more from its SM prediction. Such large values of C cb;τ τ V L could be sensitive [23] to direct searches at the HL-LHC (High Luminosity LHC). On the other hand, if R(D ( * ) ) becomes SM like in future, B → µν and B → τ ν decays would provide sensitive probe for |y L 33 | when y R 12 and y R 13 respectively are non-vanishing. The ratios R µ/e D = B(B → Dµν)/B(B → Deν) and R e/µ D * = B(B → D * eν)/B(B → D * µν) are measured by Belle and found to be 0.995±0.022±0.039 [25] and 1.04± 0.05 ± 0.01 [26] respectively. Moreover, Ref. [27] found lepton flavor universality violation between b → ceν and b → cµν transitions could still be ∼ 2% 2 . Assuming b → ceν transition to be SM like, and allowing ∼ 2% deviation in the b → cµν transitions, we find |y L 32 | 0.95 is excluded for m S1 = 1.2, which is larger than the ranges plotted in Figs. 3. This not surprising since the contribution from y L 32 is suppressed compared to SM by factor √ 2/(8G F m 2 S1 ) (see Eq. (13)), which is about ∼ 0.01 for m S1 = 1.2 and induces per-mille to sub-percent level effects in b → cµν transitions for the plotted range of |y L 32 |. This is also well within the 2σ allowed ranges of Belle R µ/e D and R e/µ D * measurements [25,26].

C. Muon anomalous magnetic moment
The muon anomalous magnetic moment a µ is defined via the coupling (e/4m µ ) a µμ σ αβ µ F αβ . The S 1 leptoquark can generate ∆a µ radiatively [9,29,30] via where Q q c = −2/3, Q S = 1/3 for the leptoquark S 1 , x t = m 2 t /m 2 S1 , and the functions f q c (x), f S (x), g q c (x) and g S (x) can be found in Ref. [29]. This will constrain y L 32 , regardless of the value of y R 32 . The current experimental world average [2] and the SM predicted [31] values show some deviation, corresponding to a long standing 3.7σ discrepancy [31] that could be due to New Physics. However, for the plotted ranges in Figs. 3 and 4, the contributions from y L 32 turn out to be negligible. 2 Ref. [27] utilized 2014 PDG [28] fit for the combined B ± /B 0 results to estimate the lepton flavor universality violation between b → ceν and b → cµν transitions. We used this value as yardstick to determine the constraint on |y L 32 |.

D. τ → µγ decay
The branching ratio for τ → µγ is given by [9,30] where Γ τ is the τ width, and The current limits are B(τ → µγ) < 4.5 × 10 −8 from Belle [32] and 4.4 × 10 −8 from BABAR [33], both at 90% C.L. Belle II may improve the limit by a factor of 100 [8], which would provide some constraint on the parameter space via A R τ µ , where the product y L * 32 y L 33 is proportional to m τ . However, we find that the present constraints from Belle and Babar are again weaker than the range plotted in Figs. 3 and 4.

E. EDM measurements
The ACME experiment has put stringent [19] constraints on electron EDM, d e , which prompted us to set y L(R) i1 to zero. The neutron EDM, d n , imparts some constraint on the parameter space for B(B → ν) decays.
The effective Hamiltonian can be written as [11,34] where the dimension-6 O T and dimension-5 O γ, g operators can be found in Ref. [11]. At one loop, the leptoquark S 1 will contribute to the neutron EDM with τ and µ running inside the loop. The contribution arising from the τ loop to the Wilson coefficients at the high scale can be written as C γ = − m τ |V ub | |y L * 33 | |y R 13 |e −iφττ 96π 2 m 2 S1 4 + 3 log(µ 2 /m 2 S1 ) , whereas the muon loop is suppressed by m µ . Note that V ub enters here through the first term of Eq. (2). Neutron EDM depends on finite CPV phase. The contribution to neutron EDM can be expressed as [35] d n /e = −(0.44 ± 0.06) Im C γ − (1.10 ± 0.56) Im C g , where C γ and C g are evaluated at 1 GeV, while C T does not contribute. We follow Ref. [11] for the RGE evolution of the Wilson coefficients from the m S1 scale.
The current 95% C.L. upper limit of neutron EDM, viz. |d n /e| < 3.6 × 10 −26 cm [36] (see also Ref. [37]) sets strong constraint on the parameter space in Fig. 4(left) for sin φ τ τ = 0. As illustration, we use Eq. (24) and find the orange shaded excluded region for φ τ τ = π/2, i.e. purely imaginary Wilson coefficient in Fig. 2(left). Future measurements are expected to push the upper limit to |d n /e| < 10 −28 cm [38], which is displayed as the thick orange line. This illustrates that future d n measurements can exclude the whole parameter space of S 1 that supports the current Belle central value for B(B → τν), if φ τ τ = π/2, i.e. the phase of C τ τ SL is near maximal. Note that the constraint vanishes for φ τ τ = 0 or π, hence it should not be confused that the φ τ τ = 0 or π contours for B(B → τν) are excluded. The parameter space of B → µν decay is less constrained due to m µ suppression.
We have mainly focused on the neutron EDM. Impact of other EDMs such as mercury, proton, deuteron and can be found in more detail in Refs. [11,34] for B → τν.

F. Direct searches
The scalar leptoquark S 1 can be singly or pair produced at the LHC in pp collisions and subsequently decay into u i − j and d i ν j final states (conjugate processes are always implied), depending on the values of y R ij and y L ij . Several searches by ATLAS (e.g. Refs. [12,16]) and CMS (e.g. Refs. [15,39]) set strong limits on leptoquark mass and branching ratios. At the current collision energy of √ s = 13 TeV, S 1 pair production via gluon fusion [16] is the dominant mechanism, while qg initiated single leptoquark production is subdominant. For the range of y R/L ij couplings in Figs. 3 and 4, we find that the most relevant constraints arise from Refs. [15,16].
The strongest constraint comes from the ATLAS search [16] for SLQs at √ s = 13 TeV with 36.1 fb −1 , with final states containing two or more jets, one muon or electron and missing energy, or two or more jets with two electrons or muons. The search gives 95% C.L. upper limits for branching ratios of leptoquark decaying into an electron and a quark, or a muon and a quark, for different values of leptoquark masses. As the final state jets are not tagged, the constraint on S 1 parameters will be modulated by B(S 1 → uµ) if y R 12 is nonzero. Using the 95% C.L. upper limit as B(S 1 → qµ ± ) [16] for leptoquark mass of 1.2 TeV, we find the purple excluded regions as displayed in Fig. 3.
The CMS search [15] sets limit on the mass vs branching ratio to tµ (and tτ ) to third generation leptoquarks. Although only bν type of SLQ couplings enter B → ν, there are corresponding t couplings. With our assumptions discussed above, for B → µν decay of left (right) panel of Fig. 3, S 1 decays to uµ, bν µ (bν τ ) and tµ (tτ ), while for B → τν decay of left (right) panel of Fig. 4, S 1 decays to uτ , bν τ (bν µ ) and tτ (tµ). With the assumed couplings, one has e.g. B(S 1 → tµ) ≈ 0.3 (0.5) for y R 12 0.2 (0.05) and y L 32 0.2 (0.2). For an SLQ with mass of 1.2 TeV, these branching ratios are below the observed [15] 95% C.L. upper limits at 0.56 and 0.64 for tµ and tτ decays, respectively. Similarly, constraints from CMS upper limits are also weaker than the parameter ranges given in Fig. 4. The sensitivity of the HL-LHC in probing S 1 → tµ is discussed in Ref. [40], where y L 32 0.4 is expected to be excluded at 95% C.L. Note that we have neglected CKM suppressed decays such as S 1 → cµ, while finding the limit on |y R 12 |. However, the limit will weaken if S 1 decays to other final states such as tµ or tτ , but the HL-LHC may be sensitive [40] (see also Refs. [41,42]) to B(S 1 → t ) ∼ 0.5 (with = µ, τ ). The impact of direct searches with full HL-LHC dataset on the parameter space of B → ν is worthy of further scrutiny, and will be studied elsewhere.
We remark that the couplings y R 12 and, y R 13 receive constraints from heavy resonance searches in the dilepton final states such as pp → µµ (dimuon), pp → τ τ (ditau) and, pp → µτ via t channel S 1 exchange as discussed in Ref. [43,44]. Utilizing the search for heavy resonances in the dimuon final states of Ref. [45], Ref. [47] find |y R 12 | 0.5 are excluded for m S1 ∼ 1 TeV, whereas future dimuon searches with full HL-LHC dataset can exclude |y R 12 | 0.1 [46]. The coupling y L 32 also receives constraint from such search, however, the limit is rather weak due to suppression from associated CKM elements. The search for heavy resonance in the ditau final state can constrain y R 13 . We utilize ATLAS √ s = 13 TeV 36.1 fb −1 ditau search result [48] to constrain |y R 13 |. We closely follow the procedure outlined in Refs. [41,44] in our analysis for extraction of the upper limit.
The ATLAS search [48] is divided the into two categories, based on τ had τ had and τ had τ lep final states. In the τ had τ had category events with two hadronically decaying τ s are selected, however, in the latter case events are selected such that it contains one leptonically and one hadronically decaying tauons. The search provides distributions for m top T [48] (see also Refs. [41,44] for definition) in different bins in both the final state categories, which can be found in HEPData repository [49]. In pp collision non-zero, y R 13 will induce uū → τ + τ − process via t-channel S 1 exchange and contribute abundantly in both τ had τ had and τ had τ lep categories. As the search results in Ref. [48] does not veto additional activity, we also include contributions from ug → S 1 τ + → τ + τ − u and gg → S 1 S * 1 → τ − uτ +ū . To determine the constraint on y R 13 , we generated these processes at Leading Order (LO) in pp collision at √ s = 13 TeV utilizing Monte Carlo event generator MadGraph5 aMC@NLO [50] with the parton distribution function (PDF) set NN23LO1 [51]. The event samples are then interfaced with PYTHIA 6.4 [52] for showering and hadronization, and finally fed into fast detector simulator Delphes 3.4.0 [53] to incorporate detector effects. Here we adopted the default ATLAS based detector card available in Delphes. We adopted MLM matching scheme [54] for matrix element and parton shower merging, and utilized the Feyn-Rules [55] model available in Ref. [56]. We have defined test statistic as [41]: where N i T = N i expec. + N i S1 , ∆N i = N i obs with N i expec. , N i S1 and N i obs are the expected number of events, events from S 1 LQ and observed number of events in the i-th bin of the m top T distribution 3 . The ∆χ 2 = χ 2 − χ 2 min = 2 corresponds to 2σ range, with χ 2 min is the minimum value of χ 2 for m S1 = 1.2 TeV for some value y R 13 = y R 13 min ≥ 0. We have found |y R 13 | 0.25 is excluded at 2σ. The 2σ excluded regions are shown by light gray shaded region in Figs. 4. We remark that |y L 33 | can also be constrained by search in Ref. [48], however the limit is rather weak due to suppression from CKM element V cb .
The search for high-mass resonances in pp → τ ν [58,59] can constrain |y L 33 | due to the presence of cτ S 1 and bνS 1 couplings, as discussed in Ref. [23]. The limit is not stringent due to weak c and b PDF, however, the excluded regions lie beyond the plotted ranges in Fig. 3 (right) and Fig. 4 (left). Such a search can also constrain coupling products |y R * 13 y L 33 |. Recasting the 2σ upper bound of the Wilson coefficient |C τ τ SL | (at m b scale) of Ref. [23], we found |y R * 13 y L 33 | 0.28 is excluded, whereas, the full Run-2 dataset can exclude |y R * 13 y L 33 | 0.17. As the neutrino flavor is not measured similar upper bound can also be set on |y R * 13 y L 32 |. In addition, search for heavy resonance decaying into µν final state (denoted as pp → µν search) [60] can constrain |y R 32 | as well as coupling products |y R * 12 y L 32 | and |y R * 12 y L 33 |. However, we find these constraints to be weaker and excluded regions lie just outside the plotted ranges of Fig. 3. As discussed earlier, current 1D global fit of b → cτ ν observables favors rather large |y L 33 |, which means parameter space of B → τν in Fig. 4 (left) (B → µν in Fig. 3 (right)) where |y R * 13 y L 33 | (|y R * 12 y L 33 |) product is somewhat large could be excluded by pp → τ ν (pp → µν) search with full Run-2 or early Run-3 data. Furthermore, search for heavy resonances decaying into the τ µ final state [61] can potentially constrain coupling products such as |y R * 12 y L 33 | and |y R * 13 y L 32 | (see Eq. (1)), however, we find the limits to be not yet relevant for the coupling ranges considered in this paper (see also Ref. [62] for similar discussion).

V. DISCUSSION AND SUMMARY
We offer a few brief remarks in passing. The decays of Z and W bosons can constrain the parameter space for y L ij . For example, Z → τ τ and W → τ ν exclude y L 33 1 at 2σ for m S1 ≈ 1 TeV [11]. Such constraints are, however, in general weaker than the ranges plotted in Figs. 3 and 4. Note that the effect of S 1 on B → τν has been previously discussed [11], and our discussion is only for comparison with B → µν. In principle, other SLQs such as S 3 , R 2 ,R 2 (see Ref. [9] for definition), as well as vector leptoquarks U 1 , U 3 , V 2 ,Ṽ 2 can all potentially affect B → ν decays. We leave a detailed study of these for the future.
In some sense, it is remarkable that somewhat large leptoquark Yukawa couplings as displayed in Figs. 3 and 4 remain unexplored. We have seen that, when light jets are involved, ATLAS data [16] provide strong constraints. Moreover, direct searches such as pp → µµ and pp → τ τ provides meaningful constraints on the available parameter space for both B → µν or B → τν decays. But this also illustrates the relative arbitrariness of the S 1 scalar leptoquark, where putting the mass above TeV scale on one hand escapes LHC detection, but on the other hand demands the rather large Yukawa couplings (the bottom Yukawa coupling is ∼ 0.02 in SM) to have an effect on purely leptonic B − decays. In contrast, the H + effect of the general 2HDM that allows for extra Yukawa couplings is much more nuanced [5]. The charged Higgs could be sub-TeV, with rather weak extra Yukawa couplings, but could still enhance B → µν (less so for B → τν) within the Belle allowed range.
In summary, we have explored the constraints placed by current Belle results on the Wilson coefficients that can affect B → µν, τν decays, and then interpreted in terms of the Yukawa couplings of the S 1 scalar leptoquark. With m S1 set at 1.2 TeV, rather sizable Yukawa couplings are needed for enhancing the purely leptonic B − decays to the 2σ upper reach of Belle measurements. As one awaits eventual Belle II observation of B → µν and improved measurements of B → τν, we find that neutron EDM can probe the CP violating phases of S 1 Yukawa coupling, while a large part of the rather large leptoquark Yukawa coupling range remains to be explored at hadron colliders.