Long-lived Sterile Neutrinos at the LHC in Effective Field Theory

We study the prospects of a displaced-vertex search of sterile neutrinos at the Large Hadron Collider (LHC) in the framework of the neutrino-extended Standard Model Effective Field Theory ($\nu$SMEFT). The production and decay of sterile neutrinos can proceed via the standard active-sterile neutrino mixing in the weak current, as well as through higher-dimensional operators arising from decoupled new physics. If sterile neutrinos are long-lived, their decay can lead to displaced vertices which can be reconstructed. We investigate the search sensitivities for the ATLAS/CMS detector, the future far-detector experiments: AL3X, ANUBIS, CODEX-b, FASER, MATHUSLA, and MoEDAL-MAPP, and at the proposed fixed-target experiment SHiP. We study scenarios where sterile neutrinos are predominantly produced via rare charm and bottom mesons decays through minimal mixing and/or dimension-six operators in the $\nu$SMEFT Lagrangian. We perform simulations to determine the potential reach of high-luminosity LHC experiments in probing the EFT operators, finding that these experiments are very competitive with other searches.


Introduction
Neutrino oscillations have proven without doubt that neutrinos are massive particles. The Standard Model of particle physics (SM) contains no right-handed neutrino fields. This forbids the generation of a neutrino mass via the Higgs mechanism, which generates the masses of the other elementary particles. This situation can be remedied by adding a sterile neutrino field to the SM [1][2][3][4][5]. The sterile neutrino, also called heavy neutral lepton (HNL), is a right-handed gaugesinglet spin-1/2 field and couples to left-handed neutrinos and the Higgs field through Yukawa interactions. This generates a Dirac neutrino mass after electroweak symmetry breaking.
In general, nothing forbids an additional Majorana mass term for the right-handed neutrino field, leading to Majorana mass eigenstates and lepton number violation (LNV). However, lepton number can be an (approximate) symmetry of extension beyond the SM (BSM), such that lowenergy LNV signals, e.g. neutrinoless double beta decay (0νββ), is suppressed. Sterile neutrinos may not only account for neutrino masses, but have also been linked to explanations of other problems of the SM. Light sterile neutrinos can account for dark matter [6][7][8][9], while sterile neutrinos with a broad range of masses can account for the baryon asymmetry of the Universe through leptogenesis [10]. Sterile neutrinos are thus a well-motivated solution to a number of major outstanding issues in particle physics and cosmology.
While the observation of neutrino masses provides a hint for the existence of sterile neutrinos, it does not specify their mass scale. They might very well be light and accessible in presentday and near future experiments. A large number of experimental and theoretical works have gone into the search for sterile neutrinos in so-called minimal scenarios, where sterile neutrinos only interact with SM fields through renormalizable Yukawa interactions (see Ref. [11] for a review). Here we take a more general approach. In broad classes of BSM models, sterile neutrinos appear sterile at lower energies, but interact at higher energies through the exchange of heavy BSM fields. Examples are left-right symmetric models [12][13][14], grand unified theories [15], or leptoquark models [16], which contain new fields that are heavy compared to the electroweak scale. Independent of the details of these models, at low energies the sterile neutrinos can be described in terms of local effective operators in the framework of the neutrino-extended Standard Model effective field theory (νSMEFT) [17,18].
In this work, we study relatively light GeV-scale sterile neutrinos (see e.g. Refs [19][20][21] for LHC searches for somewhat heavier neutrinos). Sterile neutrinos can then be efficiently produced via decays of mesons that are copiously produced at the LHC interaction points and if they are relatively long-lived, their decays lead to displaced vertices that can be reconstructed in LHC detectors. We consider a broad range of (proposed) LHC experiments: ATLAS [22]/CMS [23], CODEX-b [24], FASER [25,26], MATHUSLA [27][28][29], AL3X [30], ANUBIS [31], MoEDAL-MAPP [32,33], as well as the proposed CERN SPS experiment SHiP [34][35][36], and discuss their potential in probing νSMEFT operators. We calculate νSMEFT corrections to sterile neutrino production and decay processes and perform simulations for the various detectors, to estimate their search sensitivities. Our simulations show that the experimental reach is strong, probing dimension-six operators associated to BSM scales up to a hundred TeV. In a simple 3 + 1 model, adding just one sterile neutrino field, we compare our results to existing 0νββ decay limits, showing that these experiments are complementary.
This paper is organized as follows. In Sec. 2 we introduce the model framework of νSMEFT, followed by Secs. 3 and 4 detailing the calculation of production cross sections and decay widths of the sterile neutrinos in both the minimal model and from higher-dimensional operators. Sec. 5 shows the theoretical scenarios considered by numerical study in this work, and Sec. 6 goes through the different experiments we study in detail and briefly introduces the Monte-Carlo simulation procedure. In Sec. 7 we present the numerical results for both the minimal scenario and a number of flavor benchmarks. These results are compared with a number of other experimental probes including 0νββ decay in Sec. 8. In a set of appendices, we present the details of the production and decay rate computations, as well as the physical parameters, decay constants, and form factor input we employ. We conclude and provide an outlook in Sec. 9.

The Effective Neutrino Lagrangian
We are interested in the production and decay of sterile neutrinos at the LHC. In particular, we investigate the production of sterile neutrinos in the decay of mesons containing a single b or c quark, as these are copiously produced and are sufficiently massive to produce GeV sterile neutrinos. This is an interesting mass range that appears in scenarios of low-scale leptogenesis [37][38][39][40][41][42][43]. The sterile neutrinos are assumed to be singlets under the SM gauge group and, at the renormalizable level, only interact with SM fields via a Higgs Yukawa coupling to the lepton doublet. The renormalizable part of the Lagrangian is given by L SM denotes the SM Lagrangian, L is the lepton doublet, and H is the SM complex Higgs doublet field withH = iτ 2 H * . We work in the unitary gauge where v = 246 GeV is the Higgs vacuum expectation value of the Higgs real scalar h. ν R is an n × 1 column vector of n right-handed gauge-singlet neutrinos. Y ν is a 3 × n matrix of Yukawa couplings.M R is a complex symmetric n × n mass matrix. In general we can work in a basis where the charged leptons and all quarks are in their mass eigenstates, except the d i L , i = 1, 2, 3, for which we have d i L = V ij d j, mass L with V the CKM matrix. Ψ c is the charge conjugate field of Ψ with Ψ c = CΨ T and C is the charge conjugation matrix, C = −iγ 2 γ 0 , which satisfies the relation C = −C −1 = −C T = −C † . We define Ψ c L,R = (Ψ L,R ) c = CΨ L,R T = P R,L Ψ c , in terms of the projectors P R,L = (1 ± γ 5 )/2.  Table 1: νSMEFT dim-6 operators [18] involving one sterile neutrino field.
The Lagrangian in Eq. (1) can account for the observed active neutrino masses and mixing angles if n ≥ 2 (in case of n = 2 the lightest neutrino is massless [44]). As lepton number is explicitly violated by the Majorana masses,M R , Eq. (1) in general leads to Majorana neutrino mass eigenstates and thus to 0νββ decay and other LNV processes, unless additional structure is imposed on the matricesM R and Y ν .
In various popular extensions of the SM, right-handed neutrinos appear naturally, but are not completely sterile. For instance, in left-right symmetric models right-handed neutrinos are charged under a right-handed SU (2) R gauge group and interact with right-handed gauge bosons and new scalar fields. If such bosons exist, they must be heavy, with masses well above the electroweak scale, to avoid experimental constraints. The right-handed neutrinos on the other hand, can remain light. From this point of view, the scale separation suggests the use of an EFT framework where the degrees of freedom are the usual SM fields, as well as a set of n neutrinos, which are singlets under the SM gauge groups. The interactions in Eq. (1) form the dimensionfour and lower part of a more general Lagrangian containing higher-dimensional operators that is often referred to as the νSMEFT [17,45,46]. We begin by introducing the higher-dimensional operators at a scale Λ v, where Λ denotes the scale where we match the microscopic UV theory to the νSMEFT.
The first operators have dimension 5 At lower energies, after electroweak symmetry breaking these operators contribute to, respectively, Majorana mass terms for active and sterile neutrinos. The first operator, the famous Weinberg operator [47], can for instance be induced by integrating out sterile neutrinos with masses of O(Λ) usually referred to as a type-I seesaw mechanism. The second operator, for our purposes, can simply be absorbed inM R in Eq. (1). For n ≥ 2 there appears a dim-5 transition dipole operator, but we will not consider it here. We are mainly interested in the operators that appear at dimension-6 [18,48]. We focus on operators that involve a single right-handed neutrino 1 and limit the set of effective operators to those that lead to hadronic processes at tree level. A more general set of interactions is left for future work. The operators are presented in Table 1. For our purposes, the operator O (6) LνH can be absorbed in a shift in Y ν in Eq. (1). The related Higgs phenomenology was discussed in Ref. [49]. This leaves us with the remaining six operators that, in general, have arbitrary flavor indices.
We evolve the operators from Λ to the electroweak scale using one-loop QCD anomalous dimensions. The operators O (6) Hνe , O (6) νW , and O (6) duνe do not evolve under QCD at one loop. The remaining three operators evolve simply as [46] dC where C (6) S and C (6) T are defined as the linear combinations LνQd , C and where C F = (N 2 c − 1)/(2N c ) and N c = 3, the number of colors. At the electroweak scale, we integrate out the heavy SM fields (W-, Z-, and Higgs-boson and top quark) and match to a SU (3) c × U (1) em -invariant EFT. The tree-level matching relations for νSMEFT operators up to dimension-7 were given in Ref. [46] and here we use a subset of these results. We obtain ∆L=0 , where L SM contains dimension-four and lower operators involving light SM fields. L ∆L=0 includes dim-6 operators that conserve lepton number (∆L = 0) and is given by For the operators in Table 1 the Lagrangians L ∆L=2 and L ∆L=0 only contain a single term each where In these expressions we have suppressed flavor indices on the Wilson coefficients. Each Wilson coefficient carries indices ijkl where i, j = {1, 2, 3} indicate the generation of the involved up-type and down-type quarks, respectively, k = {1, 2, 3} the generation of the charged lepton (we will often use the labels e, µ, τ instead for clarity) and l = {1, 2, 3} for c The explicit matching relations are given by Ref. [46] for the mass terms in Eq. (7), and for the remaining operators. Here V denotes the CKM matrix, 1 the 3 × 3 identity matrix in lepton flavor space, and M e = diag(m e , m µ , m τ ) is the diagonal matrix of charged lepton masses. The first term in the expression for c (6) VL denotes the contribution from the SM weak interaction. All other entries arise from the dimension-6 operators, in some cases with additional insertions of leptonic mass matrices. The operators with Wilson coefficientsC (6) VL andc (7) VL are only induced by the dimension-6 operator O (6) νW . This operator involves a derivative acting on a charged W ± field which, after integrating out the W ± bosons, leads to operators involving an explicit derivative (c (7) VL ) or an insertion of a lepton mass by using the equations of motion. C (6) νW is strictly constrained because it generates neutrino dipole moments at one-loop [49,50]. Additionally, if these constraints are avoided N would decay relatively fast into two body final states via N → νγ [51,52]. To ensure that N is long-lived, we suppress O (6) νW in the following. This effectively implies we do not consider the effects ofC (6) VL andc (7) VL . Besides the charged currents listed in Eq. (8), we also include the effects of the SM weak neutral currents that contribute to decay processes of sterile neutrinos where i, j are the flavor indices of active neutrinos and θ w is the Weinberg angle.

