A neutrino window to scalar leptoquarks: from low energy to colliders

Leptoquarks are theorized particles of either scalar or vector nature that couple simultaneously to quarks and leptons. Motivated by recent measurements of coherent elastic neutrino-nucleus scattering, we consider the impact of scalar leptoquarks coupling to neutrinos on a few complementary processes, from low energy to colliders. In particular, we set competitive constraints on the typical mass and coupling of scalar leptoquarks by analyzing recent COHERENT data. We compare these constraints with bounds from atomic parity violation experiments, deep inelastic neutrino-nucleon scattering and LHC data. Our results highlight a strong complementarity between different facilities and demonstrate the compelling power of coherent elastic neutrino-nucleus scattering experiments to probe leptoquark masses in the MeV-GeV range. Finally, we also present prospects for improving current bounds with future upgrades of the COHERENT detectors and the planned European Spallation Source.


I. INTRODUCTION
Leptoquarks (LQs) are hypothetical particles that carry both lepton and baryon numbers, and can arise in many extensions of the Standard Model (SM) unifying matter [1][2][3][4] with the unique property of connecting leptons and quarks.This peculiar property could induce rapid proton decay.However, as they arise in many Grand Unified Theories (GUT), their mass is expected to be close to the GUT scale preventing this process from happening [5].On the other hand, there are models where the operator that gives rise to proton decay, the diquark operator (qqlq), is suppressed or even forbidden by a symmetry, allowing LQs to have low masses [6][7][8][9][10].LQ properties and signatures have been extensively studied in the literature.We refer the reader to [11] for a comprehensive review on the LQ phenomenology at precision experiments and colliders.Moreover, LQs coupling to third-generation fermions have received much attention lately as likely candidates to account for flavor anomalies, see for instance [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30].
At present, LQ models have been studied through a variety of processes involving different ranges of energy.Depending on certain assumptions, different observables can test particular regions of masses and coupling strengths in the parameter space.On the one hand, Atomic Parity Violation (APV) in cesium nuclei has allowed to test the effects of LQs at low energy.Particle colliders, on the other hand, like LEP and LHC, have allowed to probe LQs in relatively large energy ranges through processes like Drell Yan, as well as through single and double pair production.Another important test of LQ interactions are electroweak precision observables.The presence of LQs could induce effects on the self-energy of the Z and W bosons constrained by the oblique parameters, S, T and U .However, such corrections mostly depend on the difference of masses between states [11].Since we are interested in LQ multiplets with degenerate states, we will not take these constraints into account in our study.A more complete investigation constraints [130,134,135] using Monte Carlo tools in order to obtain the most up-to-date limits from their searches.We finally mention that during the completion of this work a study on LQs has been presented in Ref. [136].Their analysis focuses on LQ interactions with third-generation leptons and second-generation quarks, and on their testability using neutrino telescopes, compared to other constraints including colliders.Their study complements the results presented here.
This paper is organised as follows.In Sec.II we introduce the general LQ picture, particularly focusing on scalar LQs that couple to the first and second generation of leptons and to the first generation of quarks, which are hence accessible at CEνNS experiments.In Sec.III we briefly describe the CEνNS process, the associated LQ effects, and the COHERENT experiment.In the same section, we also detail our statistical analysis, that allows us to set stringent constraints on LQ properties.We then describe in Sec.IV the effect of scalar LQs on other low-energy observables like APV, deep inelastic neutrino-nucleon scattering, and collider processes.We provide a summary of all results and current constraints in Sec.V. Next, we go back to CEνNS and present in Sec.VI the expected sensitivities at future upgrades of COHERENT and at the ESS.Finally, we summarize and present our conclusions in Sec.VII.

