Atmospheric neutrino oscillations from upward throughgoing muon multiple scattering in MACRO

The energy of atmospheric neutrinos detected by MACRO was estimated using multiple coulomb scattering of upward throughgoing muons. This analysis allows a test of atmospheric neutrino oscillations, relying on the distortion of the muon energy distribution. These results have been combined with those coming from the upward throughgoing muon angular distribution only. Both analyses are independent of the neutrino flux normalization and provide strong evidence, above the 4 sigma level, in favour of neutrino oscillations.


Introduction
The results obtained by experiments looking for neutrinos coming from natural or artificial far sources gave the first indication of physics beyond the standard model.Neutrino flavor changing is in fact the most straightforward explanation for electron and muon neutrino disappearance and for the evidence of non electron neutrinos from the sun.In this scenario, the study of atmospheric neutrinos plays an important role to improve our understanding of the neutrino oscillation mechanism.
An atmospheric muon neutrino deficit was found in all of the three categories and an angular distribution distortion was found in category (1).Such results are explained by the neutrino oscillation hypothesis with parameters ∆m 2 = 2.5 × 10 −3 eV 2 and sin 2 2θ=1, in good agree-0 a Also INFN Milano, 20133   ment with the Super-Kamiokande results [4].
A detailed study of the upward throughgoing muon angular distribution by MACRO [3] and by Super-Kamiokande [5] allowed the exclusion at the 99%C.L. of the muon neutrino into sterile neutrino oscillation mechanism, compared to the muon neutrino into the tau neutrino.Considering a two flavor neutrino oscillation, the probability for a ν µ to oscillate to ν τ is given by: where L ν (km) is the distance between the neutrino production and interaction points and E ν (GeV) is the neutrino energy.The different categories of events quoted above have median energy from 4 GeV to 50 GeV, which provides evidence of the dependence on neutrino energies, as required by the oscillation hypothesis.In this paper we address for the first time the estimate of the upward throughgoing muon energy by using multiple scattering.The implementation of this method to the MACRO data is discussed in Section 2. Two different analysis were performed.The first one, using the tracking system in digital mode, is described in Section 3. Section 4 shows the results obtained with the streamer tube in drift mode, where electronic readout of the hit time was used to improve the position resolution and thereby the range of muon energies that could be estimated.

The muon energy estimate
The MACRO detector was extensively described in [7].The detector consisted of a lower half (the lower detector) and an upper half (the "Attico").The lower detector had 10 horizontal planes of streamer tubes separated by either crushed rock absorbers (60 g/cm 2 ) or (at its top and bottom) a plane of scintillator counters.The Attico had 4 horizontal planes of streamer tubes separated into two groups with a plane of scintillator counters between them.
Upward throughgoing muons are mainly produced in neutrino deep inelastic scattering (DIS) in the rock below the detector.The recoil hadrons are lost and the muon energy is degraded in the propagation to the detector.Nevertheless, Monte Carlo simulations show a linear relation between the parent neutrino energy E ν and the muon residual energy E µ at detector level.
The momentum resolution obtainable from multiple coulomb scattering (MCS) measurements is the result of two different contributions: the number of sampling planes N and the space resolution of the detector σ 0 .In MACRO the number of tracking planes interleaved with rock absorbers gives N = 8.The resolution of the streamer tube system used in digital mode is σ 0 ≃1 cm and was improved to σ 0 ≃0.3 cm by using the drift time.The other six horizontal tracking planes are separated by a negligible amount of material, and do not contribute to the MCS measurements.
In a tracking detector with equispaced tracking planes, separated by a given absorber and neglecting the energy losses, the characteristic momentum scale can be written as [8]: where ∆ is the distance between tracking planes and X 0 is the material radiation length.For p<p MCS , the main limitation to the momentum reconstruction comes from the number of sampling planes while for p>p MCS the space resolution dominates.Under the conditions quoted above the relative momentum error [8] can be written as: for p<<p MCS , and for p>>p MCS .As shown in [8] the 1/p distribution approaches a Gaussian only for N≥30; therefore the momentum error for MACRO is not expected to have a Gaussian behaviour.The eight horizontal streamer tube planes interleaved with rock absorbers (planes 2 to 9, from the bottom) are equispaced and have the same space resolution, hence p MACRO MCS =2.2 GeV/c.Since about 90% of upward throughgoing muons have p > 2.2 GeV/c, the intrinsic chamber resolution dominates.The MCS-based momentum estimates scale therefore as the square of the space resolution σ o .
The r.m.s. of the lateral displacement of a relativistic muon crossing a layer of material with depth X is proportional to the inverse of muon momentum p µ [9]: In MACRO, X≃25X o /cos Θ, giving on the vertical σ MCS proj ≃10 cm/E µ (GeV).Therefore, a saturation point above 10 GeV requires a space resolution better than 1 cm.An improvement in the space resolution increases the maximum energy value (saturation point) above which the MCS method is not effective.
Consequently we decided to improve the streamer tube resolution: in the second analysis, described in Section 4, we took advantage of the QTP-TDC electronics, developed for magnetic monopole searches [7].In this case, we obtained an improved space resolution, σ x ∼0.3 cm [11], which allowed us to estimate the muon energy up to E µ ≃40 GeV.
For both analyses, we used the whole sample of upward throughgoing muon events collected with the complete apparatus in a period of data taking equivalent to 5.5 years of live time.We studied upward muon events selected by the timeof-flight measured by planes of scintillators combined with the standard MACRO tracking algorithm.The original upward throughgoing muon data set has been described previously [3].To make a comparison between real data and expectation we performed a Monte Carlo simulation using the Bartol neutrino flux [12] and the GRV94 DIS parton distributions [13] for deep inelastic scattering.For low energy neutrino interactions we used the cross sections given in [14].The energy loss for muons propagating through rock is taken from Lohmann et al. [15] while the muon simulation inside the detector was performed with GMACRO (the GEANT 3.21 based detector simulation).For the second analysis, where the streamer system is used in drift mode, another simulation chain has been implemented.In this case neutrino interactions are randomly distributed in a rock semi-sphere below the detector.Muons are then transported to the detector using the FLUKA99 package [6].A Monte Carlo statistics corresponding to a live time of ≃2700 years was produced (500 equivalent experiments).This simulation was compared with that used in [1,3], obtaining a satisfactory result.

