The see-saw portal at future Higgs Factories

We consider an extension of the Standard Model with two right-handed singlet fermions with mass at the electroweak scale that induce neutrino masses, plus a generic new physics sector at a higher scale $\Lambda$. We focus on the effective operators of lowest dimension $d=5$, which induce new production and decay modes for the singlet fermions. We assess the sensitivity of future Higgs Factories, such as FCC-ee, CLIC-380, ILC and CEPC, to the coefficients of these operators for various center of mass energies. We show that future lepton colliders can test the cut-off of the theory up to $\Lambda \simeq 500 - 1000\;$TeV, surpassing the reach of future indirect measurements of the Higgs and $Z$ boson widths. We also comment on the possibility of determining the underlying model flavor structure should a New Physics signal be observed, and on the impact of higher dimensional $d=6$ operators on the experimental signatures.


Introduction
The discovery of neutrino masses and oscillations is one of the most striking evidences for the need of new physics (NP) beyond the Standard Model (SM). Arguably, the simplest extension of the SM consists in extending its field content with the right-handed (RH) counterparts of the left-handed SM neutrinos, N . At the renormalizable level this allows for Yukawa type interactions between the new states and the SM leptons, providing a Dirac mass term for the neutrinos. However, being electroweak (EW) and color singlets, the new states can also have a Majorana mass terms, M N . As it is well known, this allows to explain the lightness of the observed neutrino masses through a large hierarchy between the EW scale and the mass scale of the RH neutrinos where y is the strength of the Yukawa interaction and v the Higgs vacuum expectation value (VEV). The relation of Eq. (1.1) defines the see-saw mechanism [1][2][3][4]. For a natural choice of the Yukawa interactions, y = O(1), the lightness of neutrino masses requires RH neutrinos at around the Grand Unification scale. However see-saw models with EW-scale RH neutrinos have recently received increasing attention. On the one side they offer a compelling alternative for the generation of the matter-antimatter asymmetry via neutrino oscillations [5,6], while on the other side they can be searched for at colliders and in beam-dump experiments [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].
At the renormalizable level these extra states can only be produced or decay via their mixing with the active neutrinos. This mixing, that controls their charged-and neutralcurrent interactions, is given by For EW scale RH neutrinos the naive see-saw scaling of Eq. (1.1) requires a tiny value for the Yukawa coupling connecting the SM and the beyond the SM (BSM) sectors. This implies a tiny mixing of the RH neutrino, resulting in a small production cross-section and a small decay width. The latter can give rise to striking signatures, such as displaced decays.
Interestingly, both the properties of production via mixing and displaced decays of the RH neutrinos can be challenged. The naive see-saw scaling relation can be broken when more than one RH state is present by specific Yukawa and Majorana mass textures that ensure an approximate lepton number symmetry [24,25]. The mixing can be much larger than the one implied by the see-saw relation, thus modifying the lifetime of the BSM states. RH neutrinos can then feature a prompt, displaced or detector stable behaviour.
These predictions can also be altered by the presence of additional NP states at a scale Λ v, M . At low energy their effects can be described in the language of effective field theories (EFT) by a tower of higher dimensional operators O d Λ 4−d with d > 4, built out from the SM and RH neutrinos fields: the νSMEFT. At the lowest dimension, d = 5, two new operators intervene to induce new RH neutrinos production modes: an operator triggering a new Higgs decay channel into a pair of RH neutrinos and a dipole operator connecting the RH neutrino tensor current with the hypercharge gauge boson [26,27]. At d = 6 many more operators are present [8,28,29], which can induce new production as well as new decay channels.
Various theoretical studies have investigated the signatures at the Large Hadron Collider (LHC) of a subset of these higher dimensional operators, see e.g. [8,20,26,27,[30][31][32]. The search for EW scale RH neutrinos is however one the primary goal of future e + e − colliders, thanks to the clean detector environment and the tipically lower SM backgrounds with respect to an hadronic machine, which can help to overcome the generally small production cross-sections of SM singlet states. Various future prototypes has been designed for the post LHC era: both circular colliders, as the Future Circular Collider [33][34][35][36] (FCC-ee) and the Compact electron-positron collider [37,38] (CEPC), and linear ones, such as the International Linear Collider [39][40][41] (ILC) and the Compact Linear Collider [42,43] (CLIC). It is then the purpose of this paper to investigate the phenomenology of the νSMEFT at these future machines and study their sensitivity on the d > 4 operators inducing new RH neutrinos production and decay modes in all the possible regimes of the N decay lifetimes, for RH masses in the range 1 -60 GeV.
The paper is organized as follows. In Sec. 2 we introduce our notation and describe the active-sterile mixing formalism while in Sec. 3 we discuss the properties of the various e + e − colliders under analysis. In Sec. 4 and in Sec. 5 we describe the analysis details and show the projected reach of the collider prototypes on the d = 5 operators involving RH neutrino fields, while in Sec. 6 we discuss the possible impact of d = 6 operators inducing extra N production and decay modes. We then conclude in Sec. 7.

Theoretical framework
We work in the framework of the νSMEFT, which is described by the following Lagrangian where N is vector describing N flavors of gauge singlet RH neutrino fields with N c = CN T and C = iγ 2 γ 0 , 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 neutrino fields and O n the Lorentz and gauge invariant operators built out from the SM and the RH neutrino fields. The νSMEFT has been constructed up to d = 7 in [8,[26][27][28][29]. At dimension five only three operators exist

Neutrino mixing formalism
Without loss of generality it is possible to go from Eq. (2.1) to a basis where the matrix M N and the charged lepton mass matrix are diagonal with non negative entries. Working at d = 5, the operator O N H contributes to the neutrino mass matrix. By defining n = (ν L , N c ) and using H = 174 GeV, the mass Lagrangian in the neutrino sector can be written as where the ν L − ν L block receives a contribution only from the d = 5 Weinberg operator while the N − N one has both d = 4 and d = 5 contributions. This mass matrix can be perturbatively diagonalized in the regime in which the entries of the ν L − N block are smaller than the ones in the N − N one. For our purposes we assume that the see-saw contribution to the active neutrino masses dominates over the other ones. Under this approximation we obtain ν is diagonal with non negative entries and U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [45,46]. From Eq. (2.4) one can obtain ν and √ µ and √ M N indicate, respectively, µ 1/2 and M 1/2 N . The usefulness of this parametrization is that it allows to write in a compact way the expressions for the various matrices involved. We now restrict our analysis to the case of two RH neutrinos, thus fixing N = 2. Without loss of generality the matrix √ µ can be written using the so-called Casas-Ibarra parametrization [47] as where √ m is a 3 × 2 matrix containing the physical neutrino masses m i , while R is a complex orthogonal 2 × 2 matrix, R T R = 1. With two RH neutrinos one has m ν 1 = 0 and m ν 3 > m ν 2 in the normal hierarchy (NH) case, while m ν 3 = 0 and m ν 2 > m ν 1 in the inverted hierarchy (IH) one 1 . More in detail, for NH and inverted IH we have while we parametrize the orthogonal matrix R in terms of the complex angle z = α + iγ as For both hierarchies we can thus write where m = m NH or m IH , and obtain a compact expression for the active-sterile mixing angle It's crucial that the angle z can be taken in general as a complex number. 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 This "naive see-saw scaling" relation is drastically modified by the imaginary part of z, that gives an exponential enhancement. In the limit γ 1 12) and the relation of Eq. (2.11) is modified to The same enhancement is inherited by the active-sterile mixing, that now reads (2.14) In the previous expression we have α = e, µ, τ and i = 1, 2. As anticipated in the Introduction and as we will see more in detail below, 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 huge implications for search strategies at future colliders. 1 For the NH we take mν 2 = 8.6 × 10 −3 eV and mν 3 = 5.1 × 10 −2 eV while for the IH we take mν 1 = 4.9 × 10 −2 eV and mν 2 = 5.0 × 10 −2 eV. 2 We have assumed NH and fixed mν = mν 3 . The expression holds also for the IH case modulo order one factors.

Final state
Channel Mediator Table 1. Possible decay channels for the RH neutrino N . Here α, β, i and j are flavor indices and we do not specify the charge of the charged lepton = e, µ, τ nor the nature of the (anti)neutrino.

Heavy neutrinos decay modes
In the mass range of our interest and at the renormalizable level the RH neutrinos can only decay through charged-and neutral-currents via an off-shell W or Z boson 3 . The RH neutrino decay mode is thus completely fixed once the W and Z decay channels are specified. The various final states from N decay are reported in Tab. 1 where α, β, i and j are flavor indices and, for simplicity, we do not specify neither the charge of the charged lepton = e, µ, τ nor the nature of the (anti)neutrino. In this table we have grouped together the three decay modes giving rise to the ν final state since, given that both the W and Z boson will be non resonant for the RH neutrinos mass range of our interest, these processes will not be distinguishable. Moreover, the ν α β β process with α = β receives contributions from both neutral-and charged-current interactions, which interfere among themselves.
For computing the partial widths into the final state of Tab. 1 we use the results of [48]. For m N Λ QCD , the decay rates involving quark pairs are physical quantities, otherwise decays into hadrons should instead be considered. Following again [48] we have implemented three-loop QCD corrections through which the full hadronic width can be computed from the decay width into free quarks. Altogether the effect on the total width is found to be around 30% for m N 1 GeV, decreasing down to 10% for m N 5 GeV. By fixing the phases of the PMNS matrix δ = φ 1 = 0 4 and also γ = α = 0 we obtain the branching ratios (BRs) shown in the left panel of Fig. 1 for the case of the lightest RH neutrino N 1 , while similar rates are obtained for N 2 .
The lepton flavor composition of the final states depends on the active-sterile mixing matrix, which in turn depends on i) the hierarchy and the squared mass differences of the active neutrinos, ii) the phases of the PMNS matrix δ and φ 1 , and iii) the α and γ parameters entering in the Casas-Ibarra parametrization of Eq. (2.8). Analytical approximations for the various mixings can be derived when γ 1, see e.g. [49]. In this regime 3 A decay into an off-shell Higgs boson is generally suppressed by the smallness of the SM Yukawa couplings. 4 We remind the reader that with N = 2 RH states only two phases are present in the PMNS matrix.
We denote by δ the so-called Dirac phase, and by φ1 the so-called Majorana phase. the Casas-Ibarra parameter α has a little impact on the normalized squared mixing which is mainly determined by the PMSN phases δ and φ 1 . By varying them between [0,2π] we obtain for r 2 α1 the allowed ranges shown in the right panel of Fig. 1. The red and blue regions correspond to NH and IH, respectively, and r 2 τ 1 = 1 − r 2 e1 − r 2 µ1 , see also [20]. We also show the benchmark points that will be used in the following analysis.