II. LEPTOQUARK FORMALISM
As their name suggests, LQs are particles than can simultaneously couple to a lepton and a quark and, in general, they can be of either scalar or vector nature.These two properties allow us to add different terms to the SM Lagrangian that couple LQs with the SM fields while preserving the gauge structure of the SM.A list of all possible ways by which LQs can give rise to a lepton-quark interaction is long and has been studied, for instance, in Ref. [11].However, since we are interested in signatures at CEνNS experiments, we will mainly focus on those interactions that can connect neutrinos and quarks.Respecting the SM symmetries, these interactions are possible through either a scalar or a vector LQ as long as they do not allow proton decay at tree level.Starting from the completely general LQ list given in [11], we extract all the relevant operators involving neutrinos and we summarize them in Table I, indicating in each case the LQ quantum numbers under the SM gauge group.For each scenario, we follow the notation of [11] and we denote the SM lepton and quark doublets by L and Q, respectively.We now briefly discuss the associated Lagrangian for each of the different scenarios.As seen from Table I, there are four possible scalar LQs giving rise to operators relevant for CEνNS (our main topic of interest): one SU (2) L singlet, two doublets, and one triplet.For simplicity, we will assume that LQs only interact with first-generation quarks and first-and second-generation leptons: Since this is a singlet under SU (2) L , there is only one component under this symmetry, precisely denoted by S 1 , which carries a charge Q = −1/3.Then, the corresponding Lagrangian that adds to the SM reads where λ ij is in general a complex matrix, τ 2 is the indicated Pauli matrix, and i, j = 1, 2, 3 are flavor indices.The right-hand side of the previous equation results from expanding the doublet terms in the Lagrangian.P L denotes the left-handed chirality operator and we just keep the λ 1j (i = 1) term to explicitly remark that we are interested in LQs coupling only to first-generation quarks, with u 1 = u and d 1 = d.To illustrate how such an interaction in Eq. ( 1) can contribute to a neutral-current process as CEνNS, we can see how the matrix element depends on the fermionic spinors involving neutrino reactions, Notice that to get the right-hand side of the previous equation we have performed a Fierz transformation (see Appendix A), giving the desired neutral current shape for the interaction.An important feature arising from Eq. ( 2) is that, when interacting with neutrinos, the scalar LQ S 1 couples only to down quarks.
Being a doublet under SU (2) L , in this case there are two LQ components, which we denote 2 ) T , where the superscript indicates the corresponding electric charge Q.From Table I, we see that R 2 couples only to u-type quarks.Then, the relevant Lagrangian for neutrino interactions is given by where, again, λ ij is a complex matrix, and the right-hand side is obtained by expanding the different doublets.Assuming a degeneracy in the associated masses of the two states, the matrix element involved in neutrino interactions with matter is proportional to where we have Fierz-transformed again the operator in order to get the expression on the right side of the equation.
As in the previous case, here we have a doublet under SU (2) L , with components denoted by R2 = ( R2/3 Then, this LQ can only couple to d-type quarks, and the relevant Lagrangian for neutrino interactions reads Again we assume a degeneracy in mass between both states of the multiplet, so the matrix element reads By comparing with Eq. ( 2), we can see that the relevant matrix element has the same structure as S 1 , and R2 couples only to down quarks.
To finish our scalar LQ list, we now have a triplet under SU (2) L , whose components are denoted as S 3 = (S ).Then, being τ = (τ 1 , τ 2 , τ 3 ) the standard Pauli matrices, the associated Lagrangian involving neutrino interactions and respecting the SM symmetries reads As in the previous cases we can see that the matrix element reads Interestingly, among all cases presented till now, only S 3 couples to both up-and downtype quarks.However, the strength of the corresponding coupling is not the same as there is a difference of a factor 1/2 between them.As we will see, this will result in an enhancement of the CEνNS cross section associated with the contribution of this type of LQ.
Notice that in Table I we have also listed operators involving vector-type LQs.In principle, these operators can also be studied, giving their corresponding contribution to the CEνNS cross section.However, if we perform a Fierz transformation the contribution turns out to have the same shape as the scalar cases listed above.To illustrate this we can take, for instance, the operator associated to the LQ denoted as U 1 .Then, if we apply a Fierz transformation to the matrix element of neutrino interactions we have that is the same Lorentz structure as R 2 with a factor 1/2 of difference.Given this similarity among Lorentz structures, in the following we will focus only on scalar LQs.
In the rest of the paper we will focus on scalar LQ interactions between first-generation quarks and first-and second-generation leptons.Hence the matrix structure of the parameters λ ij appearing in the Lagrangians of the different LQs considered here can be written as, where g is the strength of the interaction and we assume it to be the same for both lepton flavors.During the rest of this manuscript, we will consider this flavor structure and we will assume same strength of the coupling for the different LQ models under study.

III. CURRENT CEνNS DATA: COHERENT
In this section, we investigate the potential of CEνNS to study LQs and we obtain constraints from COHERENT data in the parameter space of LQ masses and couplings.Within the SM, the CEνNS differential cross section, in terms of the nuclear recoil energy E nr , is given by [33] where G F is the Fermi constant, E ν denotes the incoming neutrino energy, ℓ indicates the neutrino flavor, and m N refers to the nuclear mass.Notice that at tree level, the CEνNS cross section is flavor independent, with small radiative corrections that are not relevant for present experimental sensitivities [137].The SM weak charge Q SM W is defined as where g p,n V are the proton and neutron couplings, g p V = 1/2(1 − 4 sin2 θ W ) and g n V = −1/2.The weak charge is the term which eventually encodes the typical N 2 dependence of the CEνNS cross section, and which gives rise to the relevant enhancement with respect to other neutrino processes.Notice that the proton contribution carries the dependence on the weak mixing angle, being this contribution subdominant due to an accidental cancellation generated by the SM value of the weak mixing angle at low energy. 1uclear-physics effects are encoded in the nuclear form factor F W (|q| 2 ) appearing in Eq. (11).We adopt the Klein-Nystrand parametrization, which reads where J 1 is the spherical Bessel function of order one, |q| ≈ √ 2m N E nr stands for the threemomentum transfer, R A = 1.23A 1/3 fm is the nuclear radius and a k = 0.7 fm is the Yukawa potential range.
The weak charge term may be modified in the presence of new physics.In the specific case of the LQ scenarios of interest in this paper, we will compute the CEνNS events based on Eq. ( 11) by changing where the first and second indexes in Q ij,LQ denote quark and lepton family, respectively, and LQ = S 1 , R 2 , R2 , S 3 stands for the LQ type.For simplicity, we assume that LQs couple with the same strength to electrons and muons, and to u and d quarks, with vanishing coupling to τ neutrinos and to the second-and third-generation quark families.Then, we denote g 2 ≡ λ 1i λ 1j (i, j = e, µ) and for the different models studied in Sec.II we have being m LQ the corresponding LQ mass.As it is clear from Eqs. ( 15) and ( 17), the impact of S 1 and R2 on CEνNS is expected to be exactly the same.For the last case (S 3 ), we have assumed that the two LQs arising from the SU (2) L triplet have the same mass.In that case, one of the states couples to the u-type quarks, while the other one couples to the down quarks.However,  given the parametrization used for the Lagrangian in Eq. ( 8), they do not couple with the same strength.
To set constraints on the LQ scenarios using CEνNS data, we rely on the most recent measurements of the COHERENT experiment, which were performed by using CsI [35] and LAr [139] detectors, whose specifications are summarized in the first two lines of Table II.We perform a thorough analysis by including both energy and timing information together with all relevant systematic effects, for each detector, following Ref.[58]. 2 The neutrino flux at COHERENT comes in three components, from π-DARs produced at the Spallation Neutron Source (SNS): where m µ , m π denote the muon and pion masses, while η = rN POT /4πL 2 is a normalization factor which depends on the number of neutrinos per flavor (r) produced for each proton on target (POT).We assume r = 0.0848 (0.009) and N POT = 3.198 (1.38) × 10 23 for the CsI (LAr) detector.Notice that the three different neutrino flux components come with different timing, the ν µ being prompt and the other two components delayed.Next, we proceed to evaluate the expected number of events.We assume a detector mass m det = 14.6 ( 24) kg located at a distance L = 19.3(27.5) m from the SNS source, for the CsI (LAr) detector.The expected number of events, on a nuclear target N , per neutrino flavor, ν ℓ , and in each nuclear recoil energy bin i can be written as [58,101] where N target = N A m det /M target is the number of target atoms in the detector, with M target the molar mass of the detector material, and N A the Avogadro's constant.Kinematically, the integration limits in Eq. ( 20) are found to be the maximum incoming neutrino energy, which for SNS neutrinos is ≈ 52.8 MeV.Finally, the energy resolution function R(E nr , E ′ nr ) appearing in Eq. ( 20) associates the true nuclear recoil energy (E ′ nr ) with the reconstructed one (E nr ) and the ϵ E (E nr ) is the energy-dependent detector efficiency.We refer the reader to Refs.[35,58,101] for more details.
In order to take into account the neutrino-flux timing information in our analysis, we distribute the predicted N CEνNS i,ν ℓ (N ) in each time bin j.At this scope, we rely on the time distributions P ν ℓ T (t rec ) provided in [35,140], and we normalize them to 6 µs [58,101].The predicted event number, per observed nuclear recoil energy and time bins i, j is finally obtained as where ϵ T (t rec ) is the time-dependent efficiency [35,58,101].(We include an additional nuisance parameter on the beam timing, α 6 , see [58,101].) To proceed with the statistical analysis of the COHERENT CsI data set we consider a Poissonian least-squares function [58,101], expressed as The predicted number of events, which includes both SM and LQ CEνNS events, as well as backgrounds, depends on several nuisance parameters (α i ) and reads The nuisances come together with their associated uncertainties σ i [58,101]: σ 0 = 11% (efficiency and flux uncertainties), σ 1 = 25% (Beam Related Neutrons (BRN)), σ 2 = 35% (Neutrino Induced Neutrons (NIN)) and σ 3 = 2.1% (Steady State Background (SSB)), σ 5 = 3.8% (QF).The predicted number of events N th ij also depends on three nuisance parameters: α 4 , which enters the nuclear form factor through the nuclear radius in Eq. ( 13), via R A = 1.23A 1/3 (1 + α 4 ), with σ 4 = 5%; α 6 which accounts for the uncertainty in beam timing with no prior assigned; and α 7 which allows for deviations of the uncertainty in the CEνNS efficiency.
For the statistical analysis of the COHERENT-LAr data set we instead adopt the following Gaussian least-squares approach based on [58,80,101] Here, the theoretical number of events is defined as and the experimental uncertainty is . The expected number of events depend on several nuisance parameters, dubbed β 0 , β 3 , β 4 and β 8 , which account for the normalization uncertainties of CEνNS, SS, prompt BRN (pBRN) and delayed BRN (dBRN) background rates respectively, with uncertainties {σ 0 , σ 3 , σ 4 , σ 8 } = {0.13,0.0079, 0.32, 1.0} [34].Let us notice that β 0 encodes multiple uncertainties, namely the flux (10%), efficiency (3.6%), energy calibration (0.8%), the calibration of the pulse-shape discrimination parameter F 90 (7.8%),QF (1%), and nuclear form factor (2%) [34].The additional nuisance parameters β 1 , β 2 , β 5 , β 6 and β 7 account for systematic effects affecting the shape uncertainties of the CEνNS and pBRN rates, namely the uncertainty on the CEνNS shape due to existing systematic uncertainties on the ±1σ energy distributions of the F 90 parameter (∆ CEνNS ), due to the mean time to trigger distribution (∆ t trig CEνNS ) or the pBRN shape uncertainty due to the corresponding uncertainty on the ±1σ energy, time and trigger width distributions ).These distributions are defined as departures from the central value (CV) ones [139] , with λ = {CEνNS, pBRN} and ξ λ referring to the different source uncertainties affecting the CEνNS or pBRN shapes.

