Bilinear R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document} parity violation at the ILC: neutrino physics at colliders

Supersymmetry (SUSY) with bilinearly broken R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document} parity (bRPV) offers an attractive possibility to explain the origin of neutrino masses and mixings. In such scenarios, the study of neutralino decays at colliders gives access to neutrino sector parameters. The ILC offers a very clean environment to study the neutralino properties as well as its subsequent decays, which typically involve a W\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W$$\end{document} or Z\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z$$\end{document} boson and a lepton. This study is based on ILC beam parameters according to the Technical Design Report for a center of mass energy of 500GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$500~\hbox {GeV}$$\end{document}. A full detector simulation of the International Large Detector (ILD) has been performed for all Standard Model backgrounds and for neutralino pair production within a simplified model. The bRPV parameters are fixed according to current neutrino data. In this scenario, the χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_1$$\end{document} mass can be reconstructed with an uncertainty of δmχ~10fit=(40(stat.)⊕50(syst.))MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta m^{{\mathrm {fit}}}_{\widetilde{\chi }^{0}_{1}}=(40\text {(stat.)} \oplus 50\text {(syst.)})~\hbox {MeV}$$\end{document} for an integrated luminosity of 500fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$500\,\hbox {fb}^{-1}$$\end{document} from direct χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_1$$\end{document} pair production, thus, to a large extent independently of the rest of the SUSY spectrum. The achievable precision on the atmospheric neutrino mixing angle sin2θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2 \theta _{23}$$\end{document} from measuring the neutralino branching fractions BR(χ~10→Wμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }_1^0\rightarrow W\mu $$\end{document}) and BR(χ~10→Wτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }_1^0\rightarrow W \tau $$\end{document}) at the ILC is in the same range than current uncertainties from neutrino experiments. Thus, the ILC could have the opportunity to unveil the mechanism of neutrino mass generation.