Rotating to the Neutrino Mass Basis
After electroweak symmetry breaking the neutrino masses can be written as where M ν is an ×n symmetric matrix withn = 3 + n and N = (ν L , ν c R ) T . We use an ×n unitary matrix, U, to diagonalize the mass matrix and define N = U N m . In absence of sterile neutrinos, U is the usual PMNS matrix. We write the Majorana mass eigenstates as ν ≡ N m + N c m = ν c that appear in the Lagrangian as We introduce 3 ×n and n ×n projector matrices to express the relation between the neutrinos in the flavor and mass basis as In the mass basis, the operators in Eqs. (8)-(9) become TRR ν VLR i where C (6) VLL = c (6) VL P U +C (6) VL P s U , C VLR =c (6) VL P s U * , C (6) VRR =c (6) VR P s U * , C SRR =c (6) SR P s U * , C SLR =c (6) SL P s U * , C TRR =c T P s U * , C VLR =c Each Wilson coefficient again carries four flavor indices ijkl where i, j, k = {1, 2, 3} indicate the generation of the involved up quark, down quark, and charged lepton, respectively. l now denotes the particular neutrino mass eigenstate and runs from {1, . . . ,n}. Finally, we evolve the operators down to the bottom or charm mass. The vector currents do not evolve while the scalar and tensor dimension-6 and -7 couplings evolve in the same way as the scalar and tensor currents in Eq. (4), with TRR .
In what follows we only consider the dimension-six terms in Eq. (18) and neglect the dimensionseven operator proportional to C VLR whose low-energy effects are suppressed by m b,c /v. Furthermore, as discussed below Eq. (11), C VLR is only induced by the νSMEFT operator C (6) νW which is strongly constrained by other probes.

Production of Sterile Neutrinos
In this section we discuss the production of sterile neutrinos at collider and fixed-target experiments. For concreteness we consider the case where the sterile neutrinos are Majorana particles. We consider the production through the decay of mesons, produced at the interaction points, containing a single charm or bottom quark. We neglect the subdominant contribution from B c , J/Ψ, and Υ mesons although they would allow to probe a larger neutrino mass range. Production via the decay of pseudoscalar mesons dominates over the contribution from vector mesons of the same quark composition, due to the much shorter lifetime of the latter. Sterile neutrino production via direct decays of W -, Z-, and Higgs bosons is subdominant for O(GeV) neutrinos, mainly because of their smaller production cross sections [28,53,54].

