The see-saw portal at future Higgs factories: the role of dimension six operators

We study an extension of the Standard Model with electroweak scale right-handed singlet fermions $N$ that induces neutrino masses, plus a generic new physics sector at a higher scale $\Lambda$. The latter is parametrized in terms of effective operators in the language of the $\nu$SMEFT. We study its phenomenology considering operators up to $d=6$, where additional production and decay modes for $N$ are present in addition to those arising from the mixing with the active neutrinos. We focus on the production with four-Fermi operators and we identify the most relevant additional decay modes to be $N\to \nu \gamma$ and $N\to 3f$. We assess the sensitivity of future Higgs factories on the $\nu$SMEFT in regions of the parameter space where the new states decay promptly, displaced or are stable on detector lengths. We show that new physics scale up to $5-60\;$TeV can be explored, depending on the collider considered.


Introduction
Neutrino masses and mixing can be explained by adding to the Standard Model (SM) a new Weyl fermion N , total singlet under the SM gauge group, which acts as the right-handed (RH) counterpart of the left-handed (LH) SM neutrino. The lightness of the neutrino masses can be explained by the see-saw mechanism [1][2][3][4] where v is the Higgs vacuum expectation value (VEV), y the strength of the Dirac type interaction and M N the Majorana mass of the RH neutrino. While there is no indications on the energy scale at which this mechanism takes place, there is nowadays a strong interest in models where RH neutrinos have a mass at the EW scale. From one side they are in fact an interesting alternative, in that they can generate the observed matter-antimatter asymmetry via neutrino oscillations [5,6], while they can be searched for at colliders and at beam-dump experiments [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Moreover, if lighter than O(100) MeV, they can be relevant for the solution of longstanding anomalies reported in the neutral-current [26][27][28] and charged-current [29][30][31][32][33][34][35][36][37] semileptonic decay of B mesons [38][39][40]. Their phenomenology is driven by the mixing θ with the active neutrinos which drives their production rates and and their decay width and hence their lifetime.
The naive see-saw scaling of Eq. (1.2) can be modified if multiple RH neutrinos are present with specific Yukawa and Majorana mass textures that ensure an approximate lepton number symmetry [41,42]. Consequently scenarios with a much larger mixing angles can be realized, thus altering the RH neutrinos phenomenology. It's also interesting to speculate on the possibile presence of additional NP states at a scale Λ v, M N , whose effects can be parametrized in the language of effective field theories in the so called νSMEFT, where a tower of higher-dimensional operators O d Λ 4−d built out with the SM fields and the RH neutrinos is added to the renormalizable lagrangian. At the lowest possible dimension, d = 5, there are two genuine, i.e. that contains at least a RH neutrino field, νSMEFT operators: one that triggers the decay of the SM Higgs boson into a pair of RH neutrinos and a dipole operator with the hypercharge gauge boson [43,44]. Already at d = 6 many more operators are present [9,45,46] with interesting phenomenological consequences, since they also can induce new production and decay channels.
Many of these operators have been subject of theoretical studies, especially for what concerns their phenomenology at the Large Hadron Collider (LHC), see e.g. [9,21,43,44,[47][48][49][50][51][52][53]. However, RH neutrinos with a mass at the EW scale are one of the primary goals of future lepton colliders, since the generally small production cross section proper of EW singlet states can be overcome, thanks to the clean detector environments and the typically lower SM backgrounds with respect to hadronic machines. For the post LHC era many future lepton colliders have been proposed. These includes e + e − facilities, both circular ones as the Future Circular Collider [54][55][56][57] (FCC-ee) and the Compact electron-positron collider [58,59] (CEPC), and linear ones as the International Linear Collider [60][61][62] (ILC) and the Compact Linear Collider [63,64] (CLIC). Finally, a great attention has recently arose for multi TeV µ + µ − colliders [65], which could provide a great handle to test higherdimensional operators whose effect grows with energy.
In a recent paper [66], we have investigated the prospects of these machines, commonly denoted as Higgs factories, in testing the two genuine d = 5 operators of the νSMEFT through Higgs and Z boson physics and focusing on RH neutrinos with masses in the [1 GeV − m h,Z ] range. There we have shown that future lepton colliders can test exotic branching ratios (BRs) for the Higgs and Z boson down to ∼ 10 −3 and 10 −9 respectively, greatly surpassing the reach of future indirect measurements of the Higgs and Z boson width. In this paper we extend our previous work by studying the phenomenology of the νSMEFT operators that arise at d = 6. Since these are typically generated by different ultraviolet (UV) completions than the d = 5 ones, the bounds on the cut off-scale Λ derived in [66] do not necessarily direct apply 1 .
We focus on EW scale RH neutrinos with masses in the [1 GeV − m W ] range and study the additional production and decay channels induced by the d = 6 operators. We distinguish two main decay channels: a two body decay into a SM neutrino and a photon, N → νγ, and a three body decay into a SM leptons and a fermion pair which can proceed either as N → νff or N → ff , where = e, µ, τ . In the three body decay cases the final state fermions could involve either a pair of quarks or leptons. For what concerns the production, we identify the most relevant channels as single-production and pair-production of RH neutrinos induced by four-Fermi d = 6 operators, since they induce amplitudes that grow with the energy of the process.
The paper is organized as follows. In Sec. 2 we set our notation and review the νSMEFT framework, while in Sec. 3 we present the properties of the future colliders which are under analysis. Then in Sec. 4 we study the main decay channels induced by the d = 6 operators and present the expressions for the various partial widths. We then show under which conditions these additional decay modes can dominate with respect to the one already present at renormalizable level and induced by the active-sterile mixing. We further quantify the lifetime of the RH neutrinos once these operators are switched on. In Sec. 5 we discuss the additional production modes relevant for studies at future Higgs factories. We present our results in Sec. 6, Sec. 7 and Sec. 8 for prompt, displaced and detector stable RH neutrinos. We finally conclude in Sec. 9. We then report in App. A the expressions for the spin averaged matrix elements squared relevant for the N three-body decay via an off-shell SM boson induced by d > 4 operators.

Theoretical framework
The νSMEFT is described by the following Lagrangian Here N is a vector describing N flavors of RH neutrino fields, singlet under the SM gauge group, and N c = CN T , with C = iγ 2 γ 0 . Furthermore L L is the SM lepton doublet, Y ν is the 3×N Yukawa matrix of the neutrino sector withH = iσ 2 H * , M N is a N ×N 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 Λ indicating the EFT cut-off. In [9,[43][44][45][46] the νSMEFT has been built up to d = 7 and at d = 5 only three operators exists. The first is the well-know Weinberg operator [69] while the two

four-Fermi scalar operators
Other four-Fermi operators where B µν is the SM hypercharge field strength tensor and 2σ µν = i[γ µ , γ ν ], which have been recently investigated in [66]. At d = 6 many more operators are present. They are reported in Tab. 1, where we split them between operators that involve the Higgs boson and four-Fermi operators which do not 2 .

Neutrino mixing formalism
We summarize here the properties of the neutrino sector of the νSMEFT, and we refer the reader to [66] for a more detailed discussion. Under the approximation in which the contribution to the active neutrino masses dominates over the ones induced by the effective operators the active neutrino mass matrix takes the standard form where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [70,71] and m d ν is the diagonal matrix of neutrino masses. Eq. (2.2) can be solved for the Yukawa matrix of the neutrino sector. In the Casa-Ibarra parametrization [72] one obtains where √ m is a 3×N matrix containing the physical neutrino masses m i and R is a complex orthogonal N × N matrix. We restrict now our study to the case N = 2 where for the normal hierarchy (NH) and inverted hierarchy (IH) one has for the matrix while we parametrize the orthogonal matrix R in terms of the complex angle z = β + iγ as The active-sterile mixing angle is given by It's crucial that the angle z is, in general, a complex parameter. In fact, in the limit in which z is a real number, by taking U and R with entries of order unity and by assuming an equal value for the diagonal entries of the Majorana mass term for the two RH neutrino This relation is drastically modified by the imaginary part of z, that gives an exponential enhancement. In the limit γ 1 one has (2.9) Clearly the same enhancement is inherited by the active-sterile mixing, that now reads (2.10) We use α = 1, 2, 3 for the active neutrino flavor and i = 1, 2 for the RH neutrino flavor. This deviation from the naive see-saw scaling has a crucial impact on the RH neutrinos phenomenology, especially for what concerns their decay width and consequently their lifetime, with drastic implications for search strategies at future colliders as recently shown in [66,73]. 3 When needed, for our numerical estimate we take mν 2 = 8.6 × 10 −3 eV and mν 3 = 5.1 × 10 −2 eV for the NH while we take mν 1 = 4.9 × 10 −2 eV and mν 2 = 5.0 × 10 −2 eV for the IH. 4 We have assumed NH and fixed mν = mν 3 . The expression holds also for the IH case modulo order one factors.

Future Higgs factories
In this work we study the phenomenology of the νSMEFT at future Higgs factories, both at their low energy runs, relevant for physics at the Higgs-strahlung threshold and at the Z pole, together with high-energy multi TeV runs, which can greatly enhance the sensitivity on higher-dimensional operators which induce a quadratic grow with the energy, as for the case of four-Fermi operators. For what concerns the low energy runs, various e + e − prototypes, presently at different stages of their design, have been proposed. These include circular ones, as the Future Circular Collider (FCC-ee) [54][55][56][57] and the Circular Electron Positron Collider (CEPC) [58,59], and linear ones, as the International Linear Collider (ILC) [60][61][62] and the Compact Linear Collider (CLIC) [63,64]. Regarding colliders in the multi TeV regime, prototypes include CLIC with a center of mass energy of 3 TeV [63,64] and a µ + µ − colliders with various center of mass and luminosity options [65]. We report in Tab. 2 the main parameters of these colliders prototypes. For concreteness and clarity of presentation, in this work we will present our results only for the case of FCC-ee at √ s = m Z and √ s = 240 GeV, for a µµ collider at √ s = 3 TeV and for CLIC at 3 TeV, since the ILC and CEPC prototypes will have an overall similar behavior to the considered options.

Decay channels for RH neutrinos
At the renormalizable level, the RH neutrinos only decay thanks to the mixing with their SM active counterpart, a pattern which is not altered by the inclusion of d = 5 operators except in the case of a sufficiently large mass splitting m N 2 − m N 1 (with m N 2 > m N 1 ) in which the O 5 N B operator can trigger a non-negligible N 2 → N 1 γ decay rate [66]. The inclusion of d = 6 operators can dramatically alter this behavior, leading to new decay patterns, see also [74,75]. For example, the four-Fermi operators reported in Tab. 1 can induce the decay of a RH neutrino into three SM fermions, N → 3f . Depending on the operator, the rate for this decay may or may not be suppressed by the active-sterile mixing Operator Decay Mixing Loop Table 3. Decay modes for the RH neutrinos induced by higher-dimensional operators and renormalizable mixing. We highlight whether the corresponding rates are mixing and/or loop suppressed. Neutral and charged indicate the four-Fermi operators with two and one RH neutrino respectively.
angle. In particular, it is suppressed in the case of four-Fermi operators which contain two RH neutrino fields, while unsuppressed otherwise. On the other side, the operators involving the Higgs boson can induce, after EW symmetry breaking, the decay into a final state fermion and a massive SM boson, B = Z, W ± , h. Given the range of RH neutrino masses on which we are interested in, the SM boson turns out to be off shell and the resulting decays are thus N → νB * and N → B * , with the subsequent decay B * → ff . Also in this case the rate could be, or not, suppressed by the active-sterile mixing and again the final state is composed by three SM fermions, as for the case of the four-Fermi operators. However the kinematic and the flavor composition is generally different. Finally, the SM boson could be a massless photon, and the decay is thus simply N → νγ. We now discuss the various operators and the decay that they mediate in turn, summarizing their main properties in Tab. 3, where we highlight whether the decays that they generate are suppressed by mixing and/or loop effects. We then report in App. A the spin averaged amplitude squares for the considered three body decays.

Operators that induce N → νγ
This decay mode is induced at the d = 5 level by the O 5 N B operator and at the d = 6 level by the O 6 LN B,LN W operators.
This operator gives the N i → ν α γ decay only after a mixing insertion 5 . The rate for this decay reads [44] where c ω is the cosine of the Weinberg angle and where we have explicitly introduced a loop suppression factor, since in any weakly coupled UV completion this operator arises at loop level [76,77].
where, again, we have explicitly written the loop suppression factor. This decay is not suppressed by the active-sterile mixing.

Operators that induce N → 3f
This decay mode has contribution from both operators involving the Higgs boson as well as four-Fermi operators. Also the decay induced at the renormalizable level by the activesterile mixing produces the same final state.
where the factor N c = 3 is present if the final state is a quark-antiquark pair. In our numerical analysis we use the full expression for the decay and sum over the relevant ff final states for any N mass. For the decay into a bb final state, which is the relevant one for m N > 10 GeV, one has where the Wilson coefficient has been fixed to one.

N H
This operator induces the decay N i → ν α Z * and the rate is suppressed by a mixing insertion. For a generic ff final state arising from the Z * decay one has, in the limit m f = 0 with t 3 = ±1/2 and where q is the electric charge of the final state fermion pair. For example for the decay into a final state bottom pair, where t 3 = −1/2 and q = −1/3 one has where θ schematically indicates the relevant mixing angle.

Decay from O 6 N eH
This operator induces the mixing unsuppressed decay N i → α W * . Working again in the limit where all the final state fermions are massless, the decay for one of the two charge conjugate modes reads where the estimate is with N c = 3 and an O(1) Wilson coefficient.
The combination of these two operators orthogonal to the one that induces N → νγ gives again a N → νZ * decay. In addition, the operator with the W boson produces a N → W * decay. Both these rates are not suppressed by the active-sterile mixing. In the massless limit, the neutral decay width is while for the charged case we obtain where, again, the rate is for one of the two charged conjugated modes.
Decay from four-Fermi vector operators: The first operators are of the form (N R γ µ N R )(f L/R γ µ f L/R ). They mediate the decay N → νff which is suppressed by the active-sterile mixing angle. For simplicity we assume a diagonal flavor structure for the SM fermion pair final state. In the limit of massless final states the decay rate reads where α 6 is the Wilson coefficient of the four-Fermi operator and the factor 2 comes from summing over ν andν, since also the SM neutrino is Majorana. The charged operator O 6 N edu triggers a decay N → − ud + +ū d, which has a rate (4.12)

Decay from four-Fermi scalar operators
These operators induce the decay N i → α ff , where could be a charged or neutral lepton. Each decay proceeds with a rate where, once more, α 6 denotes the generic Wilson coefficient of the four-Fermi operator.

Which decay dominates?
We can now compare the decay rates computed in Secs. 4.1 and 4.2 to see which one dominates in the different regions of the νSMEFT parameter space. This is essentially determined by three parameters 6 : the mass of the RH neutrino m N , the active-sterile mixing θ and the EFT cut-off scale Λ. We take the latter to be the same for all the considered operators. Clearly, different UV completions will generate at low energy different operators, in general suppressed by different mass scales. We will comment in Sec. 4.4 and Sec. 4.5 on the independent limits on the scale Λ that can be set for the most relevant ones. For simplicity however, in performing our main analysis, we will assume that only four fermi operators and the dipole ones triggering the N → νγ decay are active, and that they are associated to a unique scale Λ.
The first question we want to address is in which region of the parameter space the decay induced at the renormalizable level by the active-sterile mixing dominates over the decay generated by higher-dimensional operators, taking into account that current constraints forces the squared mixing angle to be 10 −6 [79][80][81]. In order to do this we need to make some assumptions on the number of four-Fermi operators that are active, since each one can contribute with a multiplicity due to the flavor structure of the operator itself. To be practical, we parametrize this with a coefficient ξ which takes into account how many channels from four-Fermi operators contribute to the decay of a RH neutrino, for example for a decay into a pair of final state quarks ξ = N c = 3. Clearly, the most important four-Fermi operators for N decay are the ones that do not pay a mixing suppression, i.e. O 6 N edu and all the scalar ones. On the other side, for the operators that contribute to the 3f final state via an off-shell h, Z and W , we can consider all possible decays by summing on their decay modes, since those are fixed by the SM symmetries. In these calculation we retain the full expressions for the various decay rates. We then show in Fig. 1 the region of parameter space where the decay induced by the higher-dimensional operators dominate over the one induced by the mixing. We illustrate them for the degenerate RH neutrino masses m N = 1 GeV (left) and m N = 50 GeV (right) 7 . Above the black solid line the decay pattern is thus the one analyzed in [66], while effects from higher-dimensional operators become relevant in the lower part of the plot. The dashed gray lines indicate the experimental bound on θ α = i=1,2 |θ αi | 2 reported in [79][80][81]. We show only the bound on 6 While we already stated that we work in the limit where the two RH neutrinos are almost degenerate and the various entries of the active-sterile mixing matrix θ are determined by the choice of the NH or IH once the RH neutrino mass has been fixed, each higher-dimensional operator can in principle have a different Wilson coefficient. For concreteness, we work under the assumption that they are all equal. We also consider the NH scenario. Results in the IH case are almost identical. 7 The dependence on ξ turns out to be completely negligible for the N mass range of our interest up to |θ µ | 2 which turns out to be the most stringent one. Finally the gray shaded area represents the see-saw limit, below which the lightness of the neutrino masses cannot be explained by the see-saw mechanism. As we see, for small enough Λ the dominant decay mode of the RH neutrino can be induced by the higher-dimensional operators of the νSMEFT while retaining compatibility with existent active-sterile mixing bounds. As previously discussed, in this region two decay modes compete: N → 3f , which produces the same final state as the decay via mixing albeit with different kinematics, and N → νγ. For |θ e | 2 10 −6 one has that the ratio Γ N →νγ Γ N →3f is almost independent on Λ and that the νγ decay dominates over the N → 3f decay for m N 15 GeV . In this region the decay is driven by the d = 6 operators O 6 LN B,LN W , since the d = 5 operator O N B gives a rate which is mixing suppressed. For larger masses, the operator that dominates the N → 3f decays is O 6 N eH , which is again not mixing suppressed. Given that we are interested in the phenomenology of the d > 4 operators in the following we will focus in the region where the decay is dominated by higher-dimensional operators and work under the assumption of negligible active-sterile mixing.

Bounds from theoretical considerations
The computation of the neutrino properties outlined in Sec. 2 rests on the assumption that the d = 4 masses and Yukawa couplings dominate over the higher dimensional contributions. In order for this to be true, the NP scale will have to satisfy some conditions. Before enumerating them, it is useful to point out that, unlike what happens in the SMEFT, the νSMEFT is characterized by two expansion parameters: the active-sterile mixing θ and the cutoff scale Λ −1 . As previously discussed and shown in Fig. 1, the phenomenology will strongly depend on the interplay between the two. In order to understand the stability of the d = 4 parameters against the additional contributions, we will consider only those effects that solely depend on Λ, neglecting possible effects that are doubly suppressed by some power of θ and of 1/Λ. where the reference value |θ| 2 ∼ 10 −6 is the approximate experimental upper bound on the mixing angle for the RH neutrino mass range of our interest. As we can see, the theory bound on the scale of O 6 LN H is pretty strong, while the one on O 6 LN B /O 6 LN W is rather weak, at least for values of the mixing close to the allowed upper bound. In order for the bounds on the scale of these operators to be of the TeV order we would need |θ| ∼ 10 −14 , which is below the see-saw limit, see Fig. 1.

Bounds from precision measurements
The operators of Tab. 1 involving the Higgs bosons will also trigger additional decay of the SM bosons, which are constrained by precision measurements from LEP and LHC data. By asking that these additional decay modes do not exceed the absolute uncertainty on the measurement of the Z and W boson width, and that they contribute less than 10% to the SM Higgs boson width, one obtains that the strongest limit arises from the constraints on h → N N decay given by O 6 LN H that reads a result compatible with the one reported in [50]. This is due to the small total width of the SM Higgs boson, which compensates for the lower absolute precision on its determination with respect to the Z and W cases. For the latter we obtain a bound of Λ 0.8 TeV 8 We have explicitly checked that the contribution from the operators O 6 LN qu , O 6 LN qd and O 6 LdqN give weaker bounds with respect to the ones shown. and 0.6 TeV, respectively. While for O 6 LN H the theoretical bound discussed in the previous section is stronger, for the dipole operators the experimental bounds are stronger.
The interplay between the O 6 LN B,LN W operator and the active-sterile mixing also generates a magnetic moment d µν σ αβ νF αβ for the SM neutrinos which can be estimated as This is another example of effect which is suppressed by both θ and powers of 1/Λ. The value of the active-sterile dipole moments constrained by reactor, accelerator and solar neutrino data [82,83] which give Λ 4 × 10 −2 |θ| 2 10 −10 TeV.
that is weaker than Λ 1 TeV for the allowed mixing angles range.

Lifetime of RH neutrinos
After having discussed the main RH neutrinos decay modes, it is important to determine the lifetime of these state, to assess whether their decay happen with a prompt or displaced behaviour or if instead they are stable on collider lengths. We quantify the three behavior as follows: Prompt decay We consider a RH neutrino to decay promptly if its decay happens within ∼ 0.1 cm from the primary vertex. At the renormalizable level, prompt RH neutrino decays require a large breaking of the naive see-saw scaling. In the notation of Sec. 2.1, this is parametrized by a large value of the γ parameter, see Eq. (2.10). Large mixing angles are however constrained by a variety of experimental searches, and too large values of γ are thus ruled out.
Displaced decay A particle is considered to decay displaced if it decays away from the primary vertex but within the detector environment. The precise distance for defining a vertex to be displaced clearly depends on the specific detector geometry. Given that our study focuses on future proposed e + e − and µ + µ − colliders, for which detailed detector characteristics have not yet been settled, we consider as displaced particles decaying between 0.1 cm and 1 m from the primary vertex. Taking into account the preliminary nature of our study, we also consider the detector to have a spherical symmetry, instead of a cylindrical one.
Decays outside the detector Also in this case, the precise value of the decay length of the RH neutrinos needed for it to be considered detector stable depends on the specific geometry of the detector. We then consider as detector stable, RH neutrinos which decay more than 5 m away from the primary vertex. The decay length in the laboratory frame βγcτ can be readily obtained for the two dominant N production modes that will be discussed in Sec. 5, i.e. pair-production and singleproduction from four-Fermi operators. The βγ factor is fixed by the kinematic of the process and reads (4.20) As discussed in the previous section, in the region where the RH decay width is dominated by the d > 4 operators, two decays compete: N → νγ and N → 3f . As an example, we show in Fig. 2 the isocontours of βγcτ for the case of exclusive νγ (left) and 3f (right) decay, fixing √ s = 240 GeV and 3 TeV and considering the pair-production case. These lifetimes are dominated by mixing unsuppressed operators and thus do not strongly depend on the mixing angle. As in Sec. 4.3, the dependence on ξ is extremely mild. The case of single-production is qualitatively similar, with more pronounced differences appearing for large m N in the case √ s = m Z . From the figures we see that the RH neutrino can have, for both final states, a prompt, displaced and stable behaviour, depending on the values of m N and Λ considered, although a detector stable N can only arise for m N 20 GeV for Λ 100 TeV. Clearly, if one considers only the decay induced by mixing suppressed operators these will in general give larger values for the proper cτ decay length, which are compatible with a displaced or stable behavior for N and that can be of the same order of magnitude as the one induced by the active-sterile mixing.

Production modes for RH neutrinos
At the renormalizable level, RH neutrinos are produced only via their mixing with the active neutrinos, while at d = 5 two different production mechanisms arise: one from an exotic decay of the Higgs boson and one from the exotic decay of the Z boson. These have been studied in [66], where the N were considered to decay only via mixing, being this the dominant mechanism for d ≤ 5. The inclusion of d = 6 operators brings new production modes for RH neutrinos. The main mechanisms can be divided in two categories.
i ) Single-and pair-production of N via four-Fermi operators, ii ) N production via Z, W and h decay from d = 6 operators involving the Higgs boson.
In this work we focus on production via four-Fermi operators while we leave the analysis of the production from SM boson decay for future work.

Single and pair-production of N via four-Fermi operators
At lepton colliders there are three four-Fermi operators that can produce RH neutrinos. The O 6 N e and O 6 N L operators generate the process + − → N i N j with a rate where the numerical approximation is valid in the massless limit. In both cases we have set to unity the Wilson coefficient of the operator inducing the process and assumed fixed flavors. Appropriate multiciplicity factors must be included to compute the inclusive crosssections in all flavors.
As a preliminary indication, we can ask what is the maximum scale Λ that can be tested by requiring the production of at least one signal event before enforcing any BR factor and selection acceptance. As mentioned in Sec. 3, we take as benchmark colliders the FCC-ee at √ s = 240 GeV, the FCC-ee at the Z pole, a µµ collider with √ s = 3 TeV and CLIC at √ s = 3 TeV. For all these options, the considered integrated luminosities are reported in Tab. 2. The maximum scales that can be tested are show in Fig. 3, where the left and right panel are for N pair-and single-production respectively. By comparing this result with Fig. 1 we see that, for light N in the majority of the allowed parameter space that can be tested, the decay of the RH neutrino will proceed via higher-dimensional operators while for heavier N the decay might also proceed via active-sterile mixing.
Even if produced via a four-Fermi operator, the heavy N can nevertheless decay into a γν final state. For instance, four-Fermi operators of the form (N γ µ N )(f γ µ f ) will induce an unsuppressed pair-production cross-section e + e − → N N and a decay N → ν ¯ which, being mixing-suppressed, will typically be subdominant. In addition, this will also always In order to be concrete, we thus analyze the two possible signatures in turn separately, assuming a 100% exclusive decay for each mode and separately considering the possibility of a prompt, displaced and collider stable behavior.

N prompt decay
As shown in Fig. 2 the RH neutrino can promptly decay into a νγ and 3f final state in all the N mass range of our interest if Λ is sufficiently small. We start by considering the exclusive N → νγ decay, moving then to the N → 3f one for both N single-and pair-production.

Decay N → νγ
When the dominant decay mode is the one into a SM neutrino and a photon we consider the following processes for pair-and single-production of N  4. 95% CL exclusion limit for the prompt decay into νγ for N pair-production (left) and single-production (right) for various collider options. Also indicated is the region where the decay cannot be prompt so that the described analysis doesn't apply. See text for more details.
In the case of N pair-production, Eq. (6.1), the final state consists in a pair of γ and / E T . Two operators can mediate the N pair-production: O 6 N e,N L whose cross section is reported, for each process, in Eq. (5.1). For simplicity, and being conservative, we assume that only one of the two operators is present and only one pair of RH neutrino is produced. When the RH neutrino is singly produced, Eq. (6.2), the final state consist of a single photon and / E T . Only one operator can mediated this process, O 6 N eN L , whose cross-section is reported in Eq. (5.2). We have implemented 9 the relevant higher-dimensional operators in the Feynrules package [84] and exported it under the UFO format [85] in order to generate parton level signal events with MadGraph5 aMCNLO [86]. Events has been then analysed with the MadAnalysis5 package [87][88][89]. The irreducible SM backgrounds + − → γγ / E T and + − → γ / E T have been generated with the same prescription. At the analysis level, we require the photon to be reconstructed with |η γ | < 2.44 and, for the pair-production case, that they are separated as ∆R(γγ) > 0.1. We enforce the following cut on the photon(s): in the pair-production case we apply p γ T > 80 GeV, 20 GeV and 300 GeV for the FCC-ee at √ s = 240 GeV, the FCC-ee at √ s = m Z and CLIC and the µµ collider at 3 TeV respectively. In the single-production case we apply instead p γ T > 20 GeV, 20 GeV and 300 GeV for the same three collider options. The statistical significance is evaluated in units of standard deviations as S/ √ S + B where S and B are the final number of signal and background events respectively.
We then show in Fig. 4 the the 95% confidence level (CL) exclusion contours for the four collider options for the pair-production (left) and single-production (right) cases respectively. In the figures the gray shaded area is the region with βγcτ > 0.1 cm, that is where the RH neutrinos do not decay promptly and the analysis doesn't apply. This region is conservatively shown for √ s = 3 TeV and is smaller for lower collider energies, see Fig. 2. In the pair-production case we observe that the FCC-ee running at the Z mass has a higher sensitivity to this scenario with respect to the FCC-ee running at √ s = 240 GeV, thanks to the higher integrated luminosity of the first option. In the region where the prompt analysis applies, the bound reaches its maximum at around m N ∼ 30 GeV, then depleting at the mass threshold for N pair-production, where the 240 GeV run of FCC-ee will retain a sensitivity up to Λ ∼ 5 TeV. Note that the bound on Λ from Higgs precision measurements, see Eq. (4.18), partially covers these regions if the O 6 LN H operator is switched on. On the other side a µµ collider running at √ s = 3 TeV will be able to test in principle up to Λ ∼ 20 TeV, while CLIC at the same center of mass energy will be able to test scales up to Λ ∼ 25 TeV. However only lower scales will be effectively tested by this analysis since for higher values of Λ the RH neutrinos will not decay promptly. We also note that the reach is dramatically reduced with respect to the maximal one, left panel of Fig. 3, due to the non negligible SM background for this process. On this respect the limits obtained in Fig. 4 can however be considered as conservative and can be improved by dedicated background treatment and reduction, thus increasing the overall reach on Λ in a realistic analysis. The results in the single-production case are qualitatively similar, albeit slightly weaker, with respect to the pair-production scenario, due to the higher rate for the SM background.

Decay N → 3f
When the dominant decay mode is the one into three SM fermions, we consider the following processes for pair-and single-production of N : where the fermion final state could also include quarks. These final state are similar to the one that arises by singly or pair-produced N that decay via mixing, albeit with a different kinematics 10 . For the pair-production case we focus on the following process with a pair of same-sign (SS) leptons, which is expected to be particularly clean where the four quarks arise from the virtual W decay and can be in any flavor combination. As for the SM background, we follow the same procedure of [66] and compute the SM background + − → + − 4q, correcting it for a (flat) lepton charge misidentification probability factor of misID = 10 −3 [90], e.g. we compute the background 10 In practice, we consider a scenario where the decay is triggered by the O 6 N eH operator, which mediate N → W * . Not being mixing nor loop suppressed, this decay is the dominant one even when the O 6

LN Le
operator that mediate single-production and that can trigger N → ν is switched on. Figure 5. 95% CL exclusion limit for the prompt decay into 3f for N pair-production for various collider options. Also indicated is the region where the decay cannot be prompt so that the described analysis doesn't apply. See text for more details.
yield as σ + − 4q × 2 × misID (1 − misID ). At the analysis level, we require p T > 2.5 GeV, p j T > 5 GeV, |η | < 2.44, |η j | < 2.4 and ∆R > 0.1 11 between the two leptons and a lepton and jet pair 12 . We furthermore consider the correct mass dependent SS branching ratio from the N decay induced by O 6 N eH . We thus obtain the 95% CL exclusion limit shown in the left panel of Fig. 5, where we see that the FCC-ee will be able to test roughly Λ ∼ 5 TeV in the whole considered N mass range for both runs at the Z pole mass and at √ s = 240 GeV, while the high-energy colliders will be able to test up to Λ ∼ 20 − 25 TeV, although only in a smaller region the RH neutrino will decay promptly. For the singleproduction case, whose results are shown in the right panel of Fig. 5, we study the single lepton channel 6) and the corresponding SM irreducible background. Other than the same basic selection cuts imposed in the pair-production case, we further impose a requirement on the missing transverse energy / E miss T > √ s/3. This is motivated by the fact that in the signal case the light active neutrino carries away ∼ 50% of the available center of mass energy, while this is not the case for the background processes, for which the / E T distribution peaks at lower values.