Introduction
Supersymmetry (SUSY) [1,2] is a very appealing extension of the Standard Model (SM). It provides an elegant solution for the Higgs hierarchy problem, makes gauge unification possible, and, apart from that, SUSY is the only non-trivial extension of the Lorentz algebra [3]. In the most general renormalizable Lagrangian of the Minimal Supersymmetric Standard Model (MSSM) trilinear and bilinear terms appear, a e-mail: benedikt.vormwald@desy.de which violate the conservation of baryon number B and lepton number L.
The presence of all these terms would lead to proton decay, which is experimentally not observed. A common way to circumvent this problem is to introduce a discrete Z 2 symmetry assigned to each field in order to suppress these terms. This quantum number, called R parity, has the form where B is the baryon number, L the lepton number and S the spin of the field. Hence, SM particles always carry R = +1 and SUSY particles R = −1. The conservation of this quantum number has the consequence that all B and L breaking terms in the SUSY Lagrangian are forbidden and the proton remains stable. However, proton decay only appears if B and L violation is present at the same time. So, breaking either B or L is well consistent with proton stability. Thus, R parity violating (RPV) SUSY scenarios are also viable alternatives to the widely studied R parity conserving (RPC) scenarios.
We will focus in the following on bilinear R parity violation (bRPV), which has the interesting feature to be able to introduce neutrino masses and mixings. The phenomenology of this mechanism has already been discussed in detail in the literature [4][5][6][7][8][9]. Our aim is to investigate the performance of the International Linear Collider and one of its proposed detector concepts for measuring the atmospheric neutrino mixing angle in such a bRPV SUSY scenario, based on full detector simulation and current beam parameters [10,11].
In this section, we briefly summarize the basic concept of bilinear R parity violation and its connection to collider physics. The superpotential and the corresponding soft SUSY breaking terms in bRPV SUSY have the form where i = {e, μ, τ } is the generation index. H u indicates the SU (2) doublet of the Higgs superfield and L i the SU (2) doublet of the lepton superfield. i and B i are bRPV parameters. In addition to that, the three sneutrinos acquire a vacuum expectation value (VEV) ν i = v i . Because of three additional tadpole equations one ends up with six free parameters for bilinear R parity violation. These parameters can be fixed by fitting them to neutrino observables, like neutrino mass differences and mixing angles. The introduction of lepton number violation allows the neutrinos to mix with the other neutral fermions of the model, i.e. the gauginos and higgsinos. Thus, in the basis of neutral fermions 0T = B, W 0 , H 0 d , H 0 u , ν e , ν μ , ν τ the corresponding mass term in the Lagrangian looks like where the mass matrix M N has additional off-diagonal entries due to bilinear R parity breaking. Diagonalizing M N generates one neutrino mass at tree level as well as two neutrino mixing angles. The atmospheric neutrino mixing angle, for instance, reads tan(θ 23 where i = μv i + v d i are so called alignment parameters. Herein, μ is the MSSM higgsino mass parameter and v d represents the VEV of the down-type Higgs. It has been show in [7] that the remaining neutrino mixing angle and neutrino masses can be derived on 1-loop level. A very interesting feature of this model is that the lefthanded part of the χ 0 1 − W − l i coupling is approximately proportional to the alignment parameters where f is a function of the soft SUSY breaking parameters. The full expression of O χ 0 1 Wl i can be found in [8]. Combining Eq. (5) with Eq. (6) makes it clear that neutrino mixing can be determined from measuring branching rations of the neutralino decays: The exact relation only holds at tree level. Reference [8] shows that via loop contributions additional SUSY parameters enter into the mass matrix and thus into the relation between branching ratios and neutrino mixing angles. For this study, we use SPheno3.2.4beta [12,13] to fit the alignment parameters to current neutrino data and to extract the physical observables, in particular the neutralino branching ratios, at loop level. This procedure is repeated for every considered point in SUSY parameter space. It is worth mentioning that for bRPV SUSY there is always a connection between LSP decays and neutrino physics independently of the type of the LSP, which is shown in [9]. The remainder of the paper is organized as follows: In Sect. 2 we summarize the status of RPV searches at the LHC. In the following Sect. 3 we focus on the study at the ILC. Here we introduce the studied simplified scenario, give some details of the used Monte-Carlo samples, define the event selection and comment on sources for systematic uncertainties. The results are presented in Sect. 4, where we discuss the LSP mass measurement, the expected signal significance in the parameter space of the simplified model and the precision in measuring the atmospheric mixing angle from the ratio of two LSP branching ratios. We finally conclude by summarizing our results obtained in Sect. 5.

Status at the LHC
In pp collisions the dominant SUSY production mode is via squark and gluino production. Those colored particles then decay via cascades down to the lightest SUSY particle (LSP), which in the RPV case then decays into Standard Model particles. The main difference with studies assuming R parity conservation is that the cut on missing energy is relaxed significantly, since the LSP does not escape undetected anymore.
The ATLAS collaboration performed a dedicated bRPV SUSY search in the framework of the CMSSM, where the RPV parameters have been fitted to neutrino data and have not been taken as free [14]. In this study contributions from all possible production modes have been taken into account. This study excludes a wide range of the CMSSM parameter plane reaching up to m 1/2 ≈ 600 GeV or m 0 ≈ 1.2 TeV. However, most of the exclusions of the parameter space result from the limits on colored particles for the specific parameter points.
Except for this, various RPV SUSY searches have been performed in the simplified model framework [15][16][17][18][19][20][21][22]. Many of these studies assume strong production, which is dominant for not too high squark or gluino masses. So, the derived limits are again predominantly limits on the colored sector of the model and the electroweak sector remains untested.
The LHC now starts to become sensitive to direct electroweakino production and is able to set limits on electroweakino masses in the RPC [23][24][25] as well as in the RPV case [26]. However, the cross section for direct χ 0 1 χ 0 depend on the assumed RPV coupling. Under the most optimistic assumptions, chargino masses up to 750 GeV have been probed. Due to the strong dependency on the assumed type and strength of the RPV couplings, it is not possible to directly re-interpret these limits in a bRPV scenario which accounts for neutrino data.

Model definition
At the ILC the situation is complementary to the LHC: Here, direct electroweakino production is dominant and the electroweak sector can be probed directly, including LSP pair production.
In our scenario we assume the lightest neutralino to be a bino, which leaves only the t/u-channel production for direct χ 0 1 pair production (see Fig. 1). In presence of bRPV couplings, not only selectrons are possible as exchange particle, but due to the additional terms in principle all other charged scalars could contribute. However, these contributions are strongly suppressed by the small RPV couplings. We define a simplified model in which we set all masses of the SUSY particles to the multi-TeV scale, except for m χ 0 1 and m e R . Thus, those two parameters fix the production cross section.
In the case of a wino LSP, the cross section would drop due to the missing coupling to the right selectron. For a light left selectron, however, the situation is comparable to the bino case. A higgsino-like LSP would allow s-channel associate production of χ 0 1 and χ 0 2 , predominantly via a Z boson. The cross section for this production process is about 100−200 fb [27]. In a light higgsino scenario χ 0 1 and χ 0 2 are usually close in mass and the decay products of χ 0 2 to χ 0 1 are rather soft. Thus, experimentally the situation is comparable to direct χ 0 1 pair production.
At the ILC, the electron beam is polarized to 80 % and the positron beam to up to 60 % [11]. The advantages of beam polarization have been discussed in many studies. A comprehensive overview can be found, for instance, in [28]. In the case of t/u-channel production with selectron exchange different combinations of beam polarization influence the production cross section significantly. Figure 2 shows the cross section in theẽ R -χ 0 1 mass plane for unpolarized beams (left) and for the baseline polarization of P(e + , e − ) = (−30 %, +80 %) (right), which enhances the cross section considerably.
Since the RPV couplings are very small if they are used to describe neutrino data correctly, light LSPs can become rather long-lived with a decay length of meters up to kilometers and escape the detector. Therefore, scenarios with very light LSPs would behave very similar to RPC SUSY scenarios, where the LSPs are stable. However, those escaping pair-produced LSPs could still be detected via radiative LSP production for cross sections of O(10 fb) as demonstrated in [29]. In this study, we focus on scenarios where the on-shell W -decay channel is available, i.e. on LSP masses larger than the W mass. Thus, the decay length of the LSP lies between  10 cm and 100 µm as depicted in Fig. 3 (left). The right-hand panel of Fig. 3 shows the impact parameter resolution of the ILD detector concept as determined from full detector simulation [10]. It indicates that the arising displaced vertices in this scenarios are well detectable for ILC detectors.
For the ILC study the following example point has been used: The mass of the neutralino has been selected as a worst case scenario, where m χ 0 1 m W/Z . This is most challenging, since in this case the LSP signal is expected to overlay significantly with SM background and the involved leptons from the LSP decay to W become relatively soft. For higher LSP masses the study is almost SM background free. The production cross section for the polarization P(e + , e − ) = (−30 %, +80 %) and a center of mass energy of √ s = 500 GeV amounts to 344 fb, driven by the choice of m e R . As can be seen from Fig. 2, the production cross section is still ∼100 fb when m e R is twice as large. The branching ratios for the decay modes that are relevant for measuring the neutrino atmospheric mixing angle [see Eq.
The remaining fraction comes mainly from on-shell Z decay modes χ 0 1 → ν i Z , where the relative fraction of ν i Z vs. l i W decays depends only on the χ 0 1 mass. Decays to electrons have a negligibly small branching fraction due to the smallness of the reactor neutrino mixing angle θ 13 , since BR(χ 0 1 → W e) ∝ e ∝ tan 2 (θ 13 ) [8]. Three-body decays are negligible in our case. They become sizable only for small τχ 0 1 mass differences, which would allow direct τ production at the ILC and thus offer a wealth of additional observables which we do not discuss in this paper. 1 These branching ratios correspond to bRPV couplings O χ 0 1 Wl2,3 in the order of 10 −7 .
For an integrated luminosity of 100 fb −1 one expects 34,400 produced neutralino pairs and among them 6,361 events ending up in the μμ-channel, 7,599 events in the τ τ -channel, and 13,904 events in the mixed μτ -channel. As the ILC is expected to deliver an integrated luminosity of 250 fb −1 /year at 500 GeV, this amount of data is going to be collected within approximately five months of operation at design luminosity.

Data samples
For the given example point a full detector simulation of the International Large Detector (ILD) based on the recently published detector description [10] has been performed. The ILD concept is one of two proposed detector concepts at the ILC. Its design is optimized for Particle Flow reconstruction, which aims at reconstructing each individual particle with the most precise detector component. To this end, properties of charged particles are only measured by the tracking system and the highly segmented calorimeters are only used for the measurement of neutral particles. The main part of the proposed tracking system at the ILD is a time-projection chamber, which is complemented by silicon strip and pixel detectors. This system is expected to obtain a tracking resolution of up to σ 1/ p t = 2 · 10 −5 GeV −1 for high momenta. Due to the more benign radiation environment at the ILC, the inner detectors can be built with very little material, amounting to only 10 % of a radiation length in front of the electromagnetic calorimeter in the barrel region. The electromagnetic calorimeter is foreseen as a 30-layer silicontungsten sampling calorimeter with cell sizes of 5 × 5 mm 2 . In testbeam operation of a prototype [30], an energy resolution of E/E = (16.6 ± 0.1) %/ √ E ⊕ (1.1 ± 0.1) % has been achieved. The highly segmented hadronic calorimeter with cell sizes of 3 × 3 cm 2 in connection with the Particle Flow Concept allows for a jet energy resolution of The full detector simulation has been performed for an integrated luminosity of 100 fb −1 and at a center of mass energy of √ s = 500 GeV. For the SM background, samples produced for the benchmarking of the ILD detector for the Technical Design Report have been used [10].
In the case of bRPV events the program Sarah [32] has been used to generate model files for the event generator. As for the SM samples, Whizard [33,34] has been used as event generator of the hard process and Pythia [35] for fragmentation and hadronization. These events have been passed through Mokka [36], the full Geant4-based [37,38] simulation of the ILD detector and finally reconstructed with MarlinReco [39,40]. For the event generation, realistic beam parameters have been taken into account, in particular the ILC specific beam energy spectrum at 500 GeV [11].
With the instantaneous luminosity foreseen at the ILC, on average N = 1.2 interactions of photons leading to the production of low p t hadrons are expected per bunch-crossing [10]. This takes into account contributions from real photons accompanying the electron beam due to bremsstrahlung and synchrotron radiation as well as from virtual photons radiated off the primary beam electrons. Therefore, each hard-interaction event (from SUSY or SM background) has been overlaid with a Poissonian distributed number of such γ γ → hadrons events before the reconstruction step.
At the time of the MC production a too large number of overlay events per hard-interaction had been assumed ( N = 1.7), which results in a conservative estimate of the γ γ → hadrons background in this study. However, as will be demonstrated in the next section, even this larger background could be removed very efficiently.

Event selection
The signal events have a rather clear signature: The produced LSPs decay into either a μ or τ plus a W boson. In the following we restrict ourselves to the hadronic W decay mode.
Thereby, the event is-except for some missing energy from a potential τ decay-fully reconstructable with six visible objects in the final state.

Event preparation and preselection
Before the actual event selection is performed, the γ γ → hadrons background is removed from the event by the following procedure: Since we expect to have six final state objects, an exclusive k T jet clustering algorithm (R = 1.3) [41,42] which is forced to find six jets is applied to the reconstructed objects in the events. This algorithm builds up six jets which are assigned to the hard interaction, as well as two very forward directed beam jets which are treated as beam background. Removing those beam jets and using only the objects which end up clustered in the main jets for the further analysis, recovers very well the bare event without background overlay. Figure 4 shows the impact of overlaid γ γ events on the visible energy (left) and the ability to remove this background with the described method (right).
Two hadronically decaying W bosons imply a relatively high particle multiplicity N objects in the event. Due to the fact that there is no major source for missing energy, the visible energy in the event is close to the center of mass energy of 500 GeV. So, the following preselection cuts have been used: This preselection on the one hand reduces Standard Model background and on the other hand cuts away some of the leptonic W contribution in the LSP decay, which is considered as background in this analysis. 95 % of the signal events pass this preselection. In order to be able to measure the ratio of different branching ratios of the LSP decay [see Eq. (7)] we have to define a selection to distinguish between the different event classes.

μμ class
Here we look for an event with at least two reconstructed muons. The muon identification is provided by the Pandora Particle Flow Algorithms [31], seeded by a minimumionizing signature in the calorimeters and the instrumented flux return yoke. The two most energetic muons are removed from the event. No further lepton isolation criterion is required. The remaining objects of the event are clustered by the Durham jets clustering algorithm into four jets. All pair-permutations of the four jets are used to find the best W candidates, where the objects with the smallest are used for the proceeding analysis. Herein, σ res is an estimated resolution factor of 5 GeV. For all combinations of W candidates and muons the invariant mass is determined and the reconstructed objets from the pair with smallest are considered as LSP candidates. If the event fulfills the condition the event is counted as μμ event, else it is tested against the μτ class.

μτ class
In this class at least one reconstructed muon in the event is required. The most energetic muon is removed and the rest of the event is forced into five jets by an inclusive Durham jet clustering algorithm. Since the muon requirement is relaxed, this class is more prone to background. Therefore we test in addition whether the chosen jet configuration describes the event well. We employ a cut on the Durham jet algorithm parameters y i−1,i and y i,i+1 [41]. y i−1,i is a measure of the distance in energy-momentum space between the two of the i jets which would be merged if a i − 1 configuration was required. y i,i+1 gives the analogous distance measure for the two jets which have been merged in the last step when clustering from the i + 1 to the final i configuration. For events well fitting to the 5-jet configuration, y 4,5 − y 5,6 gets maximal, so we use the following cut: The jet with the smallest number of constituents is considered as τ candidate if N const. < 10. Additionally it is required that the jet does not contain a muon. The further approach to find W candidates and finally LSP candidates is identical to the μμ class. If condition (14) is satisfied this event is counted as μτ event, otherwise the τ τ class is tested.

τ τ class
The event is forced into six jets by an inclusive Durham jet clustering algorithm. In analogy to the μτ class, it is required that y 5,6 − y 6,7 > 5 · 10 −5 . The two jets with the smallest numbers of constituents are considered as τ candidates if each fulfills N const. < 10 and the muon veto. The further approach to finding W candidates and finally LSP candidates is again identical to the μμ and μτ class. If condition (14) is satisfied this event is counted as τ τ event. Otherwise, the event is rejected. After this selection Fig. 5 shows the reconstructed neutralino mass in each of the event classes for 100 fb −1 of fully simulated events. The remaining background is dominated  by SM W pair production, whereas the contribution from other LSP decays is very small.

Mass reconstruction
For the reconstruction of the LSP mass the two most important sources of uncertainties arise from the reconstruction of the μ momentum as well as on the reconstruction of the momentum and energy of the W candidates. In our scenario, the muons originating from an LSP decay carry a momentum of up to 80 GeV, while the jets originating from the W bosons from the signal decay have an energy of up to 200 GeV. The muon momentum scale can be calibrated using Z boson decays. The unpolarized cross section of the process e + e − → Z → μ + μ − at the ILC500 amounts to 2.5 pb. For an assumed integrated luminosity of Ldt = 500 fb −1 this yields in total N μμ = 12.5 · 10 5 muon pairs. Assuming that ultimatively the precision of the scale calibration is limited by the available statistics for the calibration process, we estimate that a precision of 1/ N μμ = 0.09 % can be reached.
The jet energy calibration of hadronically decaying W bosons can be derived analogously from hadronic Z decays. Since the cross section for e + e − → Z → qq is ten times larger compared to the muonic decay channel, which means that the from the point of view of available control sample statistics, the jet energy scale uncertainty could reach 0.03 %, assuming that the calibration is sufficiently stable over time. The resulting impact on the neutralino mass determination is 11 MeV. Thus it could be a factor 4 larger before it becomes comparable to the contribution from the momentum scale. Alternatively, exploiting kinematic fits [43] in conjunction with the well-known beam energy, at the ILC foreseen to be controlled to 10 −4 [10], could significantly reduce the dependence of the reconstructed neutralino mass on the jet energy scale, leading to similar final precision estimates.
Systematic errors on luminosity, beam energy and beam polarization do not enter into the mass measurement, which depends solely on the reconstructed detector signals. The selection efficiency does not show any dependency on the reconstructed LSP mass within the available Monte-Carlo statistics. Thus we conclude that any potential bias due to the selection efficiency plays an insignificant role in the mass measurement.
For LSP masses below ∼105 GeV, the SM background has a steeply falling invariant mass distribution. Therefore any uncertainty related to the modeling of this slope, e.g. assumptions on hadronization, color-reconnection etc., could enter into the LSP mass determination. However, the ILC itself will offer numerous opportunities for SM precision measurements beyond today's knowledge. In particular the cross section for the dominating background process, W pair production, is in the several pb range, thus providing ample possibilities to tune the modeling of this process. We thus estimate that the residual effect on the LSP mass determination via subtraction of the SM background is not larger than 15 MeV. Table 1 summarizes the main sources for the systematic uncertainty on the LSP mass reconstruction and its propagation to the LSP mass reconstruction. The total systematic error is estimated to be about 50 MeV.

Measurement of ratio of branching ratios
The measurement of the ratio of branching ratios is a measurement of the ratio of the number of events reconstructed in the different event classes. For this reason, all systematic uncertainties which factorize with the number of events cancel. The same is true for reconstruction effects which affect all event classes simultaneously, like systematic uncertainties on jet energy scales, for instance.
The main source of systematic error is expected to arise from the determination of the selection efficiencies and purities of the different event classes from Monte-Carlo. Therefore, validation of the Monte Carlo simulation with data is very important.
The process e + e − → Z Z → l + l − qq with l = {μ, τ } offers the possibility to study the Monte-Carlo description of τ + jets and μ + jets events under comparable experimental conditions as for the signal decay. The unpolarized cross section per process at a center of mass energy of 500 GeV is 225 fb. Thus, for Ldt = 500 fb −1 the expected precision on Monte-Carlo and data comparison amounts to 0.3 %. As a second process, e + e − → W + W − → lν l qq with l = {μ, τ } can also be studied in order to validate Monte-Carlo. This process has a significantly larger cross section of 2.5 pb per process, which leads in the end to a statistical uncertainty of 0.09 % for Ldt = 500 fb −1 .
As for the mass measurement, we allow some modeling uncertainty for the main SM background, W W → hadrons, of 0.05 %.
Taking these considerations into account, a conservative estimate on the systematic uncertainties on Monte-Carlo is O(0.5 %) (cf. Table 2).

Mass measurement and resolution
The μμ channel shows a clear signal peak, which can be used for measuring the LSP mass accurately. Figure 6 shows the reconstructed LSP mass spectrum for an integrated luminosity of Ldt = 100 fb −1 after subtracting the SM back- which is in very good agreement with the ILD design goals [44]. Scaling the statistical uncertainty to an integrated luminosity of Ldt = 500 fb −1 and combining the systematic uncertainties as discussed in Sect. 3.4, the total uncertainty on the LSP mass measurement becomes δm fit

Signal significance
The precision measurement of the LSP mass in the μμ channel can be used to define a signal region m fit Thereby, N true LSPBG counts events in which at least one LSP decays differently than the targeted two-body decay χ 0 1 → Wl. As expected, the τ τ channel is by far the one with the lowest purity.
The signal over Standard Model background ratio for the different classes can be derived. Herein, every selected event originating from an LSP decay is counted as signal. The described event selection is still significant on a 5σ level for very large selectron masses well above 1.5 TeV and a large range of m χ 0 1 , as depicted in Fig. 7. Positron beam polarization further enhances the production cross section and, thus, increases the sensitivity of the analysis to selectron masses of almost 2 TeV for P(e + ) = −60 %.
It has already been pointed out that this studied parameter point is a worst case scenario with respect to the neutralino mass. It is clearly visible in Fig. 5 that for another parameter point with a higher LSP mass the signal peak would shift into an almost completely background-free region. Though, the LSPs in this model become rather long-lived at lower masses (compare Fig. 3), where the Standard Model background is large. Adding this information to the analysis and requiring from the reconstructed objects not to point to the primary vertex, would reduce the Standard Model background drastically. However, this would on the other hand introduce a strong model dependency to the analysis. Therefore the exploitation of the lifetime information is left as a future option for improvements.
Likewise, the requirement for same-sign leptons in the event classes is an option for further improving the analysis. This restriction would reduce the number of signal events by a factor of 2, but could heavily suppress the remaining SM background.

Branching ratio measurement
We assume in the following that the average number of Standard Model background events can be estimated from Monte Carlo with a precision of 0.05 % (cf. Sect. 3.4). The LSP non-signal background consists mainly of events in which one of the two LSPs decayed non-signal like into Z ν. As soon as the LSP mass is known, the relative fraction BR( χ 0 1 → Z ν)/ BR( χ 0 1 → Wl) is determined and, thus, also the number of LSP background events can be predicted. Under this assumption we can subtract the backgrounds and build a 3 × 3 efficiency matrix E, which is defined like where i, j = μμ, μτ, τ τ . The error on the entries of the efficiency matrix is dominated by the assumed systematic uncertainty of 0.5 % on the Monte-Carlo prediction on the migrations between the signal classes (cf. Sect. 3.4). The vector of selected events becomes Herein, δ N MC BG,ll with ll = {μμ, μτ, τ τ } are the systematic errors on the SM background estimation in the event classes. They are negligible compared to the statistical fluctuations of the reconstructed events per event class.
The efficiency matrix can then be inverted and used to unfold the different event classes obtaining the number of reconstructed events in the event classes The ratio of the two branching rations can be extracted in different ways since for the expected number of events the following relations hold: where N χ 0 Ldt is the number of produced LSP pairs. Because of the low selection purity in the τ τ channel, we have chosen the relation involving only the μμ and μτ channel for the further analysis. Using all three relations in Eq. (23) as input for a constrained fit can improve the precision by a factor of 2, but this is not pursued in the following.
In order to estimate the uncertainty on the resulting ratio of event numbers, an error propagation has been performed. The uncertainty is depicted in Fig. 8 in dependence of the integrated luminosity at the ILC500. For the studied parameter point the achievable precision for Ldt = 100 fb −1 is about 9 %. This scales down to roughly 4 % for 500 fb −1 , which is the desired integrated luminosity at ILC500 in a first stage. The uncertainty contains a systematic error of 0.85 % arising from the propagation of the assumed systematic uncertainties on the efficiencies Table 2 through the unfolding procedure.

Neutrino interpretation
The measured ratio of branching ratios can now be translated into the atmospheric mixing angle following Eq. (7). As mentioned earlier, the given relation is only valid on tree level and there are additional parametric uncertainties coming from residual SUSY parameter dependencies [8]. For this reason we define two scenarios: In the first scenario we assume that the LSP is the only accessible SUSY particle at ILC500.  Fig. 9 Random scan of the SUSY parameter space [45] in order to estimate the parametric uncertainty on the correlation between the atmospheric mixing angle and the ratio of branching ratios BR( χ 0 All other supersymmetric particles are randomly chosen to be heavier than 300 GeV. The result of a random parameter scan consisting of 6,000 scan points [45] is depicted in Fig. 9. We find that for 95 % of all found viable SUSY parameter points the deviation of the correlation between atmospheric mixing angle and the ratio of the branching ratios is below 17 %. Assuming that the remaining three electroweakinos are measureable at ILC1000, which is an optional upgrade of the ILC to a center of mass energy of 1 TeV, the correlation uncertainty reduces to 7 %. The derived precision of the measurement of the atmospheric neutrino mixing angle for different assumed parametric uncertainties is shown in Fig. 10 (left). One can now compare this precision with the uncertainty of current neutrino experiments [46], which is done in Fig. 10 (right). The middle red line indicates the best fit value of the atmospheric neutrino mixing angle and the upper and lower dashed red lines indicate the 1σ uncertainty. An agreement between the collider and neutrino experiment data would clearly establish bRPV as origin of neutrino masses. Improvements from future neutrino experiments or a reduction of the parametric uncertainty by observation of additional SUSY particles at the ILC or the LHC would strengthen this conclusion even further.

Conclusions
We have presented a full ILD detector simulation of a bRPV SUSY model, which is an attractive possibility to explain neutrino mass generation and mixing. A highly detailed ILD We have developed a model-independent selection strategy to disentangle the different event classes involving the two decay modes of the LSP χ 0 1 → μ ± W ∓ and χ 0 1 → τ ± W ∓ . It has been demonstrated that in the μμ event class a very accurate mass measurement with an uncertainty of δm fit χ 0 1 = (40(stat.) ⊕ 50(syst.)) MeV is possible for an integrated luminosity of 500 fb −1 . With the described selection, a signal to background ratio of 130 in the μμ event class and 60 in the μτ event class has been achieved. Even for very large selectron masses of up to 1.5 TeV a 5σ discovery is possible for a large range of m χ 0 1 . The μμ and μτ event class have been used to determine the ratio of the two branching ratios BR( χ 0 1 → μ ± W ∓ )/ BR( χ 0 1 → τ ± W ∓ ), which is related to the atmospheric neutrino mixing angle sin 2 θ 23 . For an integrated luminosity of 500 fb −1 the total uncertainty on this ratio, including statistical and systematic uncertainties, has been determined to 4 %.
Finally, we have shown that the precision in measuring the atmospheric neutrino mixing angle is in the same range than measurements from neutrino oscillation experiments, even when taking parametric uncertainties due to the unknown parts of the SUSY spectrum into account. Therefore, the International Linear Collider is highly capable to test bRPV SUSY as origin of neutrino masses and mixings.