IV. OTHER CONSTRAINTS
In this section we proceed to discuss further current constraints on the LQ scenarios presented in Sec.II.Following an increasing energy scale, we will start with APV, then proceed with deep inelastic neutrino-nucleon scattering, and finally move to collider searches.

A. Atomic Parity Violation
One very accurate determination of the weak mixing angle currently available in the low-energy regime comes from APV -or parity nonconservation -experiments on cesium atoms [110,111,141].It has been shown [55,60,142] that such a measurement can provide complementary information to CEνNS, also regarding nuclear physics parameters besides the weak mixing angle.Moreover, stringent APV bounds on LQ coupling to first-generation fermions have been obtained in the literature under the assumption of effective four-fermion interactions and that only one contribution (from u or d quarks) is present at a given time [11,114,143,144].
Here we want to exploit the low-energy measurement of the weak charge Q W of 133 Cs from APV experiments to constrain the LQ scenarios proposed in Sec.II, including model S 3 which simultaneously encodes couplings to both u and d quarks, and taking into account the explicit dependence on the LQ mass.In this subsection, we hence derive APV constraints on LQ through their effect on the weak charge.Including radiative corrections in the MS scheme, the APV weak charge in the SM reads [55,142,145,146] where the couplings of electrons to nucleons are g ep AV = 2g eu AV + g ed AV = −0.0357and g en AV = g eu AV + 2g ed AV = −0.495,α is the fine-structure constant and Z = 55, N = 78 for cesium.The theoretical expected value is therefore [55,142].As anticipated, the presence of LQs will affect the value of the weak charge as following The LQ charge Q ee,LQ encodes the dependence on the free parameters g and m LQ , and on the momentum transfer, which for APV is |q| 2 ≃ (2.4 MeV) 2 .The experimental value of Q W for cesium is extracted by measuring the ratio of the parity violating amplitude E PNC to the Stark vector transition polarizability, and by calculating theoretically E PNC as a function of Q W [146]. Taking into account small uncertainties associated with the atomic wave function calculations, most recent computations of the parity non-conserving amplitude combined with the measurements [110,111] We evaluate the APV bound on the LQ scenarios by minimizing the following least-square function: where σ APV = 0.42 is the total (experimental + theoretical) uncertainty.
B. Deep inelastic neutrino-nucleon scattering (NuTeV) The neutrino scattering off nuclei is a very accurate process by which the nature of the weak currents can be tested.Usually, neutrino scattering experiments make use of neutrino beams originating directly from colliders.This is the case for the NuTeV experiment that benefits from the Sign Selected Quadrupole Train (SSQT) beamline at the Fermilab Tevatron collider to obtain well-controlled muon neutrino beams and test the neutrino-nucleon cross section with iron targets.NuTeV measured the deep inelastic neutrino-nucleon scattering cross section with high accuracy [115], improving upon its predecessors CDHS [147] and CHARM [148].Provided the target of neutrino experiments is isoscalar, contributions to the cross section from neutral and charged currents can be written as [116] where the coupling constants are (g ) 2 and r is defined as [116] By measuring these ratios one can measure the coupling of neutrinos to quarks.In the presence of new physics, these ratios will show a deviation from the SM predictions.One relevant example is the case of neutrino NSI.In order to parameterize the presence of new physics giving rise to NSI, one can define the low-energy effective Lagrangian as [65] where P = L, R are the chiral projectors and ε f P αβ are the coefficients that parameterize the NSI, where α, β run over the different neutrino flavors, and f are the charged fermions of the SM.Still focusing on NuTeV, we can relate the NSI coefficients to the couplings (g L,R µ ) in the following way [116] where we set α = β = µ.
The introduction of LQs may alter the interaction between neutrinos and quarks, since they couple directly to them, for this reason the presence of LQs may induce non-zero NSI coefficients.For instance, in the case of S 1 , for large masses we can integrate out the LQ degrees of freedom, thus obtaining where we have performed a Fierz transformation to get the right-hand side of the equation.Analogously, for the other LQ scenarios we have Then, we can compare Eq. ( 31) with Eqs. ( 33), ( 34), (35), and (36) to obtain the relations between NSI coefficients and LQ parameters [69,104,113] It is important to note that given the Lagrangian of R2 , the coefficient for this state is the same as the one for S 1 .
NuTeV data have been studied in Ref. [116] where a detailed analysis was done to translate the experimental measurements into the NSI parameter space.We have recast their results to obtain limits on the LQ mass-coupling plane through the NSI coefficients.It is important to note that due to the fact that the NSI parametrization is given in terms of an effective interaction, NuTeV bounds apply for LQ masses larger than 10 GeV, where the effective theory holds.