Final states from pair produced heavy neutrinos
Both the operators O N H and O N B will mediate the production of a pair of RH neutrinos, which in turn will decay producing a six-body final state. These final states can be categorized into fully-leptonic, fully-hadronic, semi-leptonic and invisible channels and are reported in Tab. 2. Here we differentiate between = e, µ and τ , since the latter particle decays within the detector thus producing a different and more complex final state. With an abuse of notation we nevertheless use the nomenclature fully-leptonic and semi-leptonic also for the final states involving τ leptons. Particularly interesting are the final states that can present a pair of same-sign (SS) leptons = e, µ, a signature which generally has a low SM background. Final states with ≥ 3 will clearly have a SS pair. However, also the final state with exactly 2 can produce a SS signature, due to the Majorana nature of the RH neutrinos. The branching ratios of N 1 and N 2 into the various final states of Tab. 2 will depend upon the choice of the normalized squared mixings r 2 αi which, as explained above, depend on the PMNS parameter δ and φ 1 and can span the ranges illustrated on the right panel of Fig. 1.
For what concerns the detection and reconstruction of the final state, e and µ can be considered in first approximation on the same footage, while one needs to distinguish them with respect to τ leptons. We then choose to perform our analysis and illustrate our results for some representative points in the allowed range for the r 2 αi shown in Fig. 1 for both the NH and IH cases. In particular for each mass hierarchy we choose two points, one with a large and one with a small mixing with the third generation leptons. After having fixed r 2 τ i we choose to maximize the mixing with the electron, which does not affect the reconstruction of the final state, in the approximation of similar e and µ detection efficiencies, but does intervene in the production of RH neutrinos via mixing, e + e − → νN . The reason why we choose to maximize r 2 eα is because, conservatively, we want to analyze a configuration with the largest possible production cross-section via mixing, to see whether the additional production modes arising from the O N H and O N B operators can still dominate over it. More concretely, we choose the following two benchmark points for the NH case  Table 3. Decay rates from N N production for BP1 NH (left) and BP2 NH (right). The rates are obtained by summing on all charges and flavor configurations. Here = e, µ. The checkmarks correspond to channels that can produce a SS lepton signal. Table 4. Decay rates from N N production for BP1 IH (left) and BP2 IH (right). The rates are obtained by summing on all charges and flavor configurations. Here = e, µ. The checkmarks correspond to channels that can produce a SS lepton signal.

