NA62 sensitivity to heavy neutral leptons in the low scale seesaw model

The sensitivity of beam dump experiments to heavy neutral leptons depends on the relative strength of their couplings to individual lepton flavours in the Standard Model. We study the impact of present neutrino oscillation data on these couplings in the minimal type I seesaw model and find that it significantly constrains the allowed heavy neutrino flavour mixing patterns. We estimate the effect that the DUNE experiment will have on these predictions. We then discuss implication that this has for the sensitivity of the NA62 experiment when operated in the beam dump mode and provide sensitivity estimates for different benchmark scenarios. We find that the sensitivity can vary by almost two orders of magnitude for general choices of the model parameters, but depends only weakly on the flavour mixing pattern within the parameter range that is preferred by neutrino oscillation data.


Introduction
All fermions in the Standard Model (SM) of particle physics with the exception of neutrinos are known to exist with both, left handed and right handed chirality. A particularly strong motivation for the existence of right handed neutrinos ν R comes from the fact that they can explain the light neutrino flavour oscillations via the type I seesaw mechanism [1][2][3][4][5][6]. If right handed neutrinos exist, they could, in addition to the Dirac mass term ν L m D ν R that is generated by the Higgs mechanism, also have a Majorana mass ν R M M ν c R , which is forbidden for the SM fields by gauge invariance. Here m D = vF , where v is the Higgs vacuum expectation value and F is a matrix of Yukawa couplings between the right handed neutrinos ν R , the Higgs doublet φ and the charged lepton doublet L = (ν L , e L ) T . The scale of the eigenvalues M i of M M is entirely unknown; different choices can have a wide range of implications for particle physics, astrophysics and cosmology, see e.g. [7] for an overview. Also the number n of right handed neutrino states ν Ri is unknown. The minimal number n required to explain the data from neutrino oscillation experiments is n = 2, as two light neutrinos are known to have non-zero masses. The light neutrinos obtain their masses via a quantum mechanical mixing of the states ν R with the left handed neutrinos ν L . This mixing can be characterised by the entries of a matrix θ = m D M −1 M . It does not only give a small mass term m ν = −m D M −1 M m T D = −θM M θ T to the light neutrinos, but also a θ-suppressed weak interaction to the heavy mass eigenstates N i ν Ri + θ ai ν c La + c.c., defined in (2). The N i are a type of heavy neutral lepton (HNL) with a mass M i that can be searched for experimentally.
Another motivation for the existence of the ν R comes from cosmology. The Yukawa couplings F ai generally violate CP, and the interactions of the ν R in the early universe can potentially generate a matter-antimatter asymmetry in the primordial plasma. At temperatures above T sph = 130 GeV [8] this asymmetry can be converted into a net baryon number by weak sphalerons [9]. This process called leptogenesis can either occur during the freeze out and decay of the ν R [10] ("freeze out scenario") or during their production [11][12][13] ("freeze in scenario"). It is one of the most promising explanations for the baryon asymmetry of the universe (BAU), which is believed to be the origin of baryonic matter in the present day universe, see [14] for a discussion. The freeze in scenario is e.g. realised in the Neutrino Minimal Standard Model (νMSM) [12]. It is particularly interesting from a phenomenological viewpoint because it is feasible for masses M i as low as 10 MeV [15], which are well within reach of present day experiments.
The requirement to explain the light neutrino masses imposes a lower bound on combinations of the matrix elements θ ai because the differences between the eigenvalues m 2 i of m ν m † ν must coincide with the neutrino mass square differences obtained from neutrino oscillation experiments. This implies that the heavy neutrinos considered here necessarily come into thermal equilibrium in the early universe [16,17]. To avoid tension with the observed abundances of light elements in the intergalactic medium, they should be sufficiently short lived that their decay does not disturb the primordial nucleosynthesis [18]. Together with experimental constraints [19,20] this imposes a lower bound of roughly M i > 100 MeV on those N i that participate in neutrino mass generation. 1 The NA62 experiment [28] can probe precisely this mass range. The experiment can be operated in two different modes, the kaon mode and the beam dump mode, as will be discussed in Section 3.1. It is sensitive to heavy neutrinos that are produced in weak decays [29,30] of mesons or tauons [31] in both modes.
In the kaon mode the SPS 400 GeV proton beam hits a fixed target and produces a secondary positively charged hadron beam with a momentum of 75 GeV, which comprises to ∼ 6 % of charged kaons. Charged kaons are then used to search for the rare kaon decay K + → π + νν, which is the primary purpose of the experiment. Among the much more frequent decays K + → + a ν a , a small fraction (O θ 2 ) of the neutrinos appears as a heavy neutrino. The resulting decay K + → + a N i can be detected with a search for a peak in the missing mass distribution [32], even if the N i itself is not seen. Such peak searches have been performed and are known to reach deep into the region [33] where leptogenesis is feasible in the νMSM [15,34,35]. An interesting point in this context is the extremely good mass resolution, of order of few MeV, which would allow to test the leptogenesis hypothesis [36] in part of the parameter space [37][38][39] if any HNLs are discovered. The disadvantage is, however, that only N i that are lighter than kaons can be produced. In this mass range there exist strong constraints from past experiments, in particular PS191 [40], CHARM [41], NuTeV [42,43], E949 [44], PIENU [45], TRIUMF-248 [46] and NA3 [47].
In the dump mode the target is removed and the beam is directly sent into a dump, producing large numbers of mesons and leptons of all sorts, which can decay into final states involving N i . The N i travel in the direction of the beam with a certain angular distribution, and a fraction of their decays into charged particles can be seen in the detector. The disadvantage is that the detection relies on the N i decay, resulting in an event rate that is of O θ 4 in the regime where NA62 is sensitive. 2 On the other hand, one is not restricted to N i below the kaon mass because the N i can be produced in the decay of any meson. This in principle extends the range of M i that can be probed up to the B-meson mass, 3 though the number of produced B-mesons is so small that the D-meson mass is practically the upper limit for the range of M i that can be probed with NA62. In this paper we study the sensitivity of the NA62 experiment when operated in dump mode to heavy neutrinos, and how it depends on assumptions about the heavy neutrino flavour mixing patter in view of the most recent neutrino oscillation data. 2 The number of heavy neutrinos that are produced along with a lepton of flavour a can be estimated as ∼ σν U 2 a , where σν is the production cross section for light neutrinos. Using definitions (6), the number of events with lepton flavour b in the final state that can be seen in a detector can then be estimated as i is the HNL decay rate (see [48] for recent more precise computations) and L is the length of the detector. For LΓ 1 one can expand the exponential and the event rate is of order O θ 4 .
3 The IceCube experiment can cover a similar mass range as NA62 with searches for double cascade events [49], though with a somewhat lower sensitivity. A significant improvement in the entire mass range below the B-meson mass could be achieved with the proposed SHiP experiment [50]. Part of this mass range could also be probed by a FASER [51], CODEX-b [52] or a MATHUSLA-type [53,54] upgrade to the LHC in the future. Above the B-meson mass, displaced vertex searches at high energy hadron [55][56][57][58][59] or lepton [39,57] colliders would be more sensitive, cf. [60] for a summary. Neutrinos that are heavy enough to decay promptly can leave distinct lepton number and flavour violating signatures in high energy collisions, see [61,62] for a recent review.