C. Collider data
As LQs couple to both quarks and leptons, they are very likely to be produced in lepton or hadron colliders.Furthermore, given the nature of their couplings, LQs may give rise to a very interesting set of signatures [11,149].Depending on the collider nature, the production of LQs may differ and also the signatures obtained at the detectors.Since we focus on LQ interactions with neutrinos, it is reasonable to expect that they will also interact with charged leptons.This interaction could be mediated by the same LQ state, as it happens for S 1 , by the other charged state from the LQ multiplet, as for R 2 and R2 , or eventually by both at the same time, as it occurs in the S 3 scenario.Given that we are considering the different states from the multiplets to be mass-degenerate, we will translate the bounds on the masses to the state that couples to neutrinos.Moreover, since the couplings are the same before decomposing the multiplets, limits on the LQ couplings apply to all states.In this subsection we will focus on different kind of colliders regarding the nature of their collisions.First, we will recast data from HERA (an e − p collider), LEP (e + e − ), and both SPS and Tevatron (pp).Then, we will compute bounds from the LHC proton-proton collider.
HERA was an electron-proton collider that operated at center-of mass energies up to 320 GeV in the regime of deep inelastic scattering.In this regime, a better understanding of the nature of the proton was reached.However, other results were obtained using the deep inelastic scattering of electrons and protons.For example, at those energies LQs can contribute to the total electron jet cross section in both s and t channel, as we can see in the first two diagrams of Fig. 1.It is important to notice that LQs are only resonantly produced in ep colliders.For that reason, the experiments H1 and ZEUS performed searches for LQs in electron plus jet final states [117][118][119].To set limits in the mass-coupling plane, we recast the bounds from Ref. [119].This search focused on the production of LQs of first generation leading to an electron and a jet signature at the ZEUS experiment with a luminosity of 498 pb −1 .Using this search, the ZEUS collaboration could set constraints on the LQ production in the 150 GeV to 1 TeV mass range.To recast these results we have taken into account the coupling structure of the models presented in Sec.II and weighted the data with the corresponding branching ratios.
In electron-positron colliders it is also possible to measure the presence of LQs.As we can see in the last two diagrams of Fig. 1, LQs can contribute through t and u channels to the total dijet cross section.The L3 and OPAL experiments from LEP performed searches in the dijet cross section looking for new physics [122,124].We have recasted the search of OPAL [122] that looked for constraints on new physics in the dijet cross section.This search can be translated into limits on LQ masses from 100 to 400 GeV.We have weighted their data according to the characteristics of our models in order to recast the search.Furthermore, experiments at LEP have also looked for LQs produced in decays of an on-shell Z boson [120,121].Since no positive results were found, this can be translated into a lower limit of m LQ > 45 GeV.Finally, OPAL at LEP has also searched for charged long-lived LQ production [123].In order to recast this search we have to verify the range of validity where the LQs are long-lived.The decay length of a long-lived particle is given by L = γβcτ , where γβ is the boost factor and τ is the proper lifetime of that particle.For scalar LQs, assuming that the decay products have smaller masses than the initial particle, we find where g is the LQ coupling and m LQ is its mass.As the proper lifetime of a particle is given by τ = ℏ/Γ we can then obtain an expression for g where the long-lived regime starts (assuming that a long-lived particle can be identified when it has a decay length of 0.1 cm [150]) Using the information provided in [123] to compute the γ and β factors for different masses and center-of-mass energies, this gives us a range of validity for the long-lived LQ search of g < 2.0 × 10 −7 .On the other hand, proton-antiproton colliders have also been used to look for LQ signals.In particular, the UA2 experiment at SPS has performed a search for double production of LQs that decay into a charged lepton plus jet or into a neutrino plus jet [125].Double LQ production in a proton-antiproton collider does not depend on the specific LQ scenario, however the identification of the different LQ decays relies on the different branching ratios of the produced LQ, and hence from the model under scrutiny.We have recast the results from [125] according to our specific benchmark models, obtaining a mass limit of m LQ < 50 GeV.Other LQ searches have been done at Tevatron, where double pair production is the leading production mechanism that allows to constrain the LQ mass.Both CDF and DØ experiments at Tevatron have performed such a search [125][126][127][128][129] for different luminosities, being the production of LQs independent of the LQ coupling, as in the case of the UA2 bounds.However, the limits derived by the collaboration are dependent on the branching ratios.We have recast such limits into our scenarios and we have found that for our branching ratios these bounds are less powerful than the one obtained by UA2 [125].In the case of Ref. [129] the experimental results do not show the results for the specific branching ratios typical of the benchmark points considered in our study (i.e.BR ≲ 25% for each channel).However, we have decided to take as a reference the limits of the most similar case (BR = 50%), even if in this case the resulting exclusion region is overestimated.Taking all that into account, we have found a window in tension with Tevatron data [129] that lies in the mass range of 150 < m LQ < 260 GeV.
Finally, LQs can also be produced in multiple ways at the LHC.For this reason there are different searches from ATLAS and CMS that look for LQs taking into account their multiple production channels and decays.One of these searches is the double production of LQs.Being LQs color triplets, they can be created in the LHC in events initiated by gluons and quarks.Fig. 2 shows all the diagrams that contribute to the LQ double production.As we can see, the first four diagrams are initiated by the strong coupling, so they scale as α 2 s , which is independent on the LQ coupling.The last diagram is the only one that depends on the LQ coupling, g 4 .Because of this, we expect this kind of searches to be independent of LQ couplings when their value is small, and the limits will only display their dependence on g when they reach values of O(1) [149].For that reason, this search can exclude LQ masses independently of the LQ coupling.ATLAS and CMS have performed several searches for double production of LQs in multiple final states, mainly the presence of two quarks and the combination of charged and neutral leptons, qqℓℓ, qqℓν and qqνν, where, depending on the LQ generation, signals may vary into different flavors of leptons.Nonetheless, given the purposes of our work, we are only interested in those searches that contain first-and second-generation leptons in the final state [130][131][132].In order to set limits from LQ double production we have recast the search from Ref. [130] using the recommendations of Refs.[149,151].To simulate the LQ double production we have made use of MadGraph5_aMC@NLO-v3.5.0 [152,153] using the codes and recommendations from Refs.[154][155][156].We have then compared the production cross sections weighted by their corresponding branching fractions with those from Ref. [130] to obtain the limits.Finally, it is important to note that even if the limits from LQ pair production are independent of the LQ coupling, g, when this is small (g ≪ 1), there is an actual limit to the constraints on this coupling.This is given by the fact that the LQs could be long-lived due to the smallness of the coupling.Hence we can use the same procedure as we used above, following Eq.( 40), in order to infer an estimate of the reliability of the double production bounds.Using this method we obtain that this search is valid for values of the coupling g > O(10 −6 − 10 −7 ).
Another imprint of LQs at LHC is through dilepton production.As LQs couple to both LQ * ℓ q ℓ q LQ ℓ/ν q ℓ/ν q g ℓ/ν q ℓ/ν q g LQ * LQ FIG. 3. Diagrams contributing to Drell Yan production and single LQ production.The first diagram corresponds to the LQ contribution to Drell Yan while the last ones are the responsible for LQ single production.
quarks and leptons they can contribute to the Drell Yan cross section as it is shown in the first diagram of Fig. 3.The LQ appears in the t-channel of a process initiated by a pair of quarks giving as a final results two leptons.The presence of LQs in this process can interfere with the SM processes mediated by γ/Z.ATLAS and CMS experiments have performed several searches studying the dilepton cross section [133,134].As we can see from the first diagram in Fig. 3, this process is coupling-dependent, contrary to what we had for the double production mechanism.To understand how strong these limits are, we simulate the LQ cross section using MadGraph5_aMC@NLO-v3.5.0 [152,153] and again, using the codes from [154][155][156] we compare it against data from Ref. [134].Concerning this last step, we follow the prescriptions of Ref. [149] to obtain the cross section with the cuts from the search and compare the results against actual data.Another important channel to consider in hadron colliders is the single LQ production.LQs can be produced together with a lepton in processes initiated by a quark and a gluon, as it is shown in the last two diagrams of Fig. 3. Once the LQ is produced it will subsequently decay into a quark and a lepton, leading to different final states containing either one jet and two charged leptons, or one jet, one charged lepton and missing transverse energy, or simply one jet and missing energy.Several searches from ATLAS and CMS have been looking for these signals: we choose to recast the CMS search for energetic jets and missing transverse energy [135] at a center-of-mass energy of √ s = 13 TeV and a luminosity of L =101 fb −1 since it proves directly the coupling of LQ to neutrinos.These limits are also computed using MadGraph5_aMC@NLO-v3.5.0 [152,153], benefiting from the codes of [154][155][156] and compared with the data obtained from Ref. [135].