Heavy neutrinos lifetime
A crucial quantity affecting the phenomenology and thus the search strategies for RH neutrinos is their lifetime τ N = 1/Γ N . We focus for simplicity on the case of (almost) degenerate RH neutrinos and we start by considering only the decay modes of Sec. 2.2 that are induced at the renormalizable level by the active-sterile mixing 5 . In particular, as we will show in Sec. 4.1, one can have RH neutrinos that decay promptly, displaced or are stable on detector length scales.

Prompt decay
We consider a RH neutrino decay as prompt if it happens within ∼ 0.1 cm from the primary vertex. The production of a pair of N gives rise to a six-body final state, including signatures with high lepton multiplicity, see Tab. 3 and Tab. 4. As we will see, in order to have promptly decaying RH neutrinos, one needs to have a large breaking of the naive see-saw scaling between the active-sterile mixing, the RH neutrino and the light neutrino masses. In the notation of Sec. 2 this breaking is parametrized by a large value of the γ parameter, see Eq. (2.14). 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 − 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. Given 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 in order 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.

Collider prototypes
Lepton colliders are ideal machines for SM precision measurements due to the cleanliness of their environment, the precise knowledge of the initial-state particles configuration and, for some prototypes, the possibility of having polarized beams that can help to enhance the signal-to-background ratio. Despite their center-of-mass energy being typically smaller than the one of hadronic machines, they however offer excellent prospects in the direct search for NP. This is mainly due to the low SM backgrounds, whose rates are generically comparable to the searched signals, as opposed to what happens in hadronic machines. Amongst the various proposals for a new generation of colliders after the HL-LHC era, leptons colliders then stand out as one of the more concrete possibility. 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) [33][34][35][36] and the Circular Electron Positron Collider (CEPC) [37,38], and linear ones, as the International Linear Collider (ILC) [39][40][41] and the Compact Linear Collider (CLIC) [42,43]. We report in Tab. 5 the center-of-mass energies and luminosities considered in this study  Table 5. Center-of-mass energies and total integrated luminosities for the various collider options considered in the analysis for Higgs runs (left) and Z pole runs (right). For the Higgs runs we report the values of the Higgs-strahlung cross-section reported in [51], while for the Z pole runs the number of expected Z bosons produced with the corresponding integrated luminosity from [52].
for various benchmark configurations. Notice that for CLIC we focus on its low-energy stage at √ s = 380 GeV, to which we refer as CLIC-380 throughout the text. We report both the parameters for Higgs physics runs as well as the ones for runs at √ s = m Z . For what concerns the Higgs-strahlung cross sections σ(e + e − → Zh), in the cases of the ILC and CLIC-380, these are reported under the assumptions of a beam polarization fraction (P e − , P e + ) of (−80%, +30%) and (−80%, 0%), respectively.

Simulation Framework
In our analysis, signal events have been simulated at parton level by MadGraph5 aMC@NLO [53]. The events have then been analysed with the MadAnalysis5 package [54][55][56]. We consider our sensitivity estimates to be preliminary, not including any irreducible or reducible backgrounds. Nevertheless, we expect them to be not too far from a realistic lepton-collider sensitivity projection. For instance, an Higgs invariant-mass selection cut on the final states featuring very high multiplicity Higgs decays, as we are considering here, should be sufficient for suppressing the irreducible SM backgrounds in the relatively clean lepton-collisions environment. Although a full simulations of the actual detector performance, when available, will certainly make the corresponding projections more robust, we are confident that, in a more realistic approach, the excellent accuracy of particle-flow reconstruction, as now under consideration for Higgs Factories detectors, complemented by advanced analysis techniques, might only moderately degrade the present sensitivity estimates.

The O N H operator and the Higgs-strahlung channel
We start our analysis by discussing the O N H operator. It can be generated by the tree-level exchange of a scalar singlet or by a fermion doublet with hypercharge ±1/2 [27]. It gives rise to a new interaction of the Higgs boson with a pair of RH neutrinos. If kinematically allowed, this interaction induces an extra decay channel for the Higgs boson, h → N N . In the following we assume for simplicity degenerate RH neutrino masses For real couplings 6 the partial width reads [26] This operator can be constrained by searches for additional untagged Higgs decay modes or invisible Higgs decays [51]. More importantly, it adds a new production mode for RH neutrinos through the Higgs-strahlung process At lepton colliders the process of Eq. (4.3) is the dominant production mode for a SM Higgs boson for center-of-mass energies √ s 400 GeV. It is crucial that, by using the recoil mass technique, this process can be tagged by reconstructing the Z decay products, without any knowledge of the particles arising from the Higgs boson decay. This property makes this channel important for all the three regimes of the RH neutrinos lifetimes described in Sec. 2.4. The Higgs-strahlung cross-sections for the various colliders are normalized as reported in Tab. 5.

