First Evaluation of Meson and $\tau$ lepton Spectra and Search for Heavy Neutral Leptons at ILC Beam Dump

A beam dump experiment can be seamlessly added to the {proposed} International Linear Collider (ILC) program because the high energy electron beam should be dumped after the collision point. The ILC beam dump experiment will provide an excellent opportunity to search for new long-lived particles. Since many of them can be produced by a rare decay of standard model particles, we evaluate spectra of the mesons and $\tau$ lepton at the decay based on the PHITS and PYTHIA8 simulations. As a motivated physics case, we study the projected sensitivity of heavy neutral leptons at the ILC beam dump experiment. The heavy neutral leptons can also be produced via deep inelastic scattering and $Z$ boson decay at the ILC main detector, which we include in the projection. With the multi-track signal, the reach would be greatly extended in mass and coupling, even compared with the other proposed searches.


Introduction
The large hadron collider (LHC) has discovered the Higgs boson and measured its properties consistent with the precision measurements of the Standard Model (SM).The discovery of the last particle of the standard model strengthened the foundation of the SM.However, overwhelming evidence and hints require physics beyond the Standard Model (BSM).While the intensive searches for new particles at the weak scale or heavier have been performed experimentally, we have not found any convincing evidence of such new particles yet.In these circumstances, there is growing attention to feebly interacting light particles from both theoretical and experimental points of view.
New light particles with mass less than tens of GeV are compelling addition to the SM as they can directly resolve the central issues of the SM, such as the strong CP problem and the finite neutrino masses, explain many flavor physics anomalies, and be a portal to the dark sector accommodating the dark matter and the baryon asymmetry of the universe.On the experimental grounds, the light particles, if they exist, need to be feebly interacting with the SM sector due to various experimental constraints.However, one should be aware that the bounds are strongly model-dependent because the signature depends on the nature of the light degrees of freedom.This situation is contrasted with heavy new particle exploration with mass beyond the reach of current collider energies, which can be probed through Effective Field Theories (EFT).Therefore, there is growing interest in various scenarios based on the current experimental data and the potential for future high-intensity experiments.The status of recent progress is well-summarized in the report of Physics Beyond Collider [1].
In this paper, we study the physics potential of the dump facilities of the International Linear Collider (ILC) to hunt new light particles.The ILC is a proposed future e + e − linear collider with beam-energy of 125 GeV, which can be upgraded to 500 GeV.The main goal of the ILC is to study events of e + e − collision to perform high precision measurements of the Higgs bosons and search for new particles produced by the electroweak interaction.In the linear collider, the beam has to go into the beam dump after the collision, and it provides an excellent opportunity for a high-intensity experiment to produce feebly interacting light particles.The number of electrons on the beam dump is enormous, N EOT = 4 × 10 21 / per year, and the electromagnetic (EM) shower also leads to a highintensity photon carrying O(10%) of the beam energy.The setup is equivalent to the other fixed target experiments if a detector is placed downstream of the beam dump.Hereafter we refer to such setup as the ILC beam dump experiment.One can search for light particles produced in the dump, fly beyond the muon shield, then decay to the SM particles or scatter the detector material.The potential of the ILC beam dump project has been investigated in Refs [2][3][4][5][6][7].
In the previous works of the electron beam dump experiments, the searches for light particles dominantly interacting with electron or photon were mainly studied since they are produced by the primary electron and positron or the secondary shower particles.In the ILC beam dump experiment, the initial beam energy is much higher than those of the past electron beam dump experiments.In addition to the EM showers, heavy mesons and τ lepton can be produced by the shower photon hitting nuclei.The decay of the produced SM particles is another promising source of light particles.This paper examines the production yield and spectrum of the mesons and τ lepton at the ILC beam dump setup for the first time.
As an application of this study, we investigate a projected sensitivity of heavy neutral leptons (HNLs, called sterile neutrinos in some literature).The HNLs mix with the SM neutrinos with an angle of U ( = e, µ, τ ), resulting in a suppressed weak interaction to the SM.Therefore, it is very natural to consider the HNL production from the meson and τ lepton decays where the weak interaction dominates.The current experimental constraints are mostly for mixing with ν e and ν µ , and the high intensity of mesons would further reach the under-explored parameter space.Also the sensitivity to the τ neutrino mixing angle U τ can be significantly improved because τ lepton is accessible.
The phenomenology of the HNLs is well-reviewed in Refs [8,9], and see references therein.We stress that feebly interacting HNLs in a GeV mass range are a motivated and well-defined physics target.The seesaw mechanism can explain the neutrino masses observed by the neutrino oscillations with at least two HNLs.If the two HNLs are almost degenerate, the sum of the mixing angles have the lower bound, approximately [8,10].Furthermore, the degenerate HNLs in the early universe can produce the baryon asymmetry via the HNL oscillations [11,12].The feeble interactions are necessary for the departure of the thermal equilibrium which is one of Sakharov's conditions for generating the baryon asymmetry.Together with the neutrino mass constraints, the interesting parameter space is in a range of 10 −11 U2 10 −6 . 1his paper is organized as follows.In Sec. 2, we briefly review the ILC beam dump experiment, and in the following section, we evaluate the spectra of mesons and τ lepton.In Sec. 4, we study the HNLs at the ILC beam dump experiment.In addition to the HNL production from the SM particle decays, we consider an HNL production via deep inelastic scattering and another production from Z decays at the interaction point.We finish with the discussion in Sec. 5.

ILC Beam Dump Experiment
The ILC main beam dump has to absorb 2.6 MW (125 GeV × 21 µA) of the e ± beam energy for 125 GeV in the initial stage and 13.6 MW (500 GeV × 27.3 µA) for 500 GeV in an upgrade stage.Following a water dump designed with the length of l dump = 11 m and the diameter of 1.8 m, a muon shield of length l sh = 70 m would be placed as proposed in [4,5], to remove the secondary muon background.The cylindrical decay volume of length l dec = 50 m and radius r det = 3 m would lie between the muon shield and the downstream detector.The schematic view is seen in Fig. 1, and the similar design of the setup can be found in [4,5].Here, we additionally assume a multi-layer tracker in the decay volume.We consider a time frame of 10-year run for both ILC-250 and ILC-1000 where the beam energy is 125 GeV and 500 GeV, respectively.
The EM shower (photons, electrons, and positions) starts in the beam dump.The ILC beam dump experiment is unique compared to the past electron beam dump experiments with regards to the higher beam energy and the high intensity, i.e., the number of electrons on the beam dump N EOT is about 4 × 10 21 per year 2 .The energetic beam creates photons at O(1-10) GeV energy scale (Fig. 9  A setup for ILC beam dump experiments.It consists of the main beam dump, a muon shield, and a decay volume.We assume a multi-layer tracker is placed in the decay volume so that the charged tracks are measured.[5]), and the secondary interaction between the photon and the nucleus can produce light mesons (π and K), heavy mesons (D, B, and even B c ) and τ lepton.They loose the energy in the beam dump and finally decay, and some of them produce the HNLs.The yield and energy spectrum of the SM particles at the decay is therefore important to study the sensitivity of the ILC beam dump experiment.We show the evaluation in Sec. 3.
The secondary muons are stopped in the muon shield at ILC-250, but their penetration behind the shield cannot be neglected for E beam = 500 GeV at ILC-1000.In this case, an additional active muon shield behind the muon shield would be necessary, and we assume that the muon shield consists of the lead shield (l lead sh = 10 m) and the active shield (l active sh = 70 m − l lead sh ).The HNL with a dominant mixing with the muon neutrino can be produced inside the muon shield by scattering of the shower muon, and we approximate that the HNLs produced behind the lead shield do not contribute to the signal events.In Appendix E, we study how different depth of the muon shield affects the sensitivity for the HNL at ILC-1000.

Meson and τ lepton Spectra
This section presents the meson and τ lepton spectrum obtained by Monte-Carlo simulation at the ILC beam dump.We use PHITS 3.25 [13] for production and transport of particles other than heavy mesons.PHITS (Particle and Heavy Ion Transport code System) is a general-purpose Monte Carlo particle transport simulation code developed under collaboration between JAEA, RIST, KEK, and several other institutes.PHITS can transport most particle species for a given geometry of the materials, and it is tested thoroughly by benchmarks studies [14,15].For heavy mesons production, we implement their differential production cross sections obtained by PYTHIA8.3 [16] into PHITS.More details are given in the following.