V. RESULTS
In this section, we summarize the results for the constraints obtained with all different processes discussed in the previous section, each covering a different range of validity in the LQ parameter space (m LQ , g).Notice that, for all our computations, we have assumed the same coupling constant g for all leptons except for taus, for which we assume a vanishing coupling.The main results for each LQ scenario, in terms of the LQ mass m LQ and coupling g, are shown in the different panels of Fig. 4, where we show the 90% C.L. exclusion limits from each case.The top-left panel in the figure shows the constraints obtained for the scalar scenario listed as S 1 in Table I.From the figure, we can distinguish the different mass ranges that each experiment is able to test, and hence, we can infer the complementarity between different observables and facilities to constrain a wide region of parameter space.Colored regions in the figure indicate new results computed in this work, while grey-shaded regions correspond to previous limits found in the literature, including those from NuTeV [115] and colliders such as ZEUS at HERA [119], OPAL at LEP [121,122] and UA2 at SPS [125], which we have recast for the LQ scenarios of our interest.
Starting from low-energy observables, APV upper limits obtained from Eq. ( 28) are indicated as magenta contours,3 while COHERENT excluded regions (see Eqs. ( 22) and ( 23)) are shown as different shades of blue.Dark blue indicates the limits from the COHERENT-LAr (2020) detector alone, light blue stands for the COHERENT-CsI data set (2021), and cyan is used for their combined analysis.We can notice that, when considered individually, the excluded region for COHERENT-LAr is not continuous, as it contains a tiny allowed band.This degeneracy is a consequence of a destructive interference between the SM and the LQ contributions to the CEνNS cross section, resulting in a combination of non-zero parameters that can mimic the SM solution.Although not visible, the situation is similar for the COHERENT-CsI analysis.However, when combining the results from the two detectors, the degeneracy is lifted, resulting in an excluded cyan region which is now continuous.We refer the reader to Appendix B for a detailed discussion on the origin and effects of this interference.Moreover, when comparing to previous results that use the older CsI data set [100], we can see that new data allow to constrain slightly lower couplings and their combination with the LAr result allows to remove the SM degeneracy, as just discussed.
Still referring to the top-left panel in the figure, and moving to the heavy mass (and energy) regime, the yellow region in the figure corresponds to LHC constraints that have been obtained including channels like single production, double production and Drell Yan (see discussion in Sec.IV C).Regarding LHC constraints, it is worth mentioning that, within their validity range, double pair production limits are mass independent for g ≲ 3 × 10 −1 .This is shown as a vertical yellow band that extends for LQ masses in the range 400 GeV ≲ m LQ ≲ 1400 GeV.On the other hand, we further show as grey-shaded regions those excluded by LEP, UA2, HERA, Tevatron and NuTeV, which we recast from pre-existing analyses (see Secs.IV B and IV C).In general, these bounds apply to heavy LQ masses, and some of them have been obtained under effectivetheory assumptions, except for the LEP search, which probes LQs from an on-shell Z boson decay, and can be extended to very low masses.For the displayed mass range, we see that LQ searches in Z decays from LEP are able to constrain masses below m LQ ≲ 40 − 45 GeV, this result being overtaken by the double production of LQ at UA2, which excludes masses up to m LQ ≲ 50 GeV.Let us remark that the LQ double production channel is actually independent on the LQ coupling.However, the LEP and UA2 bounds that we derived here depend upon the flavor structure of their couplings and therefore upon their branching ratios, being these searches less sensitive to smaller branching fractions.As we work under the assumption that LQs couple to both the first and second lepton families scenarios, the best limit we get from LEP and UA2 double production searches is m LQ ≲ 50 while leaving heavier masses unconstrained.As a result, we see that COHERENT data lead to the most stringent constraints in one small region in the parameter space that goes from m LQ ∼ 50 GeV and up to m LQ ∼ 100 − 150 GeV.In this small window, we see that COHERENT clearly overtakes former bounds from NuTeV, APV and LEP.For masses greater than m LQ ≳ 100 − 150 GeV, LEP, HERA and Tevatron become the most stringent bounds, being the latter two the most powerful constraints in terms of the LQ coupling due to its nature as an ep collider that can produce on-shell LQ for the case of HERA, and due to the double LQ production in Tevatron.As mentioned before, it is important to note that the bounds imposed by Tevatron searches are over-estimated due to different assumptions on the branching ratios, and the actual limits would lie within the area set by HERA.However, the strength of the bounds weakens for masses around m LQ ∼ 400 GeV.From this LQ mass and on, LHC constraints dominate thanks to the double pair production, which leads to coupling-independent bounds on the LQ mass.As anticipated in Sec.IV C, while covering the whole g parameter space in these panels, the LHC-excluded yellow band is expected to extend down to ∼ 10 −6 −10 −7 due to LQs lifetime considerations.The small region at g ∼ 0.8 bounded at m LQ ∼ 2 TeV comes from Drell Yan processes.All in all, we can safely conclude that for m LQ ≲ 0.1 TeV the dominant exclusion process is CEνNS, thus complementing the strong collider bounds which instead dominate above 0.1 TeV.It is important to note that existing FIG. 4. 90% C.L. excluded regions, in the (m LQ , g) plane, on different LQ scenarios, S 1 , R 2 , R2 , and S 3 .Colored contours and filled areas denote new upper bounds obtained in this work: APV (magenta line); CEνNS data from the COHERENT-CsI (2021) [35] and COHERENT-LAr (2020) [34] data sets (different shades of blue); single production, double production and Drell Yan processes at LHC (yellow region).For comparison, we also show previously obtained limits in the literature and recast them here into the LQ scenarios under scrutiny (grey-shaded regions): NuTeV [116], ZEUS at HERA [119], UA2 at SPS [125], CDF and D0 at Tevatron [129], and OPAL at LEP [121,122].See main text for more details.long-lived charged particles searches, such as [123] by LEP, constrain the LQ parameter space for couplings below g ≲ 10 −7 , however these constraints lie outside the parameter range shown in these plots.
The top-right and lower panels in Fig. 4 show the corresponding results for other LQ scenarios: R 2 (top right), R2 , (bottom left) and S 3 (bottom right).Overall, we see a similar behaviour for all cases, being CEνNS the dominant channel to constrain LQ masses in one small window that goes from 50 GeV ≲ m LQ ≲ 150 GeV for the S 3 scenario, and for two small windows that go from 50 GeV ≲ m LQ ≲ 150 GeV and 300 GeV ≲ m LQ ≲ 400 GeV in the case of R 2 and R2 .Above these mass ranges, LHC data have a major ability in setting constraints on the LQ coupling g.
When comparing the different panels, notice that the most stringent CEνNS constraint is found for S 3 .This was expected from the modified weak charge defined in Eq. ( 18), where we see that the CEνNS cross section in this case is effectively enhanced with a factor (4N ) 2 for fixed m LQ and g.On the other hand, for R 2 the CEνNS constraint is less robust because of a factor 2Z in the cross section given in Eq. ( 16) which, given the different relative sign between g p V and g n V , results into a smaller cross section, and hence a lower number of events expected in the statistical analysis.Another interesting feature is that, when coupling to neutrinos, scenarios S 1 and R2 are indistinguishable for CEνNS(see Sec.II), and in consequence, the excluded blue-shaded regions in the two left panels are the same.However, this is not the case for collider observables since, when coupling to charged leptons, S 1 couples only to up quarks while R2 couples only to down quarks.Then, because of the ratio between up and down quarks within the proton, this results in different excluded yellow regions in the top and bottom left panels of Fig. 4.