The first analysis: streamer tubes in digital mode
Because of MCS in the lower detector, with this analysis we expect a measurable deflection between the incoming and the outgoing track directions for muons with energies smaller than ∼10 GeV.This corresponds to an effective armlever of ∼4 m between the lower detector and the Attico [7].
The eight lowest streamer tube planes were used, through a track refit, to estimate the direction of the incoming muon.The five upper streamer tube planes in the lower detector and the four Attico planes were used to estimate the direction of the outgoing muon.The distance r w between the intercepts of the two tracks in the z = 0 plane, and the difference ∆Φ between the two slopes depend on the muon energy E µ .We divided the upward throughgoing muons in three sub samples, according to the values of r w and ∆Φ: sample L= Low energy, if r w >3 cm and ∆Φ > 0.3 • ; H= High energy if r w ≤3 cm and ∆Φ ≤ 0.3 • .The remaining events were classified as M=Medium Energy.These cuts were optimized using two large samples of real atmospheric muons: i) downward throughgoing (average energy of ∼ 300 GeV [16]), and ii) downward going stopping muons (average energy of ∼1 GeV).
The aforementioned cuts were applied on the upward throughgoing muons (crossing the whole apparatus, lower detector + Attico), both real and simulated.They were selected from the same analysis chain [1] and had the same data format.From the simulation we expect 430 events: 178 of Low energy, with <E ν >=11 GeV, 59 of Medium energy, with <E ν >=33 GeV and 193 of High energy, with <E ν >=72 GeV.From the real data, 316 events were selected: 101 of Low, 51 of Medium and 164 of High energy.
The L and H events were further divided according to their zenith angle Θ: events from the vertical direction −1.< cos Θ < −0.8, and events with cos Θ > −0.8.For each event topology the number of detected and expected events were determined and ordered with decreasing value of the <L ν >/ <E ν > (<L ν >= average value of neutrino path length).The relative systematic uncertainty on each of the five <L ν >/<E ν > values is 12.5%, coming from the uncertainties on the neutrino spectrum and angular shape (discussed in Sect.4.3), from the detector related effects and analysis cuts.This includes the number of planes used for the refit, the definition of the cuts, the fluctuations in the streamer tube and scintillator efficiencies, and detector acceptance uncertainties.
We evaluated χ 2 for the hypothesis of no neutrino oscillations using the five <L ν >/<E ν > bins, plus the additional point from the analysis of semi contained upward going events (IU) [2] (upgoing neutrinos with <E ν >∼4 GeV).We found that the distribution agrees with the no oscillation hypothesis with a probability lower than 2%).For two-flavor oscillations with parameters ∆m 2 = 2.5 × 10 −3 eV 2 and sin 2 2θ = 1 we get a probability of 45%.