⇡ ±
< l a t e x i t s h a 1 _ b a s e 6 4 = " b q 7 C j u F k t p z J d X q P L 3 W 9 d m Y K i R M = " > A A A B 8 H i c b V D L S g M x F M 3 4 r O O r 6 t J N s A i u S q a 2 t l 2 I B T c u K 9 i H d G r J p G k b m s y E J C O U o V / h x o U i 4 s 4 f 8 D / c i H / j d K a I i g c u H M 6 5 l 3 v v 8 S R n 2 i D 0 a S 0 s L i 2 v r G b W 7 P W N z a 3 t 7 M 5 u U w e h I r R B A h 6 o t o c 1 5 c y n D c M M p 2 2 p K B Y e p y 1 v f D 7 z W 7 d U a R b 4 V 2 Y i a V f g o c 8 G j G A T S 9 e u Z D e R K 8 W 0 l 8 2 h f L V S K h a q E O V R g h k p l B E 6 g c 5 c y Z 2 9 2 a f y 5 c O u 9 7 L v b j 8 g o a C + I R x r 3 X G Q N N 0 I K 8 M I p 1 P b D T W V m I z x k H Z i 6 m N B d T d K D p 7 C w 1 j p w 0 G g 4 v I N T N S f E x E W W k + E F 3 c K b E b 6 r z c T / / M 6 o R l U u h H z Z W i o T 9 J F g 5 B D E 8 D Z 9 7 D P F C W G T 2 K C i W L x r Z C M s M L E x B n Z a Q g J Y E r K x T m p O t 8 h N A t 5 5 z h f u k S 5 W g G k y I B 9 c A C O g A P K o A Y u Q B 0 0 A A E C 3 I E H 8 G g p 6 9 5 6 s p 7 T 1 g V r P r M H f s F 6 / Q L s / Z S A < / l a t e x i t > E beam = 125 GeV < l a t e x i t s h a 1 _ b a s e 6 4 = " v N / e / z 9 t V S 3 c e i U N 0 O v S h x i g 6 l Q = " > A A A C B H i c b V B L S w M x G M z W V 6 2 v W o + 9 h B b B g 5 T d 2 l J 7 E A o q e q x g H 9 A u S z Z N 2 9 B k d 0 m y Q l l 6 8 O J P 0 Y s H R T z q L / D k z X 9 j u i 3 i a y A w m f m G L x k 3 Y F Q q 0 / w w E g u L S 8 s r y d X U 2 v r G 5 l Z 6 O 9 O U f i g w a W C f + a L t I k k Y 9 U h D U c V I O x A E c Z e R l j s 6 n v q t K y I k 9 b 1 L N Q 6 I z d H A o 3 2 K k d K S k 8 6 e O l F X c O j q z O T I K p a 7 + / H 9 j D Q n T j p v F s w Y 8 C + x 5 i R f S 7 6 9 Z E 5 u c 3 U n / d 7 t + T j k x F O Y I S k 7 l h k o O 0 J C U c z I J N U N J Q k Q H q E B 6 W j q I U 6 k H c W f m M B d r f R g 3 x f 6 e A r G 6 v d E h L i U Y + 7 q S Y 7 U U P 7 2 p u J / X i d U / U M 7 o l 4 Q K u L h 2 a J + y K D y 4 b Q R 2 K O C Y M X G m i A s q H 4 r x E M k E F a 6 t 1 R c Q j U G n J F K a U 6 q 1 l c J z W L B O i i U L 3 Q b J T B D E m R B D u w B C 1 R A D Z y D O m g A D K 7 B H X g A j 8 a N c W 8 8 G c + z 0 Y Q x z + y A H z B e P w E V Y p p r < / l a t e x i t > K ± < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 M g y y 2 d O E n Q l K h I 1 7 3 R w F R F F 4 C c = " > A A A B 7 n i c b Z D L S g M x F I Y z 9 V b r r e r S T b A I r o a M v U y 7 s u B G c F P B X q A d S y Z N a 2 j m Q p I R y t C H c C N Y E b e + g q / h z r c x M 9 N F F X 8 I f P z n H M 7 J 7 4 a c S Y X Q t 5 F b W 9 / Y 3 M p v F 3 Z 2 9 / Y P i o d H H R l E g t A 2 C X g g e i 6 W l D O f t h V T n P Z C Q b H n c t p 1 p 1 d J v f t I h W S B f 6 d m I X U 8 P P H Z m B G s t N W 9 u Y 8 H o T c f F k v I R N V 6 z b Y h M i v l R s 2 u a 0 C o W m 5 Y 0 N K Q q H T 5 + Z J o 0 R o W v w a j g E Q e 9 R X h W M q + h U L l x F g o R j i d F w a R p C E m U z y h f Y 0 + 9 q h 0 4 v T c O T z T z g i O A 6 G f r 2 D q r k 7 E 2 J N y 5 r m 6 0 8 P q Q f 6 t J e Z / t X 6 k x n U n Z n 4 Y K e q T b N E 4 4 l A F M P k 7 H D F B i e I z D Z g I p m + F 5 A E L T J R O q J C G 0 E g F M 7 A r S 1 g J o X N h W m W z e o t K z Q r I l A c n 4 B S c A w v Y o A m u Q Q u 0 A Q F T 8 A Q W 4 N U I j W f j z X j P W n P G c u Y Y / J L x 8 Q N E x 5 U G < / l a t e x i t > K 0 L < l a t e x i t s h a 1 _ b a s e 6 4 = " + B H u I T Q r G v T A L C M m r r 6 E O t n x F 8 w = " > A A A B 7 H i c b V D L S s N A F J 3 4 a q 2 v + t i 5 G S y C G 0 u i l T a 7 g h t B F x V M W 2 x j m U w n 7 d D J J M x M l B L 6 D W 5 c K K J L P 8 i d f + M 0 K V L F A x c O 5 The kinetic energy distribution of light mesons when they decay, where π ± and K ± indicates the sum of charged pions and kaons.We consider two beam energies, E beam = 125, 500 GeV at ILC-250 and ILC-1000, respectively.

Light Mesons
The light mesons are mainly produced by the interaction of real photons in the electromagnetic shower with the nucleons in the beam dump.If the decay length of the produced light mesons is in the same order of magnitude of or greater than the mean free path in the material, the particles reduce their energy or change into different flavors by (in-)elastic scattering or multiple scattering.We use the following codes and models which are available in PHITS to simulate the electromagnetic shower and the production and transport of the light mesons.For the electromagnetic shower, the simulation is performed by EGS5 [17].For the light meson production and transport, the JAM [18] and JQMD [19] models modified for photoproduction (photonuclear interaction) are used.In addition to these models, INCL4.6 [20] is also employed to calculate the interaction of the mesons with nuclei during transport.The energy loss of the charged particles due to multiple scattering is evaluated by ATIMA [21].Fig. 2 shows the kinetic energy distribution of light mesons when they decay.The decay energy distribution is more important than the production energy distribution because the kinematics of new particles is determined by the parent particle distribution at the decay.For reference, the production energy distribution of light mesons is shown in Fig. 3.

Heavy Mesons
We use PYTHIA8 to calculate the differential cross sections of the γp(n) → B(D) + X process for the heavy mesons. 3We have checked that the sum of the direct and non-diffractve cross section of the D meson production agrees very well with the photoproduction data [22,23], see Appendix B. Therefore, we regard the sum of the two cross sections as the total cross section, σ total (γp(n)) = σ non-diff (γp(n)) +σ direct (γp(n)).The total cross section of atomic neucleus is obtained by taking < l a t e x i t s h a 1 _ b a s e 6 4 = " / j p 8 R s 6 j g l q t x t 0 g p 6 9 5 6 s p 7 T 1 g V r P r M H f s F 6 / Q L s / Z S A < / l a t e x i t > K < l a t e x i t s h a 1 _ b a s e 6 4 = " e P Q B S E G 9 X g A N h a K r U P z U x A P C t w I x y y x B q S 8 E U 7 N I K 1 k J o X 5 l W 0 S w 3 U a F e A q m y 4 A y c g 0 t g A R v U w S 1 o g B Y g g I I n 8 A J e j Q f j 2 X g z F m l r x l j N n I J f M j 5 + A H 5 l k j s = < / l a t e x i t > D < l a t e x i t s h a 1 _ b a s e 6 4 = " + d A H 0 w A K l p u h i S e 9 w a 4 e 0 9 a M s Z w 5 A L 9 g v H w B L e m S g g = = < / l a t e x i t > E beam = 125 GeV < l a t e x i t s h a 1 _ b a s e 6 4 = " v N / e / z 9 t V S 3 c e i U N 0 O v S h x i g 6 l Q < l a t e x i t s h a 1 _ b a s e 6 4 = " / j p 8 R s 6 j g l q t x t 0 e w B N 4 N u 6 N R + P F e J 2 3 L h j p z B 7 4 I e P t E + d S l i M = < / l a t e x i t > ⇡ ± < l a t e x i t s h a 1 _ b a s e 6 4 = " b q 7 C j u F k t p z J d X q P L 3 g p 6 9 5 6 s p 7 T 1 g V r P r M H f s F 6 / Q L s / Z S A < / l a t e x i t > K < l a t e x i t s h a 1 _ b a s e 6 4 = " e P Q B S E G 9 X g A N h a K r U P z U x A P C t w I x y y x B q S 8 E U 7 N I K 1 k J o X 5 l W 0 S w 3 U a F e A q m y 4 A y c g 0 t g A R v U w S 1 o g B Y g g I I n 8 A J e j Q f j 2 X g z F m l r x l j N n I J f M j 5 + A H 5 l k j s = < / l a t e x i t > D < l a t e x i t s h a 1 _ b a s e 6 4 = " + d A H 0 w A K l p u h i S e 9 w a 4 7 1 5 6 0 5 6 w 1 . The production rate per electron injection in the beam dump for mesons and τ lepton with respect to the kinetic energy at production, where the energy is normalized by the beam energy.The results for , and τ (τ + , τ − ) produced in the beam dump are shown, which represent the sum of the particles in the parenthesis.We consider two beam energies, E beam = 125, 500 GeV at ILC-250 and ILC-1000, respectively.
into account the shadowing effect [24,25], where l = 0.92 [25].The differential production cross sections for heavy meson photoproduction are obtained by PYTHIA8 [16] and implemented in PHITS.Since the decay lengths of the heavy mesons are much shorter than the mean free path in the material, the spectra at their production and decay are similar.In Fig. 3, we show the production rate per electron injection in the beam dump for mesons and τ lepton with respect to the kinetic energy at production, where the energy is normalized by the beam energy.The results for π , and τ (τ + , τ − ) produced in the beam dump are shown, which represent the sum of the particles in the parenthesis.We consider two beam energies, E beam = 125, 500 GeV at ILC-250 and ILC-1000, respectively.The overall yield of the heavy meson production increases as the beam energy gets higher, and B c becomes accessible at ILC-1000.For the sake of comparison, we also include π, K distributions at production.