Prompt decay
As anticipated in Sec. 2.1, in order to have the RH neutrinos to decay promptly, one needs a large breaking of the naive see-saw scaling, parametrized by a large value of the parameter γ entering the Casas-Ibarra parametrization of Eq. (2.8) and enhancing the active-sterile mixing angle. This mixing is however constrained by a variety of experimental searches and too large values of γ are ruled out. Using the bounds on θ α = i=1,2 |θ αi | 2 reported in [57][58][59], we show in Fig. 2 as a gray dashed line the exclusion contour in the m N − |θ 2 e | plane for the NH case. Similar results are obtained for IH. We show only the bound arising from θ µ , which turns out to be the most stringent one. Notice that |θ e | 2 inherits the bound from θ µ through its dependence on γ and m N . In the same plot we also show as black dashed lines the isocontours of proper decay length cτ of the RH neutrino N 1 , in order to identify the regions where the decay is prompt, displaced or outside of the detector 7 . We neglect the dependence of the lifetime on α, δ and φ 1 which is mild for γ 1. We also show the colored regions in which the RH neutrino pair production from the Higgs-strahlung process of Eq. (4.3) is larger than the production of a single RH neutrinos via mixing. For concreteness we consider the case of the FCC-ee collider with √ s = 240 GeV, and compare 6 For imaginary couplings there is a different dependence on the RH neutrino velocity βN due to CP properties of the matrix element see, e.g., [26]. 7 Strictly speaking, the important quantity is the laboratory frame decay length, which is larger than the proper decay length due to Lorentz time dilation. We will accurately compute this quantity in Sec. 4.2 when dealing with displaced vertices while in this section we assume to be in a region of the mN − |θ 2 e | parameter space where the Lorentz factor does not modify the behavior between prompt, displaced or stable. This is the case for cτ values away from the boundary regions indicated in the plot for not too light N . Figure 2. Regions in the m N − |θ e | 2 parameter space where the RH neutrino pair production from the Higgs-strahlung process of Eq. (4.3) is larger than the production of a single RH neutrinos via mixing, for BR(h → N N ) = 10% (red), 1% (blue) and 0.1% (green) for the FCC-ee case. The NH case is assumed. The gray dashed line represents the limit on the mixing angle arising from existing experimental searches, while the black dashed lines represent isocontour of proper decay length cτ . In the gray shaded region the lightness of the neutrino masses cannot be explained by the see-saw mechanism.

���-��� �����
the Higgs-strahlung cross section with the one of production via mixing. Since the result depends on the branching ratio of the Higgs boson into RH states, we show our results for BR(h → N N ) = 10% (red), 1% (blue) and 0.1% (green), using [17] for the normalization of the mixing cross-section. Finally, in the gray area at the bottom of the plot the lightness of the neutrino masses cannot be explained by the see-saw mechanism.
Altogether, the figure makes clear that there are large regions in parameter space in which Higgs-strahlung production may dominate over mixing production. Moreover, depending on m N , in these regions the RH neutrinos can have prompt, displaced or outside the detector decays depending on the active-sterile mixing.