The second analysis: streamer tubes in drift mode
The performance of the streamer tube system, read out in drift mode, is described in [11].Here an absolute muon energy calibration was performed at the PS-T9 and SPS-X7 beams, where a slice of the MACRO detector was reproduced.The MCS information was handled by a neural network, based on JETNET 3.0 [17], calibrated with the muon beams quoted above.The neural network(NN) output obtained using test beam data, was compared with that expected from the Monte Carlo simulation, obtaining a satisfactory agreement [11].The application of such analysis to the MACRO data resulted in an improvement of the space resolution from σ x ≃1 cm to σ x ≃0.3 cm.
Figure 2 shows the Monte Carlo prediction for NN downward throughgoing muons and downward going stopping muons compared with experimental data.A nice agreement between simulation and real data is found for both categories, representing a large energy interval.
The result of the neural network output energy calibration shows that NN output is almost linear with log 10 (E µ ), increasing with the muon energy up to E µ ≃40 GeV, where a saturation effect occurs.With respect to the approach used in [11] few details of the neural network were optimized: the neural network training and the energy calibration was performed separately for events with hits in the upper part of the detector(Attico).
Although smeared by the energy carried away by hadrons and by energy loss in the rock, the detected neutrino induced muons still carry on average ≃40% of the original neutrino energy.By using the full Monte Carlo simulation quoted in Sect.2, we calibrated the NN output as a function of log 10 (E ν ): the calibrated NN output is linear up to log 10 (E ν (GeV))≃2.15.We may write: where δ log 10 (E ν ) is the difference between the log 10 of the real and of the reconstructed neutrino energy.Taking into account that the NN output was calibrated as a function of log 10 (E ν ), the energy resolution δE ν /E ν was obtained by plotting the quantity δ log 10 (E ν )/log 10 (e).In Fig. 3 we show the resolution that could be obtained with an ideal muon energy resolution(dotted line) and that obtained with the present analysis(continuous line).The precision of the neutrino energy estimate obtained with an ideal muon energy resolution detector is δE ν /E ν ≃70%, while with the present method a resolution of δE ν /E ν ≃150% is obtained.The asymmetry in both curves comes from neutrino interactions occurring far from the detector, for which a large fraction of the muon energy is lost during the transport.

Data Selection
We used the 783 upward throughgoing muon data set collected in the full detector run started in 1994.Table 1 describes further event selection to arrive at the sample for MCS analysis.We required a single track in both the wire and the strip view.We selected hits belonging to the track and made of a single fired tube, to associate unam- biguously its QTP-TDC time information.This cut is effective for muon tracks with large zenith angles (Θ > 30 • ), while it is quite loose around the vertical: we restricted the present analysis to events with Θ≤60 • .Comparisons were performed between simulated and experimental downward going muon data, to ensure that the selection efficiency, in the used angular window, is the same in the Monte Carlo and in the real data.Spurious background hits have been avoided by requiring a time window of 2 µs around the trigger time.Finally, we selected events with at least four streamer tube planes with valid QTP-TDC data.We fitted the drift circles using the same tracking developed to analyse test beam muons.
A minimum path length of 200 cm in the lower detector is required for tracks hitting the Attico and 400 cm for tracks not hitting it.These geometric requirement ensure a minimum depth of material where muons may experience a measurable amount of multiple scattering and a lever arm long enough for a comfortable tracking.After the selection cuts 300 events survived, giving an overall efficiency of 38.3%.Data selection cuts and selection efficiencies.

Qualitative tests
We used the information provided by the neural network to separate the neutrino events into four energy subsamples, as shown in Table 2.The same selection was applied to simulated events.
Figure 4 shows the zenith angle distributions

