Probing right-handed neutrinos dipole operators

We consider the minimal see-saw extension of the Standard Model with two right-handed singlet fermions $N_{1,2}$ with mass at the GeV scale, augmented by an effective dipole operator between the sterile states. We firstly review current bounds on this effective interaction from fixed-target and collider experiments as well as from astrophysical and cosmological observations. We then highlight the prospects for testing the decay $N_2 \to N_1 \gamma$ induced by the dipole at future facilities targeting long lived particles such as ANUBIS, CODEX-b, FACET, FASER 2, MAPP and SHiP.


Introduction and framework
The see-saw mechanism [1][2][3][4][5] is arguably the simplest extension of the Standard Model (SM) that is able to explain the observed pattern of neutrino masses and oscillations.In its simplest incarnation, it consists in adding to the SM particle content a right-handed (RH) neutrino, that is a spin 1/2 fermion, singlet under the SM gauge group, which has a Yukawa interaction with SM leptons, as well as a Majorana mass term.One of the active neutrinos acquire thus a non-vanishing mass m ν and a mixing θ with the new sterile state, parametrically expressed by the relations where y ν and m N are the RH neutrino Yukawa interaction and mass respectively and v is the electroweak (EW) vacuum expectation value (VEV).Since experimental data point to at least two massive neutrinos, at least two RH states must be added to obtain a realistic phenomenology.In this case the essence of the see-saw mechanism is unaltered, with the obvious promotion of y ν and m N to matrices in flavor space, but the relations of Eq. (1.1) turn out to be modified.In particular the mixing angles can receive an exponential enhancement with respect to the naive see-saw scaling case, that may drastically modify the phenomenology.This is best seen in the Casas-Ibarra parametrization [6].From the practical point of view, this means that masses and mixings can be treated as independent parameters.Irrespective of this consideration, by fixing m ν Eq. (1.1) doesn't uniquely point to a preferred mass range for m N , which could lie all the way up the grand unification scale if y ν is an O(1) parameter.However, in recent years RH neutrinos with mass below the EW scale have gained more and more attention in that they can explain the matter-antimatter asymmetry of the Universe via neutrino oscillations [7,8] and, crucially, can be tested at present and future colliders and fixed-target experiments, see e.g. .
While the see-saw model is a full-fledged ultraviolet (UV) complete theory, at least in the same way as the SM is, in the case where RH neutrinos lie at the EW scale it is interesting to consider it as a low energy effective field theory (EFT) extended with higher dimensional operators built from the SM and the RH neutrino fields.The resulting theory is called νSMEFT and is described by the following Lagrangian where N is a vector describing N f flavors of RH neutrino fields and N c = C N T , with C = iγ 2 γ 0 .Furthermore, Y ν is the 3 × N f Yukawa matrix of the neutrino sector with H = iσ 2 H * , M N is a N f × N f Majorana mass matrix for the RH neutrinos and O n the Lorentz and gauge invariant operators with dimension n built out from the SM and the RH neutrino fields, with Λ parametrizing the Wilson coefficient of the operator.A complete and independent set of operators has been built up to dimension nine [36][37][38].
Interestingly, already at d = 5 two genuine νSMEFT operators appear 1 .The first is an operator coupling the RH neutrinos with the Higgs boson, O 5 N H = N c N H † H.This triggers a new decay mode for the Higgs, with interesting consequences for collider phenomenology, both at the Large Hadron Collider (LHC) [39][40][41][42] and future colliders [43].The second operator is a dipole with the hypercharge gauge boson 2 O 5 N B = N c σ µν N B µν , which has so far been less investigated [43][44][45][46][47][48] 3 .Among other effects, this operator generates the decay 4   N heavy → N light + γ , m N heavy > m N light . ( This interaction is the subject of our study.Our focus will be on light RH neutrinos with masses up to a few GeV.Such light states can be produced not only at high energy colliders via parton interactions, but also at fixed-target experiments, typically via meson decay.More specifically, we will analyze in detail the current bounds from colliders experiments, such as LHC, LEP and BaBar, and fixed-target experiments, such as CHARM [53], NuCal [54,55] and NA64 [56].We will then compute the predicted sensitivity to the νSMEFT parameter space of the proposed experiments ANUBIS [57] , CODEX-b [58][59][60], FACET [61], FASER 2 [62,63], MAPP [64,65] and SHiP [66,67].In addition, we will also discuss constraints from cosmology and astrophysics.Throughout this work we will consider the theory of Eq. (1.2), focusing on the d = 5 dipole operator.Given its symmetry properties, this operator is non-vanishing only for N f ≥ 2. Since we are primarily interested in probing the effect of the dipole operator, we will work under the assumption that the active-sterile mixing effects are negligible for what concerns the heavier sterile neutrinos phenomenology, in such a way that their decay proceeds only via the dipole operator under our scrutiny through the process of Eq. (1.3).As for the lightest RH neutrino N 1 , its decay pattern is completely determined by the active-sterile mixing, as in the standard see-saw case.As we are going to discuss in Sec. 2, the N 1 lifetime can be strongly constrained by cosmological observations, especially the ones related to the epoch of Big Bang Nucleosynthesis (BBN).A relatively safe scenario is the one in which N 1 mixes dominantly with the third generation of SM neutrinos, ν τ .We found that this configuration can be easily obtained by choosing N f = 3, satisfying at the same time all other relevant constraints.Interestingly, in this case the heaviest RH neutrinos N 3 can be decoupled from the spectrum without affecting the mixing pattern for N 1 , leaving only the two lightest RH neutrinos N 1,2 as dynamical states.In presenting our main findings we will thus consider a framework with only these two states living at the EW scale and interacting via the dipole operator which we normalize as where m N 2 > m N 1 , and where g Y and B µν are the U (1) Y coupling and field strength tensor respectively.The loop suppression factor is explicitly introduced since this operator only arises at loop level in any weakly coupled UV completion, see e.g.[68,69], while the hypercharge coupling is added because of the presence of B µν .Explicit UV completions include models with additional scalar and fermions or models with additional vectors and fermions, with non-vanishing hypercharge [44,70].We will comment later on possible strongly interacting UV completions.Since the Wilson coefficient can be complex, we show explicitly its phase α.In this scenario, O5 N B completely governs the RH neutrinos phenomenology 5 .In particular, it dictates the heaviest neutrino N 2 total decay width.For RH neutrinos below the Z mass the dominant decay mode is where c w is the cosine of the Weinberg angle and the last equality holds for small values of δ, which is defined as The three-body decay into an off-shell Z boson provides a subdominant contribution.As it is clear from Eq. (1.5), the relative mass splitting δ is crucial in determining the RH neutrino decay length, and hence its lifetime.This gives an indication on the type of experiments that can have a sensitivity to this scenario, depending on how far the detector is located with respect to the N 2 production points.For example the neutrinos N 2 could decay promptly, i.e. with a typical decay length smaller than O(mm).In this case they are a primary target for standard collider searches.Their lifetime could also be longer, with corresponding decay lengths in the O(1 m − 100 m), for which different strategies need to be envisaged.In the more extreme case, they can be stable with respect to the length scale of any terrestrial experiment and hence completely invisible for what concerns laboratory searches.We will comment upon all these possibilities in the following, mainly focusing however on a region of parameter space in which the heavier neutrino N 2 is a long-lived state with a macroscopic decay length.This choice has a twofold motivation.
From one side, light new states with suppressed interactions, as the one inherited from the dipole operator, have usually a long lifetime.From the other side, the study of long-lived particles is an active field which has received a lot of attention in the last years, following the philosophy of leaving no stones unturned and lighting new lampposts in the quest of new physics beyond the SM.In this area, big experimental progresses are foreseen in the mid-and short-term.
The relative mass splitting is also important in determining the photon energy arising from the decay of Eq. (1.3).From basic kinematics in the N 2 rest frame one has (1.7) Assuming the photon to be produced collinearly with the direction of N 2 in the laboratory frame, which maximizes the photon energy in this frame of reference, one has where P N 2 is the modulus of the N 2 spatial momentum and the last equality holds for m N 2 /P N 2 1 and δ 1 6 .Thus the smaller the relative mass splitting the softer the final state photons, which however should satisfy some minimal threshold requirement in order to be identified in a detector.Hence too small relative mass splittings will hardly be experimentally testable.
The rest of the paper is organized as follows.In Sec. 2 we review the existing limits on the dipole operator from cosmology, colliders and other type of experiments, while in Sec. 3 we discuss the future sensitivity of SHiP and FASER 2 on the model parameter space, wrapping up our conclusion in Sec. 4. We also add three appendices.In App.A we present results for other future LHC experiments targeting long-lived particles, in App.B we report useful formulae for computing the decay of a QCD meson into a pair of RH neutrinos via the dipole operator and in App.C we discuss possible constraints arising from searches of electrons recoil signatures in laboratory experiments.

Current limits from cosmology, colliders and other experiments
The parameter space of the simplified scenario considered in this work is spanned by the lighter neutrino mass m N 1 , the relative mass splitting δ with the heavier RH neutrino, and the Wilson coefficient of the dipole operator, parametrized by Λ and its phase α.This parameter space is already constrained by laboratory data from colliders and past beam dump experiments, as well as by astrophysical and cosmological measurements.In this section we will review the most important and stringent ones.Particular care must be taken in ensuring the validity of the EFT in the various considered processes.The dipole operator in Eq. (1.4) induces N 1 N 2 production through the exchange of a photon or a Z boson.Following Ref. [71] and assuming couplings of order one, we identify the EFT cut-off scale with Λ, and for the EFT to be valid we require where ŝ = (p N 1 + p N 2 ) 2 is the Lorentz invariant energy that enters the vertex.One important production mechanism for light N 1,2 is via heavy meson decay, that can be copiously produced in fixed-target experiments.In this case the heavy neutrino production proceeds via an s-channel γ and we have ŝ = m 2 M , with m M the meson mass.For higher masses direct production at collider can be relevant.In this case, ŝ is the center of mass energy squared of the parton pair that exchange the photon or the Z boson, e.g. e + e − for LEP or q q for the LHC.Analogous considerations apply for other production modes, as for example production via photon bremsstrahlung.

Fixed-target experiments
We start our discussion with fixed-target experiments, for which we have considered data collected at CHARM [53], NuCal [54,55] and NA64 [56].We consider the production of RH neutrinos from the decay M → N 1 N 2 , with M a vector meson produced at these experiments. 7HARM: In the CHARM experiment, a 400 GeV proton beam was dumped on a copper target.The detector, placed at a distance of 480 m from the interaction point (IP) and 5 m off the beam axis, consisted of a decay volume 35 m long and with a surface area of 9 m 2 .We have modeled the detector following [72] and recast the analysis of [53], in which an axionlike particle (ALP) decaying into a pair of photons was searched for.Since the analysis required only a single electromagnetic shower, it can be applied to the decay N 2 → N 1 γ.We compute the number of events expected at CHARM following the equations that will be described in more details in Sec. 3, see Eq. (3.1) and subsequent ones.In our analysis, we simulate the production of N 1 N 2 pairs from the decay of the mesons ρ, ω, J/Ψ and Υ using PYTHIA 8.3 [73,74], finding the following production multiplicities: N ρ = 0.58, N ω = 0.57, N J/Ψ = 4.7 × 10 −6 and N Υ = 2.2 × 10 −9 , see Sec. 3 for their definition.Then, we require the energy of the photon in the laboratory frame to satisfy E γ ≥ 1 GeV and, following [53], we take a signal acceptance of 51%.The number of protons-on-target (POT) is taken to be N POT = 2.4×10 18 .Since no signal events were observed in the search of [53], we set an upper limit at 95% confidence level (CL) of N signal = 3.The region excluded by the CHARM experiment is shown in Fig. 1 , 2 and Fig. 3.
NuCal: In the ν-calorimeter I experiment (NuCal), a 70 GeV proton beam from the U70 accelerator was dumped on an iron target.The detector consisted of a cylindrical decay volume 26 m long and with a diameter of 2.6 m, placed at 64 m from the iron target.We implement such geometry accepting N 2 events with a maximum angle of 0.014 rad from the beam axis.To set a limit on the parameter space of our scenario, we simulate N 1 N 2 production from ρ and ω decays8 using PYTHIA 8.3 obtaining the following production multiplicities: N ρ = 0.30, N ω = 0.30.Then, we follow the analysis in [75], requiring the photons produced in the N 2 → N 1 γ decay to satisfy two conditions: their energy in the laboratory frame must be E γ ≥ 3 GeV, while their angle with respect to the beam axis must satisfy θ γ < 0.05 rad.After these selection cuts, 5 events were observed, with an estimated background of 3.5 events from the simulated neutrino interactions in the detector [54].Given these numbers, assuming Poisson likelihood we set a 95% CL upper limit of N signal ∼ 7.1 [75].The region excluded by NuCal is shown in Fig. 1 , 2 and 3.
NA64: In the NA64 experiment, an electron beam of 100 GeV was dumped on a lead target.We have considered the analysis of [76], in which an ALP decaying into a pair of photons was searched for.Since the two photons are too collimated to be distinguished, the final state was reconstructed as a single photon, allowing us to reinterpret this search.In this case the N 1 N 2 pair is produced via photon bremsstrahlung.Using 2.84 × 10 11 electrons-on-target [76], we find that the number of events produced at NA64 is too small to put any bound on the parameter space of the model.

Colliders
We now analyze the bounds enforced by collider experiments, by considering searches performed at LEP, BaBar and LHC.In this case different searches apply, depending on whether the N 2 → N 1 γ decay is prompt, i.e. it happens at a distance smaller than ∼ 1 mm from the IP, displaced, i.e. it happens within ∼ 1 mm and ∼ 1 m, or else is detector-stable, i.e when the decay happens at a distances greater than ∼ 1 m.For a 2 → 2 scattering, in terms of ŝ = (p N 1 + p N 2 ) 2 , we have while the N 2 lifetime is given by τ , with the decay width of Eq. (1.5).
In the case of prompt decays, we consider two analyses: one from LEP [77] and one from BaBar [78].In the LEP analysis, data were taken at various center of mass energies around the Z peak.We employ the largest dataset, taken at √ ŝ = 91.26GeV with an integrated luminosity of 52.462 pb −1 .We have simulated our signal at the parton level by using MadGraph5 aMCNLO [79].We enforce the analysis selections by requiring a single photon with | cos θ γ | < 0.7 and considering two signal regions.In the first one a minimum energy of the photon was required E γ > 22 GeV, and no events were observed.In the second region the cut is loosened to E γ > 3 GeV, with 73 observed events and 72 ± 5 expected SM events.Therefore, we set a 95% CL upper limit on the number of signal events N signal = 3 for the first signal region and N signal ∼ 22 for the second one.For the weakly coupled normalization of the operator shown in Eq. (1.4), the stronger bound corresponds to Λ 17 − 50 GeV for δ = 0.1 − 1.However, these values of Λ lie outside the range of validity of the EFT, implemented as in Eq. (2.1), therefore we conclude that no meaningful constraints on the cutoff scale Λ can be set.
We then move to the analysis of BaBar in Ref. [78], which derived bounds on single photon events produced in association with an invisibly decaying dark photon.The selection of signal events makes use of a multivariate Boosted Decision Tree discriminant.Given the complexity of the analysis, we adopt a simplified strategy to estimate the constraint.We use MadGraph5 aMCNLO to simulate a sample of e + e − → N 1 N 2 events at the center of mass energies corresponding to the Υ(2s), Υ(3s) and Υ(4s) resonances, considering the luminosities reported in [78].We enforce the selections −0.4 < cos θ γ < 0.6 and E γ > 3 GeV, corresponding to the LowM region of Tab. 1 of [78].To extract a bound, we focus on the loose R L selection of [78], and assume an equal number of observed and background events.This allows us to set an upper limit at 95% CL of N signal ∼ 28.The excluded region is largely independent of m N 1 but depends quite strongly on δ, since for δ 1 the energy of the photon in the laboratory frame is too small to pass the 3 GeV cut, see Eq. (1.8).We find that, for δ = (0.5 − 1), the bound extends up to Λ ∼ 60 GeV, while for δ < 0.5 the bound disappears.Given the approximate nature of our computation, we do not show these results explicitly in our figures.
Turning to searches for displaced events, we have considered analyses from ALEPH, ATLAS, CDFII and DELPHI [80][81][82][83].Their reinterpretation is generally less straightforward than the ones for prompt signatures, due to the need of cutting on additional quantities as the photon time of flight and pointing variable.We decide to firstly apply a simplified strategy, by only imposing energy threshold and angular selection cuts on the final state photon.In this way the limit that we extract will be stronger than the one obtained by a full implementation of the analysis.Using this strategy we find that the strongest bounds come from the DELPHI search [83], see also [84].In this case we enforce E γ > 10 GeV and |η γ | < 4.04, which corresponds to the angular coverage between 2 • and 178 • reported in the analysis.We have once again simulated e + e − → N 1 N 2 and the subsequent N 2 → N 1 γ decay using MadGraph5 aMCNLO, with center of mass energies between 180 and 209 GeV, and with the corresponding luminosities as reported in [83].Following this approach, we obtain 95% CL limits which are in the 20 − 40 GeV ballpark, for δ = 0.1 − 1.Given these results, we avoid implementing the full selection for the displaced analysis, since the obtained limits are already to be discarded because they lie beyond the validity of the EFT.
Finally, when N 1 and N 2 are both detector-stable, we consider searches of mono-γ with missing energy, and the LEP limits on the invisible Z width.It turns out that the strongest bound comes from the latter.By requiring Γ(Z → N 1 N 2 ) < 0.56 MeV [85] we obtain Λ 9 GeV, which again lies beyond the validity of the EFT 9 .We conclude this section by observing that some of the searches mentioned above could put meaningful constraints on the parameter space of strongly coupled UV completions of the EFT dipole operator.Implicit in our identification of the EFT cut-off scale with Λ in Eq. (1.4), is the hypothesis that the dipole operator is generated perturbatively at one loop level by some heavy states.An alternative possibility could be for the dipole operator to be generated by some strong dynamics, similarly to what happens for the neutron in QCD.In this case, adopting the convenient parametrization of the dipole operator O 5 N B = 1/Λ N c σ µν N B µν , one expects the EFT cut-off scale to be of the order of Λ = Λ (16π 2 )/g Y .The bounds discussed above from the LEP search [77] valid for prompt decays are simply rescaled into Λ (8−23) TeV for δ = (0.1−1), essentially independent of m N 1 .Analogously, the simplified approach adopted for the DELPHI search [83] for displaced decays leads to the constraint Λ (9 − 18) TeV, again for δ = (0.1 − 1), essentially independent of m N 1 .Clearly in this case a more thorough reinterpretation of the analysis will be needed, with respect to the simplified approach previously described.Finally, the bound from the invisible Z decay width valid for the detector-stable case would read 4 TeV.Clearly, these constraints are probing a relevant part of the parameter space lying inside the regime of validity of the EFT, i.e. √ ŝ < Λ .It is important to notice that such bounds might be quite at odds with the range of N 1 and N 2 masses to which we are interested in.For example in a QCD-like strongly coupled scenario, we expect N 1 and N 2 to emerge as baryons, with masses of order Λ and not much lighter as it would emerge from our analysis.A possibility of having a composite state much lighter than Λ could be envisaged in a scenario where a light baryon arises in order to match the anomaly of an unbroken global chiral symmetry in the UV, along the lines of [86,87].We are not aware of any realistic model realizing such a framework.For this reason, in the remainder of the paper, we will consider only the weakly coupled scenario of Eq. (1.4).

Bounds from astrophysics and cosmology
In addition to the bounds presented above, limits from astrophysics and cosmology may be important for the scenario that we are considering.The constraint which is more relevant for us comes from BBN.Although in our simplified scenario the dipole operator O 5 N B completely governs the N 2 decays, the fate of N 1 is determined by its mixing with the active sector.Particularly dangerous is the situation in which the N 1 decays could potentially spoil the predictions of the standard BBN model [88][89][90].There are two natural ways to avoid this bound.Either N 1 is stable on cosmological scales, or N 1 decays with τ N 1 10 −2 9 Notice that in any UV completion of the dipole operator new states with masses around the EFT cut-off scale will be present, among which there will also be states with non-vanishing electroweak charges.Therefore, if their masses are small enough, additional bounds, that we are not discussing, could arise from the on-shell production of these particles.
s.In the first case, N 1 would be a dark matter candidate 10 .In the second case, it has been shown in [88] that the combination of limits from BBN and terrestrial experiments exclude m N 1 (0.4 − 0.5) GeV, for N 1 mixing dominantly with ν e or ν µ , while lighter masses can be allowed for mixing dominantly with ν τ .Since in the first case an important region of parameter space that can be tested by the experiments we consider would be excluded, we turn to the case of dominant mixing with ν τ .Can such mixing be obtained in a way which is compatible with neutrino mass generation?As shown in [89], in a scenario with only two RH neutrinos this is possible for m N 1 0.5 (0.1) GeV for normal (inverted) hierarchy.The situation becomes less constrained considering three RH neutrinos, since in this case we have explicitly checked that a dominant N 1 − ν τ mixing can be obtained, in a way compatible with the generation of neutrino masses, for m N 1 0.1 GeV, which is the range we consider.This can be obtained also by assuming a mass hierarchy m N 1 ∼ m N 2 m N 3 , i.e. in a situation in which the phenomenology is driven by N 1 and N 2 as the one we are considering by using the simplified scenario of Eq. (1.4).In what follows, we will always implicitly suppose this to be the case.For what concerns N 2 , in Figs. 1, 2, 3 we show contours of constant N 2 lifetime, in order to highlight the region of the parameter space where N 2 decays fast enough to avoid BBN bounds.
Limits derived from supernovae may also be important.Light particles produced in the interior of supernovae, which reach temperatures of several tens of MeV, can escape from the star, therefore cooling the system.From this argument, masses up to a few hundred MeV can be constrained.In the δ 1 limit, we expect our scenario to be qualitatively similar to the one studied in [85], where the Authors focus on a dipole operator constructed with a single new Dirac fermion playing the role of the dark matter.They obtain bounds up to masses of ∼ 0.1 GeV.In the mass range 1 − 100 MeV, these limits exclude 2 TeV Λ 50 TeV.The constraints disappear for larger Λ because the production inside the supernovae is suppressed, while for smaller Λ efficient scattering processes can partially trap the particles inside the system.For the case δ 1 we expect the situation to be qualitatively different from the one above.For small enough values of Λ, N 2 decays quickly, thus leaving a dominant population of N 1 inside the supernovae.The only relevant N 1 scattering process is N 1 SM → N 2 SM, which however might be not kinematically allowed for large δ, possibly making the bound disappearing at small values of Λ.This scenario has been studied in the context of inelastic dark matter with a dark photon mediator in [92], albeit by using a simplified and conservative approach.Given the complexity of performing a detailed analysis of the supernovae bounds, and the fact that we expect these limits to affect only a quite limited region of parameter space, we defer a detailed study of this problem for future work.

Projected sensitivity of the SHiP and FASER 2 experiments
In this section we study the prospects for detection of long-lived RH neutrinos at the proposed future experiments SHiP [66] and FASER 2 [62,63].These experimental facilities have the capability to probe long-lived particles in a variety of hidden sector models [63,93].SHiP is a fixed-target experiment, based on a high intensity 400 GeV proton beam dumped on a heavy target.Instead, FASER 2 is an LHC experiment, which aims at exploiting the proton collisions occurring at √ s = 14 TeV during the High-Luminosity LHC (HL-LHC) program.In this kind of experiments, RH neutrinos can be copiously produced by the decay of mesons generated by the proton collisions.The decay of long-lived N 2 particles can then show up in dedicated detectors located around the IP of these experiments.The expected number of signal events can be computed as: where N prod is the total number of N 2 produced, f dec corresponds to the probability for N 2 to decay inside the detector volume, and det accounts both for the efficiency for the reconstruction of the events, that we simply take as 100%, and selection cuts.Finally, • indicates a statistical average, that we define in the following.We consider the production of RH neutrinos from the decay of the following mesons: M = ρ, ω, J/Ψ, Υ.The number of N 2 produced then reads: where N POT is the total number of collected protons on target, N M is the average number of mesons M produced per proton interaction, and BR(M → N 1 N 2 ) is the branching ratio of the decay of the meson M into RH neutrinos, computed in App.B. Similarly, the number of RH neutrinos produced by the decay of mesons at the LHC is given by: where L = 3 ab −1 is the integrated luminosity at the HL-LHC, and σ ine = 79.5 mb is the inelastic proton-proton cross-section [94].The quantity f dec is: where L entry (L exit ) is the distance between the IP where the N 2 particle is produced, and the point at which N 2 enters (exits) the detector.Finally, L N 2 is the decay length of N 2 in the laboratory frame, given by Our calculations are based on simulations of the production of mesons performed with PYTHIA 8.3 and EPOS-LHC [95].More details will be given in the following.From these simulations we obtain a sample of mesons events, and we compute the associated multiplicities N M .Then, for each meson event in the sample, we simulate its decay in N 1 N 2 pairs.These data are used to statistically evaluate Eq. (3.1), averaging ( • ) Eq. (3.4) over all the possible kinematical configurations of the N 2 particles in our sample.Finally, we impose a minimum energy of the photon produced in the decay N 2 → N 1 + γ.We compute the efficiency of this selection cut, det , using the N 2 events in our sample, simulating the N 2 decays, and selecting the events for which the photon energy in the laboratory frame is larger than a threshold E cut .The procedure explained in this section has been used to compute the number of signal events also at the CHARM and NuCal experiments described in Sec. 2.

SHiP
The SHiP fixed-target experiment aims at accumulating N POT = 2 × 10 20 protons on a target composed by Molybdenum and Tungsten in 5 years of operation.A description of the experiment can be found in [67].The decay volume of the detector has a length of 50 m and it is located at ∼45 m from the proton target.A spectrometer and a particle identification system with a rectangular acceptance of 5 × 10 m 2 are placed behind the decay volume.The rectangular face of the decay volume closer to the IP has a size of 1.5 × 4.3 m 2 .Following these specifics, we approximate the detector as a cylinder with an opening angle of 31.8 mrad.
The production of the different mesons at SHiP is based on simulations of protonproton collisions performed with PYTHIA 8.3.For the ρ and ω mesons we obtain the production multiplicities N ρ = 0.58 and N ω = 0.57, in good agreement with previous results present in the literature [91,[96][97][98][99].We assume the same production rate for proton interactions in the target material of SHiP.In principle a dependence on nuclear target is expected, however detailed simulations or measurements are needed to fully capture these effects.Instead for the J/Ψ, we normalize our simulation in order to reproduce the total number of mesons predicted by the SHiP collaboration in [100]: N J/Ψ = 2 × X cc × f (q → J/Ψ)×f cascade .The cc production fraction is X cc = 1.7×10 −3 , the J/Ψ production fraction is f (q → J/Ψ) = 0.01 and the enhancement from cascade events is f cascade = 2.3.Finally, the production rate of the Υ mesons is directly obtained using PYTHIA 8.3, as for the case of the ρ and ω, since detailed simulations of the production at SHiP are not available.We find 11 Currently, there are no studies of the background rates at SHiP for the single photon signature arising in our scenario.By assuming that the backgrounds can be reduced at a negligible level as it happens for other searches, see e.g.[100], the 95% CL upper limit on the number of signal events is N signal = 3.For a more conservative approach we follow [101], which estimated ∼ 1000 background events after rescaling the background events observed at the NOMAD detector by the number of POT in the two experiments.This number will likely be reduced by vetos, as noticed in [101].Assuming Poisson likelihood, we set a 95% CL upper limit on the number of signal events of N signal ∼ 63.8.Finally, we impose a minimum energy of the photon E cut = 1 GeV, which is a reasonable threshold for SHiP [102].Summarizing, in Fig. 1 we show sensitivity contours for the two choices of signal events discussed here.This corresponds to a range between an optimistic and a conservative assumption of the background level.

FASER 2
The FASER collaboration has proposed to build a suite of forward detectors to be placed within the LHC environment along the beam axis, nearby the ATLAS experiment [63].These experiments are dedicated to study several interesting topics, as the properties of neutrinos, QCD in the forward regime and new physics beyond the SM, including dark sectors.Several small size pilot detectors have already been constructed, which are FASER [62,103], FASERν [104,105] and SND@LHC [106].However, to fully exploit the potential of the HL-LHC, a dedicated facility to host larger detectors is under study.In particular, we focus on the proposed FASER 2 detector, which is dedicated to the study of long-lived particles.According to current design, it will be placed at 620 m from the IP, and it will have a cylindrical shape, with a radius of 1 m and a length of 10 m [63].
To simulate the production of the ρ and ω mesons at the LHC, we use EPOS-LHC, which has been tuned to forward LHC data.For the J/Ψ and Υ mesons we again use PYTHIA 8.3, rescaling its rates to match the production cross-sections as measured by the LHCb experiment [107,108].In addition, following [109], we modify the production rate as a function of the transverse momenta p T with respect to the default setup of PYTHIA 8.3.Employing this procedure, we obtain a good agreement with the measured p T distributions [107,108].The resulting production multiplicities in one hemisphere are N ρ = 2.3, N ω = 2.2, N J/Ψ = 5.0 × 10 −4 and N Υ = 6.1 × 10 −6 .In addition to the production of RH neutrinos from the decay of mesons, we include Drell-Yan processes q q → γ/Z → N 1 N 2 .More details are provided in App. A. This production mechanism is however in most cases subdominant with respect the one from mesons decay.FASER 2 will be sensitive to photon signals and, at the same time, it will have the capability to strongly reduce the relevant backgrounds, see the discussion in [110] where single photon signals have been studied in the context of a model of sterile neutrinos coupled to active neutrinos via a dipole operator 12 .Since a thorough simulation of the relevant backgrounds has not been performed yet, we decide to follow the strategy of [110], presenting our results as isocontours of N signal = 3 and N signal = 30 events.A cut on the energy of the photon E cut > 100 GeV has been employed in our analysis, again inspired by [110].

Results
Our main findings are presented in Fig. 1, where we fix the phase to α = π/2 to maximize the production rate from meson decays, see App.B. The results remain qualitatively the same for other choices.We show two different slices of the parameter space: either we fix the mass splitting δ and we explore the plane m N 1 − Λ or we fix the mass of N 1 and we project the results on the δ − Λ plane.In both cases we consider three benchmark scenarios, namely δ = 0.01, 0.1, 1 in the first case and m N 1 = 0.3, 0.6, 1 GeV in the second case.As explained in Sec. 3, the sensitivities of SHiP and FASER 2 are computed for two numbers of signal events, corresponding to different choices of the background rate at these experiments.The strategy followed to compute the regions excluded by CHARM and  NuCal is detailed in Sec. 2. When the line associated to a specific experiment is missing in our plots, this means that the corresponding experiment has not enough sensitivity to probe the parameter space.
As evident in Fig. 1, for a mass splitting δ = 0.1, both SHiP and FASER 2 will be able to extend the current limits from CHARM and NuCal, and probe an uncharted region of the parameter space.In particular the sensitivity of SHiP reaches N 1 masses around the kinematical threshold for production from the decay of the J/Ψ meson, i.e. m N 1 ∼ 1.5 GeV.A more modest sensitivity is obtained for FASER 2. It is worth recalling that despite small values of Λ are formally not excluded in the EFT of Eq. (1.4), weakly coupled UV completions with Λ 100 GeV are likely already ruled out from direct searches of additional EW charged states.The constraints from BBN on N 2 decays are of the order τ N 2 = O(10 −2 − 1) s, see Sec. 2.3.Looking at the isocontour of τ N 2 in Fig. 1, one can notice that these bounds are not overlapping with the sensitivities of SHiP and FASER 2. For a larger mass splitting, δ = 1, the region probed by SHiP tends to shift to larger values of Λ, while FASER 2 can not probe this slice of the parameter space.This behaviour can be understood by recalling that increasing δ tends to reduce the lifetime of N 2 , see Eq. (1.5).This can be compensated by increasing Λ at the price, however, of reducing the production rate of N 1 N 2 pairs.The correlation between δ and Λ can be appreciated in the plots with m N 1 fixed.The different experiments that we have studied are probing proper decay length c τ N 2 ∼ 10 −2 − 10 3 m.
In the case of a smaller mass splitting, as in the case of δ = 0.01, the threshold on the energy of the photon plays an important role.Small δ reduces the energy of the photon, see Eq. (1.8).This implies that at FASER 2 most of the events do not satisfy the cut E cut > 100 GeV and therefore no sensitivity is obtained.To highlight the role of the energy threshold, in Fig. 2 we show the isocontours of N signal = 3 for different values of E cut , namely E cut = 0.1, 0.5, 1, 2, 10 GeV for SHiP and E cut = 10, 50, 100, 200 GeV for FASER 2. While for δ = 1 the sensitivities are almost unchanged, for smaller δ the energy threshold has a significant impact.In particular, for E cut ∼ 10 GeV and provided that background can be kept negligible, FASER 2 will be able to test up to Λ ∼ 400 GeV for δ = 0.01, to be compared with a zero sensitivity scenario with E cut ∼ 100 GeV, shown in Fig. 1.
In addition to FASER 2, other proposed LHC detectors targeting long-lived particles, as ANUBIS, CODEX-b, FACET and MAPP, might potentially be sensitive to the monophoton signature.We estimate the number of signal events expected at these experiments in App. A. In App.C we discuss potential constraints arising from electron recoil searches in fixed-target experiments, finding much weaker sensitivities than those presented in this section.
Finally, before concluding, we shall mention that the RH dipole operator might also be tested by the currently operating e + e − collider experiment Belle II [112], and by the future neutrino experiment DUNE [113].Dedicated analyses are in order to investigate their sensitivities.

Conclusions
In this work we have studied the phenomenological consequences of a dipole operator between RH neutrino fields.This is described by the νSMEFT d = 5 operator N2 σ µν N 1 B µν and triggers the decay N 2 → N 1 γ, which is the subject of our study.Motivated by the current experimental and theoretical interest, we have focused on RH neutrino masses in the GeV range and considered the regime in which N 2 is long-lived, with a proper decay length of O(10 −2 − 10 3 m), while N 1 is considered to be stable on these length scales.
More in details, we have firstly considered the existing bounds on this scenario, arising from terrestrial experiments like CHARM, NuCal and colliders, as well as constraints from cosmological considerations, in particular in relation to the Big Bang Nucleosynthesis epoch.
We have subsequently investigated the sensitivity of the future proposed experiments FASER 2 and SHiP.In these facilities the RH neutrinos are produced in N 1 N 2 pairs through the dipole operator, either via meson decay or via direct production.Then, RH neutrinos give rise to single−γ events through N 2 → N 1 γ decays, which can be detected by these experiments in a background controlled environment.
Our main results are summarized in Fig. 1 where we show that SHiP will be able to probe ample regions of the parameter space not yet excluded by current data, testing Wilson coefficients up to Λ ∼ 10 5 GeV, while the sensitivity of FASER 2 is more limited.Given the early design stage at which these experiments are, and the preliminary nature of the background estimates for the scenario under consideration, we have then studied how different cuts on the photon energy enforced at the analysis level affect the sensitivity reach.Our results are shown in Fig. 2. We found that relaxing the cut on the photon energy has a limited impact on the sensitivities predicted for SHiP, while for FASER 2 ampler regions of parameter space can be reached, provided that the background can be maintained at a negligible level.Finally, in Fig. 3 in the appendix, we also present results on the number of signal events expected at other future LHC experiments, namely ANUBIS, CODEXb, FACET, MAPP.
In conclusion, our work provides a first realistic estimate on the reach of experiments targeting long-lived particles on the lowest dimensional effective dipole operator that appears in the minimal see-saw extension of the Standard Model.A Projected sensitivity of other future LHC experiments In addition to FASER, several other LHC detectors dedicated to search for long-lived particles have been proposed in recent years: MATHUSLA [114,115], CODEX-b [58][59][60], AL3X [116], MAPP [64,65], ANUBIS [57] and FACET [61].These facilities, to be placed around the LHC IP points, could potentially probe the radiative decay of the RH neutrinos that we are considering.Concretely, we focus on CODEX-b, ANUBIS, MAPP and FACET, that, in principle, can have the potential to reconstruct photons (for CODEX-b this assumes an extension of the baseline design, according to [59]) 13 .For the single photon signature under scrutiny, the background rates at these experiments have not been computed, and a detailed discussion of their capability to reduce the relevant backgrounds is not currently available in the literature.Given this limited information, the estimate of their sensitivity reach is quite uncertain.We show the region of the parameter space where the following numbers of signal events at these experiments are obtained: N signal = 3, 10, 100, 1000.These results are intended to give an idea of potential sensitivity reach of these proposals, if the backgrounds are reduced to the appropriate rates.
For the calculation of the signal rate we follow the same procedure described in Sec. 3. The geometry of the different experiments and the cut on the photon energy E cut are the same adopted in [117].We consider two mechanisms for the production of the RH neutrinos: meson decays, see Sec. 3.2 for details, and Drell-Yan processes.For the latter, we employ MadGraph5 aMCNLO for our simulations.A couple of comments are in order.As discussed in Sec. 2, we should require that the energy scale of this process is smaller than the cut-off of the EFT.Assuming a weakly coupled extension of the SM and couplings of O(1), we impose √ ŝ < Λ, where ŝ is the Mandelstam variable associated to the process q q → γ/Z → N 1 N 2 .In practice, from our MadGraph5 aMCNLO simulations, we select only the events satisfying this condition.In addition, we impose √ ŝ > 2 GeV, in order to work in the regime of perturbative QCD.
The results are shown in Fig. 3 for two different slices of the parameter space: fixing the mass splitting to δ = 0.1 or fixing the mass of the lightest RH neutrino to m N 1 = 0.3 GeV.For some experiments, FACET and ANUBIS, N signal > 10 2 − 10 3 can be obtained in some parts of the parameter space, while more modest signal rates are obtained in other cases, as for CODEX-b and MAPP.
In general, we find that for the forward detector FACET and FASER 2, the production from meson decays is more relevant than the one from Drell-Yan processes.The opposite situation happens for the off-axis detectors ANUBIS and CODEX-b.In some cases the sensitivity disappears at small Λ, see the bottom left panel of Fig. 3 or the flattening of the curves in the other panels.This is due to the requirement √ ŝ < Λ that we impose in our simulation: for small enough Λ most of the events in the simulation are rejected.

B Mesons decay into RH neutrinos
For the decay V → N 1 N 2 we need the following matrix element: where m V is the vector meson V mass, µ (p) its polarization vector and explicit expressions for the coefficients f q V can be found in Appendix A of [98].Given the range of masses to which we are interested, in our computation we will consider only photon exchange.The explicit expression for the decay width is given by where c w is the cosine of the weak angle, Q q the electric charge of quark q, in units of the electron's electric charge e.To produce our plots we set α = π/2 to maximize the number of events, although the results remain qualitatively the same for other choices of the phase.

C Electron recoil searches
The dipole operator of Eq. (1.4) induces an inelastic scattering processes between RH neutrinos and electrons, namely • N 21 : an N 2 particle produced by mesons decays produces an e-recoil signal in the detector through N 2 e − → N 1 e − scattering; • N 212 : an N 1 particle produced by N 2 decay produces an e-recoil signal in the detector through N 1 e − → N 2 e − scattering.
Each term has been evaluated through a Montecarlo simulation of the process.The N 1,2 neutrinos have been assumed to be produced from meson decays, and the fluxes of mesons have been simulated with PYTHIA 8.3, as explained in the main text.The electron number density of the detectors has been obtained from the weight and the material of each experimental apparatus.
We have computed the differential cross section with respect to the recoil energy of the inelastic scattering processes, assuming initial electrons at rest.Explicitly We find the number of recoil events to be negligible for values of Λ above 1 GeV, independently of the values of m N 1 and δ.We then conclude that searches through electron recoils do not impose relevant constraint on the νSMEFT parameter space.In the limiting case δ = 0, we reproduce existing results available in the literature, which have focused on the case of elastic scattering processes [91].

Figure 1 .
Figure 1. Green and blue lines are the sensitivity reach of the SHiP and FASER 2 experiments.For SHiP we show isocontours of N signal = 3 and N signal = 63.8.For FASER 2 we show isocontours of N signal = 3 and N signal = 63.8, see Sec. 3 for more details.The orange and magenta shaded regions are excluded by the CHARM and NuCal experiments.The dashed lines are contours of constant N 2 lifetime or proper decay length.We fix α = π/2.

N 2 e
− → N 1 e − (C.1) and N 1 e − → N 2 e − .(C.2) Note that the latter process is kinematically open only if the center of mass energy is sufficiently large, given that m N 2 ≥ m N 1 .Such processes can give rise to a signal in experiments sensitive to e-recoils.We have considered experimental searches at SHiP, CHARM II and DUNE and estimated present and future constraints.These are fixedtarget experiments, whose number of POT, detection efficiencies and cuts enforced in the analysis are recollected in Tab. 1.The expected number of signal events is the sum of three contributions N sig = N 12 + N 21 + N 212 , (C.3) which are • N 12 : an N 1 particle produced by mesons decays produces an e-recoil signal in the detector through N 1 e − → N 2 e − scattering;

Table 1 .
[66,118]ments POT(10 20) E R ∈[1,20]GeV, θ R ∈[10,20]mrad[66,118]Summary of the main characteristics of the experiments that we considered.Here E R and θ R are the recoil energy and the recoil angle with respect to the incoming neutrino's momentum of the scattered electrons, while eff is the detection efficiency of the signal.