Projected sensitivities: e and µ mixing
We start by considering the benchmark points which minimize the mixing with the third generation leptons, namely BP1 NH and BP1 IH of Eq. (2.16) and Eq. (2.17). For these two cases we focus for simplicity on the final state with the highest rate, which according to Tab. 3 and Tab. 4 is the 2 4q one. Furthermore, we require this final state to contain a pair of SS leptons, thus halving the rates reported in the tables, and focus on the case in which the Z boson decays leptonically. In computing our limits we sum on both e and µ flavor combinations for the Z and the N pair decay modes, and on both the RH neutrinos Figure 3. Feynman diagram for 2 4q production in the SS leptons final state through Higgsstrahlung production.
N 1 and N 2 . All together the process we analyze is where the first bracket indicates the Z boson decay products, while the second bracket the Higgs boson ones and α, β, γ = 1, 2 are flavor indices. This is shown in Fig. 3. For our analysis we require the leptons to have p T > 2.5 GeV and |η| < 2.44, while jets should satisfy p T > 5 GeV and |η| < 2.4. Leptons are required to be separated by ∆R > 0.15 among themselves and with respect to the selected jets. In order to tag the Higgs-strahlung topology, we require two same-flavor opposite-sign leptons with an invariant mass |m + − − m Z | < 10 GeV. If more than one pair that satisfies this condition is present, we choose the pair with an invariant mass closer to the Z mass. Furthermore this pair is also required to have a recoil mass m rec within 10 GeV of the true Higgs mass m H , where and E + − is the energy of the leptons pair. We then ask for exactly two other SS leptons in the event, while no requirements on the number of jets and missing energy is imposed. Three final state flavor configurations from the N N decay contribute to the selected final state: e + e + , e + µ + and µ + µ + plus the charge conjugated processes. The parton level acceptances as a function of the RH neutrino mass 8 for these final states are shown in Fig. 4. As it can be seen the acceptances increase with the increase of the RH neutrino mass. This is to be expected, since for light RH neutrinos their decay products turn out to be more collimated, thus making it harder to pass the ∆R > 0.15 isolation criteria.
Regarding the irreducible SM backgrounds to the process of Eq. (4.4), we expect it to be negligible. In principle, since we are not requiring any jet in the final state, any SM process matching a signature (Z → + Lepton-number conservation in the SM however implies that at least two extra neutrinos should be present in the final state. Electric-charge conservation moreover implies that there should also be at least two further quark pairs, arising e.g. from two off-shell W 's. A typical SM irreducible background to the process under consideration is therefore where the final state arising from the Higgs decay involves multiple off-shell gauge bosons, such a h → 4W * decay. We then expect it to be totally negligible. Due to the very large multiplicity of the chosen signal, we also expect other channels that match the same final state without occurring via Z and Higgs resonances to be negligible. As for other reducible backgrounds arising from the limited efficiencies in the detector performances, such as e.g. mis-identification of the final particles, we are confident that they can be efficiently removed thanks to the strong kinematical characterization of the (Z → + α − α )(h → + β + γ + . . . ) process. We thus assume the process in Eq. (4.4) to be background free. According to Poisson statistics, in case no signal event is observed, one can then set a 95% confidencelevel (CL) exclusion limit corresponding to a maximally allowed number N s = 3 of new physics events [60].
We show our results for BP1 NH in the left panel of Fig. 5. Similar results are obtained for the BP1 IH benchmark. The colored regions represent the 95% CL sensitivity of the analysis on the exotic Higgs BR described above in the m N − BR(h → N N ) plane for the various collider options. This reach has to be compared with the limits coming from untagged Higgs decay [51], represented by horizontal dashed colored lines. For simplicity we do not show the limit that might be obtained at CEPC, which are comparable to the ones from FCC-ee, due to their similar integrated luminosity, see Tab. 5. We also show in dashed black the isocontours of the scale Λ expressed in TeV, fixing α ii N H = 1. Finally, we highlight in gray the region where cτ > 1 mm in which the RH neutrino cannot decay promptly without being excluded by experimental searches, see Fig. 2.
All together we see that future e + e − colliders will be able to set a bound on the exotic Higgs BR that goes from ∼ 5×10 −3 for the case of CLIC-380 down to ∼ 7×10 −4 for the case of FCC-ee/CEPC, and that these bounds are significantly stronger than the corresponding ones arising from untagged Higgs decay measurements. In terms of NP scale Λ these limits translate in a bound which, in the most favorable case, is Λ 500 TeV.

Projected sensitivities: τ mixing
We now consider the case where the active-sterile mixing with the third generation of leptons is maximized, as is the case for the benchmarks BP2 NH and BP2 IH in Eq. (2.16) and Eq. (2.17). Again we choose for simplicity the channel with the highest BR for both benchmark points, i.e. the 2τ 4q channel. In this case however the τ leptons will promptly decay into a ν τ and an off-shell W boson. For our analysis we adopt the following strategy. We first consider the τ leptons as stable particles, and apply the same parton-level selection efficiencies computed in the e and µ case of Sec. 4.1.1. Then we focus on the hadronic decay modes τ − → π − π 0 ν τ , τ − → π − ν τ and τ − → π − π 0 π 0 ν τ (which account, respectively, for approximately 25.5%, 10.8% and 9.3% of the total τ branching ratio) and apply, following [61], a flat 90% reconstruction efficiency for each τ lepton. We consider the inclusive 2τ 4q final state with all charge combinations, τ ± τ ± and τ ± τ ∓ , and, again, neglect the SM background, although in this case we expect a larger contamination with respect to the 2 4q SS case 9 . Under these assumptions we obtain the results shown in the right panel of Fig. 5 for BP2 NH . Similar, albeit slightly weaker, results are obtained for the BP2 IH case. We see that also in the case in which N mixes dominantly with ν τ , both FCC-ee and ILC can probe values of BR(h → N N ) that go well beyond the limits arising from untagged Higgs decay, while CLIC-380 can marginally surpass this reach.

Determination of the flavor structure
Should a RH neutrino signal be detected, a crucial question to be asked is with which accuracy it will be possible to determine the flavor structure of the underlying BSM theory. This is directly related to the number of signal events which in turn depends on the physics of the underlying theory, namely BR(h → N N ), and the detector performances in reconstructing the signal. To quantify this we adopt the following strategy. We choose the 2 4q final state with a pair of SS leptons, a representative RH neutrino mass of 30 GeV and the benchmark points BP1 NH and BP1 IH . We then build a Poisson distribution p(n obs |n th ) = 1 n obs ! e −n th n n obs th where n obs correspond to the expected number of observed events in the case of the presence of a BSM signal, and n th is the theoretical prediction for the number of events in the chosen final state, which is a function of r 2 e1 and r 2 µ1 . For fixed values of BR(h → N N ), which then fixes the final event yield, we thus compute the 68% and 95% CL confidence interval around the chosen benchmark points. The results 9 For instance we estimated the cross section for the irreducible background e + e − → (Z → µ + µ − )(h → τ + τ − udūd) to be about 5.4 × 10 −6 pb, corresponding to about 27 events at FCC-ee. We expect that these can be efficiently reduced by exploiting the kinematical features of the h → N N decay. Table 6. Range of parameters that can be probed in case of detection of a RH neutrino decay for the representative m N mass and BR(h → N N ) considered in Fig. 6. The range for φ 1 corresponds to the red lines in the figure, while the range for δ to the blue ones. In both case the other phase is kept fixed to the benchmark point value.
are illustrated in Fig. 6 for the case of BR(h → N N ) = 1% (left) and 0.3% (right) for the FCC-ee collider option. We see that with a 1% Higgs exotic BR, which correspond to ∼ 37 total signal events in the final state, the normalized squared mixings r e1 and r µ1 can be determined with an absolute error of ∼ 0.1, which rapidly degrades down to ∼ 0.3 with a 0.3% Higgs BR into a pair of RH neutrinos, for which one has ∼ 11 signal events. It is interesting to interpret these results in terms of the phases appearing in the PMNS matrix, δ and φ 1 . The benchmark points we choose correspond to fixed values for both phases: (δ, φ 1 ) = (0.76, 4.59) for the BP1 NH benchmark, and (δ, φ 1 ) = (3.25, 1.60) for the BP1 IH benchmark. By keeping δ fixed and allowing φ 1 to vary we obtain the blue lines in Fig. 6, while by keeping φ 1 fixed and allowing δ to vary we obtain the red lines. The interval that δ and φ 1 can span inside the 95% CL confidence intervals around the benchmark points values are reported in Tab. 6.

Displaced decay
We now study the sensitivity for RH neutrinos decaying with a displacement which, as discussed in Sec. 2.4, we take to be between 1 cm and 100 cm from the primary vertex, see Fig. 2. We consider decays into first and second generation leptons and we focus again on the 2 4q final state, which is the one that maximizes the decay rate of the N N pair. The final state therefore consists of two prompt same flavor opposite sign leptons from the Z boson decay and the 2 4q system. The final event yield is given by where Zh is the acceptance for reconstructing the Z boson and the Higgs recoil mass from the two same flavor opposite sign prompt leptons. The parameter 2 P ∆L represents instead the acceptance for having both neutrinos decaying within a certain displacement from the primary vertex. This probability can be computed from the exponential decay law, taking into account the time dilation factor obtained by boosting the events from the RH neutrino rest frame to the laboratory frame. In practice we have computed, for each event and for each RH neutrino, the Lorentz β and γ factors  Figure 7. 95% CL exclusion for BP1 NH in the m N − |θ e | 2 plane from displaced decay searches in the Higgs-strahlung channel assuming BR(h → N N ) = 1% (left) and 0.3% (right) for the case of FCC-ee. The shaded red area and the dashed red lines represent the sensitivities for different assumptions on the efficiencies for reconstructing the displaced vertices disp. . The gray dashed line represents the limit on the mixing angle arising from existing experimental searches. In the gray shaded region the lightness of the neutrino masses cannot be explained by the see-saw mechanism. and assigned a probability for having the RH neutrino decaying at a distance ∆x = (4.9) We have then accepted events where both RH neutrino decays happened between 1 cm and 100 cm from the primary vertex. Finally, with disp. we parametrize the acceptance for reconstructing the displaced decay, including various detector inefficiencies, which will 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 lenghts. We show in Fig. 7 the FCC-ee reach for the BP1 IH benchmark in the m N − |θ e | 2 plane. In the left panel we fix BR(h → N N ) = 1%, while in the right panel we fix BR(h → N N ) = 0.3%. We show our results for two different efficiencies disp. as reported in the plots. In both figures the gray shaded area represents the see-saw limit, while the values of the mixing angles above the gray dashed line are excluded by experimental searches, see the discussion in Sec. 4.1. Our results show that the search for RH neutrinos arising from Higgs decay and decaying displaced can offer a great handle in testing the active-sterile mixing angle. This is clearly due to the fact that, unlike in the case of production via mixing, the Higgs-strahlung cross section does not depend on the active-sterile mixing angle, which only enters in the lifetime determination, but only on the additional Higgs decay rate into an N N pair. As a consequence, searches at FCC-ee could test values of |θ e | 2 down to the see-saw limit. We do not show for simplicity the limits arising from ILC and CEPC, which are comparable to the ones shown, and the one from CLIC-380, which turns out to be slightly weaker.
In Fig. 8 we show instead the reach of the same search projected on the plane BR(h → N N ) versus m h 2m N cτ , which is roughly the laboratory-frame decay length of the RH neutrinos. We do this again for two different choices of disp. : 50% (left) and 100% (right). In the plots we also show the limits arising from untagged Higgs decays searches for the various collider options. All together we see that FCC-ee will be able to test values of the Higgs exotic BR down to 0.2% for disp. = 50%, largely surpassing the indirect limits from Higgs untagged decays, while ILC and CLIC-380 have a slightly weaker reach.

Detector stable
The last case we study is the one in which the RH neutrinos lifetime is large enough that they will decay outside the detector. In this case the decay will contribute to the invisible Higgs width. This quantity can be strongly constrained at future lepton colliders, which will set a 95% CL bound on BR(h → inv.) of 0.22% (FCC-ee), 0.28% (CEPC), 0.26% (ILC) and 0.63% (CLIC-380) [51]. These limits can be directly translated on a bound on the NP scale Λ through Eq. (4.1). We obtain Λ 360 TeV (FCC-ee), 320 TeV (CEPC), Λ 330 TeV (ILC) and 210 TeV (CLIC-380) fixing m N = 10 GeV. The limits degrade by roughly a factor 20% for m N = 35 GeV due to the reduced phase space. Above this mass threshold their decay will instead happen inside the detector, see Fig. 2. Figure 9. Isocontour of equal partial widths for N 2 → N 1 γ and N 2 decaying via mixing for various choices of BR(Z → N 1 N 2 ). We fix the relative mass splitting between N 2 and N 1 to 0.1% (left) and 1% (right). Above the lines, for that specific BR(Z → N 1 N 2 ), the decay via mixing dominates.
In the red shaded area Γ BSM Z is larger than the current experimental uncertainty on the Z boson total width. The gray dashed line represents the limit on the mixing angle arising from existing experimental searches. In the gray shaded region the lightness of the neutrino masses cannot be explained by the see-saw mechanism.

The O N B operator and the s-channel Z production
We turn now to the study of the O N B dipole operator. It can be generated only a looplevel by a scalar-fermion of vector-fermion pair with opposite hypercharges [27]. Due to the presence of σ µν , the flavor structure of this operator is antisymmetric, i.e. only different RH neutrinos can participate in the interaction. Since the operator induces a new interaction between the RH neutrinos and the SM neutral EW gauge bosons, it may provide an additional production channel for a pair of RH neutrinos through an intermediate photon or Z boson. In particular, if kinematically allowed, the Z boson can decay into a N 1 N 2 pair with a rate [27] Γ Z→N 1 N 2 = 2 3π where s w is the sine of the Weinberg angle, λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2bc − 2ac and In the following we will fix φ 12 = 0, i.e we assume real Wilson coefficients.
Future colliders with an operating stage at √ s = m Z , as is the case of FCC-ee and CEPC, will produce a large number of Z bosons, see Tab. 5, and can thus probe the operator responsible for the decay of Eq. (5.1) with high precision. We focus in particular on the case of FCC-ee, where one expects to have 6.5 × 10 12 Z candidates produced, after a total integrated luminosity of ∼ 150 ab −1 , while CEPC might have a luminosity smaller by roughly one order of magnitude, and weaker results are generally expected.
In addition to new production modes, the O N B operator can also trigger new decay channels for the heavier RH neutrino N 2 . In particular we have [27] while we do not consider the possibility of N 2 → N 1 Z decay which is outside of the mass range of interest here. Crucially, the decay mode of Eq. (5.3) can compete, and even dominate, with the one induced by the mixing between the active and sterile sector. It is thus important to assess in which region of parameter space the N 2 → N 1 γ decay rate can dominate over the decay induced by active-sterile mixing. In Fig. 9 we show such regions in the m N 1 − |θ e | 2 plane considering the NH case. As in previous plots, we show the limits coming from θ µ as gray dashed line. The black continuous lines correspond to Γ(N 2 → N 1 γ) = Γ(N 2 ) mix , where the latter is the total N 2 decay width due to mixing, i.e. in the channels listed in Tab. 1, for the specific value of BR(Z → N 1 N 2 ) reported. Above the lines the decays induced by mixing dominate over the decay induced by the O N B operator. Since this operator involves different neutrinos, different mass splittings may give different physical situations. We quantify this by defining the relative mass splitting r = (m N 2 − m N 1 )/m N 2 , which is fixed to 0.1% in the left plot and to 1% in the right plot. We also show in red the region in which the total BSM width Γ BSM Z is larger than the than the current experimental uncertainty on the Z boson total width, δΓ exp Z = 2.3 MeV, and the region in which the see-saw mechanism is not able to reproduce the observed neutrino masses. As we see, in both cases for reasonable values of the branching ratio there are large regions in parameter space in which the decay via mixing dominates over N 2 → N 1 γ.
As a consequence, in this Section we will focus on the signatures of N 2 decaying via mixing. Notice that, barring the contribution of d = 6 operators on which we will comment in Sec. 6, this is always true also for N 1 . As it happened in Sec. 4, we find that for mixing angles compatible with experimental bounds, the N 1 N 2 production via Z decay can easily dominate with respect to the N ν production via mixing. In the following Sections we will thus analyze the sensitivity of future experiments for the case in which the N 1 N 2 production occurs via Z decay and N 1,2 decay via mixing. As already discussed, also in this case the decay can be prompt, displaced or outside the detector. For simplicity we fix we study the case of almost degenerate RH neutrinos.

Prompt decay
We again focus on the final state with the highest rate, i.e. the 2 4q channel, restricting our analysis to the final state with a pair of SS leptons. This process is shown in Fig. 10. We require exactly two SS leptons with p T > 2.5 GeV and |η| < 2.44, while jets should satisfy p T > 5 GeV and |η| < 2.4. Leptons are required to be separated by ∆R > 0.15 among themselves and with respect to the selected jets. With these selections, the parton level acceptances for the e + e + flavor combination are reported in Fig. 11. Regarding the background analysis, in the case of the O N H operator we assumed that the physical background mimicking the Higgs-stralung production and subsequent six-body decay in a sample of around 10 6 Higgs bosons would be well under control and negligible in first approximation. On the other hand, in case of the FCC-ee production of O(10 12 ) Z bosons at the Z peak, quite a number of reducible backgrounds are expected to limit the statistical reach of the Z boson sample, depending both on the limited detector performances and on the collider characteristics. We then expect that a realistic background analysis for the Z → + + 4q channel might degrade the ideal background free estimate in a non negligible way.
As an example of how different reducible backgrounds can degrade the Z sample sensitivity, we will consider below the possible background arising from the limited leptoncharge identification power of a LHC-like detector. Indeed, conservation of lepton number implies that no genuine irreducible Z → + + 4q background arises in the SM. As regarding the Z → + + 2ν4q background, for which we expect the cross-section to be relatively small, this might also be efficiently reduced by asking for limited missing energy in the events. On the other hand, a realistic analysis of the reducible SM background Z → + − 4q, where a lepton charge is misidentified, or an hadron is misidentified as a lepton, would require a dedicated full simulation at the experimental level, which is beyond the scope of the current work. For the lepton-charge misidentification effect, assuming LHC performances of lepton charge identification, we provide an approximate estimate of the corresponding background by applying a (flat) mis-identification probability factor of misID = 10 −3 [62] to the partonic cross-section e + e − 4q and µ + µ − 4q after the signal selections described in the text 10 . With this procedure we obtain a background yield of 130 fb× misID (1 − misID ) 0.26 fb, which we use as estimate for our analysis.
We show our results for BP1 NH in Fig. 12 Figure 13. Left: 95% CL exclusion for BP1 NH in the m N − |θ e | 2 plane from displaced decay searches assuming BR(Z → N 1 N 2 ) = 10 −10 for the case of FCC-ee. The gray dashed line represents the limit on the mixing angle arising from existing experimental searches. In the gray shaded region the lightness of the neutrino masses cannot be explained by the see-saw mechanism. Right: 95% CL exclusion for BP1 NH in the m Z 2m N cτ − BR(h → N 1 N 2 ) plane from displaced decay searches. The horizontal red dashed line represents the limit arising from the FCC-ee Z boson total width measurement with an uncertainty of 100 keV.

��(�→��)=�� -��
following the procedure outline above. The red region is the one that will be probed by the FCC-ee, while the horizontal dashed red line represent the expected sensitivity on the Z boson width of 100 keV [33]. From the figure it is clear that, when considering the presence of the lepton charge mis-identification background, FCC-ee will be able to exclude values of BR(Z → N 1 N 2 ) down to roughly 10 −8 . This corresponds to a reach on Λ 10 3 TeV. Given that this operator can only be generated at loop-level, this bounds is rescaled by a factor 16π 2 when mapped into the physical masses and couplings of an ultraviolet complete model. We stress that the present estimate might be further degraded by other background sources, such as other particle identification inefficiencies which will be strongly dependent on the actual detector performance.

Displaced decay
For the case in which the RH neutrinos decay displaced from the primary vertex we follow again the strategy outlined in Sec. 4.2. As opposed to the case of prompt decay, in this case we expect negligible irreducible SM background. We show in the left panel of Fig. 13 the FCC-ee reach for the BP1 IH benchmark, projected in the m N − |θ e | 2 . We fix BR(Z → N 1 N 2 ) = 10 −10 and assume different efficiencies displ , as reported in the plots. The gray shaded are represents the see-saw limit, while the values of the mixing angles above the gray dashed line are excluded by experimental searches. As for the case of the O N H operator our results show that the search for RH neutrinos arising from Z decay and decaying displaced from the primary vertex can offer a great handle in testing the active-sterile mixing angle. Again this is due to the fact that the production cross-section does not depend on this quantity, but only on the Z decay rate into the NP final state. In the right panel of the same Figure we show instead the FCC-ee reach projected in the m Z 2m N cτ − BR(h → N 1 N 2 ) plane. We do this again for two different choices of displ as reported in the plot. All together we see that FCC-ee will be able to test values of the Z boson exotic branching ratio down to 10 −9 for disp. = 20%, largely surpassing the limit arising from the FCC-ee Z boson total width measurement with an uncertainty of 100 keV, which is represented by the horizontal red dashed line.

Detector stable
The last case we study is again the possibility that the RH neutrinos lifetime is large enough to cause them to decay outside the detector. Similarly to the O N H operator, in this case the decay will contribute to the invisible Z width. FCC-ee will measure the ratio R ν = Γ Z→inv /Γ Z→ at the level of 0.27 × 10 −3 [51] which, under the SM hypothesis, correspond to an additional contribution to invisible decay width of the Z boson smaller than ∼ 135 KeV. This limit can be translated in a bound on the NP scale Λ which will be constrained by this measurement to be Λ 16 TeV for m N ∼ 10 GeV while this limits degrade down to Λ 9 TeV for m N = 35 GeV due to phase space effect. Above this mass threshold their decay will happen inside the detector, see Fig. 2.

The impact of d = operators
We have so far considered the phenomenology induced by the presence of d = 5 operators in the νSMEFT, discussing how they can induce additional RH neutrinos production and decay modes and showing how they can be efficiently tested at future Higgs Factories. It is important however to notice that at d = 6 many more operators are present, and they might give observable signatures. For example, among the operators involving the Higgs field, of particular interest are (6.1) They will trigger the decay to an active neutrino and a photon with a rate where, with an abuse of notation, we also indicate with α LN B and α LN W the relevant entries of the corresponding Wilson coefficient matrices. Also this decay mode can in principle dominate over the one induced via the active-sterile mixing in some regions of the parameter space. A detailed analysis would proceed in a very similar fashion to the one described in Sec. 5 for the dipole operator O N B , with however an additional suppression due to their higher dimensionality. These two operators have also been recently studied in the context of LHC in ref. [30] where a bound on Λ 2.2 TeV has been obtained via the pp → h → νN γ with subsequent decay N → νγ. Notice that in deriving this bound the Authors of [30] have considered also the Higgs decay to be triggered by the d = 6 operators of Eq. (6.1). At d = 6 it is also interesting to notice the presence of four fermions operators. In the case of future Higgs Factories of particular relevance is which triggers a direct production channel for the RH neutrinos pair with a rate σ(ee → N N ) where β = 1 − for m N = 10 GeV. However, unless the coefficients of the d = 5 operator O N H and the one of the d = 6 operator O N e have a different scaling, for equal NP scale Λ the production induced by the former is always dominant with respect to the latter for regions where the additional Higgs decay width is compatible with current constraints, Λ 50 TeV. The situation can be drastically different in the case where the underlying theory has a particular symmetry, as in the case of Minimal Flavor violation recently analyzed in [50]. In this case the Higgs operator receives an extra suppression of a factor ∼ m N /Λ, while four fermions operators as O N e will not be affected by the insertion of any spurion. In this case the main production mechanism, especially for low neutrino masses and large center of mass energies, will be via the one through the d = 6 operators.

Conclusions
The observed pattern of neutrino masses and oscillation parameters require extending the Standard Model. One of the simplest possibilities is to add to the SM particle content two or more RH neutrinos. In this framework active neutrino masses compatible with current experimental measurements are generated via the see-saw mechanism, through an interplay of the active-sterile Yukawa coupling and the RH neutrinos Majorana mass. Naturalness consideration and the observation of a large baryon asymmetry in the Universe, motivates the study of scenarios where RH neutrinos have a mass M N at around the EW scale v, hence testable at current and future collider experiments. In general, such RH neutrino states are assumed to be produced through mixing with the SM sector, which is also responsible for their decay. The presence of additional New Physics at a scale Λ v, M N can drastically modify their phenomenology and have thus a huge impact for present and future experimental search strategies. The impact of these extra deformations can be parametrized at low energy as an effective field theory with d > 4 operators built out from SM and RH neutrino fields. In this work we have focused on the effect of the two new d = 5 operators containing both the SM and the RH neutrinos, O N H and O N B .
They induce additional production modes for RH neutrino pairs through the decay of the Higgs and the Z boson respectively. We have then studied the phenomenology of these two operators at future Higgs Factories, such as FCC-ee, CEPC, ILC and CLIC-380, for different regimes of RH neutrinos lifetimes. In particular we have considered RH neutrinos with prompt, displaced or outside the detector decays. For both operators we have shown that, in favorable scenarios for detector performance and systematics, future Higgs Factories have a great potential in testing the NP scale Λ, well above the limit that can be set by indirect probes such as the search for additional untagged Higgs decay or the measurement of the Z boson decay width. In particular, in the case of prompt decays the Higgs and Z branching ratios into N N pairs can be tested at the level of 10 −3 and 10 −8 , respectively, in the most favorable scenarios, while search for displaced decays can probe active-sterile mixing angles down to the see-saw limit of |θ e | 2 ∼ 10 −13 . We have moreover discussed the possibility of disentangling the underlying flavor structure should a signal be observed in future experiments, and commented on possible additional signatures that can be triggered by d = 6 operators, and that deserve a separate and dedicated study.