τ lepton
As discussed in the previous literature [26], the primary source of τ lepton is the D s decay with approximately 5% branching ratio.PHITS simulates D s meson propagation and decay, which accounts for the τ lepton production.The sub-dominant source of τ lepton is the τ lepton pair production, γ +nucleus/nucleon → τ + +τ − +X.We implement a complete differential cross section calculated with the Born approximation in QED in the PHITS code and generate events for the process [27,28].The form factors for coherent (nucleus elastic), quasi-elastic and inelastic interactions are included.The spectrum of τ lepton is shown in Fig. 3.
We find that the number of τ leptons from the pair production is about 20 times smaller than those from the decay of D s .However, the pair production process becomes dominant in the high-energy region where the kinetic energy of τ lepton is above 65% of the beam energy.So the pair production process will be necessary when considering the physics of high-energy τ leptons or τ neutrinos at the ILC beam dump.

Heavy Neutral Leptons
If gauge singlet fermions N exist in the BSM sector, a renormalizable interaction with the SM sector is possible where L is the SM lepton doublet of flavor = e, µ, τ , and I is the index of N .If M I λ iI v, the standard seesaw mechanism make the SM neutrino mass light, m ν, ∼ I (λ I λ I )v 2 /M I .For O(1) Yukawa coupling, the singlet fermion N has to be extremely heavy to satisfy the bounds on the neutrino mass, m ν 0.05 eV [29].On the other hand, if the size of Yukawa coupling λ is small, M can be in MeV-GeV mass scale to satisfy the same condition.Such particles can be searched directly in laboratories, and they are often called heavy neutral leptons (HNLs) or sterile neutrinos.The active neutrinos ν from L doublet and the HNLs N are almost in the mass eigenstates up to a small admixture between ν i and N I characterized by the mixing angle Since the mixing and mass determine the HNL interactions with the SM particles, we use the mixing parameter U I instead of the Yukawa couplings to discuss the HNL phenomenology.The HNLs in the GeV mass range have another exciting aspect other than explaining the mass of active neutrinos and testability.They can be responsible for the baryon asymmetry by leptogenesis via HNL oscillation [11,12].The effective leptogenesis occurs for fast N I oscillation, and therefore, two degenerate HNLs are an excellent benchmark model to investigate.
In the minimalistic scenario with two quasi-degenerate HNLs, the target parameter space is welldefined by the baryon asymmetry and the neutrino mass.We schematically show it in Fig. 4. Note that the vertical axis is the sum of the mixing angles over both the active neutrinos and the HNLs, We also include a bound on U 2 from Big Bang Nucleosynthesis (BBN).The small mixing angle is disfavored as the HNLs are sufficiently long-lived to decay during or slightly before BBN and affect the ratio of the neutron and proton number densities [30].The target parameter space in a scenario of two degenerate HNLs with respect to the HNL mass and the sum of the mixing squared defined in Eq. (4.3).The region between dashed (dotted) green lines is favored as the HNL can generate the baryon asymmetry of the universe [11,12], and the lines are adopted from Fig. 4.17 of [8].The bottom shaded region cannot explain the neutrino oscillation data in Type-I seesaw mechanism with the two degenerate HNLs.Another shaded region with the BBN label is excluded by that the long-lived HNLs affect the successful big bang nucleosynthesis.The bound is obtained with respect to U e , U µ , or U τ in [30], and we take U 2 < min =e,µ,τ [|U (BBN) | 2 ] for this plot.
It is essential to probe all flavor mixings to cover the target parameter space for baryon asymmetry characterized by U 2 .The sensitivity to U eI and U µI of the current and proposed experiments is high, but the τ neutrino mixing is poorly constrained.In the ILC beam dump experiment, the relevant region of U τ I can be probed because τ leptons can be copiously produced from D s meson decay thanks to the higher beam energy.
In many phenomenological studies of HNLs, one assumes a single HNL (say N 1 ) in the low energy because having two HNLs will have little impact on the search sensitivity as long as the two HNLs are degenerate.In the following, we deal with one HNL for simplicity, and thus the index I is omitted.Furthermore we turn on only one U , a mixing with one of the active neutrinos in the flavor eigenstate (ν e , ν µ , ν τ ), at a time, which helps us to understand which underlying process matters.Under these assumptions, the phenomenology is well-described by the HNL mass and the single mixing in each benchmark model.Also, these benchmarks are commonly used in the literature, which allows us to compare our results with the previous works.