VI. FUTURE SENSITIVITIES
After having analyzed the current picture of LQ constraints in the parameter space (m LQ , g), we now turn our attention to sensitivities that can be reached at future CEνNS experiments.We consider upcoming upgrades of both the CsI and LAr detectors planned by the COHERENT collaboration, as well as two of the different detectors from a proposal at the ESS discussed in [50].We discuss these prospects in the following.
A. CEνNS data (COH-CsI-700 and COH-LAr-750) The intense experimental program of the COHERENT collaboration envisages, among others, upgrades of current detectors, namely a 700-kg cryogenic CsI scintillator and a tonne-scale LAr time-projection chamber detector [106,107].Moreover, planned up-scales of the SNS proton beam foresee an upgrade of the proton energy E p = 0.984 → 1.3 GeV and of the beam power P = 1.4 → 2 MW.By assuming a data-taking time of 5000 hr per year, this leads to N POT = 5.18 × 10 23 (for three years) [60,106,107], and to a predicted number of neutrinos per flavor produced for each POT r = 0.0848 → 0.13 [157].We estimate the future sensitivities for the COH-LAr-750 and COH-CsI-700 updates of the COHERENT detectors, assuming the technical upgrades summarized in Table II, and a detector mass of 750 and 700 kg, respectively.We perform a statistical analysis in energy and time following that done for current data and previously detailed in Sec.III.In the case of COH-CsI-700 we take into account the expected improvement in energy sensitivity, by considering a threshold of 1.4 keV nr [60,107] while keeping the shape of the energy efficiency unaltered.Pragmatically, we add in the statistical analysis an extra bin in energy, from 1 to 4 PE.Concerning backgrounds, keeping in mind that the collaboration anticipates that the cryogenic technology will allow to reduce them (in particular to remove the Cherenkov radiation background), we choose to be conservative and re-scale the current BRN, NIN and SSB backgrounds to the new detector's mass.Moreover, again assuming a conservative approach, we fix the numbers of background events in the first, new energy bin ([1-4] photoelectrons) to be exactly the same as in the second one ([4-8] photoelectrons).
The expected sensitivities for these two future COHERENT detectors are shown in Fig. 5, where we give the results for the S 1 scenario.The left and right panels in the figure correspond to the COH-CsI-700 and COH-LAr-750 upgrades, respectively.The region within the colored blue lines in each panel indicates the expected excluded values for masses and coupling constants.Similarly to what already observed for current (individual) sensitivites in Fig. 4, we see the presence of an allowed band within each of the excluded regions, which corresponds to combinations of LQ masses and parameters that allow for a destructive interference with the  The corresponding results for scenarios R 2 and S 3 are shown in Fig. 6 and Fig. 7, respectively, showing a similar qualitative behaviour to S 1 .(We recall that model R2 is equivalent to S 1 from the CEνNS point of view.)Overall, we notice that a future LAr detector is expected to enhance current constraints of up to around 50% when compared to the constraints obtained through its (current) predecessor, while CsI will be able to improve by almost one order of magnitude in some regions of the parameter space.