N displaced decay
We now study the sensitivity for RH neutrinos decaying with a displacement which, as discussed in Sec. 4.6, we take to be between 1 cm and 100 cm from the primary vertex. The final event yield for having reconstructed displaced events is parametrized as where σ prod is the pair-production or single-production cross section for N and L denotes the total integrated luminosity. P ∆L represents the acceptance for having a RH neutrino decaying within a certain displacement from the primary vertex. This can be computed from the exponential decay law, taking into account the Lorentz time dilation factor. We then assign a probability for having the RH neutrino decaying at a distance ∆x = x f − x i which reads where the βγcτ factors are reported in Eq. (4.20) for pair-production and single-production cases, for which the parameter n in Eq. (7.1) takes the value of 2 and 1 respectively. This means that for the pair-production case we ask to reconstruct both RH neutrinos as decaying displaced. With disp we instead parametrize the acceptance for reconstructing the displaced decaying neutrino, which depend on the actual detector design and performances, and which therefore we assume as a free extra parameter in the analysis. The irreducible SM background is expected to be negligible on the considered decay lengths and we thus work in the zero background hypothesis. We then show the expected 95% CL exclusion limits, now obtained by requiring N s > 3, in Fig. 6 and Fig. 7 for the pair-production and single-production cased under the assumption of exclusive νγ and 3f decay respectively. The solid and dashed lines correspond to the choice disp = 1 and 0.3 respectively, while the different colors represents the different collider options. From the results we observe that a displaced analysis at the FCC-ee running at √ s = 240 GeV can be sensitive to O(10 TeV) NP scale with a 30% efficiency on the reconstruction of the displaced for m N 10 GeV in the pair-production case, while a higher reach can be attained in the single-production scenario. The FCC-ee running at the Z pole mass can slightly increase these reach due to the large integrated luminosity, while the 3 TeV collider prototypes can reach up to Λ ∼ 50 − 60 TeV for m N ∼ 40 GeV.