Heavy Neutrinos in the type-I Seesaw Model
The most general renormalisable addition to the SM Lagrangian that can be constructed from SM fields and the ν R reads Here ε is the antisymmetric SU(2) tensor and we have suppressed all SU(2) indices. We work in this minimal model, which adds no other New Physics to the SM. This is literally the case in the νMSM, to which our results can directly be applied. It is also valid in models where all other New Physics a) involves only new particles that are much heavier than 27 GeV (the centre of mass collision energy of the NA62 experiment) and b) does not significantly contribute to neutrino mass generation. After electroweak symmetry breaking the Higgs field obtains an expectation value v = 174 GeV, which generates the Dirac mass term ν La (m D ) ai ν Ri from the term F ai a εφν Ri . Due to the seesaw hierarchy M M m D we expand all quantities to second oder in the θ ai in what follows. The light and heavy mass eigenstates after electroweak symmetry breaking are then respectively. Here V ν = 1 − 1 2 θθ † U ν is the light neutrino mixing matrix that diagonalises the light neutrino mass matrix and U ν is its unitary part, known as Pontecorvo-Maki-Nakagawa-Sakata matrix. V N = 1 − 1 2 θ T θ * U N is the equivalent matrix that diagonalises the heavy neutrino mass matrix M N = M + 1 2 θ † θM + M T θ T θ * after electroweak symmetry breaking, and Θ = θU * N . In practice we can set V N 1 and Θ θ. The mixing gives a θ-suppressed weak interaction to the heavy mass eigenstates (2). The coupling of the mass eigenstates N i to the SM can then be described by the Lagrangian contribution The first two lines are the couplings of the N i to the weak currents due to the mixing Θ, and the last line is the Yukawa coupling to the physical Higgs field h in the unitary gauge expressed using the definition of θ and the relation m W = 1 2 vg, with g the weak gauge coupling constant. This Yukawa term is not relevant for NA62 because the branching ratios for decays mediated by virtual Higgs bosons are suppressed due to the small Yukawa couplings of the kinematically accessible electron and muon final states. Due to the interactions (4) the heavy neutrinos can participate in all processes that involve ordinary neutrinos if this is kinematically allowed, but  Table 1: Definitions of the neutrino masses and their differences, for "normal ordering" (NO) and "inverted ordering" (IO). The smaller mass difference ("solar mass difference") is given by ∆m 2 sol = m 2 2 − m 2 1 for both hierarchies, while the larger mass difference is hierarchy depended. For historical reasons one also defines the "atmospheric mass difference" ∆m 2 atm = m 2 3 − m 2 1 . One can roughly identify ∆m 2 atm and ∆m 2 sol as the larger and smaller mass splitting, respectively. The values for m 2 atm differ slightly for the two hierarchies and the difference is order ∆m 2 sol /∆m 2 atm . For n = 2 the lightest neutrino is massless (m lightest = 0), which allows in both cases to identify ∆m 2 atm with one of the neutrino masses.
with amplitudes suppressed by Θ ai . The event rates are proportional to combinations of the quantities In the symmetry protected scenario defined further below the heavy neutrinos tend to have quasi-degenerate masses. If the mass splitting is smaller than the experimental resolution, it is not possible to resolve them individually. It is convenient to introduce the notation If the ν Ri are the sole origin of the light neutrino masses, then the properties of the N i are strongly constrained by relation (3). This connection has been studied by various authors, see e.g. [16,19,20,31,37,38,[63][64][65]. The constraints depend on the number n of heavy neutrinos. In general, the number of parameters in the seesaw Lagrangian (1) in addition to those in the SM is 7n − 3. n of them can be identified with the masses M i , the remaining ones are mixing angles and phases. We discuss the cases n = 2 and n = 3 below. The n = 2 scenario is the minimal model that allows to explain the two observed neutrino mass differences and predicts the lightest neutrino to be massless (m lightest = 0). It also effectively describes the neutrino mass generation in the νMSM, which in principle contains three heavy neutrinos. The lightest N i in the νMSM is a Dark Matter candidate, and the strong observational constraints on its interactions imply that it has no significant effect on the light neutrino masses, see e.g. the review [26] and references therein. The n = 3 scenario is the minimal model consistent with neutrino oscillation data if m lightest > 0; in many models with extended gauge sectors it is also motivated by the requirement of anomaly freedom and the fact that there are three generations in the SM.

Connection to light neutrino oscillation data
The smallness of the light neutrino masses m i can be explained by the seesaw relation (3) in different ways. One possibility is that the N i are superheavy. In this case the smallness of  Table 2: Best fit values of neutrino mass differences and mixing angles from the NuFIT 3.1 release by the ν-fit collaboration [77,78], for "normal ordering" (NO) and "inverted ordering" (IO). The mass differences are defined in Table 1.
This conventional seesaw mechanism cannot work for the masses M i < v that are accessible to NA62. If one simply "downscales" the M i to experimentally accessible values, one would roughly expect |F ai | 2 ∼ M i m a /v 2 and U 2 ai ∼ m a /M i , where the symbol "∼" indicates that we have neglected light neutrino mixing, i.e., we neglect the difference between mass and flavour eigenstates and assume that the flavour eigenstate ν La has the mass m a of the mass eigenstate that it is mostly composed of. If we would further assume that all N i have the same mass M , then we can define the parameters to quantify the deviation of the Yukawa couplings and mixings from the "naive seesaw expectation". Here ∆m 2 atm is the larger of the two observed neutrino mass splittings, with ∆m 2 sol being the smaller one, cf. Table 1. This would be very discouraging for experimentalists because the production rates for HNLs with U 2 U 2 0 is tiny even if the seesaw scale M 1 is as low as 100 MeV.
Another possibility is that the m i are small as the result of a slightly broken symmetry. Comparably low values of the M i are technically natural because the B − L symmetry of the SM (which is in general broken by the M i ) is restored in the limit of all M i → 0. An approximate B − L symmetry can be realised even if M i = 0 if the ν Ri come in pairs with equal mass M i = M j and couplings F ai = iF aj that can be represented by Dirac spinors ν Ri + ν c Rj [63,66,67]. This symmetry protected scenario provides a theoretical motivation for a low scale seesaw and allows for experimentally accessible U 2 ai U 2 0 . Specific examples that motivate this limit include in "inverse seesaw" type scenarios [68][69][70][71], a "linear seesaw" [72,73], scale invariant models [74], some technicolour-type models [75,76] or the νMSM [66].
A convenient way to connect the heavy neutrino couplings to neutrino oscillation data is provided by the Casas-Ibarra parametrisation [79] Here We use the common parameterisation with U ±δ = diag(1, e ∓ i δ/2 , e ± i δ/2 ), and non-vanishing entries of V (ab) for a = e, µ, τ are Here θ ab are the light neutrino mixing angles. For the purpose of fixing the parameters in U ν , we approximate V ν U ν and use the data set NuFIT 3.1 from the ν-fit collaboration [77,78], cf. Table 2. The NuFIT 3.2 update released shortly after our work was finalised. We checked that our results (in particular the benchmark scenarios listed in Table 3) are not significantly affected by the update, cf. Figures 3 and 4.
In addition to the neutrino oscillation data, N i with masses below the electroweak scale are constrained by a number of direct and indirect searches, see e.g. [80,81] and references therein. With the exception of neutrinoless double β decay [20,37,[82][83][84][85][86], the indirect observables are mostly irrelevant in the mass range that can be tested by NA62.

The Minimal Model with n = 2
For n = 2 the matrix R in the parameterisation (8) can be expressed in terms of a single complex parameter ω, where ξ = ±1 and "NO" and "IO" refer to "normal ordering" and "inverted ordering" amongst the light neutrino masses m i . Moreover, only one combination of the Majorana phases is physical. For normal ordering this is α 2 , while for inverted ordering it is α 2 − α 1 . We can therefore without loss of generality set α 2 ≡ α and α 1 = 0. The symmetry protected regime in this parametrisation (and for ξ = 1) corresponds to small values of the symmetry breaking parameters In terms of these parameters one can express the total interaction strength U 2 of the heavy neutrinos as where ∆m = 1 2 While µ can in principle be set exactly to zero, there are theoretical lower limits on from the requirements to remain perturbative (|F ai | 4π) and justify the expansion in θ (U 2 1). Practically the lower limit on from the existing experimental constraints in the mass range accessible to NA62 is stronger than these theoretical considerations. The relations between the individual U 2 ai and the parameters in U ν are given in Appendix A. We define the symmetric limit by setting µ = 0 (i.e.   Table 3a. particles for heavy neutrinos necessarily require 1. 4 This requires cancellations in the light neutrino mixing matrix (3), which should be considered "fine tuned" unless one also takes the limit µ 1, in which case a generalisation of the global B − L symmetry in the SM is approximately respected by the Lagrangian (1). Hence, the condition µ 1 that is required for leptogenesis with n = 2 [12] is automatically fulfilled in a (technically) natural way in the symmetric limit. The limit µ, 1 yields A remarkable feature of the limit → 0 is that the ratios U 2 a /U 2 become independent of the seesaw scale M and the unknown parameter Re ω irrespective of the choice of µ. That is, they can be expressed in terms of the low energy parameters in U ν alone. On one hand, this in principle allows one to determine the Majorana phase α from measurements of the U 2 a [87]. On the other hand this means that existing neutrino oscillation data can be used to make predictions for the allowed range of U 2 a /U 2 in the minimal seesaw model. A simple estimate can be made if one fixes the light neutrino mass splittings and mixing angles to the best fit values given in Table 2 and freely varies all other parameters, see 10 -10 (a) Normal ordering. 10 -10 (b) Inverted ordering. Figure 2: The applicability of the symmetric limit can be illustrated by plotting the allowed These regions correspond to the interval predicted in the symmetric limit.
the entire hashed regions are allowed, NA62 can only find N i with mixing angles which are large enough that the symmetry protected limit can be applied, cf. Figure 2. Hence, if NA62 finds HNLs with U 2 a /U 2 outside the filled regions, this would clearly disfavour the minimal model with n = 2.
The present status of neutrino oscillation experiments allows to do a more quantitative analysis. One can use the statistical information about the light neutrino parameters gathered in various neutrino oscillation experiments to obtain a probability distribution for the U 2 a /U 2 . As input we use the results of the NuFIT 3.1 release for this purpose, from which one can first determine likelihoods for ∆m 2 31 , ∆m 2 32 and all parameters in U ν with the exception of the Majorana phase α. Due to the large deviations from the Gaussian limit, for δ and θ 23 we use the two-dimensional ∆χ 2 projection, while for all other parameters we use the one-dimensional projections. The total ∆χ 2 for a set of parameters is then given by the sum of the individual ∆χ 2 .
To identify how favoured each choice of parameters is according to the current neutrino oscillation data we generate a sample of points using a simple implementation of the Metropolis-Hastings algorithm [88,89]. We start with a random choice of the low-energy parameters x 0 = (θ 12 0 , θ 13 0 , θ 23 0 , ∆m 2 210 , ∆m 2 3l0 , δ 0 ). Each following point in the Markov chain is chosen by generating a candidate point x from a multivariate normal distribution centred around the last point in the chain x i . The candidate point is then accepted with the probability       If the point is accepted, it becomes the next point in the Markov chain, x i+1 = x , otherwise we keep x i+1 = x i . In the limit of large i, the density of points reflects the probability distribution in the parameter space in view of the experimental data. Since α is experimentally unconstrained we have to pick an a priori distribution "by hand" from which we chose the samples. We generate the samples of α in two different ways. One choice of "prior" is a flat distribution in α between 0 and 4π. In the second approach we use a distribution that is "flat" in U 2 e /U 2 . This is done by choosing α from a flat distribution in sin(α/2 + δ) for normal hierarchy and sin(α/2) for inverted hierarchy according to relations (21) and (23). This in principle brings in a dependence on the assumptions (or prior) that one imposes on the value of α. In Figures 3 and 4 we show the likelihoods for the ratios U 2 a /U 2 for the two different choices of prior. The fact that the allowed regions are rather independent of the prior imposed on the unknown parameter α shows that they reflect actual experimental uncertainties (rather than theoretical prejudice on the model). We compare the results based on the NuFIT 3.1 and NuFIT 3.2 releases. Though the 3.2 update leads to a visible change in the likelihoods for some of the parameters in U ν (in particular θ 23 and δ), it does not lead to a significant change in the preferred choice of benchmark scenarios.
We finally check whether the above conclusion may change if constraints from searches for neutrinoless double β decay are taken into account. For µ, 1, the effective Majorana mass m ββ that governs the rate of the decay in the n = 2 model can be approximated by [20,37,86]  in the decay. Λ 2 is the momentum exchange in the decay, specific values can e.g. be found in [90]. For µ = 0 these expressions reduce to m ββ = 1 − f A (M ) m ν ββ , which is always smaller than m ν ββ [82]. Hence, the present non-observation of neutrinoless double β decay cannot directly constrain the U 2 a /U 2 . For finite µ the correction to m ββ is not necessarily small because it scales as ∝ µ/ . This implies that for 1, small values of µ are favoured by neutrinoless double β decay. However, since we previously established that the range of allowed U 2 ai /U 2 i for U 2 i within reach of NA62 is in good approximation identical to the range that is allowed for µ = 0, present constraints from neutrinoless double β decay do not further reduce the allowed regions in Figure 1.
In summary, the model with n = 2 is very predictive as far as the flavour mixing pattern is concerned. Out of the 11 new parameters in addition to the SM, five are fixed by neutrino oscillation data (cf. Table 2). In the symmetric limit, the mass splitting ∆M vanishes and the phase Re ω does not affect the U 2 a at leading order in . By forming the ratios U 2 a /U 2 , also the dependence on M can be eliminated, so that these ratios can be determined from fixing the Dirac phase δ and Majorana phase α in U ν alone. Since δ is constrained by light    Figure 3 with DUNE [91]. Here we assume Gaussian errors around the parameter δ = −π/2 ± π/9 used as benchmarks in [92]. For all other parameters we take the one dimensional χ 2 projections from NuFIT 3.1 (a) and NuFIT 3.2 (b). We assume that all values of α are equally valid. In reality we expect even stronger constraints as all low-energy parameters will be measured with a higher precision, we here have only taken into account the expected improvement in δ.
neutrino oscillation data, there is only one truly unknown parameter α. Figures 3 and 4 show that the predictions for the allowed range of U 2 a /U 2 are rather independent of the theoretical prejudice about the value of α. Hence, they can be used to define benchmark scenarios that mark the corners of the allowed flavour mixing patterns for n = 2. Figure 5 shows how these constraints are expected to improve with the DUNE experiment. 5

The Model with n = 3
The model with n = 2 is highly predictive because of its minimality. For n = 3 the flavour mixing patterns are far less restricted. There is no lower bound on the individual U 2 ai from neutrino oscillation data [65,81]. It is, for instance, possible to set F e1 = F µ1 = 0 by fixing the entirely unconstrained Majorana phases α 1 and α 2 for arbitrary choices of all other parameters, including m lightest . 6 Of course, in this case N 2 and N 3 have to mix with the first two SM generations in order to explain the mixing of light electron and muon neutrinos. However, M 2 and M 3 can be outside the reach of NA62, so that the only observable heavy neutrino can 5 NA62 is expected to probe the range of U 2 where the symmetric limit can be applied to determine the allowed values for U 2 a /U 2 . It is, however, worthwhile noting that the excellent mass resolution of the experiment allows to resolve the HNL masses for µ < 10 −2 [33]. A measurement of the individual U 2 ai would, in combination with a measurement of δ in neutrino oscillation experiments, in principle allow to extract all parameters in the Lagrangian (1) for n = 2, making this a fully testable model of neutrino masses and baryogenesis [38]. 6 For m lightest = 0 it is even possible to set Fe1 = Fµ1 = Fτ1 = 0, which effectively reduces the model with n = 3 to the model with n = 2. exclusively mix with the third generation, which is impossible for n = 2. This implies that it is not possible to make any reliable predictions for the values of the U 2 ai in general. However, observably large U 2 ai still require cancellations in the seesaw relations (3). This is only possible without significant tunings if the m i are protected by a symmetry. The implementation of the same generalised B − L symmetry discussed in Section 2.3 requires that two of the N i (which we may choose to call N 2 and N 3 ) effectively behave like the two heavy neutrinos in the n = 2 model, while the third one has much smaller mixings ∼ m lightest /M 1 . This is precisely the behaviour that is observed in the νMSM, where the feebly coupled neutrino is a DM candidate. Since only N 2 and N 3 can realistically be observed at NA62, the considerations from the previous Section 2.3 also apply to the n = 3 case. The only difference is that, while these are unavoidable predictions in the case n = 2, they can be circumvent by "tuning" in the parameters in the scenario with n = 3 if one chooses to explain the small m i by accidental cancellations in relation (3).

The NA62 Experiment
The main goal of the NA62 experiment [28] which is currently taking data at the CERN SPS is to measure the Branching Ratio (BR) of the K + → π + νν decay with a precision of at least 10 %. In order to achieve this goal the experiment needs to collect about 10 13 kaon decays of which O few × 10 12 have already been collected in the current run [93].
In its normal operation mode, the kaon mode, the primary 400 GeV proton beam impinges on a 400 mm long cylindrical beryllium target with a diameter of 2 mm which is used to produce a secondary positively charged hadron beam with a momentum of 75 GeV. 100 m downstream of the target the secondary beam reaches the 120 m long evacuated decay volume which has a diameter of 2 m. About 6 % of the hadron beam are kaons, which are identified and timestamped by a N 2 filled Cherenkov counter located along the beam line. Three silicon pixel stations measure momentum and time of all the particles in the beam at a rate of 750 MHz. A guard ring detector tags hadronic interactions in the last pixel station at the entrance of the decay volume. Large angle electromagnetic calorimeters made of lead glass blocks surrounding the decay vessel are used to veto particles up to 50 mrad. A magnetic spectrometer made of straw tubes in vacuum measures the momentum of the charged particles. A Ring-Imaging Cherenkov (RICH) counter filled with Neon separates π, µ and e for momenta up to 40 GeV. The time of flight for charged particles is measured both by the RICH and by the scintillator hodoscopes placed downstream of the RICH. An electromagnetic calorimeter covers the forward region and complements the RICH for the particle identification. The hadronic calorimeter provides further separation between π and µ based on hadronic energy and a fast scintillator array identifies muons with sub-nanosecond time resolution.
At the nominal beam intensity of 3 × 10 12 protons per pulse, with pulses of 4.8 s, the NA62 experiment can collect up to 3 × 10 18 protons on target (POT) per year.
When the experiment is operated in the dump mode, the target is pulled up and the primary proton beam is send directly onto the Cu-Fe based collimators that act as a hadron stopper (or dump) located 20 m downstream of the target. In this configuration, about 2 × 10 15 D-mesons and ∼ 10 11 b-hadrons are produced from the 10 18 POT, which correspond to about 80 days of data taking at the nominal NA62 beam intensity. This dataset will be collected during Run 3 (2021Run 3 ( -2023 and is assumed to produce the sensitivity plots discussed in Section 3.   Table 3: We consider six benchmark scenarios A)-F) listed in Panel (a). For comparison we also list the benchmark scenarios used in [50] for the analysis for the SHiP experiment in Panel (b).

Benchmark scenarios
As for any fixed target experiment, the sensitivity of NA62 can only be computed for fixed ratios U 2 ei : U 2 µi : U 2 τ i . In the symmetric limit this can be used exchangeably with U 2 e : U 2 µ : U 2 τ , cf. relation (14). We consider the benchmark scenarios listed in Table 3a. The scenarios A)-D) extremise ratios of U 2 a in the symmetric limit. Scenario A) minimises the mixing with the first generation and maximises it with the second generation. In the n = 2 model it can be realised for (α, δ, ξ) = (−π, π, 1) with NO. Scenario B) maximises the mixing with the first generation and minimises it with the second generation for NO. In the n = 2 model it can be realised for (α, δ, ξ) = (π, π, 1). Scenario C) minimises the mixing with the first generation for IO. In the n = 2 model it can be realised for (α, δ, ξ) = (−π, π, 1). The scenarios are marked by stars in Figure 1. We have not included the other "corners" of the allowed regions because they either require δ = 0 for n = 2, which is disfavoured by current neutrino oscillation data, or are phenomenologically very similar to one of our scenarios. We practically use scenario D) to model the maximal mixing with the first generation for IO.
The scenarios A)-D) may be compared to the choices 1)-5) listed in Table 3b and used in the estimates for the SHiP experiment [94,95] in reference [50], which were introduced in references [31,34]. Scenarios 1)-3) maximise U 2 e , U 2 µ and U 2 τ for NO in view of the constrains from neutrino oscillation experiments at the time when they were proposed. However, such large U 2 τ are disfavoured by more recent data, while a stronger contributions from U 2 e and U 2 µ now seem to be allowed. Scenario 4), which at the time was motivated by IO, is very similar to scenario 1), while scenario 5) is similar to our scenario C). Hence, the scenarios A)-D) can be seen as updates to the scenarios 1)-5) in view of the most recent neutrino oscillation data (in particular the measurement of θ 13 , which was unknown at the time when the articles in references [31,34] were written). The scenarios D)-F) are the extreme cases in which a N i couples only to one SM generation. Scenarios E) and F) are not allowed in the n = 2 model for U 2 ai within reach of NA62, and in the n = 3 model they can only be realised with significant tuning. It is nevertheless instructive to include them in order to understand the most optimistic and most pessimistic predictions for the NA62 sensitivity without theoretical prejudice. As mentioned before, scenario D) can almost be realised for IO and may therefore be classified as "allowed".