Sensitivities at ILC beam dump
In this subsection, we evaluate the sensitivity of the ILC beam dump experiment to the HNLs.We consider the following two production mechanisms of HNL: (i) Productions from meson and τ lepton decays; (ii) Direct productions from electrons and muons in EM shower interacting with nucleons.
In both cases, we consider the HNLs decaying inside the decay volume as the signal.We adopt the decay widths of HNL to the SM particles based on Ref. [26], and most decay patterns of the HNL would leave multiple tracks, which is distinguishable from the background.We then require that the HNL decays inside the decay volume with two or more charged tracks, i.e., the decay modes of HNL → νν ν, and π 0 ν are excluded.
We estimate that the background events are limited.One of the possible backgrounds is hardons produced by high-energy neutrinos, such as K S,L decaying to the charged pions.This type of background was estimated with GIANT 4 in the SHiP setup [31,32].They found that the relevant hadrons are always produced at the edge of the muon shield and can be significantly reduced by requiring the charged tracks consistent with the HNL kinematics.Comparing the number of produced neutrinos between SHiP and the ILC beam dump experiment, we can infer from the number of expected background events.The SHiP experiment expects 7 × 10 17 neutrinos per 2 × 10 20 protons on target, which leads to ∼ 10 7 neutrino interactions and ∼ 10 4 events with two tracks of opposite charge.99.4% of the two-track events can be discarded by the simple topology cuts, which leads to 50 remaining events.Furthermore, the veto system of SHiP can reduce the background events by 10 −4 to the level of 0.07 events.At the ILC 250 (1000), we estimate 8 × 10 16 (3 × 10 17 ) neutrinos for the 10-year run, and, assuming the similar reduction factors as in SHiP, we expect about 5 (20) background events after the two-track requirement and the simple topology cuts.Although the veto system could further reduce the background, the reduction factor strongly depends on the detector setup.Therefore, we conservatively take a pessimistic scenario such that 5 and 20 background events are expected at ILC-250 and 1000 for 10-year statistics.The corresponding upper bounds on the signal events are respectively 5.5 and 9.1.
As another source of background, cosmic muons may produce beam-unrelated events.However, the deep underground location of the detector (100 m from the ground surface), direction of tracks, and coincidence time window based on the pulsed electron beam significantly reduce the backgrounds.The detector setup and event selections to suppress the background further are beyond the scope of this paper.
In the HNL production process (i), evaluating the position and 4-momentum of the HNL is not straightforward because it involves various intermediate particles with possible higher multiplicity.Therefore, it is suitable to estimate the sensitivity using Monte-Carlo simulation.We include the SM particle decays to an HNL which are calculated in Refs.[26,[33][34][35][36][37] and summarized in Appendix F and Ref. [26].
For the other production (ii), associated physical processes are tractable without simulation.Therefore, we evaluate its sensitivity by the numerical integration and provide relevant formulae for it.To help the quantitative understanding, we also provide an approximate formula of the signal rate of the process (i) that can be integrated numerically.The signal rate of the numerical integration will be compared with the one obtained by the Monte-Carlo simulation.
In the following, we describe the detail of each method.

Monte-Carlo method
We simulate particle production and transport by PHITS with the help of PYTHIA8 for electron injection as described in Sec. 3. We also modify the decay tables of mesons and τ lepton of PHITS to include decay modes to an HNL as discussed in Appendix F. It is essential to perform the Monte-Carlo simulation to track how the system evolves from the incident electron since the intermediate steps are very involved in the production (i).The result of the Monte-Carlo simulation also becomes the guiding post for the coarse-grained integration method which is described in Sec.4.1.2.However, a naive Monte-Carlo simulation of the long-lived particle would suffer from technical challenges in obtaining sufficient statistics, in particular in the following cases: (a) Small cross section of new physics process or photoproduction of mesons.
(b) Small decay branching ratio from a SM particle to the new particle.
(c) The decay length of the new particle much shorter than the shield length.
(d) The decay length of the new particle much longer than the length of the experimental setup.
In cases (a) and (b), the processes resulting in the signal process are so rare that it is difficult to obtain a sufficient number of events.The issue can be solved in PHITS-based simulation using the biasing technique.In this technique, PHITS provides a biasing parameter for photoproduction.In this technique, the biased process occurs more frequently according to the biasing parameter, and an appropriate weight of the produced particle lower than that of the incident particle is assigned.We obtain the correct expected value of physical quantities by adding up the weights rather than adding up the number of particles produced.
In case (c), the new particles generated by the simulation predominantly decay inside the shield, so sufficient signal statistics are not easily obtained.This issue can be avoided by using the importance sampling technique which PHITS supports.This technique allows us to assign an importance to different regions of the simulation geometry.When a particle passes through the boundary between different regions with increasing importance, several copies of the particle are created in an event, and their weights are reduced depending on the importance ratio between the regions.We divide the region of the shield in the direction of the beam axis and increase their importance value exponentially so that a short lifetime particle passes through the shield more efficiently without exponential loss.
The importance technique is also useful in the opposite case (d) as we can set a large importance value in the decay volume to increase the event sampling in the relevant region.In addition, forceddecay technique is more useful when the decay length of the long-lived particle X is extremely longer than the length of the shield and decay volume.In this technique, we introduce a maximum decay length, l max .When the decay length of X, l X , is sufficiently longer than the typical length of experiment, l exp (∼ l dump + l sh + l dec ), the differential decay rate of X with respect to the flight distance z is This is a very small value, so it's hard to get enough Monte-Carlo statistics for the decay of X in the decay volume.To deal with this problem, we multiply the decay probability by b = l X /l max ( > 1) and the weight of X by 1/b.The rescaled probability is This corresponds to the decay rate when the lifetime of X is multiplied by 1/b.In summary, when l max < l X , multiplying the lifetime and weight of X by 1/b forces a higher probability of decay within the decay volume and reduces the computational cost by a factor b.This approach is an approximation, and the higher order correction to weight is O(l 2 X /l 2 max ).We typically use l max ∼ 10 3 m, so that we can safely ignore the error.We count the number of the HNL decays to the visible SM particles in the decay region taking into account all the weight factors appropriately.

Coarse-grained integration method
In this subsection, we describe the approximate estimation for the HNL productions (i) and (ii).Here, we apply suitable approximations to obtain simplified expressions involving multiple integrals for the signal yield, and then perform the integrals numerically.We refer this calculation scheme to as the coarse-grained integration (CGI) method.The same scheme was used in Ref. [4,5] in which the ALPs and dark photon productions by the EM shower were studied.
For the production mechanism (i), we consider processes where a photon from the EM shower interacts with a nucleus in the beam dump and produces the SM particle k.The number of signal events is schematically given by where l γ is the photon track length in the beam dump, n nucleus is the number density of beam dump nucleus and the production cross section of particle k is denoted by σ γ nucleus→kX .The branching ratio Br(k → HNL X) and Br vis = Br(HNL → visible SM), where the visible SM denotes the decay modes except for HNL → νν ν, and π 0 ν, are summarized in Ref. [26] and also Appendix F. Acc (i,ii) (HNL) denotes the detector acceptance depending on the production process.The distribution of photon track length l γ in the beam dump can be obtained by the Monte-Carlo simulation, and the result is given in Ref. [5] and also Appendix A.
In the production mechanism (ii), an incoming lepton = e, µ from the EM shower interacts with a nucleon and would produce an HNL.The number of signal events is schematically expressed as where the track length l ± of the lepton in the beam dump is provided in Ref. [5] and Appendix A, and n dump/shield N is the number density of nucleon in the beam dump or muon shield.For the incoming muon, we ignore the HNL production in the beam dump.The production cross section of HNL is denoted by σ ± nucleon→HNL X , which is summarized in Appendix C.
We elaborate the schematic formulae Eqs.(4.6) and (4.7) in the following.Given that the differential distribution of mesons and τ lepton were obtained by PHITS and PYTHIA8, as the results are shown in Sec. 3, we fit the differential distribution of the number of SM particle N k as the function of E k .This quantity is expressed by the more fundamental quantity as follows, This corresponds to the schematical terms in the parentheses of Eq. (4.6).Note that dl γ /dE γ and σ γ nucleus→kX is absorbed in dN k /dE k .By using dN k /dE k the production rate approximately expressed as follows, where we assume momentum of the SM particle k is aligned to the electron beam direction.The HNL production via the SM particle decay is the leading production channel at the ILC beam dump experiment for the most accessible parameter space of mass and mixing angle.As we will see later, this production mode even contributes to the sensitivity of U τ .Additionally, the direct production initiated by the charged leptons of the EM shower takes an essential role to explore the higher mass region of the HNL in the presence of U e or U µ .The direct production from the electron and positron occurs inside the beam dump, while the one from the muon occurs dominantly in the muon shield.Several direct production processes contribute to the HNL productions [38], such as quasi-elastic scattering, resonance production, and deep inelastic scattering (DIS).Among those, DIS would be a significant contribution to the HNL productions around E e ± ,µ ± ≥ 10 GeV which is favored by the acceptance of the HNL production angle.Therefore, we focus on the DIS process, and the formulae of the number of signal events are expressed as follows, The energy spectrum of the muon yield per incident electron behind the beam dump dY µ0 /dE µ0 where E µ ± 0 ≥ E µ ± is evaluated in Ref. [28], also see Eq.(A.3) of Ref. [4] and Appendix A, and the track length of muons l µ ± in the muon shield is provided in Ref. [4] ) where, as in the Monte-Carlo simulation, all the HNL decays inside the decay volume are assumed to be the visible signal except for the decay modes of HNL → νν ν, and π 0 ν.Here z is the decay position of HNL in the decay volume, and P dec denotes the decay probability of HNL given by with p HNL being the momentum, m HNL the mass, and Γ HNL the total decay width of HNL evaluated in [26].The angular distribution of the HNLs is denoted as dP ang /dθ dec HNL , where θ dec HNL is the decay angle of HNL from the particle k with respect to the direction of k, and we use an approximated formula as follows, where m k is the mass of the SM particle k, and p k denotes the momentum of k in the lab-frame.Then, the energy of the HNL is approximately given by The distance that a muon passes through the muon shield before emitting the HNL is the function of E µ0 and E µ [4], with the energy independent stopping power: dE/dx Lead = 0.02 GeV/cm.For E beam = 500 GeV, the additional active shield behind the muon shield may be necessary due to their strong penetrating power, and the muon's trajectory becomes non-trivial after the active shield.Then we account for only the HNLs produced before the active shield, which leads to the Heaviside step function in Eq. (4.14).
Regarding the impact of varying the active shield length, we study how different depth of the muon shield affects the sensitivity for the HNL at ILC-1000 in Appendix E. The angular acceptance is expressed by the Heaviside step function associated with the transverse radius r ⊥ (z) in Eqs.(4.12), (4.13), and (4.14).The typical deviation of the visible SM particles emitted from HNL from the beam axis is estimated as where θ shower e (θ shower µ ) is the angle of shower electrons and positrons (muons) with respect to the beam axis, and θ prod HNL is the production angle of HNL.For the DIS process, the production angle is expressed as where p e,µ is the momentum of incoming electron and muon, and