Detector stable N
Finally, we discuss the possibility of detector stable RH neutrinos, e.g. the case in which the decay happens more than 500 cm away from the interaction vertex. In this case, both pair-production and single-production give rise to a totally invisible final state. This process can be targeted through the emission of an initial state photon, producing a monoγ signature, + − → γ / E T , which has as SM background + − → νν / E T . In [91] exclusion prospects for various four-Fermi operators producing a weakly interacting massive particle dark matter candidate were given using a full detector simulation of the International Linear Detector prototype for the International Linear Collider. Moreover, rescaling factors for different collider energies, luminosities and beam polarizations where provided. Based on these results at the FCC-ee with 5 ab −1 of integrated luminosities, cutoff scales up to Λ ∼ 1.5 TeV can be tested in the pair-production case. In the single-production case the cross-section is larger than in the pair-production case but the photon spectrum is expected to be more similar to the SM due to the presence of only one heavy particle in the final state. Overall we thus expect the exclusion reach on Λ to be similar to the one of the pair-production case. However for such low scale the RH neutrino N → γν decay happens inside the detector, see Fig. 2, unless there is a cancellation among the α LN B and α LN W Wilson coefficient, see Eq. (4.2). If the dominant decay is N → 3f instead, the RH neutrino can be stable on detector lengths if Λ > 750 GeV and m N < 2 GeV, so that the derived limit of 1.5 TeV applies.
For higher center of mass energies we can again use as a guidance the results of [91]. Here the derived reach of CLIC at √ s = 3 TeV with 1 ab −1 of integrated luminosity is Λ ∼ 10 TeV. For CLIC and the 3 TeV µµ collider at the same center of mass energy we expect a reach in the same ballpark, although a dedicated study is required for a quantitative assessment. By again a comparison with Fig. 2 we see that a reach of 10 TeV on Λ will be able to prove detector stable RH neutrinos up to 5 GeV and 10 GeV if the only available decay mode is the one into νγ and 3f respectively.