NA62 sensitivity computation
The computation of the NA62 sensitivity for U 2 ia (a = e, µ, τ ) in different scenarios has been performed using a toy Monte Carlo in which all the kinematics of the N i production and decay processes have been implemented and the geometrical acceptance for the decay products has been properly evaluated using the geometry of the experiment as described in [28].
Heavy neutrinos can be produced in hadron decays where a SM neutrino ν a is replaced by a N i through the mixing θ ai . At the SPS energies this happens via production of strange, charm and beauty hadrons. Above the kaon mass, the main production processes are the decays of charm and beauty hadrons which are kinematically allowed. The number n N of N i produced in a beam-dump experiment can be quantified by (17) where N POT is the total number of protons on target, χ c and χ b are the ratios of the production cross sections of c and b quarks with respect to the total pp cross section for a thick target 7 , f D j and f B k are the production fractions of charm or beauty mesons and finally BR (D j → XN i ) and BR (B k → XN i ) are the branching fractions of charm or beauty mesons into heavy neutrinos 8 , which depend on their mass M i and on the coupling parameter U 2 i . In the sensitivity computation we assume 10 18 POT, as discussed in Section 3.1. The cand b-hadrons can be created in the hadronisation process originated from primary protons, and from all secondary products of the hadronic shower in the dump, such as protons, neutrons, and pions. Given the ratio χ c /χ b ∼ 10 4 , the N i production via charm decays is the dominant process up to the D-meson masses, while the contribution from b-hadrons decays at the NA62 intensity is almost negligible. The composition of the shower and the kinematics of the produced c-and b-hadrons have been studied by simulating the 400 GeV proton beam on a thick (∼ 11 interaction lengths) high-Z target with Pythia 6.4 [96].
Once produced, the N i decay to SM particles via their θ-suppressed weak interactions. These massive neutrino states can decay to a variety of final states through charged-and neutral-current processes. The main decay channels of N i with masses below a few GeV are where = e, µ, τ . The branching fractions as a function of the N i mass are shown for scenario A) in Figure 6. The relative weight of the branching fractions depends on the scenario under consideration. The analytical formulae can be found in [31]. The NA62 detector is able to reconstruct all final states with two charged tracks. The number of events reconstructed in the NA62 detector is given by where n N,I is the number of produced N i for a given production process I, via decays of charm and beauty hadrons as shown in (17),  for a N i produced in the process I of a given mass M i and coupling U 2 e, µ, τ which decays into a final state with two charged tracks (f + f − ) and other decay products X as, for example, photons and neutrinos. The efficiency ε(f + f − X) is the product of the trigger, reconstruction and selection efficiencies for a given final state, and is currently assumed to be 100 %. The curves shown in Figure 7 can therefore slightly change when all the experimental effects are taken into account. 10 The position of the decay vertex of the N i is required to be inside the 75 m long fiducial volume, located 80 m downstream of the dump, and within a circle with 1 m radius in the transverse plane. The decay products are required to be within acceptance of the charged hodoscope. The geometrical acceptance is given by the convolution of the probability that the N i decays inside the fiducial volume with the probability that the two charged tracks in the final state are reconstructed in the magnetic spectrometer.
Monte Carlo techniques are used to assess the sensitivity contours in the mass-coupling parameter space. This is done by generating heavy neutrinos with different mass and coupling values, letting them decay in the kinematically accessible final states, and computing the corresponding geometrical acceptance.
during a ∼ 10 hours long run in the dump mode using 40 % of the nominal beam intensity about 2 × 10 15 POT were recorded and 18 kHz of muons have been measured within the NA62 acceptance. A preliminary study of the background rates and topologies has been performed.
The analysis of this dataset shows that, within the current available statistics, NA62 can reduce the background to zero for all hidden sector particles decaying into fully reconstructed final states with two charged tracks. In order to control the backgrounds of partially reconstructed final states or final states containing photons an upgrade of the current apparatus is required. A detailed discussion of the background observed in NA62 when operated in the beam dump mode can be found in [103]. Up to date, the dataset collected by NA62 in the beam-dump mode is about a factor 10 larger and the analysis is under progress.