Sterile Neutrino Production in Minimal Models
We begin by discussing sterile neutrino production in minimal models where sterile neutrinos interact with SM fields only via mixing. For simplicity, we consider a single sterile neutrino with mass m N and set n = 1 (n = 4). The production then arises solely from the first term in Eq. (18) with (C where V is the CKM matrix and U the lepton mixing matrix. A broad range of processes are relevant. Naively one might think that leptonic meson decays M ± ij → N + l ± k would dominate because of phase space suppression associated to semi-leptonic decays, but CKM factors and powers of meson/neutrino masses in the amplitude expressions change this picture. To calculate the production rate, we require the number of mesons produced at the various experiments and the branching ratio to final states including a sterile neutrino. The former is discussed below, while here we calculate the latter. For minimal models, these  branching ratios have been calculated in the literature, see Refs. [55,56] for recent discussions, and here we confirm (most of) these results. We consider leptonic and semi-leptonic decays of D ± , D 0 , D s , B ± , B 0 , and B s mesons. For semi-leptonic decays, we consider final-state pseudoscalar and vector mesons. The decay rate formulae and the associated decay constants and form factors are given in the Appendix. In the left and right panels of Fig. 1 we depict a selection of branching ratios for decay processes of D − , D s , and B − , B s mesons, respectively. Branching ratios for analogous decays of neutral D 0 or B 0 are similar and not shown to not clutter the plots too much. For these examples we considered a final-state electron and set U e4 = 1. All branching ratios are in excellent agreement with Ref. [55]. From the plots it is clear that, depending on the mass m N , both leptonic (solid lines) and semi-leptonic processes (dashed, dotted, and dot-dashed lines) must be included and the latter involve both final-state pseudoscalar and vector mesons.

Sterile Neutrino Production from Higher-Dimensional Operators
For higher-dimensional operators the quark flavor structure of the Wilson coefficients is unknown in contrast to the minimal case where the CKM matrix provides the relation between processes involving different quarks. As such, each flavor structure is independent unless model assumptions are used. This leads to a large number of possible cases corresponding to several flavor structures for each effective operator in Eq. (18). All branching ratios can be calculated from the expressions given in the appendices.
Here we discuss a few cases only. We consider the operators with Wilson coefficients C SRR , C TRR , and C (6) VLR . We consider the flavor structures {ijkl} = 13e4 and {ijkl} = 21e4 that allow for leptonic decays B → N + e and D → N + e, respectively, if the Lorentz structure permits this. These choices also allow for semi-leptonic decays of the form B → N + e + X and D → N + e + X where X is a pseudoscalar or vector meson consisting of just up, down, and strange quarks (strange quarks only if the decaying meson contains a strange quark as is the case for B s and D s mesons). For the plots in this section, we assume the weak interaction is turned off and consider only one non-zero EFT operator at a time. The results are depicted in Fig. 2.
The three chosen operators correspond to quark bilinears with different Lorentz structures (scalar, tensor, and vector, respectively). In the scalar case, leptonic decays are allowed and these dominate over semi-leptonic decay modes for all considered values of the sterile neutrino mass. For the tensor operator, however, the leptonic decay mode is forbidden and a final-state meson must be produced. In these cases, the dominant decay modes are those with a final-

CTRR=0.025
VLR vector operator has a similar Lorentz structure as the SM charged weak current, but with different flavor structure. As was the case in Fig. 1, depending on the sterile neutrino mass, either leptonic (solid lines) or semi-leptonic processes (dashed, dotted, and dot-dashed lines) can dominate the production of sterile neutrinos, and must all be included.

Sterile Neutrino Decays in Minimal Models
We begin by considering the minimal scenario, where we assume that the only non-zero term in Eq. (18) is (C For concreteness we consider decays of Majorana sterile neutrinos into final-state electrons and set U e4 = 1 and U µ4 = U τ 4 = 0. In addition, we consider the SM weak neutral current (see Eq. (12)), which leads to N → ν +f +f decays where f denotes any SM fermion that is kinematically allowed (in case of quarks, we consider a final-state neutral meson). These decay rates have all been calculated in the literature, see e.g. Refs. [55][56][57][58][59][60]. Most results agree with each other and with our findings given in the Appendix, with the exception for decay processes into final-state neutral mesons where some differences appear. For these cases, our results agree with Ref. [56].
We consider N → leptons through both charged and neutral weak currents. The latter leads to the invisible three-neutrino decay mode. We include decays into a single pseudoscalar (π, K, η, η , D, D s , η c ) and vector meson (ρ, ω, K * , φ, D * , D * s , J/Ψ). This effectively also takes into account decay modes into two pions through the intermediate decays of ρ mesons [55], assuming m N > m ρ . For heavier sterile neutrinos other multi-meson final states become relevant and summing exclusive channels becomes impractical. Instead we follow Refs. [55] and estimate the total hadronic decay width by calculating the decay width to spectator quarks times appropriate loop corrections. The loop corrections are taken from a comparison to hadronic τ decays where D denotes a d or s quark and where dots denote higher-order corrections. This gives a good description of the inclusive hadronic τ decay rate and we assume this to hold for sterile neutrino decays in the minimal scenario as well. That is, we use to calculate the inclusive hadronic sterile neutrino decay rate through both charged and neutral weak currents. We find that single meson channels dominate for m N 1 GeV, while the decay to quarks become relevant for larger masses, indicating that multi-meson final states become significant. We demonstrate this in the left plot of Fig. 3, which shows branching ratios to individual mesons, compared to the sum of all single-meson final states, and compared to quarks, for m N 1 GeV. At m N = 5 GeV, the single meson final states make up roughly 20% of the hadronic decay rate. Our results are in good agreement with Ref. [55], apart from decays to neutral vector mesons, which only play a small role.
We write the total decay rate as In the right panel of Fig. 3, we show a plot of the scaled proper decay length, U 2 e4 cτ N , as a function of m N in the minimal scenario, where c is the speed of light and τ N is the proper lifetime of N . The branching ratios to individual mesons, leptons, and three neutrinos (invisible) for m N < 1 GeV are shown in the left panel of Fig. 4 while the branching ratios to quarks, leptons, and invisible for m N > 1 GeV are shown in the right panel of the same figure.

Sterile Neutrino Decays from Higher-Dimensional Operators
The EFT operators in Eq. (18) involve two quarks and a charged lepton and induce new channels for sterile neutrinos to decay into hadrons. Depending on the flavor structure of the operators, typically only one or two single-meson decay modes are relevant. For instance for an operator (C VLR ) ije4 with ij = 11 we consider decays into a single charged pion or ρ meson, while for (C (6) SRR ) 11e4 only pions are relevant. As is the case for the minimal scenario, one can imagine multi-meson states to become relevant for larger sterile neutrino masses (for this flavor choice, such states would be three-or more pions). We do not include such states here, although we find that decays to quarks become dominant around m N 2 GeV for scalar operators, as we do not have the benchmark of hadronic τ decays to verify our results for non-SM currents. We only consider decays into individual mesons. This leads to a potential underestimate of the sterile neutrino decay width and consequently renders our sensitivity limits for future experiments conservative in the large decay length regime. As all our results below are given on log scales, we do not expect significant deviations from our findings. In Fig. 5 we plot the proper decay length of N , cτ N , against its mass for various choices of Wilson coefficients and flavor assignments. We consider one effective coupling at the time, turn off minimal mixing, and the Wilson coefficients are set to unity. (C SRR (6) ) 11 e4 (C VLR (6) ) 21 e4

Theoretical Scenarios
In this work we focus on the production of sterile neutrinos through the decays of B and D mesons. We consider three classes of scenarios that are representative of the effective Lagrangian in Eq. (18). We always consider the case of a single sterile neutrino which is mixed with the active electron neutrino. The three classes of scenarios are listed below: 1. Here we consider the minimal scenario without higher-dimensional operators and 1 sterile neutrino (a 3 + 1 model). Because of minimal sterile-active mixing, the sterile neutrino only interacts through charged and neutral SM weak interactions. In the 3 + 1 minimal seesaw model, the mixing angle is related to the ratio of active and sterile neutrino masses |U e4 | 2 ∼ m ν /m N , but we treat the mixing angle as a free parameter. We stress that the minimal 3+1 model leads to two massless active neutrinos and is thus ruled out by neutrino oscillation experiments, but with its simplicity it provides a useful benchmark.
2. In this scenario we extend the minimal 3 + 1 model by interactions generated by the exchange of leptoquarks. All possible representations of LQs are summarized in Ref. [16]. We focus on the representationR (3, 2, 1/6), which can couple to (sterile) neutrinos through the Lagrangian where a, b are SU (2) indices and i, j, k, l are flavor indices, respectively. Note that l = 1 in the 3 + 1 model. LHC constraints force the leptoquark mass, m LQ , to be above a few TeV and for low-energy purposes we can integrate it out. At tree level this leads to the effective operator where We read from Eq. (11)c LdQν .
Finally, going to the neutrino mass basis and focusing on the couplings to electrons and sterile neutrinos, we obtain the matching contributions to the effective operators in Eq. (18) where i and j denote the up-and down-quark generation, respectively, and V ij elements of the CKM matrix. In this scenario, we have contributions from minimal mixing proportional to the mixing angle U e4 and from leptoquark interactions proportional to U * 44 . We will use the canonical see-saw relations and set m ν = 0.05 eV as representative for the active neutrino masses. For a specific quark flavor choice of i and j, this reduces the effective number of free parameters to two: the sterile neutrino mass m N and the combination of the LQ couplings and mass 3. The final scenario we consider is inspired by models, such as left-right symmetric models, with right-handed charged gauge bosons, which can mix with W ± . Instead of implementing the full left-right symmetric model, we take a simplified scenario and consider the effects of a nonzero C VLR in Eq. (18) in combination with the SM left-handed weak interactions. That is, we consider In left-right symmetric scenarios nonzero C VRR and C VRL are generally induced as well. The resulting phenomenology is very similar to C (6) VLL and C (6) VLR and for simplicity we do not consider these effective operators here. They can be easily added to the analysis if so required. For the active-sterile neutrino mixing parameter, we follow the leptoquark scenario and set U e4 = m ν /m N . In minimal left-right symmetric models with exact P or C symmetry, the dependence of the effective operators on the quark flavor indices ij can be calculated. We do not consider this here and consider one particular choice of ij at a time. It is straightforward to generalize this choice.

Collider and Fixed-target Analysis
We proceed to explore the search possibilities for the above scenarios at the LHC, considering both existing and proposed experiments, as well as the proposed fixed-target experiment SHiP at the CERN SPS [34][35][36]. Currently at the LHC we have the operational experiments ALICE [61,62], ATLAS [22], CMS [23], and LHCb [63,64]. Of these we shall focus on the search sensitivity for longlived sterile neutrinos at ATLAS, which is the largest experiment in detector volume, and can thus in principle explore the largest decay lengths. Beyond this, a series of new experiments has recently been proposed at various locations near the LHC interaction points (IPs), namely in alphabetical order: AL3X [30], ANUBIS [31], CODEX-b [24], FASER and FASER2 [25,26], MATHUSLA [27][28][29], and MoEDAL-MAPP1 and MoEDAL-MAPP2 of the MoEDAL collaboration [32,33]. These latter experiments, including SHiP, are all designed to specifically look for neutral long-lived particles (LLPs). We shall discuss the search potential of all the above-mentioned experiments specifically for sterile neutrinos.
In this section, we review briefly the setup of the beam and detector geometries for each experiment and introduce the Monte Carlo (MC) simulation procedure for the event simulation. We focus on the production via the rare decay of B-and D-mesons. This can proceed via either purely leptonic two-body decays, or semi-leptonic three-body decays, as discussed in Sec. 3. Similarly, for the displaced decay of sterile neutrinos, we consider both the two-body decay into a charged lepton and a charged meson, and decays into multiple hadronic states plus a lepton, as explained in Sec. 4. It is worth mentioning that the heavy meson M 1 that is to decay into a sterile neutrino, is produced in the process pp → M 1 + X and decays promptly into e.g. a sterile neutrino (X denotes the remaining decay products). Such signal events can be observed in the near detector such as ATLAS with e.g. the prompt lepton accompanying the production of M 1 . The long-lived sterile neutrinos decay at a macroscopic distance, where the displaced vertex (DV) is reconstructed if at least two tracks are observed stemming from the same DV inside the detector.
We do not study the production of the sterile neutrinos from the decay of lighter mesons such as π ± and kaons as these are only relevant for very light sterile neutrinos, and their simulation in Pythia 8 [65,66] is insufficiently validated in the forward direction, which is relevant for the FASER experiments. In fact, even for D ± and B ± mesons, we will use FONLL [67][68][69][70] to correct the behavior of Pythia 8 in the very large pseudorapidity regime. Finally, we ignore the vector mesons decays into sterile neutrinos, as their decay width is typically many orders of magnitude larger than that of pseudoscalar mesons leading to tiny branching ratios for sterile neutrino production.
Instead of performing a detailed study considering different components of the detectors separately, for simplicity we will take the whole detector as the fiducial volume, and make a comparison between the various experiments. Since the ATLAS detector was not designed to look for neutral LLPs, we shall add an estimate for the efficiency factor based on existing neutral LLP searches, which, however, search for heavier candidates than considered here. One potential issue relates to possible background events which, depending on the placement of the detector, may consist of long-lived SM hadron decays, cosmic rays, hadronic interaction with the detector material, etc. All the experiments considered in this study, except ATLAS, employ a far detector with a distance 5 − 500 meters away from the IP. The space between the IP and the far detector is usually sufficiently large to allow for the installation of veto and shielding segments, as argued in Refs. [24,25,27,28,30,31,33,71,72]. For MATHUSLA the rock and shielding below ground should remove the SM background. As for cosmic rays, directional cuts will be applied. To assess and compare the sensitivities of different experiments, we show 3-event isocurves which correspond to 95% confidence level (C.L.) with zero background. This is not appropriate for the ATLAS detector which has an almost 4π coverage immediately around the IP, and a large irreducible SM background is expected. Depending on the signal type, the number  Table 2: Summary of integrated luminosities for the various experiments. "POT" for SHiP stands for "Protons on Target". We also list the simple geometric coverage for each experiment.
of such background events may vary. Instead of performing an estimate of background events for each scenario, we shall implement an estimate for the efficiency, which we discuss below.
In addition, since the detector efficiency of the future experiments is unknown, in order to make a fair comparison, we assume a 100% reconstruction efficiency for all the experiments except ATLAS.
The search potential of these experiments, but also of other fixed-target experiments, has been investigated for example for neutralinos [73][74][75][76][77]. For various other theoretical scenarios, see Ref. [78] for a recent review.
We start with a brief introduction of the fixed-target experiment SHiP, as its beam setup is different from the other experiments, and then discuss the other experiments, which are all associated to either the ATLAS, CMS, or LHCb IP and the LHC accelerator 2 . Since these experiments differ in the projected integrated luminosity, we summarize them in Table 2.

SHiP
The SHiP facility was proposed to make use of the high-intensity CERN SPS beam of 400 GeV protons incident on a fixed target made of e.g. a hybrid material composed of (solid) molybdenum alloy and pure tungsten [34][35][36]. It has not been approved yet. With a center-of-mass energy of approximately 27 GeV, large production rates of D-and B-mesons are expected. The SHiP experiment is proposed to have a cylindrical detector downstream at roughly 70 m away from the IP. The experiment is specifically designed for detecting long-lived neutral particles, which are produced from e.g. charm or bottom meson rare decays, fly an extended length, and then decay inside the detector chamber downstream, especially if the lab-frame decay length of the LLP lies within the SHiP sensitivity range. For the lifetime of the SHiP project, a total of 2 × 10 20 protons on target (POT) are planned.
At SHiP, the initial meson M 1 of sterile neutrinos is produced in a hadronic collision between the beam protons and the target material. The differential production cross section is strongly forward peaked, and the M 1 will have a significant forward boost. This will be passed onto the decay products, including N , the sterile neutrino. The active decay chamber is 68.8 m downstream of the target. It has a cyclindrical shape with a length of 60 m, where the first 5 m are to be used for placing background suppression vetoes. The front surface has an elliptical shape with semi-axes of 5 m and 2.5 m. The optimal sensitivity is for particles with where β z N , γ N are the relativistic speed along the beam axis and the Lorentz boost factor of N , respectively.
In order to study sterile neutrinos produced from the decays of D-and B-mesons, we need to know the total number of these mesons expected at SHiP in its 5 year lifetime: N D ± , N D 0 , N Ds , N B ± , N B 0 , and N Bs , respectively. These numbers can be estimated by following Ref. [55]. See also the earlier work in Ref. [75]. The cc (bb) production rate is 1.7 × 10 −3 (1.6 × 10 −7 ) per collision. After the fragmentation factors are taken into account, the numbers of the heavy-flavor mesons can be estimated and reproduced below from Ref. [55]: We note that B c mesons are in principle also produced and may extend the upper reach in m N . However, given the much smaller production cross section, we do not take it into account. Similarly for the other LHC experiments, as discussed below.

Experiments at the LHC
For ATLAS and extended programs at the ATLAS/CMS/LHCb sites, we consider pp collisions at √ s = 14 TeV with equal beam energies. These experiments benefit from a beam energy that is orders of magnitude higher compared to SHiP. As a result, the mesons and the therefrom produced sterile neutrinos are much more boosted, leading to good sensitivities even for detectors that are to be installed hundreds of meters away from the IP.
To perform the sensitivity estimate for experiments at ATLAS/CMS/LHCb, it is necessary to know the inclusive production rate of the heavy-flavor mesons at the HL-LHC with up to 3 ab −1 integrated luminosity. To achieve this, we follow the procedure in Refs. [75,76,79]. The LHCb collaborations reported D-meson and B-meson production cross sections for a 13 TeV pp-collider, for certain kinematic range. Using the simulation tools FONLL [67][68][69][70] and Pythia 8 we can extrapolate these cross sections to the whole kinematic range. We find that for 3 ab −1 integrated luminosity over the full 4π solid angle the following list of numbers of the produced charm and bottom mesons: For the estimate of N HL-LHC D ± , the decay branching ratio of D * ± into D ± is included, and for the number of B-mesons the corresponding fragmentation factors determined by Pythia 8 are used. These numbers will be used not only for the evaluation at the ATLAS experiment, but also for all the extended programs discussed below with a possible overall re-scaling by the integrated luminosity.

FASER
At the ATLAS IP, the differential cross section for the production of GeV-scale mesons is strongly peaked at large pseudorapidities, i.e. close to the beam pipe in both directions. The production rate of light neutral LLPs, resulting from the decay of the mesons, should then also be peaked in the large pseudorapidity regime. A far detector known as FASER (Forward Search ExpeRiment) [25,26] has been proposed, which is designed to specifically make use of this feature. It has been officially approved by CERN and should be collecting data during Run 3 of the LHC. 3 It is placed in the existing TI12 tunnel at a distance of 480 m from the ATLAS IP and has a cylindrical shape exactly aligned with the beam collision axis, but slightly off of the beam due to the curvature of the accelerator. The FASER detector has a cylindrical radius of 10 cm and a length of 1.5 m, and the expected integrated luminosity is 150 fb −1 . The corresponding angular coverage is η ∈ [9.17, +∞] in pseudorapidity and full 2π in the azimuthal angle.
After Run 3 is finished, FASER is currently planned to be upgraded to a larger version known as FASER2, to be under operation during the HL-LHC era, collecting up to 3 ab −1 of data [26], and located in the same place. Its geometry is specified as a cylinder with 1 m radius and 5 m length, which is 300x larger by volume, than FASER. The pseudorapidity coverage is correspondingly enlarged to η ∈ [6.86, +∞] while the azimuthal angle remains fully covered. In this work, we will study the search potential of both FASER and FASER2 for sterile neutrinos as neutral LLPs.
As mentioned earlier and discussed in Ref. [76], Pythia 8 is not well validated in the large pseudorapidity regime for the differential production cross section of charm and bottom mesons. To solve this issue, we re-scale the meson production cross section in Pythia 8 at different ranges of the transverse momentum and pseudorapidity by using the more reliable results given by FONLL.

MATHUSLA
A much larger experiment called "MATHUSLA" (MAssive Timing Hodoscope for Ultra Stable neutraL pArticles) has been suggested to be constructed at the surface above the CMS experiment [27][28][29]. 4 It is proposed to be built for the HL-LHC phase with 3 ab −1 integrated luminosity. With a box shape, it has a base area of 100 m× 100 m and a height of 25 m. Relative to the CMS IP, the front edge of the detector should be horizontally shifted along the beam axis by 68 m, and vertically upwards by 60 m. Despite its huge size, the large distance of MATHUSLA from the IP still leads to a small geometric coverage. It corresponds to a solid angle or geometric coverage of about 3.8%. 5 Nevertheless, it has been shown to be one of the far detectors at the LHC that may have the strongest reach in the very small active-sterile neutrino mixing at |U e4 | 2 ∼ 10 −9 , when only the weak interaction is considered [29].

ANUBIS
It has recently been proposed [31] to construct a detector, named ANUBIS (AN Underground Belayed In-Shaft search experiment), in one of the service shafts just above the ATLAS or CMS IP. Consequently, a total of 3 ab −1 integrated luminosity from the LHC is projected. Similarly to the detectors discussed above, it also has a cylindrical shape however with the axis oriented in the vertical direction. The cyclinder diameter is 18 m and the length is 56 m. The axis of the cylinder runs from the middle point of the bottom: (x b , y b , z b ) = (0, 24 m, 14 m) to the top (x t , y t , z t ) = (0, 80 m, 14 m), where z is along the beam axis, x is horizontally transverse, and y is vertically transverse, and the IP is at (x, y, z) = (0, 0, 0). Please refer to Refs. [31,54] for a sketch. MC integration finds the angular coverage of ANUBIS to be 1.79% [79].

CODEX-b
An external detector extension has also been proposed at the IP8 of the LHCb experiment. CODEX-b (COmpact Detector for EXotics at LHCb) [24] is designed as a cubic box of dimensions 10 m × 10 m × 10 m. Occupying an empty space with a distance ∼25 m from the LHCb IP, it covers the pseudorapidity range of η ∈ [0.2, 0.6] with the azimuthal angle coverage of ∼ 6.4%. The total geometric coverage of the solid-angle is about 1%.

MoEDAL-MAPP
MAPP (MoEDAL's Apparatus for Penetrating Particles) is proposed as one sub-detector of the MoEDAL experiment [32,33] at the IP8 of LHCb, and designed for searching for neutral LLPs. Similar to FASER, the MAPP experiment will have a significant upgrade in a second version. MAPP1 is a rather small detector of volume ∼ 130 m 3 , currently under deployment and expected to collect data during LHC Run 3 with up to 30 fb −1 integrated luminosity. It is roughly 55 m from the IP8 at a polar angle 5 • from the beam. During the HL-LHC era, the MAPP2 program is planned to be under operation at IP8 until the end of Run 5, accumulating up to 300 fb −1 of data. It is designed to occupy almost the whole of the UGC8 gallery in the LHC tunnel complex, taking up a volume of about 430 m 3 . With a larger integrated luminosity and a bigger volume, MAPP2 is predicted to have higher sensitivity. The two detectors cover about 0.17% and 0.68% of the total solid angle [79], respectively.

AL3X
AL3X (A Laboratory for Long-Lived eXotics) was proposed in Ref. [30] to be built at the LHC IP2 where the ALICE experiment sits. Placed at a horizontal distance along the beam line of 5.25 m from the IP, the detector has a cylindrical shape aligned with and surrounding the beam line, corresponding to a full azimuthal coverage. The proposed inner (outer) radius is 0.85 m (5 m) with the length 12 m. This gives a pseudo-rapidity coverage of η ∈ [0.9, 3.7]. The authors of Ref. [30] proposed two values of the integrated luminosity, 100 fb −1 and 250 fb −1 , in order to accommodate practical concerns including the move of the IP, beam quality, and investigation of background events. Here we focus on the benchmark value of 250 fb −1 integrated luminosity.

ATLAS
In the previous sections we have discussed proposed new experiments which are specifically designed to look for long-lived particles. They are typically placed some distance away from the various IPs at the LHC or, in the case of SHiP, designed as a fixed-target experiment with a long decay path. In e.g. Ref. [75] some of us considered the search for light long-lived particles at ATLAS. 6 ATLAS can study decays upto a length of about 11 2 + (43/2) 2 = 24 m, where 11 m is the cylindrical radius and 24 m is half the detector length. Over this length, as we discuss below in more detail, the fraction 1 − exp[24 m/(βγcτ )] of the LLP's decay where τ and γ denote the LLP lifetime and boost factor, respectively, and β is the relativistic speed of the particle. However, ATLAS offers almost 4π in angular coverage, which is significantly larger than the other detectors at the LHC. In Ref. [75], we found that for 100% signal efficiency, an integrated luminosity of 250 fb −1 , and assuming zero background, ATLAS is competitive with or slightly better than SHiP for the LLPs produced from B-mesons decays, and was somewhat worse in the case of D-mesons.
Here we describe in more detail the geometrical parameters we use for the fiducial volume of the ATLAS detector. It has a cylindrical shape with inner (outer) radius of 0.0505 m (11 m) corresponding to the beginning of the inner detector (the end of the muon spectrometer), and a length of 43 m. In principle, even if a DV is inside the beam pipe, as long as its distance from the IP is larger than the detector spatial resolution and consists of displaced tracks, it can be reconstructed. However, to be more conservative, we choose to include only DVs that are located inside the detector volume. We take 3 ab −1 as the benchmark value for the integrated luminosity over the lifetime fo the HL-LHC.
As we have mentioned, all the above proposed new experiments are specifically designed to look for light long-lived particles. They are thus further away from the LHC IPs and shielded from many SM backgrounds generated at the IP. All the same they will all have separate background issues, depending on where they are located. MATHUSLA is to be installed on the ground and thus highly susceptible to cosmic rays. ANUBIS is in a shaft which has minimal overhead shielding. FASER is far down along the tunnel, but still close to the beam pipe. In order to compare these experiments we have for now assumed that they can tackle the issues concerning background, and assumed zero background for all. The precise design of these detectors is also not yet well-established, except for the first phase of FASER. We thus furthermore assume 100% signal efficiency.
All this does not hold for ATLAS (or CMS, of course). ATLAS is a well-established experiment, operating in the immediate vicinity of the IP. There is purposely no shielding, as a priori all events are of interest. With respect to light long-lived particles there are thus large backgrounds which must be dealt with. Any cuts which are imposed to reduce the background, will also affect the signal efficiency, most likely in a significant manner. In order to compare a search at ATLAS with the other experiments we must impose a realistic signal efficiency, to take this into account. To our knowledge at the moment, there is no dedicated study at ATLAS or CMS for the scenarios we are considering here. We discuss several related analyses and estimate a signal efficiency based on this.
In Ref. [81], ATLAS considered the following scenario proposed in Ref. [82]. A Higgs boson decays to two dark fermions, which in turn decay to one or two dark photons plus an invisible hidden lightest stable particle. The dark photons are the long-lived particles and decay either to a pair of muons, or a pair of electrons/pions (light hadrons). Dark photon masses between 0.4 and 3.5 GeV are considered, i.e. similar to our sterile neutrino mass range. However the Higgs boson masses are 125 or 800 GeV and thus the dark photons would have a much higher boost than our sterile neutrino's. For the lighter Higgs mass and the purely hadronic decay of the dark photon ATLAS finds a signal efficiency of at best 5·10 −5 for cτ ≈ 35 mm, and dropping off rapidly for smaller or larger cτ . In Ref. [83], ATLAS considered the following scenario. A Higgs boson decays to a pair of neutral long-lived scalars, which in turn each decay to 2 jets, one in the inner detector and one in the muon spectrometer, thus utilizing different parts of the ATLAS detector. For the lightest mass scenario the Higgs boson is 125 GeV and the scalar is 8 GeV. A detector efficiency of 3·10 −5 is found, similar to the other analysis. In Ref. [84] ATLAS considered a related scenario with scalar masses down to 5 GeV. For a Higgs boson of 125 GeV, and a scalar of 5 GeV with a decay length of 75 cm they find a signal efficiency of 7·10 −4 , somewhat better than in the other studies.
In all these analyses, ATLAS searches for pairs of produced particles. This reduces background but also efficiency. Overall we have applied a from the ATLAS point of view somewhat optimistic flat signal efficiency factor of 10 −3 to all the ATLAS searches. In addition we assume that the corresponding cuts reduces the background to zero.

Monte-Carlo Simulation
We perform the MC simulation with the tool Pythia 8.243 [65,66], in order to extract the kinematics and to estimate the number of signal events. We express the number of sterile neutrinos, N , produced as where M i is summed over all mesons that can decay to N in a given scenario. "Br" stands for decay branching ratio. N M i is the number of initial mesons, M i , produced in the initial collisions at the LHC or for SHiP. The number of sterile neutrino decays in a detector volume N dec N can then be estimated by taking into account the boost factor and traveling direction of N , geometries of the detector, and the decay branching ratio of sterile neutrinos into the signal final-state particles. For the latter, we consider all the decay channels except for the fully invisible state, which contains solely three neutrinos, and is mediated by the SM weak interaction. We use the expressions where P [N in f.v.] denotes the average probability of all the simulated sterile neutrinos to decay inside the fiducial volume ("f.v.") and P [N i in f.v.] is the individual probability of the i−th simulated sterile neutrino to decay in the f.v., discussed in more detail below. N MC N is the total number of sterile neutrinos N generated in the simulation. For all the experiments except MAPP1 and MAPP2, we use formulas for calculating the individual decay probability with the exponential decay distribution, extracted from existing references [54,[75][76][77]. As input we use the boosted decay length and the traveling direction of each simulated sterile neutrino, denoted by λ i = β i γ i c τ N , where β i is the speed and γ i the boost factor.
For the MoEDAL-MAPP detectors, following Ref. [79] we implement a code which determines whether each simulated sterile neutrino travels in the direction pointing towards the detector and if so returns L 1i and L 2i , where L 1i denotes the distance from the IP to the position where the i−th sterile neutrino would enter the detector, and L i2 the distance the i−th sterile neutrino would travel across the detector, if it leaves the detector without having decayed. If the travel direction of the long-lived sterile neutrino points towards the detector, we compute P [N i in f.v.] for the MoEDAL-MAPP detectors through We use the modules "HardQCD:hardccbar" and "HardQCD:hardbbbar" of Pythia 8 to simulate the production of the charm and bottom mesons, respectively, including the processes qq, gg → cc/bb. For each benchmark scenario, we simulate 10 6 events of the corresponding process, and fix the initial-state mesons to decay to the various channels, mediated by both the SM weak interaction and the EFT operators, with the corresponding decay branching ratios. From the MC simulation with Pythia 8, we obtain the average decay probability from which we calculate N dec N in Eq. (38). We emphasize that for the partial decay widths of the heavy mesons into N and of the N into light mesons, we take into account the weak interaction via active-sterile neutrino mixing, the EFT operator(s), and also the interference between them. The computation for the production and decay processes of sterile neutrinos are presented in Secs. 3 and 4 and in the appendices.

Numerical Results
To present the results of this collider study in a compact and representative manner we consider a subset of possible EFT operators and focus on the three scenarios described in Sec. 5. For scenario 1, the minimal scenario, there is no EFT operator involved, and we consider the full standard model charged-and neutral-currents mediated via the active-sterile neutrino mixing. For scenarios 2 and 3, which involve EFT operators, we must specify their flavor. In these scenarios we consider 5 different flavor assignments. Each assignment involves two EFT operators of different quark flavors. All EFT operators are dimension-6, and we drop the (6) superscript in the following. First, we fix the production mode via the "production operator" (C P ) ij , with associated Wilson coefficients. We shall consider two cases for the indices ij, which indicate the up-and down-type quark generations considered in a flavor benchmark, respectively.
Here X indicates a potential final state meson. See examples below. Second, we simultaneously turn on another EFT operator 7 the "decay operator" (C D ) ij , which leads to the decay of sterile neutrinos via semi-leptonic processes. Here, the ij indicate the flavor content of the final state meson. We shall consider several decay modes Sterile Neutrino Decay Modes : 7 It is possible to have only one non-vanishing operator e.g. C 11 SRR = 4C 11 TRR or C 11 VRR , leading to ρ ± decaying to N + e ± and N decaying to π ± + e ∓ . However, such scenarios probe only a small sterile neutrino mass range and in addition require the decaying meson to be a vector particle. This results in a too small sensitivity reach because of the small production rate and large decay width of the vector mesons. Both the production and decay will be induced by various operators with associated Wilson coefficients, which we discuss in detail below. For all scenarios we include the production and decay of sterile neutrinos via the SM weak interaction (see Sect. 5 for details) and the corresponding interference terms with the EFT operators. For the theoretical scenarios 2 and 3 we impose the type-I Seesaw relation to include weak interaction contributions through minimal mixing, see Sect. 5.
As indicated in Eq. (41), we only include sterile neutrino production via B and D decays, also for lighter sterile neutrinos with masses m N < m K , m π , the kaon and pion masses, respectively. We do not include possible production via kaon or pion decays for the following reasons. Despite the larger production rates of these mesons, the simulation of soft pions is not well validated in quick MC simulation tools. Furthermore, the kaons are long-lived, leading to further complications. Finally, sterile neutrinos produced from pions and kaons are necessarily light, resulting in limited sensitivity reach in m N . To summarize, while we show results for light sterile neutrinos in the following, these must be taken as conservative as we underestimate the number of the produced sterile neutrinos from pion and kaon decays both via the SM weak interaction as well as via the "decay operator" C D .

The Minimal Scenario
In the minimal scenario, the interactions are purely mediated by the W -and Z-bosons via the active-sterile neutrino mixing. Using the analytical expressions given in the previous sections for the production and decay of the sterile neutrino, we estimate the sensitivity reach of the experiments discussed in Sec. 6. We present the results in Fig. 6, shown in the plane |U e4 | 2 vs. m N . Here, we lift the requirement of the type-I seesaw relation |U e4 | 2 m νe /m N , and treat the mixing angle and the sterile neutrino mass as two independent free parameters. The light gray area shows the present bounds obtained by various experiments including the searches from CHARM [85], PS191 [86], JINR [87], and DELPHI [88]. The dark gray area corresponds to the part excluded by big bang nucleosynthesis (BBN) [89,90]. We also show a brown band of "Type-I Seesaw target region" for m νe between 0.05 eV and 0.12 eV with the relation |U e4 | 2 m νe /m N . These two limits are derived from neutrino oscillation and cosmological observations, respectively. The former finds that there is at least one active neutrino mass eigenstate of mass at least 0.05 eV [91] while the latter imposed an upper limit of 0.12 eV for the sum of the active neutrino masses [92].
The sensitivity reaches of the various experiments have been determined in the literature [26,28,53,54,77,93], except for the MAPP1 and MAPP2 experiments, for which we present the estimate for the sensitivity reach for the first time. Our estimate for the ATLAS experiment considers an integrated luminosity of 3 ab −1 and takes 3000 signal events before taking into account an universal efficiency factor 10 −3 as the 95% C.L. exclusion limit. To the best of our knowledge, a similar estimate for sterile neutrinos in the minimal scenario produced from heavy-flavor mesons rare decays has not been conducted for ATLAS. See Ref. [94] for a related study within the minimal model, with the sterile neutrino produced from W -decays, however with promptly decaying sterile neutrinos. Furthermore, in Ref. [95] ATLAS investigated a sterile neutrino mixing with ν µ and with a delayed decay, i.e. a displaced vertex, as well as a promptly decaying sterile neutrino, which mixes with ν e however only for m N 6 GeV. For the other experiments, we assume zero background and 100% detector efficiency. Hence we take 3 signal events as the 95% C.L. exclusion limit. Comparing the sensitivity reach of the experiments shown in Fig. 6 with those from the literature, we find a good agreement in most cases. For the ANUBIS exclusion limits, our results shown in Fig. 6 are inferior by a factor ∼ 3.5, to those given in Ref. [54]. This difference is due to a corrected meson-production-rate and sterile-neutrinodecay-width calculation.
Given the general agreement with the existing results in the literature, we proceed to evaluate the sensitivities of these experiments to a set of benchmark scenarios, where the sterile neutrino interactions with the SM particles are enhanced by heavy new physics, encoded by EFT operators. There are a large number of possibilities for the flavor structure of the EFT operators. Here we consider a few representative flavor choices, to get an understanding of the general features and the sensitivity reach of various experiments. The calculations can easily be repeated for other flavor choices, for instance those inspired from specific UV-complete scenarios.
We note that in the following EFT flavor benchmarks, with the choice of the canonical type-I Seesaw relation, m N 1 GeV appears to be disfavored by BBN considerations. However, the inclusion of the EFT operator C D can reduce the lifetime of the sterile neutrinos, possibly circumventing BBN constraint and a potential a lower bound on the EFT Wilson coefficients. Different flavor assignments result in different final-state particles, affecting the primordial helium and deuterium abundances to different extents. A detailed study is necessary to investigate the limits that can be set from BBN consideration on EFT operators and we do not present BBN exclusion bands below.

Flavor Benchmark 1
We consider the theoretical scenarios 2 and 3 of Sec. 5 for a specific flavor choice. We focus on sterile neutrino production through decay of D-mesons, and thus consider the production operator (C P ) 21 . For scenario 2, the leptoquark scenario, for C P and C D we consider scalar and tensor operators that are related through C SRR = 4C TRR . For scenario 3 corresponding to anomalous vector interactions, the production and decay operators are both C VLR . The following Flavor benchmark 1. production process via C P D ± /D 0 /D s → N + e ± (+X) . decay process via C D N → π ± + e ∓ , ρ ± + e ∓ Table 3: Summary of flavor benchmark 1 for the theoretical scenarios 2 and 3 of Sec. 5, respectively. X denotes any additional final state particles.
decays then produce sterile neutrinos We are sensitive to sterile neutrino masses m N < m D − m e . For the decay operator we choose (C D ) 11 resulting in the decays N → e ± + π ∓ , and N → e ± + ρ ∓ .
The essential features of this benchmark are summarized in Table 3. In addition, we include production and decay modes via minimal mixing which extends both the upper and lower reach in m N . Fig. 7 presents the sensitivity reach of all considered experiments for the flavor benchmark 1. The left panels correspond to theory scenario 2 (leptoquark) and the right panels to scenario 3 (anomalous vector interactions). In the upper row we plot the value of the decay operator C D vs. the production operator C P for a fixed sterile neutrino mass m N = 1 GeV. In the bottom row we have fixed C D = C P and show the dependence of the sensitivity of the experiments on C D = C P and the sterile neutrino mass. Both the top and bottom panels show that the sensitivity reach in scenarios 2 and 3 are rather similar, indicating that the specific Lorentz structure of the EFT interactions does not greatly affect the overall sensitivity. In the upper-left and lower-right part of the top plots the curves become horizontal and vertical, respectively. In this part of the parameter space either the production (horizontal, upper left) or decay (vertical, lower right) of sterile neutrinos through EFT operators becomes sub-leading with respect to the contributions from minimal mixing. This roughly happens for couplings C P,D ≤ 10 −5 , indicating that EFT operators can dominate over minimal interactions for a new physics scale of Λ ∼ v 2 /C P,D = O(80) TeV. This scale does not include possible small dimensionless couplings or loop suppressions of the EFT operators and is thus only indicative of the sensitivity range.
For some experiments (CODEX-b, MAPP1, and FASER), there is no sensitivity in the upper left corner (small C P and large C D ). This is caused by the rather weak detector acceptance of these experiments for the light sterile neutrinos produced from D-mesons decays.
In the lower set of plots, we assume equal C P = C D , and vary the sterile neutrino mass m N and jointly the Wilson coefficients. The plots for scenario 2 and 3 look rather similar although in scenario 2, the sensitivity to smaller sterile neutrino masses is a bit better. The most sensitive experiment would clearly be SHiP, which reaches roughly C P = C D ∼ 2 · 10 −5 in the range 0.5 GeV < m N < 1.8 GeV. For couplings at this level, both the production and decay of sterile neutrinos are still dominated by the EFT operators and minimal mixing plays a sub-leading role. The sensitivity then depends a lot on the experimental setup under consideration. FASER reaches couplings at the 10 −3 level (corresponding to scales of roughly 8 TeV in Λ). Next in sensitivity is MAPP1 at 3 · 10 −4 , and then FASER2, MAPP2, CODEX-b, and ATLAS at roughly the 10 −4 level. MATHUSLA, ANUBIS, and AL3X, should be sensitive to couplings down to around 5·10 −5 , and finally SHiP at the aforementioned C P,D ≤ 2 · 10 −5 level. The hierarchy in sensitivity reach shown by the various experiments is essentially the same in scenarios 2 and 3, and is very similar to the hierarchy in the minimal scenario (see Fig. 6) for masses m N < m D . Again, the sensitivity reach in m N goes beyond the kinematical thresholds set by the pion and D-meson masses. For m N < m π , sterile neutrinos can still decay leptonically via the weak interaction. Thus larger C P values can still lead to detectable rates of sterile neutrino production. We stress again that for m N < m π , we underestimate the production of sterile neutrinos by omitting production via pions and kaons. For m N > m D , sterile neutrinos for this benchmark can still be produced from the B-meson decays via the weak current. If C P and C D are large enough, sufficiently many sterile neutrinos are produced. Furthermore, specifically for AL3X, and SHiP, and to a lesser extent MATHUSLA, the boosted decay lengths of these sterile neutrinos can fall into the respective geometric sensitivity ranges. This corresponds to the extended sensitive parameter regions, as shown on the right-hand side of the two lower plots of Fig. 7 production process via C P D ± /D 0 /D s → N + e ± (+X) . decay process via C D N → K ± + e ∓ , K * ± + e ∓

Flavor Benchmark 2
In flavor benchmark 2 we choose a different flavor-structure for the decay Wilson coefficient. For the production operator we take again (C P ) 21 but for the decay now set (C D ) 12 . This leads to sterile neutrino decay processes N → e ± + K ∓ , and N → e ± + K * ∓ .  The format is the same as in Fig. 7.
The sensitivity limits for this scenario are shown in Fig. 8 with the same format as in Fig. 7. On the left we consider C P = C 21 SRR = 4C 21 TRR , and on the right C P = C 21 VLR . Similarly for the decay we have on the left C D = C 12 SRR = 4C 12 TRR and on the right C D = C 12 VLR . In the upper Flavor benchmark 3.2 Flavor benchmark 3.3 production operator: C P C 13 SRR = 4C 13 TRR C 13 VLR decay operator: C D C 11 SRR = 4C 11 TRR C 11 VLR production process via C P B ± /B 0 /B s → N + e ± (+X) . decay process via C D N → π ± + e ∓ , ρ ± + e ∓ Table 5: Summary of flavor benchmark 3.
row, we show plots in the plane C D vs. C P with m N at 1.2 GeV. Similar features as in the previous scenario are observed and there seems hence to be little sensitivity in the event rates to the specific final-state meson.
The hierarchy in sensitivity of the different experiments is also very similar. In the lower panels we see some differences compared to flavor benchmark 1. The sensitivity to lighter m N is reduced due to the need to produce a heavier kaon in the final state. For m N < m K the EFT operators no longer contribute to the decay rate and the SM weak interaction becomes the only mechanism for sterile neutrinos to decay and be detected. This leads to a further reduction in sensitivity. We stress again that our results in this regime are conservative as we have not considered sterile neutrino production via kaon decays.

Flavor Benchmark 3
We proceed to study a scenario where sterile neutrinos are mainly produced through decays of B-mesons. Compared to flavor benchmark 1, we keep the same flavor structure for the decay Wilson coefficient, but turn on C 13 P . This leads to the sterile neutrino production via the decay processes B ± → e ± + N, B ± → π 0 + e ± + N, The relevant information is summarized in Table 5.
Our results for the sensitivity reach for this benchmark are shown in Fig. 9. The two top panels show results for the leptoquark (left) and C VLR (right) scenarios in the C P -C D plane for fixed m N = 2.6 GeV. The resulting curves are rather different from the scenarios, where sterile neutrinos are produced via D-meson decays. In the earlier flavor benchmarks, C P can be turned off and sufficient sterile neutrinos will be produced via minimal mixing to still detect sterile neutrinos, as long as C D is sufficiently large to ensure sterile neutrinos decay in the respective detector volumes. This feature has disappeared in this benchmark scenario and for C P < 10 −7 no detection is possible in any of the experiments, even for large C D . This lack of sensitivity is explained by the fact that the production rates of B-mesons are smaller than that of D-meson by roughly a factor ∼ 20 at 14 TeV pp collisions and by a factor ∼ 3000 at SHiP.
The two lower panels in Fig. 9 assume C P = C D and show sensitivity curves as a function of the sterile neutrino mass. In the previous flavor benchmarks SHiP showed the strongest sensitivity, but here it performs worse than MATHUSLA, ANUBIS, and AL3X. The reason is twofold. First, the ratio between the number of B-mesons produced and that of D-mesons is much smaller at SHiP than at the 14 TeV pp−collision experiments. Second, given their larger masses, the B-mesons, are less boosted in the very forward direction, leading to weakened acceptance of the SHiP detector for the long-lived sterile neutrinos. The hierarchies in sensitivity of the other experiments is the same as in the minimal scenario and the other flavor benchmarks. However, the overall reach is increased over the previous flavor benchmarks with the MATHUSLA sensitivity to couplings at the impressive 5 · 10 −6 level, corresponding to scales of O(100) TeV.

Flavor Benchmark 4
Flavor benchmark 4.2 Flavor benchmark 4.3 production operator: For flavor benchmark 4 the production primarily proceeds via B-meson decay, as for the previous benchmark, but here sterile neutrinos decay to a kaon through C 12 D . The relevant information is summarized in Table 6. Fig. 10 shows the numerical results for this benchmark. In the two top panels we show plots in the C P -C D plane for fixed m N = 2.8 GeV. In general these plots show very similar features as their counterparts in Fig. 9, except for an overall small reduction in the reach of C P because of the choice for a slightly larger mass of the sterile neutrino. The two lower plots show sensitivity curves in the plane C P = C D vs. m N . Compared to the lower panels of Fig. 9, they show similar exclusion limits for m N m K . However, for lighter sterile neutrinos, since only the decay modes via the weak current and active-sterile neutrino mixing are open, the sensitivity is significantly reduced. The hierarchies of the various experiments is the same as in the previous benchmark for both EFT scenarios.

Flavor Benchmark 5
In flavor benchmark 5, we turn on the operators C 13 P and C 21 D . In this case, the decay operator also leads to production of sterile neutrinos, but the resulting sterile neutrinos are restricted to a mass range where they can only decay via minimal mixing. We summarize the benchmark features in Table 7. Fig. 11 shows the resulting sensitivity reach. In the upper row we fix the sterile neutrino mass at 3.5 GeV. In general the absolute and relative sensitivities to C P and C D are comparable to the previous flavor benchmark, but the sensitivity in the bottom panel drops a bit for m N < m D . In this case sterile neutrinos only decay via minimal mixing. The exception is SHiP for which the sensitivity grows for m N < m D , where SHiP becomes the most sensitive experiment in fact, because the production cross section difference between D-and B-mesons is Flavor benchmark 5.2 Flavor benchmark 5.3

Discussion and Comparison with other Probes
Our results indicate that proposed experiments to detect long-lived particles are very sensitive to higher-dimensional operators in the νSMEFT Lagrangian. In the mass range m K < m N < m B the sensitivity curves are rather stable with respect to different EFT operators in Eq. (18) and the particular flavor configuration, as long as the sterile neutrino can be produced via the decay of D-or B-mesons, and the sterile neutrino can, in turn, decay semi-leptonically. While each particular choice of EFT operators and flavor assignment requires a detailed study, we find that FASER is sensitive to Wilson coefficient couplings of about ∼ 10 −3 (this is extended to ∼ 10 −4 for FASER2), while experiments such as MATHUSLA, ANUBIS, AL3X, and SHiP can reach down to coupling strengths of 5· ∼ 10 −6 . From the matching relations in Eq. (11), we see that such limits can be used to constrain the νSMEFT operators C Hνe , C duνe , C LνQd , C LdQν , and C QuνL . Assuming a scaling of these Wilson coefficients as ∼ v 2 /Λ 2 , the sensitivities range from Λ ∼ 8 TeV for FASER up to Λ ∼ 100 TeV for the larger experiments.
The νSMEFT operators we consider here can also be probed in other experiments. These include meson and tau decays, elastic coherent neutrino-nucleon scattering (ECνNS), missing transverse energy searches, etc. Depending on the probe, the relevant sterile neutrino mass range and the flavor assignment can differ from the cases considered in this work. For instance, limits from pion decay or from neutron or nuclear beta decay require sterile neutrino masses below the pion mass or the respective Q value of β decay, considerably lower than the GeV-scale sterile neutrinos considered in this work. Here we briefly give an overview of the literature.
Refs. [96,97] investigated limits from pion decays, tau decays, and singular leptons with missing transverse energy. The most restrictive bounds on the new physics scale are obtained from pion decays with Λ 36 TeV. However, tau decays allow for a neutrino mass range m N more comparable to our studies, while searches for l + E T are largely independent of sterile neutrino masses. The latter investigations set the new physics scale to Λ 2 − 5 TeV. We did not explicitly consider processes involving τ leptons, but there should be good sensitivity in the appropriate mass range m τ + m π < m N < m B − m τ . More quantitative statements require a detailed study that includes sterile neutrino production via τ decays and an efficiency factor for reconstructing decays of τ mesons in the final states.
Refs. [98,99] consider a larger set of pseudoscalar meson decays corresponding to several flavor configurations. Additionally, the effects of νSMEFT operators on lepton flavor universality (LFU), CKM unitarity, and β-decays are examined. The most stringent bounds are on the operators C (6) LνQd , C LdQν , and C (6) QuνL involving an up quark, a down (strange) quark and an electron using LFU constraints. The new physics scale is limited by Λ 74 (110) TeV in the limit of massless sterile neutrinos and thus cannot be directly compared to results obtained here. Bounds on other operators and different flavor combinations are in the range of Λ 0.5 − 8 TeV. Similar sensitivities are found examining anomalies in the transition b → cτ ν including light sterile neutrinos (m N 100 MeV) [100]. Further, limits from ECνNS based on the COHERENT experiment [101] are considered in Refs. [98,102,103]. Sterile neutrinos considered in these works are again much lighter than the GeV-scale (m N 0.5 MeV) and the resulting bounds are at the level of Λ 1 TeV. We conclude that the sensitivities of the experiments considered here are competitive with and complementary to existing constraints.
In this work we have focused on Majorana neutrinos (although our sensitivity curves are not affected dramatically if we had considered Dirac neutrinos instead) for which strong constraints can be set from 0νββ experiments. In Ref. [46] some of us developed a framework to calculate 0νββ decay rates in the presence of light sterile neutrinos and the νSMEFT Lagrangian. In particular, we investigated the reach of current and future experiments to probe scenario 2: the 3 + 1 leptoquark model. As sterile neutrinos appear as virtual states, 0νββ experiments are sensitive to a broad range of neutrino masses with a peak sensitivity at m N 100 MeV, which drops off for larger or smaller masses. To make a comparison, we consider the case y LR 11 y RL * 11 = y LR 21 y RL * 11 = 1 and y LR 11 y RL * 11 = y LR 11 y RL * 31 and vanishing couplings for other flavors. We can then compare flavor benchmarks 1 and 3 to the sensitivity of 0νββ experiments, which only depend on sterile neutrino couplings to first-generation quarks and leptons. For these choices of couplings, we can calculate 0νββ decay rates and determine the LHC and SHiP sensitivity curves as a function of m N and m LQ (0νββ rates have a very small dependence on phases appearing in the 3 + 1 neutrino mixing matrix and we neglect this dependence here for simplicity). The results are shown in Fig. 12. We stress that the uncertainties associated with hadronic and nuclear matrix elements for 0νββ decay rates are sizable and not included in the plot, for details we refer to Refs. [46,104]. For flavor benchmark 1 (left panel of Fig. 12), the limits from 0νββ are somewhat stronger than the prospected sensitivity of FASER2 and MATHUSLA, chosen as representative experiments, in the relevant mass range. For flavor benchmark 3 the prospected MATHUSLA overtake current 0νββ limits for masses between 1 and 5 GeV.
We stress that the bounds from 0νββ decay experiments are only valid for Majorana neutrinos and final-state electrons. However, the sensitivity curves for the various LHC experiments discussed here are (roughly) valid for (pseudo-)Dirac neutrinos, and in the appropriate mass range also for couplings to muons and, to lesser extent, taus instead of electrons, and in general to a broader range of quark flavors. They are thus more general than 0νββ limits albeit in a much small sterile neutrino mass range.

Conclusions
The possibility of sterile neutrinos provides one of the main motivations for the search for longlived particles (LLPs). Sterile neutrinos provide compelling solutions to major problems in particle physics and cosmology, such as active neutrino masses and the baryon asymmetry of the Universe. Sterile neutrinos are in fact predicted in a variety of theoretical models, ranging from the minimal scenario where they interact with Standard Model (SM) particles through minimal mixing with active neutrinos, to more exotic scenarios involving new fields with masses well above the electroweak scale such as left-right symmetric models or grand unified theories.
In this work, we focused on relatively light sterile neutrinos with masses at the GeV scale, down to about 100 MeV. This mass range is interesting, as it is linked to low-scale leptogenesis and opens up the possibility of efficiently producing sterile neutrinos through the decays of pseudoscalar mesons, which are copiously produced in collider experiments. In particular at CERN, besides the ATLAS and CMS LHC collaborations, various proposed experiments are presently under discussion specifically targeting the detection of LLPs, such as the fixed-target experiment SHiP and a number of so-called far detectors at various pp-collision interaction points e.g. FASER and MATHUSLA. A large number of mesons are expected to be produced at the interaction points of these experiments, which in turn can decay to sterile neutrinos. We focused on sterile neutrinos which can be produced from bottom and charm meson decays in hadronic collisions, and investigated the sensitivity reach of present and future LHC experiments, and SHiP to detect sterile neutrinos.
To avoid theoretical bias we approached this problem in the framework of the neutrinoextended SM effective field theory (νSMEFT). This framework allows for a light gauge-singlet fermion, a sterile neutrino or heavy neutral lepton, and assumes other BSM fields to have masses M v, the SM Higgs vacuum expectation value. This framework describes effective local interactions between sterile neutrinos and SM fields in a systematic expansion. We considered dimension-6 operators that allow for sterile neutrino production via mesonic decays at tree level. Our framework is general, but for concreteness we considered specific scenarios where a subset, or only a single, EFT operator is turned on at the same time. These scenarios are motivated by UV completions like leptoquark or left-right-symmetric models. Other EFT scenarios or specific UV-complete models can be straightforwardly investigated by using the extensive formulae given in the appendix. To benchmark our calculations with the literature, we also considered the minimal scenarios where higher-dimensional operators are turned off.
We performed Monte-Carlo simulations to evaluate the sensitivity reach of the considered experiments: ATLAS, CODEX-b, FASER, MATHUSLA, AL3X, ANUBIS, MoEDAL-MAPP, and SHiP. For the minimal scenario, the obtained sensitivity curves are in agreement with existing results with minimal discrepancies, while we obtained sensitivity curves for the two MoEDAL-MAPP experiments for the first time. For the EFT scenarios we consider active-sterile neutrino mixing, the EFT operators, and their interference terms simultaneously. For each EFT scenario, we considered a series of different flavor benchmarks, where the EFT operators induce either D-or B-meson decays into sterile neutrinos, and the sterile neutrinos decay to an electron and various mesons, N → e + M ij , cf. Tables 3-7. For the D-meson benchmarks, we found that SHiP and MATHUSLA have the most extensive sensitivity reach. They are sensitive to dimensionless Wilson coefficients at the 10 −5 level, for most of the kinematically allowed sterile neutrino mass range. For such values of couplings, the production and decay of sterile neutrinos is dominated by the higherdimensional operators with minimal sterile-active mixing playing a subleading role. Apart from dimensionless couplings and potential loop suppressions, the dimensionless Wilson coefficients scale as v 2 /Λ 2 , where Λ is the high-energy scale where the νSMEFT operators are generated. The sensitivity drops at the edges of the allowed mass range, but does not disappear completely, even for sterile neutrinos with masses above m D , because of the contributions from SM weak interactions and active-sterile mixing. Assuming a v 2 /Λ 2 scaling, SHiP and MATHUSLA could probe scales around 80 TeV. This scale is lowered to 8 TeV for FASER and 25 TeV for FASER2, which are much smaller experiments.
For our B-meson benchmarks, because of its much smaller B-meson production rates and a weaker acceptance, SHiP does not show the best performance. Instead, we found that MATHUSLA and ANUBIS would be sensitive to Wilson coefficients at the 5 · 10 −6 level. The sensitivity curves appear to be fairly independent of the Lorentz structure and the exact flavor configuration of the EFT operators, as long as sterile neutrinos can be produced through D-or B-meson decays and sterile neutrinos can decay semi-leptonically.
For our ATLAS study we applied an overall flat efficiency factor of 10 −3 . Of all the experiments we have studied here only ATLAS is operational. Thus we strongly encourage our colleagues at ATLAS and CMS collaborations to perform a proper full analysis of the scenarios we have presented and investigated here. This should put our approximations on a much firmer footing.
We compared our projected sensitivities with (projected) constraints obtained from other probes of sterile neutrinos with effective interactions such as light pseudoscalar meson decays, tau decays, coherent neutrino-nucleus scattering, LHC searches for missing transverse energy, and neutrinoless double beta decay. Our results are very competitive and complementary. We conclude that searches for displaced vertices of long-lived sterile neutrinos at the LHC and SHiP are an excellent probe of νSMEFT operators and of sterile neutrinos in general.
A straightforward extension of our work is to final state muons instead of electrons. Here we expect basically the same results, except at the very lowest end of the sterile neutrino mass we have considered. A more involved extension would be to also consider the production of sterile neutrinos from the decay of the light K and π mesons. Furthermore in a future project, we shall consider the case of final state τ 's. This will restrict the kinematically viable regions, but would also require a proper investigation of the detection efficiencies.
for the pseudoscalar current. The vector and tensor currents only induce two-body final states for vector mesons. We define where |M * (q, ) denotes a vector meson M * with mass m M * , momentum q and polarization µ . Heavy-quark symmetry relates f T M f V M . All decay constants are given below in Appendix C. Armed with these decay constants, we calculate the production and decay rates of neutrino mass eigenstates starting from Eq. (18). We begin with neutrino production via the decay of pseudoscalar mesons M − ij → l − k + ν l , and the corresponding decay ν l → M + ij + l − k where ij denotes the generation of quark flavors that make up the meson (we drop these indices below for notational convenience), k the charged lepton generation, and l = {1, . . . ,n} the neutrino mass eigenstate. For the decay of the neutrino mass eigenstate we also include the decay to the charge-conjugate final state which is equally likely due to the Majorana nature of ν l .
We obtain for the summed-over-spins squared amplitudes for sterile neutrino (N ) production VRR )(C where all Wilson coefficients carry flavor indices ijkl and we defined the functions For sterile neutrino decay, we include the charge-conjugated final states leading to an additional factor 2 compensating the additional 1/2 from initial-spin averaging. We then obtain where ab = {V V, SS, V S}, i = {1, 2}, and g ab,i = −f ab,i .
Similarly, for processes involving vector mesons we find VLL (C VLR + C For sterile neutrino decay where ab = {V V, T T, V T }, i = {1, 2}, and g * where n is the appropriate spin-averaging factor, m 1 is the mass of the decaying particle, and m 2 and m 3 are the masses of the final-state particles. The phase space function is defined as

A.2 Sterile Neutrino Decays into Neutral Pseudoscalar Mesons
We follow Ref. [55] and write the neutral weak axial current as where j A 3,µ = 1 √ 2 (ūγ µ γ 5 u −dγ µ γ 5 d) , j A 8,µ = which correspond to the neutral vector mesons ρ 0 , ω, φ, and J/Ψ, respectively. We rewrite The hadronic matrix elements are defined as where a = ρ 0 , ω, φ, J/Ψ. The decay width becomes where f a and g a are listed in Table 10 and m a is the mass of M 0 a,V .

B.1 Automizing Three-body Phase Space Integral Calculations
This work involves a large amount of three-body phase-space computations, which are straightforward but tedious to evaluate. We briefly discuss here our approach to automize these computations using Mathematica. In the rest-frame of the decaying pseudoscalar meson, the decay rate is given by where p, p , p l , and p N denote the momentum of the decaying meson, outgoing meson, SM lepton, and sterile neutrino, respectively, and M is the mass of the decaying pseudoscalar meson. |M| 2 denotes the spin-averaged product of the leptonic and hadronic matrix element squared. It can be decomposed into a hadronic form factor, which only depends on q 2 , where q ≡ p − p , and a function of four-momentum invariant scalar products. The form factors are defined below and the spin-averaged matrix elements are calculated with standard techniques and checked with PackageX [106]. We convert to four-dimensional integrals and write where for notational convenience we have omitted three Heaviside step functions. We perform the p N integral by setting p N = q − p l and introduce the variable a via a factor of 1 = da δ(a − q 2 ) to obtain Γ = 1 2M 1 (2π) 5 da d 4 p d 4 p l δ(a − q 2 )δ(p 2 − m 2 )δ(p 2 l − m 2 l )δ((q − p l ) 2 − m 2 N )  f ηc 0.335 GeV [115] θ 0 -6.9° [114] θ 8 -21.2° [114] Decay constants for pseudoscalar and vector mesons are given in Tables 8. Parameters to calculate decay constants for the neutral pseudoscalar and vector mesons are given in Tables 9  and 10, respectively. meson M f a [GeV 2 ] g a ρ 0 [56] 0.171 1 − 2 sin 2 θ w ω [56] 0.155 -2 3 sin 2 θ w φ [56] 0.232 √ 2(− 1 2 + 2 3 sin 2 θ w ) J [116] 1.29 √ 2( 1 2 − 4 3 sin 2 θ w )     Table 12: Best fit parameters for the form factors in D → K and D → π transitions from Refs. [55,118].

C.2 Form Factors for D → M P
For D → π and D → K transitions, we use the methods of Ref [118]. We write and list the best-fit parameter values in Table 12.

C.3 Form Factors for B (s) and D Decaying into M V
Eq. (86) contains seven form factors which must be determined. The pseudoscalar form factor is related to f (q 2 ) and a ± (q 2 ) through Eq. (87). To better present these form factors and follow   [55,120].   [55,120].
and choose the following three-parameter formula f (q 2 ) = f (0) (1 − q 2 /M 2 pole )(1 − σ 1 q 2 /M 2 pole + σ 2 q 4 /M 4 pole ) , to describe V , T 1 and A 0 . M is the pole mass, which is M P (0 − ) for A 0 and M V (1 − ) for V and T 1 . For the remaining form factors A 1 , A 2 , T 2 and T 3 , we use the simpler form [120] f (q 2 ) = f (0) In all scenarios we consider, the transition B s → D * s is only induced by SM weak interactions and we do not list form factors associated to BSM currents. The values of the best fit parameters are given in Table 13 to Table 15.