Conclusions
In this paper we have considered the νSMEFT and studied how the RH neutrinos N production and decays may be affected by the inclusion of d = 6 operators. More specifically, we have studied the reach of future Higgs factories machines on the cutoff scale Λ at which the EFT is generated. We focused on four representative machines: the FCC-ee at two different center-of-mass energies, √ s = 90 GeV and √ s = 240 GeV, CLIC at a center of mass energy of 3 TeV and a representative muon collider with √ s = 3 TeV. The complete list of non-redundant d = 6 operators is presented in Tab. 1. At the level of production, the d = 6 operators induce either N pair-or single-production. On the other hand, at the level of decays, they induce the modes N → νγ and N → 3f , where various fermions combinations are possible. The former will dominate for RH neutrino masses m N 15 GeV, while the latter will dominate for larger masses, unless the only operators switched on induce a mixing-suppressed decay. Even more interestingly, depending on the RH neutrino mass and on the cutoff scale Λ at which the EFT is generated, the decays can be prompt, displaced or the RN neutrinos can be collider stable. The phenomenology crucially depends on their decay behavior and we have analyzed in detail all three possibilities. Our analysis is reported in Sec. 6, Sec. 7 and Sec. 8 for the three possible RH neutrinos lifetime. We then summarize the results for convenience in Fig. 8, in which, for the Higgs factories considered in this work, we show the 95% C.L. exclusion on the scale Λ as a function of m N . We consider RH neutrino masses up to 80 GeV. For larger masses, the W boson can be produced on-shell in the N decays and our analysis should be slightly modified. We postpone the analysis of such case to future work, although we do not expect major changes with respect to the results shown here. In the left panel we consider the decay channel N → νγ, while in the right panel we show the results for N → 3f . In both panels, the gray region denotes the parameter space in which the RH neutrino decay is displaced. The solid lines show the exclusion (combining pair and single-production) computed with prompt decays, an analysis valid in the white region. The dashed lines, on the contrary, show the exclusion limit considering displaced decays with an efficiency of reconstruction of 30%. In the region of validity of the prompt analysis, the FCC-ee will be able to probe scales up to Λ ∼ 7 TeV, while larger values, up to Λ ∼ 20 − 30 TeV, can be probed with a displaced analysis. These conclusions are valid for both decay channels. In the case of the colliders at 3 TeV, on the other hand, scales up to Λ ∼ 20 ÷ 30 TeV can be probed while the displaced analysis, on the other hand, allows to probe scales up to Λ ∼ 60 TeV.

A Spin averaged matrix elements for N decay
We list here the spin-averaged matrix elements |M| 2 = 1 2 spins |M| 2 for the three body decays of the RH neutrino via the d = 6 operators that proceed through an off-shell boson considered in the text. The kinematics is fixed as 1 → 2, 3, 4 and we define m 2 ij = (p i +p j ) 2 . The final state SM neutrino is always considered to be massless while, depending on the simplicity of the expressions, some of the amplitudes are reported in the limit of vanishing masses for the other final state fermions. From these amplitudes squared the partial widths are readily obtained as [92] dΓ = 1 (2π) 3