B. CEνNS data (ESS)
In addition to the COHERENT program, there are other collaborations aiming at performing new CEνNS measurements.Here we consider the particular case of the ESS, a facility that will be located at Lund, Sweden, and that at full power will become the most intense neutron beam source in the world.The physics potential of the ESS within the context of particle physics is summarized in Ref. [158].Furthermore, a proposal of measuring CEνNS at the ESS was presented in Ref. [50], and different analyses have explored its sensitivity to new physics, particularly within the context of NSI [50,159] and electromagnetic properties of neutrinos [50].
Here we explore the sensitivity of the ESS to scalar LQ models by analyzing two of the proposed detection technologies [50], namely silicon and xenon, characterized by having a very different ratio of protons to neutrons.In contrast to CsI and LAr, for these detectors we compute the expected number of events simply through Eq. ( 20), by separating the data in nuclear recoil energy bins as done in [159].
Being a spallation source, the total neutrino flux at the ESS will also have the contributions from prompt and delayed neutrinos as given in Eq. (19).At full capacity, this facility will operate at a beam energy of 2 GeV and a beam power of 5 MW, resulting on an N POT of 2.8 × 10 23 per calendar year of operations (≈ 5000 hours), with a number of released neutrinos per flavor of r = 0.3.We assume these values for our analysis.As a result, the ESS will provide larger statistics when compared to current measurements at the SNS, the upshot being a smaller beam frequency pulse that disfavors background discrimination in particular from SSB, expected to be the dominant contribution among all backgrounds.Regarding detectors' characteristics, the considered mass, baseline, and threshold for Xe and Si are given in Table II.In addition, we follow the procedure proposed in Ref. [50] and for the analysis we consider a Gaussian smearing distribution with a resolution σ = σ 0 √ T Th E nr , with 0 = 0.40 (0.60) for Xe (Si), T Th being the energy threshold.
Given the absence, at present, of timing information for this proposal, to infer LQ sensitivities for this experiment, we minimize the following Poissionan χ 2 function where the index i runs over the recoil energy bins.The predicted number of events in this case is given by where N CEνNS i stands for the expected number of CEνNS events as a function of the LQ model parameters under study, and N bckg i is the number of background events.As discussed above, the large pulse shape at the ESS makes it more difficult to discriminate background events from SSB contributions.Being this the dominant background component, we model it as expected counts per keV per kilogram per day (ckkd), as also done in Ref. [50], where it is assumed a value of 10 ckkd (1 ckkd) for Xe (Si).Going back to Eq. ( 41), N exp i is the experimental number of events, which we assume as the SM prediction.To perform the analysis, the χ 2 function in Eq. ( 41) is minimized with respect to the nuisance parameters α and β, which are associated to the predicted CEνNS and background events, respectively, each with its corresponding uncertainty taken as σ α = 10% and σ β = 1% [50].
The expected sensitivities for the described ESS detectors, at 90% C.L., are shown in Fig. 8 for all LQ models.Colored lines in the figure represent the exclusion regions for Si (red) and Xe (yellow).Notice again the presence of an allowed band within the region, whose position and width depend not only on the specific LQ scenario considered, but also on the ratio of protons to neutrons of the target material (see Appendix B).The grey-shaded regions correspond to current excluded limits from colliders, DIS, APV, and the CEνNS bounds obtained by the combination of COHERENT CsI and LAr detectors.Given the larger nuclear mass, we see that, for all models, better sensitivies are expected for Xe when compared to Si, getting a better improvement for the S 1 and S 3 cases when compared to current bounds.