��
o w N 0 g q q o h h g K 0 T 1 6 R E 8 W W A / W s / U y a p 2 w x j P r 6 B e s 1 y 8 D R p X J < / l a t e x i t > N A 6 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " V q h n r S P Y < l a t e x i t s h a 1 _ b a s e 6 4 = " j l 2 c u h Z w j b B F j < l a t e x i t s h a 1 _ b a s e 6 4 = " q U y x 8 u a j S 9 p i Z x / 9 k P b 6 C e 7 u l w 0 = < / l a t e x i t > ILC-1000 The darker grey region is from the laboratory bounds, and the lighter gray region is the BBN bound for the HNL, which is roughly τ HNL > 0.02 s [30].Sensitivity reach through 10 9 Z-decays at ILC is shown as a blue solid line.See Sec.4.2.For a comparison, dashed lines show a sensitivity reach of the DUNE experiment [40] (brown), the FASER2 experiment [41] (purple), the NA62 experiment [1] (orange), and the SHiP experiment [8] (magenta), the MATHUSLA experiment [42] (green), and 10 12 Z-decays that could be realized at the FCC-ee experiment (cyan).
Although for the production mechanism (i), the angle of shower photon and the production angle of the particle k are precisely calculated in the Monte-Carlo simulation, we omit them in the coarsegrained integration method.This is because the heavy mesons are produced by high-energy photons and the mean value of the angle of shower photons θ shower γ = 8 mrad • GeV/E γ [4,5] is suppressed.

Projected sensitivities
We give the prospects of the ILC beam dump experiment in ILC-1000 and ILC-250 for 10 year run.We assume the parameter region with more than 5.5 (9.1) events can be probed at ILC-250 (1000), corresponding to the expected 95% C.L. exclusion sensitivity.After independently evaluating the signal events by the productions (i) and (ii), we combine the plots of the sensitivity of ILC.
y u y H q 0 X 6 3 X Z m r O y m U P 4 I e v t E 7 3 6 j t 8 = < / l a t e x i t >

ILC-1000
< l a t e x i t s h a 1 _ b a s e 6 4 = " 4 j J x r c u w

HNL Kinetic energy at decay [GeV]
< l a t e x i t s h a 1 _ b a s e 6 4 = " x E e 6 g F q l q 6 i J O Z D 2 M 7 4 r g l l Y a s O k J X a 6 C s f p z I k R c y i 5 3 d C d H q i 1 / e w P x P 6 8 a q G a u H l L X D U e dominance Fig. 5 shows the sensitivity of ILC for the HNL whose mixing to active neutrinos is dominated by the U e mixing.The region above the red (black) curves is the region expected 95% C.L. exclusion sensitivity that can be searched by ILC-1000 (ILC-250) with 10-year statistics.The gray-shaded regions are constrained from past experiments which is discussed in Sec.4.3.Also, for a single benchmark point, we show the kinetic energy distribution of the decayed HNLs inside the decay volume in Fig. 6.
Let us explore the results in Fig. 5, highlighting each of the HNL production processes.The limit obtained by the Monte-Carlo simulation can be checked using the approximate formula Eq. (4.9).As explained in Appendix D, the approximate formula agrees well with the results of the Monte-Carlo simulation.For the electron beam energy E beam = 500 GeV, the number of events by the leptonic D meson decays D ± → e ± + HNL is given by Above m HNL ∼ 3 GeV, the direct production channel with the incoming electrons and positrons become significant.In particular, the production with the high-energy primary electron is essential to extend the mass reach of the U e dominant scenario.In the small |U e | 2 regions, the decay probability in Eq. (4.15) is approximately dP dump dec /dz 1/l (lab) X , and the energy of the produced HNL is approximately E lab HNL E e ± .Then, the number of events for E beam = 500 GeV is given by U µ dominance Fig. 7 summarizes the prospects of the ILC beam dump experiment in ILC-1000 and ILC-250 for the HNL mixing dominantly with ν µ .The red (black) curves show the expected sensitivity for 95% C.L. exclusion of ILC-1000 (ILC-250) with 10-year statistics given by the meson and τ lepton decays, and the DIS process for incoming muons from the EM shower.We provide approximated formulae of the number of signal events for each HNL production mode for ease of understanding.
o m q B H 9 I x e j A f j y X g 1 3 u a j K W O x k 0 O / Z L x / A e W a l J 8 = < / l a t e x i t > ILC Giga-Z < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 T + o u g a 3 a E H 9 K j 1 t X v t S X s e j 0 5 p k 5 0 N 9 E v a 6 x f w y p Z D < / l a t e x i t > N A 6 2 o w N 0 g q q o h h g K 0 T 1 6 R E 8 W W A / W s / U y a p 2 w x j P r 6 B e s 1 y 8 D R p X J < / l a t e x i t > m HNL [GeV] < l a t e x i t s h a 1 _ b a s e 6 4 = " j l 2 c u h Z w j b B F j u a j S 9 p i Z x / 9 k P b 6 C e 7 u l w 0 = < / l a t e x i t > Figure 7. Sensitivity reach of ILC beam dump experiment to HNLs mixing with the mu-neutrino in the mass and mixing plane.For the description of this plot, see the caption of Fig. 5.
u a j S 9 p i Z x / 9 k P b 6 C e 7 u l w 0 = < / l a t e x i t > ILC-1000 o u g a 3 a E H 9 K j 1 t X v t S X s e j 0 5 p k 5 0 N 9 E v a 6 x f w y p Z D < / l a t e x i t > q q J j V E M U h e g e P a I n 4 9 Z 4 M J 6 N l + H o h D H a W U O / Z L x 9 A 8 + f l z 0 = < / l a t e x i t >

(a) Meson and τ lepton decays
The production process from meson and τ lepton decays of the HNL that mix dominantly with ν µ can be calculated in parallel with U e dominant case except for minor threshold difference.The mass dependence of the dominant HNL production process is the same as that of U e dominant case with e replaced with µ.
As shown in Figs. 17  As shown in Fig. 1, incident real photons produce muon pairs in the beam dump through the electromagnetic interaction with the nucleus or nucleon [28], and the HNL can be generated by the DIS process in the muon shield.The muon pair production cross section is about (m µ /m e ) 2 10 5 times smaller than that of the electron pair production.
In the small |U µ | 2 regions, the projected sensitivity scales with the ratio of numbers of muon and electron in the shower, and the projected sensitivity via the DIS process is less significant than the U e dominant case.On the other hand in the large |U µ | 2 regions, the number of signal events is mostly determined by the muon shield length and less sensitive to the number of muon pairs because the lifetime is shorter.The shape of the upper side of contours in Fig. 7  HNL .Since the DIS process with a muon happens in the muon shield, the shorter distance from the HNL production point to the decay volume enhances the acceptance compared the DIS process in the beam dump.

U τ dominance
The projected sensitivities in the U τ dominant scenario are shown in Fig. 8 with a similar notation as in the U e,µ dominance.The red (black) curves show the expected 95% C.L. exclusion sensitivity with 10-year statistics in ILC-1000 (ILC-250) by B and D s mesons and τ lepton decay production.The direct (DIS) production is absent in this case due to very limited flux of τ leptons.
As shown in Figs. 3, 19, 20, and 23, below the threshold of m Ds −m τ 0.2 GeV, the main HNL production channel is the D s → τ + HNL decay.For 0.2 GeV m HNL ≤ m τ − m e 1.8 GeV, the production of HNL is dominated by τ lepton decay.Most of the τ leptons are produced by D s meson decays, and the τ pair production in the electromagnetic showers are subdominant, which is therefore dropped in the HNL study.The branching ration of D s → τ + ν τ is 5.48%, and an energy dependence of τ lepton spectra is determined by that of D s meson spectra in Fig. 3.The main decay channels of τ are τ → ρ + HNL, τ → µ + ν + HNL, τ → e + ν + HNL, τ → π + HNL, τ → K * + HNL, and τ → K + HNL, see Fig 23.In the small |U τ | regions with the longer lifetime of HNL, the decay probability in Eq. (4.15) is approximated by dP dump dec /dz 1/l (lab) X .For the electron beam energy E beam = 500 GeV, the number of events by the τ ± → ρ ± + HNL is approximately given by

Sensitivities from Z decays at ILC and FCC-ee
The future e + e − colliders can be seen as a Z boson factory, such as Giga-Z program of the ILC and Tera-Z program of the CERN Future e + e − Circular Collider, dubbed FCC-ee.The HNLs can leave a clear signal with displaced tracks at e + e − colliders once they are produced via Z → N ν, N ν.This type of signal was examined at DELPHI detector of the Large Electron-Positron collider (LEP) [43].We briefly study the future sensitivities of the HNL search using the displaced tracks from the Z decay, which is complementary to the sensitivities of the ILC beam dump experiment.The proposed ILC detector [44] has three layers of vertex detector (VTX) starting from 1.5 cm to 19.5 cm surrounding the beam, Time Projection Chambers (TPC) as the main tracking detector spreading from 33 cm to 170 cm, and two layers of silicon strip detectors (SIT) arranged between TPC and VTX.We assume the ILC detector is sensitive to the HNL signal with negligible background if the HNL decays between 3 cm and 170 cm from the collision point.We adopt the radius-dependent track-detection efficiency trk linearly decreasing from r = 3 cm ( trk = 100%) to r = 170 cm ( trk = 0%) [45].Then, we assume 10 9 Z decays and zero background after requiring the displaced tracks and consider all the HNL decays except for the fully invisible mode, N → νν ν.The 95% C.L. future sensitivity corresponds to three HNL decays with the multi-displaced tracks in the fiducial volume, and the result is included in Figs. 5, 7, and 8.
Similarly, we study the future projection at the FCC-ee.For the simplicity, we take the same detector setup of the ILC and rescale the statistics assuming 10 12 Z decays.This is also shown in Figs. 5, 7, and 8. Ref. [46] studied also the projected sensitivities at the FCC-ee, but our estimate is different with respect to the detector coverage and the signal efficiency.