Sample Energy cuts Median energy
Table 2 Subsamples selected according with the reconstructed neutrino energy.
of the upward throughgoing muons in the four energy subsamples, compared with the expectations of the Monte Carlo simulation, assuming no-oscillations(dotted line) and oscillations with parameters ∆m 2 = 2.5×10 −3 eV 2 and sin 2 2θ ≃1 (wall boxes).The Monte Carlo including the oscillation hypothesis, with the parameters quoted above, reasonably reproduces the real data in each subsample.We point out the strong differ- ence between the oscillations/no-oscillations hypotheses at low energies, while such difference is reduced by increasing the reconstructed neutrino energy.
Finally, we used information on the ratio L ν /E ν .The distance L ν , travelled by the neutrinos from production to the interaction points, was measured by MACRO relying on the muon zenith angle determination, with a precision ∆L ν /L ν ∼ 3%.The resolution on the ratio L ν /E ν is therefore fully dominated by the uncertainty in the neutrino energy estimate, giving a relative error of ≃ 150%.The ratio DATA/Monte Carlo as a function of log 10 (L ν /E ν ), is plotted in Fig. 5.The black circles are the real data over the Monte Carlo predictions, assuming no oscillations, the shaded regions represent the Monte Carlo predictions with oscillations, ∆m 2 = 2.5×10 −3 eV 2 and sin 2 2θ=1, divided by the Monte Carlo predictions with no oscillation.
The log 10 ( Lν Eν ) distribution of the neutrinos de- tected by MACRO spans from 0.8 to 3.5.The left-pointing arrow at low log 10 ( Lν Eν ) represents the effect of the neural network saturation.Since the maximum energy reconstructed by the NN is E ν ≃140 GeV and since the minimum L ν for the present analysis is L min ν ≃6500 km, the minimum estimated value is log 10 ( Lν Eν ) min ≃1.7.As far as the right-pointing arrow at high Lν Eν is concerned, the reconstruction of neutrino energy based on the residual muon energy, does not allow one to reconstruct the ratio L/E beyond a given limit.This is due to far neutrino interactions, originated by high energy muons, which lost a large fraction of its energy.An ideal detector with ideal muon residual energy resolution, requiring E µ <2 GeV, would select an average log 10 ( Lν Eν )=3.2.Finally, the neural network resolution in muon energy estimate results in a maximum value log 10 ( Lν Eν )=3.The dashed line at 1 is the expectation without oscillations.The last point(empty circle) is obtained from semicontained upward going muon rate [2].It is not used for the evaluation of the oscillation probabilities (Sect.4.4) nor in the allowed region plot(Fig.7).Good agreement is found with the oscillations expected with the pa-rameters quoted above.