VII. CONCLUSIONS
In this work we have explored the potential of CEνNS in probing scalar leptoquarks.We have considered four different models, each of them giving rise to a different contribution to the weak charge.First we have analyzed current COHERENT data, from the CsI (2021) and the LAr (2020) detectors.By means of a detailed statistical analysis, which took into account timing information and all experimental uncertainties, we obtained stringent constraints on the LQ mass and couplings.We further obtained upper limits on the LQ parameter space from atomic parity violation experiments, which turned out to be comparable (although slightly less stringent) to COHERENT bounds.Next, we have obtained bounds on the same LQ models from LHC data, considering different processes and production mechanisms: single production, double production and Drell Yan.These strong collider bounds lead to an exclusion region in the mass range 0.4 ≲ m LQ ≲ 1.5 TeV, independent of the LQ coupling.To complete the picture on the LQ parameter space we have also recast bounds from deep inelastic neutrino-nucleon scattering (NuTeV) and older colliders (HERA, SPS, LEP and Tevatron).Among them, UA2 at SPS and OPAL at LEP set strong constraints on m LQ ≲ 50 GeV, while HERA and Tevatron disfavor a thin region around 0.2 TeV.However, we have identified two regions in parameter space where CEνNS data may improve upon existing constraints and provide a complementary probe, at 50 GeV ≲ m LQ ≲ 150 GeV and 300 GeV ≲ m LQ ≲ 400 GeV, depending on the LQ scenario.
Additionally, we have computed sensitivities at future upgrades of the COHERENT CsI and LAr detectors, and at the European Spallation Source.We have found that these future facilities, thanks to their larger exposures and exquisitely low thresholds, will allow to improve upon current bounds by up to a factor 3 at m LQ ∼ 100 GeV.Let us mention that we analyzed CEνNS data using π−DAR neutrinos motivated by current available measurements, however there is a vast array of experiments using reactor neutrinos [36][37][38][39][40][41][42][43][44][45][46][47][48] that can also provide valuable information on LQ scenarios and are therefore worth studying in a future work.
As a last remark, it is important to note that the present Run 3 of LHC will soon allow to explore the LQ parameter space involving large masses.With increasing luminosity, searches like double pair production will be able to set coupling-independent limits on LQ masses that lie in the TeV range.Furthermore, single LQ production and Drell Yan will be able to cover even higher masses, imposing constraints on couplings of the order of O(10 −1 ).In addition, there are new searches in the literature that are specific for LQ signatures and that are not yet exploited by the experiments.One example is the single-lepton channel initiated using the lepton content in the proton [160].This search is more sensitive than others listed before, in particular being more powerful than Drell Yan up to LQ masses m LQ ∼ 4 TeV, so it could be decisive in the search for high mass LQs. with C 1 and C 2 ∈ N. The shape of the weak charge in Eq. (B2) allows for different pairs of g and m LQ to reproduce the SM cross section, and hence the SM prediction for the number of events, giving as a result a degeneracy in the parameter space, which corresponds to an allowed band.For instance, in the case of the future sensitivities studied in Sec.VI, if we want to reproduce the SM number of events (for which ∆χ 2 is minimum), under the assumption in Eq. (B2), we need which can be satisfied when Q ii,LQ = −Q SM W or, in other words, when the parameters g and m LQ are such that Notice that for m LQ ≫ 2m N E nr the needed value of g 2 is energy-independent and we have For SNS and ESS neutrinos, the condition m 2 LQ ≫ 2m N E nr can be easily satisfied with a relatively large LQ mass, m LQ .For instance, let us consider the S 1 LQ scenario, for which C 1 = 1 and C 2 = 2.Then, regardless of the target material, for a mass of m LQ = 100 GeV, we have 2m N E max nr /m 2 LQ ≈ 1 × 10 −6 and Eq.(B5) safely applies.Then, in this case it is possible to find a solution in the parameter space that reproduces the SM solution and the χ 2 reaches a minimum, as we can see in the two panels of Fig. 9, where we show the ∆χ 2 profile as a function of g for m LQ = 100 GeV.The profiles are shown for the COH-LAr-750 and COH-CsI-700 detectors in the left panel, and for ESS-Si and ESS-Xe in the right panel.These results are consistent with the allowed bands observed in Figs. 5 and 8 at m LQ = 100 GeV.The situation would be different for low LQ masses, when the terms 2m N E nr and m 2 LQ are comparable and Eq.(B5) does not hold.However, we are not interested in such low masses given the LEP and

FIG. 1 .
FIG.1.Relevant LQ diagrams in ep and ee colliders.The first two diagrams correspond to the LQ contribution to electron jet production in HERA while the last two represent the LQ contributions to the dijet production in LEP.

2 .
FIG.2.Diagrams contributing to LQ double production at the LHC.

FIG. 5 .
FIG. 5. Expected 90% C.L. sensitivities, in the (m LQ , g) plane, obtained for the COH-CsI-700 (left panel) and COH-LAr-750 (right) detectors and assuming model S 1 .These bounds apply also to model R2 .The grey-shaded regions refer to current limits previously presented in Fig. 4. In the case of current CEνNS bounds, only the constraint obtained with the corresponding target is shown.See text for more details.

FIG. 7 .
FIG. 7. Expected 90% C.L. sensitivities, in the (m LQ , g) plane, obtained for the COH-CsI-700 (left panel) and COH-LAr-750 (right) detectors and assuming model S 3 .The grey-shaded regions refer to current limits previously presented in Fig. 4. In the case of current CEνNS bounds, only the constraint obtained with the corresponding target is shown.See text for more details.

FIG. 8 .
FIG. 8. Expected 90% C.L. sensitivities, in the (m LQ , g) plane, obtained for the Si (red contour) and Xe (yellow) detectors at the ESS, assuming different LQ models.The grey-shaded regions refer to current limits previously presented in Fig. 4. In the case of current CEνNS bounds, we show the combined COHERENT CsI+LAr result.See text for more details.

FIG. 9 .
FIG. 9. Left panel: Reduced χ 2 profile as a function of the LQ coupling g for two future detectors proposed at COHERENT, CsI-700 (light blue, dashed) and LAr-750 (blue).We fix m LQ = 100 GeV and we consider the S 1 model.Right panel: Same as left panel but for two future detectors proposed at the ESS, Si (red) and Xe (yellow, dashed).

TABLE I
. Relevant LQ operators connecting neutrinos and quarks.Left/right column shows the possible operators when LQs are scalars/vectors.The numbers between parenthesis denote the quantum numbers under SU (3) c , SU (2) L , and U (1) Y , respectively.

TABLE II .
Details of the CEνNS experiments considered in this paper.