Existing constraints on HNL and projected sensitivities at other experiments
The presence of HNL can have drastic consequences on early-universe observables and has been explored extensively in the literature.The HNL is thermally produced in the early universe, and its late decay can disrupt the standard Big Bang Nucleosynthesis.The typical constraint is for τ HNL 1s, see [47] and references therein.In addition, Ref [30] studies that the pions from the HNL hadronic • Z boson decay -The search of Z boson decaying into HNL using the DELPHI data at the LEP collider constraints all the scenarios.It sets strongest limit in mass range around m HNL ∼ O(10) GeV [43].See Sec.4.2 for more detail and the projection at the future e + e − colliders.
• Large Hadron Collider -At the ATLAS and CMS detectors of LHC, the W boson mediated process efficiently produces the HNL.Depending on the final states and decay length, different analyses were performed at ATLAS [67,68] and CMS [69][70][71].In particular, the search for the displaced decays of HNL [68,71] excludes the large parameter space.
In Figs. 5, 7, and 8, the estimated sensitivity contours of ILC is better than the other proposed searches.The ILC sensitivity is very close to the one of SHiP in the low mass region (m HNL < 2 GeV).Above 2 GeV, the sensitivity of the ILC beam dump experiment is better than SHiP because the initial electron energy is high enough to produce B mesons at a higher rate.Thanks to the higher HNL mass, the HNL decay products are more clearly separated from the background so that the neutrino background would be less critical in the region.

Discussions
The ILC beam dump experiment is a seamless extension of the ILC program, which provides a unique opportunity to test feebly interacting light particles.The electron beam energy of ILC is much higher than the ones at the current and past electron beam dump experiments, and therefore, not only the light mesons but also the heavier SM particles, such as heavy flavor mesons and τ lepton, can be produced in the beam dump.Moreover, a large number of electrons on target are expected thanks to the high-intensity beam.In many BSM scenarios, new particles can be efficiently produced by decays of the SM particles.Therefore, it is essential to estimate the production rate of the SM particles at the ILC beam dump experiment.In this paper, we evaluate the light and heavy mesons and τ lepton spectra at the decay for the first time using PHITS and PYTHIA8.PHITS is responsible for producing and transporting light SM particles, and we incorporate the heavy meson productions calculated by PYTHIA8 into PHITS.The main results are in Fig. 2 for the light mesons and in Fig. 3 for the heavy mesons and τ lepton.
The spectra of SM particles can be used to estimate the yield of BSM particles from SM particle decay.As a demonstration, we studied the projected sensitivity of the heavy neutral leptons at the ILC beam dump experiment.Figs. 5, 7, and 8 show that ILC would explore the HNLs with the heavier mass and small mixing and cover the large parameter space motivated by the baryon asymmetry of the Universe.For the reader's convenience, we show the sensitivities after one year of operation in Appendix G.We use the Monte-Carlo simulation to evaluate the HNL signal from the SM particle The dotted and long dashed lines correspond to non-diffractive and direct production cross section of charm quark.The solid line is the total charm meson production cross section.They are consistent with the measured cross section of Tagged Photon Spectrometer (circle) [23] and SLAC (triangle) [22].
photoproduction cross section is measured by the Tagged Photon Spectrometer Collaboration [23] and the SLAC Hybrid Facility Photon Collaboration [22].
In Fig 9, we show the measured cross section(the dots and the triangle ) and the non-diffractive and direct production cross sections of PYTHIA8 as a function of √ s.The sum of the two cross sections (the solid line) is consistent with the measured cross sections. 5

C Production cross sections of the HNLs by deep inelastic scattering
We consider the production of the HNLs through deep inelastic scatterings (DIS) in addition to the SM particle decays.For convenience, we summarize the production cross sections used in our numerical evaluations; using Eqs.(4.10) and (4.11), one can estimate the number of signal events.The cross sections for the DIS on the unpolarized nucleon are given by five structure func-tions [29,38,[81][82][83][84][85]] where |U | is the mixing angle, E ∓ is the initial lepton energy in the nucleon rest frame, M N is the mass of nucleus, M W is the W mass, G F is the Fermi coupling constant, and m is the incoming charged lepton mass.The squared four-momentum transfer is given by Q where k and k HNL are the four-momenta of the incoming and outgoing leptons, respectively.In this work, for simplicity, we impose the Callan-Gross relations [86], and the Albright-Jarlskog relations [81] as follows: Here, the nucleon structure functions in the quark-parton model are given by [29,83] . The left (right) panel shows the curves corresponding to three signal events with 10-year statistics for the U e (U τ ) dominant scenario at ILC-1000.For simplicity, we assume the signal events are the HNL decays inside the decay volume except for HNL → νν ν.The solid (dashed) lines correspond to the limit based on the Monte-Carlo simulation (the coarse-grained integration method).Note that the curves of the coarse-grained integration method almost reproduce that of the Monte-Carlo simulation.
where x denotes the Bjorken scaling defined as x = Q 2 /2M N (E ∓ − E HNL ) with the lepton energy loss in the rest frame, and y = (E ∓ − E HNL )/E ∓ is the fraction of the lepton energy loss in the rest frame.The squared four-momentum transfer can be expressed as the x and y kinematic limits are obtained.In the numerical analysis, we adopt the MSTW parton distribution functions [87].

D Comparison between Monte-Carlo method and coarse-grained integration method
In our analyses in Sec. 4, we evaluated the number of signal events by the Monte-Carlo simulation based on the PHITS framework.To explain how to use the meson and τ lepton spectra, and perform easy consistency checks, we also evaluate the number of signal events by the coarse-grained integration method using the meson spectra in Sec. 3. The formulae of coarse-grained integration for the number of signal events are summarized in Sec. 2. In this section, we assume zero background and require three signal events that the HNL decays inside the decay volume as signal except for the invisible decay mode, HNL → νν ν.
In Fig. 10, we show the sensitivity of ILC-1000 with 10-year statistics for the U e and U τ dominant scenarios.The red solid (dotted) curves show the sensitivities based on the Monte-Carlo simulation (coarse-grained integration method).For the U e dominance, the D and B meson decays dominate the HNL productions in 0.5 GeV m HNL 2 GeV and 2 GeV m HNL 3 GeV, respectively.In 3 GeV m HNL , the DIS process becomes significant.For the U τ dominance, D and B mesons dominate the HNL productions in m HNL 0.2 GeV and 0.2 GeV m HNL , respectively.In a wide range of the parameter space, the coarse-grained integration method almost reproduces the results of the Monte-Carlo simulation.Note here that the coarse-grained integration method does not overestimate the number of signal events.

E DIS process and muon shield configuration at ILC-1000
For ILC-1000, high energy muons are produced in the beam dump and penetrate the muon lead shield.One can reduce background events by placing an active muon shield behind the muon lead shield as shown in Fig. 11.In our analyses in Sec. 4, we assumed that the HNLs produced in the muon lead shield with the length of l lead sh = 10 m contribute to the signal events.Since, however, the length of the muon lead shield highly depend on the details of the detector concepts, we study the sensitivity of The colored dotted curves correspond to three signal events with 10-year statistics for the DIS process.We assume the signal events are the HNL decays inside the decay volume except for HNL → νν ν We assume the signal events are the HNL decays inside the decay volume except for HNL → νν ν.The results for the different lengths of the muon lead shield, 10, 20, 40, 60 m are compared.The red curve shows the expected 95% C.L. exclusion sensitivity with 10-year statistics for the HNL productions from the SM particle decays.The blue curve corresponds to the sensitivity of the Giga-Z program at ILC. ILC-1000 depending on l lead sh .In this Appendix, we assume zero background and account for all the HNL decays inside the decay volume except for HNL → νν ν as signals.
In Fig. 12, the sensitivities of the DIS process for incoming muons are shown for l lead sh = 10, 20, 40, 60 m.The constraint to the large |U µ | 2 region enlarges by shortening active muon shield, because more HNLs decay in the decay volume.

F Branching ratios of meson and τ lepton decays to an HNL
For convenience, we summarize the branching ratios of the meson and τ lepton decays to produce an HNL.The details of the branching ratio calculation are found in Ref. [26].The branching ratios for the light unflavored and strange mesons are shown in Fig. 13 and 16.The relevant mesons for the HNLs productions can be: π + (u d, 139.6),K + (us, 494), K 0 S (ds, 498), and K 0 L (ds, 498).B s (bs, 5367), B c (bc, 6276).The branching ratios for the beauty mesons are shown in Fig. 15, 18, and 20.
In addition to the mesons, τ lepton can be produced by D s meson decays and τ pair productions in the electromagnetic showers.For 0.5 GeV m HNL 2 GeV, the HNL production from the τ decays is significant effect.The branching ratios for the τ lepton are shown in Fig. 21, 22, and 23.