Discussion and Conclusion
We have studied the sensitivity of the NA62 experiment in the dump mode to heavy neutral leptons in the low scale seesaw model. The NA62 sensitivity was obtained for 10 18 POT which is the data sample expected to be collected in 2021-2023. This dataset corresponds to 3 months of data taking at full intensity and will be spread over said period of three years in order to allow NA62 to complete the kaon programme. If the kaon programme will be completed earlier, more time can be dedicated to the beam dump mode, depending on the collaboration strategy. At full intensity NA62 can collect up to 3 × 10 18 POT per year.
Our main results are shown in Figure 7a. Below the kaon mass, the accessible parameter space is already strongly constrained by past kaon and pion decay experiments. Heavy neutrinos that are heavier than D-mesons can only be produced in B-meson decays at NA62, and their number is too small to give sizeable event rates. As shown in Figures 7b-d NA62 is the only existing experiment that can probe significant fractions of unexplored parameter space for heavy neutrinos that are lighter than D-mesons. Moreover, based on current estimates [38], it appears to be the only existing experiment that can enter the parameter region where leptogenesis is possible in the νMSM.
As one can see in Figures 7a and 8, the sensitivity can vary by about two orders of magnitude, depending on the "flavour mixing pattern", i.e., the way how the total mixing U 2 i = a U 2 ai of a heavy neutrino N i is distributed amongst the different flavours a = e, µ, τ . The primary reason is the tauon mass, which kinematically blocks the decay of D-mesons into tau and heavy neutrinos with masses above the kaon mass. Hence, production of N i via U 2 τ i is only possible via τ decays and the decay of N i via U 2 τ i is only possible through neutral current interactions. As a result, the sensitivity primarily depends on how U 2 i is distributed amongst U 2 τ i and U 2 ei + U 2 µi . For M i above the di-muon threshold, the sensitivity is essentially independent of U 2 ei /U 2 µi . This covers the entire mass range above the kaon mass, where NA62 can probe unexplored parameter space. Remarkably, the sensitivity of NA62 in this regime shows a strong dependence on the flavour mixing pattern only if U 2 i is dominated by U 2 τ i , cf. Figure 8.
The flavour mixing pattern in the type-I seesaw model is strongly constrained by neutrino oscillation data. These constraints depend on the number n of heavy right handed neutrinos that generate the light neutrino masses. For n = 2 and U 2 i that can be probed with NA62, global fits to current neutrino oscillation data constrain the ratios U 2 ai /U 2 i to the filled areas in Figure 1. These clearly exclude a strongly U 2 τ i dominated scenario. Hence, the sensitivity in the minimal model with n = 2, which effectively also describes the νMSM, is rather independent  Table 3a, which mark the most extreme patterns that are allowed for n = 2 (cf. stars in Figure 1). They are represented by the first four lines in Figure 7a.
In Figures 3-4 we perform a more detailed analysis to identify the likelihood for different flavour mixing patterns U 2 ai /U 2 i , based on present neutrino oscillation data. Remarkably, these likelihoods only weakly depend on theoretical assumptions, i.e., directly reflect the statistical preference due to experimental data. That is, the though the masses M i and overall coupling strengths U 2 i of the hypothetical heavy neutrinos N i are unknown, current neutrino oscillation data already allows to make robust statements on the relative size of their couplings U 2 ai to the different SM flavours. This is the second main result of the present work, in addition to the NA62 sensitivity estimates in Figure 7. In Figure 5 we estimate how the constraints can improve if the DUNE experiment measures the Dirac phase δ in the light neutrino mixing matrix. For n = 3, mixing patterns outside the filled areas in Figure 1 can be made consistent with neutrino oscillation data, including extreme U 2 τ i dominated scenarios. However, these require considerable tunings in the parameters. The preferred range of U 2 ai /U 2 i in absence of tunings still roughly coincides with the filled regions in Figure 1.
In summary, we find that NA62 is the world's most powerful existing experiment to search for heavy neutrinos with masses between those of kaons and D-mesons. In this mass range, the sensitivity within the parameter region that is preferred by neutrino oscillation data depends only weakly on the heavy neutrino flavour mixing pattern.