Experimental results
To quantify the contribution of the multiple scattering measurement in stand-alone mode, we performed a blind analysis using the Monte Carlo data, looking for the variable, based on the energy estimate, showing the maximum sensitivity to separate the oscillation from the no-oscillation hypothesis.We found that the best performance is given by the ratio: R = N low /N high (7) where N low and N high are the number of events with E rec ν <30 GeV and E rec ν >130 GeV respectively.
We considered systematic uncertainties due to the neutrino flux calculation and neutrino cross section.Due to the large uncertainty in the absolute flux [18,19,20,21], we use in this paper only the angular distribution and the ratio between different event categories selected according to the reconstructed energy.Nevertheless, the theoretical uncertainty of the existing neutrino flux, based on the old CR spectrum [12], has to be accounted for.We varied the input primary cosmic ray spectral index γ in our simulation by ∆γ=±0.05,obtaining a theoretical error on R, ∆R R f lux =±13%.Another source of systematics comes from the neutrino cross section.We looked at R varying the Monte Carlo input cross sections.The most important contribution to the error comes from the "Low" category, since the cross section at low neutrino energy is more uncertain.By comparing the cross section computed under several hypotheses, varying for instance the structure function ( [13], [23]) in the deep inelastic scattering, including or neglecting the contribution of resonant scattering, we found ∆R R σ =9 %.We estimated a total theoretical error The systematic error on R reconstruction, evaluated in 6%, includes the uncertainties on absorber density, drift velocity, streamer tube efficiency and detector acceptance. .Low energy events over high energy events as a function of ∆m 2 : the area between the solid lines includes a 17% systematics(the error is not gaussian, see text).The hatched band represents the real data.
Figure 6 shows the ratio R as a function of ∆m 2 , assuming maximal mixing, for the Monte Carlo simulation: the area between the two solid curves represents the 17% systematic error.The Monte Carlo prediction in case of no oscillations is R ′ =1.5±0.25(theo.+sys.), while for ∆m 2 =2.5×10 −3 eV 2 and sin 2 2θ =1, R=1.00±0.17(theo.+sys.).As pointed out in [3], the ratio does not have a Gaussian distribution.The errors on the ratio are therefore also not Gaussian: they are reported just to give a crude estimates of the significance.The experimental ratio is R exp =0.85±0.16(stat.).The one sided probability of measuring a value smaller than the measured one, was computed by using two different methods.In the first one we let the Monte Carlo simulation (no oscillations) to fluctuate according to statistical and systematic errors of the considered ratio.We then evaluated the fraction of events giving a value smaller than the measured one.In a more pessimistic view, rejecting the hypothesis that the lower number of events detected by MACRO with respect to the Monte Carlo expectation is due to oscillation effects, the total number of Monte Carlo events (sum of the numerator and the denominator in the ratios), were normalized to the experimentally measured sum and then they were allowed to fluctuate.
The corresponding one sided probability of measuring a value smaller than R exp according to method 1 (2) is 0.75%(1.9%)corresponding to 2.4σ(2.1 σ).Finally, we combined the information coming from the energy estimate and from the angular distribution.We considered the ratio: where N vert is the number of upward throughgoing muon events with cos θ≤-0.7 and N hor is the number of events with cos θ≥-0.4,as discussed in [3].The probability of measuring a value A ′ lower than A exp according to the method 1 (2) is P(A'<A exp )=0.001%(0.01%),corresponding to 4.3 σ(3.7σ).It is worth noting that, combining the two independent probabilities on the ratios R and A, the probability that a fluctuation of the expected values(no oscillations) generates the experimental results is P(A ′ •R ′ <A exp •R exp )=1.3×10 −6 (3.2×10 −5 ), corresponding to 4.7 σ(4 σ).
The 90% confidence level allowed regions in the oscillation parameter space have been computed according to the prescription given in [24].Figure 7 shows the 90% C.L. for the ratio R, the angular distribution [3] and for their combination.

Conclusions
Muon multiple coulomb scattering has been used to estimate the energy of the neutrino induced upward throughgoing muons detected by MACRO.The different tests performed on the data, using the estimated energy (separation in subsamples with different energies, L ν /E ν estimates, etc.) give a consistent picture, all of them supporting the neutrino oscillation hypothesis with parameters ∆m 2 = 2.5×10 −3 eV 2 and sin 2 2θ≃1.To quantify such effect, the ratio R = N low /N high was used, in stand-alone mode and in combination with the angular distribution, A = N vert /N hor .Both of them are independent of the neutrino absolute flux.The significance of the MACRO observation of the neutrino oscillations is above 4σ.Real Data and Monte Carlo ratios R=N low /N high and A=Nvert/N hor .

Figure 1 .
Figure 1.Cross section of the MACRO detector and topologies of events induced by neutrino interactions.

Figure 2 .
Figure 2. Neural network(NN) output for downthroughgoing muons: Real data (black squares) and Monte Carlo expectations (continuous line).NN output for downward going stopping muons: real data(empty circle) and Monte Carlo expectations(dotted line).

Figure 3 .
Figure 3. Neutrino energy resolution that can be obtained with an ideal residual muon energy resolution (dotted line) and with the MACRO energy estimate(continuous line).

Figure 4 .
Figure 4. Number of events versus the cosine of the zenith angle Θ for four energy ranges.Black points are the real data, dotted line is the Monte Carlo simulation, assuming no oscillation, and dotted boxes are the Monte Carlo expectation with ∆m 2 =2.5× 10 −3 eV 2 and sin 2 2θ = 1, including a 17% error.

Figure 5 .
Figure 5. Ratio (Data/MC) as a function of the estimated L ν /E ν for the MACRO upward throughgoing muon sample.The black circles are the real data over the MC(no oscillation), the solid line is the MC assuming ∆m 2 = 2.5× 10 −3 eV 2 and sin 2 2θ=1 over the MC with no oscillation.The shaded region represents the MC uncertainties.The last point(empty circle) is obtained from semicontained upward going muons.

Figure 6
Figure 6.Low energy events over high energy events as a function of ∆m 2 : the area between the solid lines includes a 17% systematics(the error is not gaussian, see text).The hatched band represents the real data.