G Projected sensitivities for one year operation
We assumed 10-year operation in the main results.Here, for the reader's convenience, we add the projected sensitivities after one year of operation.The estimated background is expected to be a factor of 10 smaller than in the 10-year case.Therefore, background of 0.5 and 2 events are assumed at ILC-250 and ILC-1000, respectively.The 95%CL sensitivities are shown as dashed lines, while the solid lines are for the 10-year operation for comparison.
Figure1.A setup for ILC beam dump experiments.It consists of the main beam dump, a muon shield, and a decay volume.We assume a multi-layer tracker is placed in the decay volume so that the charged tracks are measured.
9 z L v f d 4 E a N S m e a X s b C 4 t L y S y 6 8 W 1 t Y 3 N r e K 2 z t N G c Y C E w e H L B R t D 0 n C K C e O o o q Rd i Q I C j x G W t 7 o f O q 3 7 o m Q N O Q 3 a h w R N 0 A D T n 2 K k d K S c 3 l n 9 q 5 6 x Z J Z N l P A O W L b N a t m Q 2 u m l O p 7 x 4 O H 3 P t t o 1 f 8 7 P Z D H A e E K 8 y Q l B 3 L j J S b I K E o Z m R S 6 M a S R A i P 0 I B 0 N O U o I N J N 0 m M n 8 F A r f e i H Q h d X M F X n J x I U S D k O P N 0 Z I D W U f 7 2 p + J / X i Z V f c x P K o 1 g R j r N F f s y g C u H 0 c 9 i n g m D F x p o g L K i + F e I h E g g r n U 8 h D c F O A T N S r c y I b f 2 E 0 D w p W 6 f l s 2 u d R g V k y I N 9 c A C O g A W q o A 4 u Q A M 4 A A M K H s E z e D G4 8 W S 8 G m 9 Z 6 4 I x m 9 k F v 2 B 8 f A O 6 S 5 E / < / l a t e x i t > K 0 S < l a t e x i t s h a 1 _ b a s e 6 4 = " l d g S q 7 p i Z O k n A w 5 d u M P a w r i S V 4 0 = " > A A A B 7 H i c b V D L S s N A F J 3 U R 2 t 9 1 c f O z W A R 3 F i S 2 u e u 4 E Z w U 9 G 0 x T a W y X T a D p 1 M w s x E K a H f 4 M a F I r r 0 g 9 z 5 N 0 6 T I i o e u H A 4 5 1 7 u v c c N G J X K N D + N 1 N L y y m o 6 s 5 Z d 3 9 j c 2 s 7 t 7 L a k H w p M b O w z X 3 R c J A m j n N i K K k Y 6 g S D I c x l p u 5 O z u d + + I 0 J S n 1 + r a U A c D 4 0 4 H V K M l J b s i 1 u z f 9 X P 5 c 1 C p V a 2 i k V o F s w Y m l T q J a t c g 9 Z C y T f 2 T 0 b 3 6 b e b Z j / 3 0 R v 4 O P k g P 8 m j c + / 8 c H 4 5 v + e j C 0 6 x s 0 N e y P n z H 4 + T n S 4 = < / l a t e x i t > ⇡ ± < l a t e x i t s h a 1 _ b a s e 6 4 = " b q 7

2 NFigure 4 .
Figure 4.The target parameter space in a scenario of two degenerate HNLs with respect to the HNL mass and the sum of the mixing squared defined in Eq.(4.3).The region between dashed (dotted) green lines is favored as the HNL can generate the baryon asymmetry of the universe[11,12], and the lines are adopted from Fig.4.17 of[8].The bottom shaded region cannot explain the neutrino oscillation data in Type-I seesaw mechanism with the two degenerate HNLs.Another shaded region with the BBN label is excluded by that the long-lived HNLs affect the successful big bang nucleosynthesis.The bound is obtained with respect to U e , U µ , or U τ in[30], and we take U 2 < min =e,µ,τ [|U(BBN) | 2 ] for this plot.
denotes the lab-frame decay length of HNL given by l (lab) HNL = p HNL m HNL 1 Γ HNL (4.17) µ xy with the mass of nucleus.As provided in Refs.[4] and[5], θ shower e is estimated by the Monte-Carlo simulation, and we use the mean value θ shower e = 16 mrad • GeV/E e ± .For the incoming muons, we use θ shower µ = θ MCS .
o 0 b o 1 7 4 8 l 4 H o 5 m j N H O M v o l 4 + U L P 6 q Z E A = = < / l a t e x i t > SHiP < l a t e x i t s h a 1 _ b a s e 6 4 = " I k k u 6 b B F D V n e J 3 E S 0 8 B 1 j / 0 W m 0 g = " > A A A B 8 n i c b V D L T g J B E J z 1 i f h C P X q Z S D R 4 2 e y K C J w k 8 c I R o z w S 2 J D Z Y Y A J s 4 / M 9 B r J h m / w 5 M W D x n j l K / w E b 3 6 I d 5 d d Y p R Y S S e V q u 5 0 d 9 m + 4 A o M 4 1 N b s e j 0 5 p k 5 0 N 9 E v a 6 x f w y p Z D < / l a t e x i t > D U N E < l a t e x i t s h a 1 _ b a s e 6 4 = " L h Z C 0 a U l w 6 h I F V D 6 x c t w c 5 v M U c g = " > A A A B 8 n i c b V D L S g N B E J z 1 b X x F P X o Z D I p e l t k k m 8 S T g g q e g e P a I n 4 9 Z 4 M J 6 N l + H o h D H a W U O / Z L x 9 A 8 + f l z 0 = < / l a t e x i t > BBN < l a t e x i t s h a 1 _ b a s e 6 4 = " b K z h b p 7 T 8 P z 8 t Y + P l Q O S w 2 E x m n s = " > A A A B 8 X i c b Z D L S g M x F I Y z 9 V b r r e r S T b A o d V M y 9 d J 2 p e j G l V S w W m x L y a R p G 8 x k h u S M W I Y + g x s 3 L h R x q 0 / h I 7 j z Q d w 7 n S m i 4 g + B j / + c w z n 5 H

0 12 Z
< l a t e x i t s h a 1 _ b a s e 6 4 = " 9 A Z s 5 P t Y B i 2 4 n V p s Y J r J 9+ G + B 7 s = " > A A A B / H i c b Z B L T w I x F I U 7 + E J 8 o S z d N I L G h S E z A w K z I 3 H j E h N B I y D p l A I N n U f a O 0 Y y w b / i x o X G u P W H u P P f W A Z i 1 H i S J i f n 3 J v e f G4 o u A L T / D R S S 8 s r q 2 v p 9 c z G 5 t b 2 T n Z 3 r 6 W C S F L W p I E I 5 L V L F B P c Z 0 3 g 4 N G I j H v j y X j O W p e M x c w + / J D x 8 g n 1 e J N c < / l a t e x i t > ILC-250 < l a t e x i t s h a 1 _ b a s e 6 4 = " K S r k n / e U S u I D 4 E z L 4 i 1 6 g M s T m 2 o = " > A A A B + 3 i c b Z B L S 8 N A F I U n P m t 9 x b o U Z L A I b i x J a 6 l 1 Y 6 E b B R c V 7 A P a U i b T S T t 0 8 m D m R l p C / o o b F 0 p x 6 x 9 x 5 w 9 x b 5 o W 8 X V g 4 O O c e 5 n L s X z B F R j G u 7 a 0 v L K 6 t p 7 a S G 9 u b e / s 6 n u Z h

Figure 5 .
Figure 5. Sensitivity reach of ILC beam dump experiment to HNLs mixing with the electron neutrino in the mass and mixing plane, assuming 10 year run at ILC-250 (black solid) and ILC-1000 (red solid).The number signal events more than 5.5 (9.1) is required at ILC-250(1000), which corresponds to the 95% C.L. sensitivity.The discussion about the background is in Sec.4.1.The current exclusion bounds are shown in the gray region, see Sec. 4.3.The darker grey region is from the laboratory bounds, and the lighter gray region is the BBN bound for the HNL, which is roughly τ HNL > 0.02 s[30].Sensitivity reach through 10 9 Z-decays at ILC is shown as a blue solid line.See Sec.4.2.For a comparison, dashed lines show a sensitivity reach of the DUNE experiment[40] (brown), the FASER2 experiment[41] (purple), the NA62 experiment[1] (orange), and the SHiP experiment[8] (magenta), the MATHUSLA experiment[42] (green), and 10 12 Z-decays that could be realized at the FCC-ee experiment (cyan).

< l a t e x i t s h a 1 _ b a s e 6 4 =
X s C r 9 W C 9 W e / W x 6 w 1 Z c 1 n D s A C U u l P n n O z O w = = < / l a t e x i t > " 4 Y H P e V c T E c 9 h e s P N Y u u q g e e X b e n P e i d c w p Z 1 b J D z g f n 5 G t q s o = < / l a t e x i t > dP /dE [1/GeV] < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 y B J m g B D O 7 A I 3 g C z 9 q D 9 q K 9 a m / z 0 Z y W 7 e y B B W n v X 0 b F m V 4 = < / l a t e x i t > ILC-250< l a t e x i t s h a 1 _ b a s e 6 4 = " b S z o C N u K + O U 3 r 2 q B R 6 J U R T k 5 s b 4 = " > A A A B 7 n i c b Z D L S s N A F I Z P 6 q 3 W W 9 W l m 8 E i u L E k 1 V K 7 K 3 S j 4 K K C v U A b y m Q 6 a Y d O J m F m I p T Q h 3 D j Q h G3 P o 8 7 3 8 Z p G s T b D w M f / z m H c + b 3 I s 6 U t u 0 P K 7 e y u r a + k d 8 s b G 3 v 7 O 4 V 9 w 8 6 K o w l o W 0 S 8 l D 2 P K w o Z 4 K 2 N d O c 9 i J J c e B x 2 v W m z U W 9 e 0 + l Y q G 4 0 7 O I u g E e C + Y z g r W x u t c 3 z b n 1 a P r + p l O o 4 i y M P B 3 A I x 2 B D F e p w C U 1 o A Q U B D / A E z 9 a d 9 W i 9 W K + L 1 p y V z e z D D 1 l v n y Q w j x M = < / l a t e x i t > U e dominance < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 G J t / 8 L S k 5 9 K i Z 6 x h 2 5 z x J F E e N E = " > A A A B + H i c b Z B L S 8 N A F I U n P m t 9 N O r S z W A r u C p J V W p 3 B T c u K 5 i 2 0 I Y w m d y 0 Q y e T M D M R a u k v c e N C E b f + F H f + G 9 M 0 i K 8 D A x / n 3 M t c j p 9 w p r R l f R g r q 2 v r G 5 u l r f L 2 z u 5 e x d w / 6 K o 4 6 H F 0 x i p 1 D 9 E P G 2 y f j u J L z < / l a t e x i t > U µ dominance < l a t e x i t s h a 1 _ b a s e 6 4 = " f 3 0 1 A a w a L M / P h G w k 4 v p 3 3 1 S 6 U S s = " > A A A B / H i c b Z B L S 8 N A F I U n P m t 9 R b t 0 M 9 g K r k p S l d p d w Y 3 5 5 Z a z h e N u b J v A G D i F O B B C a e W b e a V i H 4 F + w S 6 q B U z z P f x 0G E E 0 6 E x g w p N b K t W L s Z k p p i R u b V c a J I j P A M T c g o R 4 E 4 U W 5 W H D + H J 7 k T w D C S + R M a F u 7 3 j Q x x p V L u 5 5 M c 6 a n 6 n S 3 M / 7 J R o s N L N 6 M i T j Q R e P l R m D C o I 7 h o A g Z U E q x Z m g P C k u a 3 Q j x F E m G d 9 1 U t S u g U g k t o n 5 f Q s b 9 K 6 L e a 9 l n z 4 q Z V 7 1 p l H R V w B I 7 B K b B B G 3 T B N e g B B 2 C Q g g f w B J 6 N e + P R e D F e l 6 M r R r l T A z 9 k v H 0 C V O q U 7 A = = < / l a t e x i t > U ⌧ dominance< l a t e x i t s h a 1 _ b a s e 6 4 = " s m a P g P g v v c n 9 r v 5 l w W h o v s 3 o x 5 0 = " > A A A B / X i c b Z B L S 8 N A F I U n 9 V X r K z 5 2 b g Z b w V V J q l K 7 K 7 h x W c G 0 h T a E y W T S D p 3 M h J m J U E P x r 7 h x o Y h b / 4 c 7 / 4 1 p G s T X g Y G P c + 5 l L s e P G V X a s j 6 M 0 t L y y u p a e b 2 y s b m 1 v W P u 7 n W V S C Q m D h Z M y L 6 P F G G U E 0 d T z U g / l g R F P i M 9 f 3 I 5 z 3 u 3 R C o q + I 2 e x s S N 0

Figure 6 .
Figure 6.The normalized kinetic energy distribution of HNLs when they decay inside the decay volume.The mixing angle squared is 10 −8 , and the mass is 1 GeV.The red (black) lines are for ILC-1000(250).
(a) Meson and τ lepton decays In the small |U e | 2 regions where the lifetime of HNL becomes longer, the decay probability in Eq. (4.15) is approximately dP dump dec /dz 1/l (lab) X .Depending on the mass of HNL, the HNLs are produced by the decay of mesons and τ lepton, see Figs. 13, 14, 15, and 21 in Appendix F. Below the kaon mass m K 0.5 GeV, the HNLs are mainly produced by the kaon decay.For m HNL ≤ m D,τ − m e ∼ 2 GeV, the D(D s ) meson decay dominates the HNL productions.In the region of 1 GeV m HNL 2 GeV there are many thresholds of production mode, such as D → K + e + HNL, D s → η + e + HNL, D → π + e + HNL and τ → ν τ + e + HNL, see Fig.

4 ,
(4.25)    where we used follwing approximations:Br(D ± → e ± + HNL) ∝ |U e | 2 , Br vis 1, (l lab X ) −1 ∝ |U e | 2m 4 HNL /E lab HNL .(4.26) Near the kinematic threshold, the number of events is rapidly decreased by a suppression of the branching ratio Br(D ± → e ± + HNL).For m D − m e m HNL m B − m e , the D meson decay channels are closed, and the B meson decay channels become dominant in the HNL productions up to around m HNL ∼ 3 GeV.(b) Direct production
r y X r O R p e s x c 4 + + C H r 5 R M f 8 J X H < / l a t e x i t > m HNL[GeV]    < l a t e x i t s h a 1 _ b a s e 6 4 = " j l 2 c u h Z w j b B F j+ 9 f W N 5 n o m H s E p w = " > A A A B / H i c b Z B L S 8 N A F I U n 9 V X r K 9 q l m 8 F W c F W S q t T u C i 7 s Q q S C r Y U 0 h M l 0 2 g 6 d m Y S Z i R B C / S t u X C j i 1 h / i z n 9 j m g b x d W D g 4 5 x 7 m c v x Q 0 a V t q w P o 7 C 0 v L K 6 V l w v b W x u b e + Y u 3 s 9 F U Q S k y 4 O W C D 7 P l K E U U G 6 m m p G + q E k i P u M 3 Pr T 8 3 l + e 0 e k o o G 4 0 X F I X I 7 G g o 4 o R j q 1 P L N c 5 V 4 y k B y 2 r y 5 n V e h c k J 7 r e m 0 F B U P B A 7 n 3 E s u n x 0 I r s A 0 P 4 2 J y a n p m d n E X H J + Y X F p O b W y e q 7 8 U F J W o 7 7 w 5 a V N F B P c Y z X g I N h l I B l x b c E u 7 O 7 B o L + 4 Z l J x 3 z u D X s C a L u l 4 3 O G U g I 5 a q e U G s B u I

Figure 8 .
Figure 8. Sensitivity reach of ILC beam dump experiment to HNLs mixing with the tau-neutrino in the mass and mixing plane.For the description of this plot, see the caption of Fig. 5.
and 18, for m D − m µ m HNL m B − m µ , the B meson decay gives a significant contribution to the HNL production.The zigzag curves near m HNL ∼ 3 GeV correspond to the threshold of B → D + µ + HNL.Above m HNL 3 GeV, the leptonic B meson decay such as B ± → µ ± + HNL dominates the HNL productions.(b) Direct production is determined by the exponential factor in Eq. (4.15), which is characterized by m HNL Γ HNL E lab HNL (l sh − δ µ ) ∼ const.(4.30) Combining E lab HNL E beam and Γ HNL ∝ |U µ | 2 • m 5 HNL , Eq. (4.30) becomes |U µ | 2 ∝ m −6

Figure 9 .
Figure 9.The photoproduction cross section of D mesons (D ± , D 0 and D 0 ), and the PYTHIA8 prediction.

Figure 11 .
Figure 11.A setup for the ILC-1000 beam dump experiment.It consists of the main beam dump, a muon lead shield, an active muon shield, a decay volume, and a detector.The HNLs can be produced by the DIS process for incoming muons in the muon lead shield with the length of l lead sh .

Figure 13 .
Figure 13.Dominant branching ratios of HNL production from light mesons.In this Figure, we take U e = 1, U µ = U τ = 0.The decay width to HNLs are normalized by experimental values of the total decay width of mesons.Below the kaon mass, the HNL productions are dominated by the leptonic decays of pions and kaons.

Figure 14 .
Figure 14.The same plot as Fig.13but for charmed mesons.