text:
# Are there MACHOs in the Milky Way halo?

## 1 Introduction

Massive compact halo objects (MACHOs) with mass in the range $`10^8M_{}`$ to $`10^3M_{}`$ can be detected via the temporary amplification of background stars which occurs, due to gravitational microlensing, when the MACHO passes close to the line of sight to a background star (Paczyński pac (1986)). In the early 1990s several collaborations began monitoring millions of stars in the Large and Small Magellanic Clouds (LMC and SMC) 111http://wwwmacho.mcmaster.ca/, http://eros.in2p3.fr/, in order to search for Milky Way (MW) compact halo dark matter, and indeed a number of candidate microlensing events have been observed (Alcock et al. 2000b ; Lasserre et al. eros (2000); Milsztajn & Lasserre millas (2001)). The interpretation of these microlensing events is a matter of much debate. Though the lenses may easily reside within the MW halo and constitute halo dark matter, it is also possible that the contribution to the lensing rate due to other, so far unknown non-halo populations of objects has been underestimated (see e.g. Bennett nomacho (1998); Zhao nomacho2 (1999)). In particular the lensing rate due to stars, or MACHOs, within the LMC itself (Wu wu (1994); Auborg et al. sl1 (1999); Salati et al. sl2 (1999); Evans & Kerins sl3 (2000); Gyuk, Dalal & Griest gg (2000)), or a very thick Milky Way disk (Gyuk & Gates gg01 (2001)) may be significant.

The EROS collaboration have five candidate events, resulting from three years of observation of 25 million stars in the LMC (Milsztajn & Lasserre millas (2001)). Due to the small number of events, and poor fit to microlensing of several of these events, they use these events to place constraints on the halo fraction as a function of lens mass, assuming a standard Maxwellian MW halo. In contrast, the MACHO collaboration claim that roughly 20$`\%`$ of the halo is in compact objects with $`M0.5M_{}`$, using their 13/17 candidates events which result from the 5.7 years of observation of 10.7 million stars (Alcock et al. 2000b ). While not in direct conflict with the EROS exclusion limits these MACHO best fit parameters lie just outside the EROS exclusion limits. These differences could be due to EROS covering a larger solid angle than MACHO, such that the lensing rate due to lenses in the LMC itself should be smaller, as well as the use of less crowded fields by EROS as compared to MACHO, simplifying the estimation of the event durations. See Milsztajn (mil (2001)) for an extended discussion of the current observational situation.

In this paper we revisit the analysis of the MACHO collaborations candidate microlensing events (Alcock et al. 2000b ). We contemplate whether the advocated interpretation of the microlensing observations in terms of a MACHO halo dark matter component with Maxwellian velocity distribution within a standard MW halo and with best-fit typical mass $`0.6M_{}`$ and halo dark matter fraction $`f20\%`$ actually provides, in absolute terms, a good fit to the data.

## 2 Candidate events

The MACHO collaboration apply two sets of selection criteria to their data. Criteria A is the most restrictive, only accepting events with a single highly significant bump in their lightcurve as expected for simple, single lens, events, and results in 13 events. Criteria B is looser with 17 events passing the selection criteria, being designed to also accept low signal to noise events and exotic events, such as those where the lens is binary and the lightcurve has distinctive caustic features. Marginal events which are suspected of being supernovae occurring in host galaxies behind the LMC are rejected from criteria A but kept in criteria B. Further study of event number 22 222We use the MACHO collaboration’s event numbering nomenclature (Alcock et al. 2000b ) throughout. has found that the source is extended and has emission lines uncharacteristic of a stellar object and it is hence very unlikely a microlensing event (Alcock et al. 2001a ). This extraordinarily long event will henceforth be excluded from the sample, leaving 16 candidate events for criteria B.

Known stellar populations within the MW disk and bulge, as well as within the LMC disk, are expected to account for approximately one quarter of the observed events (Alcock et al. 2000b ). Follow-up observations of fields containing each of the source stars have been carried out using the Hubble Space Telescope Wide Field Planetary Camera 2 (Alcock et al. 2001b ). In the case of event number 5 (which passes both sets of selection criteria) a faint red object has been detected close to the source star. A chance super-position is extremely unlikely and a combined analysis of the HST and microlensing data finds that this object is most likely an M dwarf located in the MW disk.

It is possible to obtain more information about the properties of a given lens if either the lens or source is within a binary, as the microlensing lightcurve then exhibits additional features. The lens responsible for event 9 is a binary, producing distinctive caustic features in the lightcurve from which the source crossing radius and hence the lens projected velocity could be estimated (Bennett et al. event9 (1996); Alcock et al. 2000a ). The low projected velocity found suggests that either the lens resides in the LMC disk or the source star is also binary. Accurate photometry of event number 14 (which passes both sets of selection criteria) revealed a periodic modulation in the lightcurve. Fits to a rotating binary source have been performed which indicate that it is most likely that the lens resides in the LMC disc, though direct spectral observations of the source are required to confirm this (Alcock et al. 2001c ).

While there is clear evidence that event number 22 is not due to microlensing and that event 5 is due to a non MW halo lens, the locations of the lenses responsible for events 9 and 14 are not unambiguously known. Furthermore no conclusions on the locations of the entire lens population can be drawn from the locations of the lenses responsible for events 9 and 14, due to the biases involved in observing exotic events (Alcock et al. 2001c ; Honma honma (1999)).

A further complication is the estimation of the event durations. In crowded stellar fields it is not possible to resolve individual stars and some fraction of the baseline flux may come from unlensed stars. This phenomena, known as blending, results in an underestimate of the event timescale. The MACHO collaboration take this into account by carrying out sophisticated Monte Carlo simulations to find the average factor by which the real timescale is underestimated, and then correcting each of the fitted durations by this factor (Alcock et al. 2000b ; Alcock et al. 2001a ).

## 3 Likelihood analysis

The MACHO collaboration defines the likelihood of a given model as the product of the Poisson probability of observing $`N_{\mathrm{obs}}`$ events when expecting $`N_{\mathrm{exp}}`$ events and the probabilities of finding the observed durations $`\widehat{t}_\mathrm{j}`$ (where $`j=1,\mathrm{}.,N_{\mathrm{obs}}`$) from the theoretical duration distribution (Alcock et al. likeref (1996); Alcock et al. likeref2 (1997)):

$$=\mathrm{exp}\left(N_{\mathrm{exp}}\right)\mathrm{\Pi }_{j=1}^{N_{\mathrm{obs}}}\mu _\mathrm{j},$$
(1)
where the expected number of events is given by

$$N_{\mathrm{exp}}=E_0^{\mathrm{}}\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}t}ϵ(t)dt,$$
(2)
and $`\mu _\mathrm{j}`$ by

$$\mu _\mathrm{j}=Eϵ(t_\mathrm{j})\frac{\mathrm{d}\mathrm{\Gamma }(t_\mathrm{j})}{\mathrm{d}t}.$$
(3)
Here $`E=6.12\times 10^7`$ star-years is the exposure, $`ϵ(t)`$ is the detection efficiency, $`t`$ is the estimated event duration statistically corrected for blending and $`\mathrm{d}\mathrm{\Gamma }/\mathrm{d}t`$ is the differential event rate. The likelihood approach can be used to compare how well different sets of parameter values within a given theoretical model fit the data. In the present context these parameters are the MACHO mass and halo fraction and the model is a standard halo population of compact dark matter. However, the likelihood approach does not tell one, in absolute terms, how good a fit really is. In particular, it is possible that though ”best-fit” model parameters have been found, the model provides a generally bad description of the observational data. In order to illustrate this we compare $`\mathrm{d}\mathrm{\Gamma }/\mathrm{d}t`$ of the best fit standard halo model with delta MACHO mass function (S-DMF) with that of a best fit (i.e. maximum likelihood) Gaussian differential event rate with variable amplitude, center and width. We note that the Gaussian has no physical motivation and is intended to simply provide some arbitrary reference to the S-DMF model. 333Note that the Gaussian model has one more fitting parameter than the S-DMF model.

The best fit Gaussian and S-DMF differential event rates which arise in a likelihood analysis assuming that all 13 events found using criteria A are due to MW halo lenses (i.e. neglecting the event rate due to background events) are plotted in Fig. 1 multiplied by the detection efficiency. The Gaussian clearly gives a better fit to the observed event distribution than the S-DMF model, and its likelihood is greater by a factor of 32. We have also carried out this comparison for all possible 10 event sub-sets (i.e. subtracted every possible combination of 3 background events) and in every case the Gaussian has a larger maximum likelihood (by a factor ranging from 5 to 4500).

## 4 Event duration distribution

Fig. 1 illustrates the reason why the artificial Gaussian event duration distribution provides a better fit to the data than the event duration distribution produced by a MACHO population within a standard halo model; namely, whereas the S-DMF predicts the occurrence of long duration events none are observed. To quantify this discrepancy we have performed a Monte Carlo simulation which compares the observed width of the distribution

$$Width=\frac{\mathrm{Max}_{\{j=1,N_{obs}^{}\}}(t_j)\mathrm{Min}_{\{j=1,N_{obs}^{}\}}(t_j)}{\mathrm{Avg}_{\{j=1,N_{obs}^{}\}}(t_j)},$$
(4)
to the $`Width`$ obtained by statistically drawing a number $`N_{obs}^{halo}`$ of events from a standard halo event duration distribution (multiplied by the relevant detection efficiency and taking into account the transverse velocity of the line-of-sight). We note that the above statistic is designed to be independent of the best-fit MACHO halo fraction and mass 444There remains a residual dependence of the $`Width`$ on the typical MACHO mass, however, due to the duration dependence of the detection efficiencies. and is thus a parameterization of the observations which is essentially orthogonal to the best-fit parameters determined by the MACHO collaboration 555This approach is similar to that of de Rujula et al. (der (1995)) who showed that the duration dispersion could in principle be used to separate background events..

In our Monte Carlo analysis we find the fraction of generated event samples that yield a $`Width`$ which is smaller than that observed. To account for background (i.e. non-halo) lenses we employ the estimates for the expected number of background events $`<N_{BG}>`$ found in Alcock et al. (2000b ). We do not distinguish between the different background populations, but simply subtract a number $`N_{BG}^{}`$ from $`N_{obs}^{}`$ to obtain $`N_{obs}^{halo}`$, the number of events which are expected to be due to MACHOs and are to be generated from the theoretical event duration distribution. We do this using two different methods: subtracting a fixed number $`N_{BG}^{}=<N_{BG}>^{}`$, and subtracting $`N_{BG}^{}`$ as drawn from a Poisson distribution with average $`<N_{BG}>^{}`$666In some of our samples $`<N_{BG}>`$ as given in Alcock et al. (2000b ) is adjusted to $`<N_{BG}>^{}`$ to account for events which have been excluded from the sample by hand, cf. Table 1. We note here that while the background events could potentially account for the longest or shortest events in the sample, we nevertheless do not change the observed $`Width`$ on subtracting background events, such that our estimated probabilities, in a conservative spirit, should generally reflect an overestimate.

The results of this procedure, in particular the fraction of simulations which have smaller $`Width`$ than observed, are shown in Table 1 for a variety of different subsamples. This fraction is generally between 5% and 10%, such that the underlying S-DMF model serving as an explanation for the observations is ruled out at the 90 to 95% confidence level. Note that this conclusion does not vary strongly if the events which are suspected of being due to background populations are excluded by hand. This would not be the case if event number 22, which has a duration of 297.8 days, was a genuine microlensing event.

The discrepancy between the observations and the theoretical expectations even worsens when the assumption of a delta function MACHO mass function (DMF) is dropped. To illustrate the effects of a MACHO mass function with finite width we take the number of MACHOs with mass within a logarithmic interval to be

$$\frac{\mathrm{d}P}{\mathrm{d}\mathrm{ln}(M/M^{})}=\frac{1}{\sqrt{2\pi }}\frac{1}{\sigma }\mathrm{exp}\left\{\frac{1}{2}\left(\frac{\mathrm{ln}\left(\frac{M}{M^{}}\right)}{\sigma }\right)^2\right\},$$
(5)
with

$$\sigma =\sqrt{\mathrm{ln}\mathrm{\hspace{0.17em}2}/2},$$
(6)
such that MACHOs with either twice or one-half the average mass $`<M>=\mathrm{exp}(\sigma ^2/2)M^{}`$ (which is assumed to take the best-fit value of the DMF) are a factor of two less common. These models are indicated in the table by GMF.

It may be noted that the probabilities are generally larger when background event subtraction is performed by employing Poisson statistics, rather than by subtracting a fixed number of events. This is due to the probabilities in the Poisson case being dominated by samples which have $`N_{BG}^{}`$ larger than $`<N_{BG}>^{}`$, and have correspondingly smaller $`N_{obs}^{halo}`$, making a comparison of theoretical expectation and observed distribution more favorable. Nevertheless, such reasoning makes the implicit and nontrivial assumption that all background events have event durations falling within the range of those of a putative MACHO population.

One may ask, making this assumption about coincidence in event durations between different populations, by how much one has to increase the number of background events, in order to achieve a ”decent” probability (say 20%) of the observed distribution indeed being drawn from the distribution expected for the standard halo model. Fig. 2 shows the probability for different samples, for the delta mass function (DMF) and that given by Eq. 5 (GMF), indicating that an increase of $`N_{BG}^{}`$ by about a factor of two would be necessary in order for this to be the case. This suggests, on statistical grounds, that the contribution of the background lens populations to the observed events has been underestimated, and consequently the MACHO halo mass fraction overestimated. Nevertheless, due to the small number of observed events, a significant contribution to the observed events by a population of MACHOs within a standard MW halo may not be ruled out. This is interesting in light of the recent detection of a population of old very cool white dwarfs which could comprise a few per-cent (Oppenheimer et al. wd (2001)) of the MW dark halo, see Richer (wd1 (2001)) for a critical discussion of these observations.

## 5 Conclusions

In this paper, we have shown that the most straightforward interpretation of the existing microlensing observations of stars within the LMC, assuming the majority of candidate events to be due to compact lenses within the MW halo, may actually be in conflict with the observational data itself. The MACHO collaboration uses a likelihood analysis to derive two ”best-fit” parameters from their microlensing observations: the typical MACHO mass and MACHO halo fraction. A naive comparison of the theoretical event rate distribution produced by these ”best-fits” assuming a standard halo model, with an ad-hoc Gaussian event rate distribution indicates that the spread in the durations of the observed events seems smaller than would be expected for a standard halo of compact dark matter. In order to quantify this suspicion we have formulated a statistic, the $`Width`$ (essentially the spread in the observed event durations relative to the average event duration), that is virtually independent of the MACHO best-fit parameter values.

We have carried out Monte Carlo simulations to calculate the expected distribution of the $`Width`$, under the assumption of a standard halo MACHO population, utilizing two different background subtraction methods (fixed number and Poissonian), for both detection criteria, for delta function and finite width MACHO mass functions and excluded various sets of events suspected of being due to background populations. The probability of observing a value of the $`Width`$ as small as observed falls, in all cases, below $`10\%`$, with values often below $`5\%`$. In other words the underlying model, in terms of the expected background event rate and a standard halo MACHO population, is excluded at between the $`90\%`$ and $`95\%`$ confidence level.

The discrepancy between the observational data and its halo dark matter interpretation may be alleviated if the number of events which are due to non-halo populations is significantly larger than expected, at the expense of correspondingly reducing the MACHO halo fraction. Such a reasoning also implicitly assumes the non-halo population to produce microlensing events with event durations being by chance very similar to those of the halo population. Possibly more natural would be an explanation in terms of ”all background”, though the nature of this background, and the reason for it producing such a narrow distribution in durations, remains to be established. Alternatively, the discrepancy may indicate that the standard halo model is a poor approximation to the actual MACHO distribution (Widrow & Dubinski sim (1998)). For instance, there could be a clump with small intrinsic velocity dispersion crossing the line of sight to the LMC, which, nevertheless, would be required to be surprisingly massive. Another possible resolution to the problem could be that the statistical correction to the event durations, due to blending, leads to a reduction in the spread of the durations i.e. that the underlying durations of some observed events are significantly longer than their estimated values. If this were the case, however, then conclusions based on these estimated timescales, not only the analysis in this paper, but also the derivation of exclusion limits on the halo fraction in massive black holes (Alcock et al. 2001a ) due to the purported absence of long-duration events, would seem unreliable.

Further clarification of the nature of the lenses may hopefully result from an enlarged sample of microlensing candidates, as expected from the upcoming analysis of the MACHO 8 year data and EROS 5 year. Evidently a longer survey duration will increase the number of candidate events, reducing the ambiguity associated with Poisson statistics, as well as increasing the efficiency for detecting longer duration events, thus making the width of the distribution a significant discriminator between models. Nevertheless, a signifcant reduction of Poisson ambiguity may only result by an increase of the number of events by an order of magnitude, such as anticipated in a next generation SUPERMACHO survey (Stubbs stubbs (1998)). This later survey is expected to also reduce uncertainties in event duration corrections due to blending as it is planned to employ the technique of difference image analysis.

###### Acknowledgements.

We are grateful to David Bennett, Alain Milsztajn and Esteban Roulet for useful comments. AMG was supported by PPARC and the Swedish Research Council, and acknowledges the use of the starlink computing facilities at QMW.


format: text

text:
# A Search for 𝜈_𝜇→𝜈_𝑒 and 𝜈̄_𝜇→𝜈̄_𝑒 Oscillations at NuTeV

## Abstract

Limits on $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations are extracted using the NuTeV detector with sign-selected $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ beams. In $`\overline{\nu }_\mu `$ mode, for the case of $`\mathrm{sin}^22\alpha =1`$, $`\mathrm{\Delta }m^2>2.6`$ $`\mathrm{eV}^2`$ is excluded, and for $`\mathrm{\Delta }m^21000`$ $`\mathrm{eV}^2`$, $`\mathrm{sin}^22\alpha >1.1\times 10^3`$. The NuTeV data exclude the high $`\mathrm{\Delta }m^2`$ end of $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations parameters favored by the LSND experiment without the need to assume that the oscillation parameters for $`\nu `$ and $`\overline{\nu }`$ are the same. We present the most stringent experimental limits for $`\nu _\mu (\overline{\nu }_\mu )\nu _e(\overline{\nu }_e)`$ oscillations in the large $`\mathrm{\Delta }m^2`$ region.

Experimental evidence for oscillations among the three neutrino generations has been recently reported. For two-generation mixing, the probability that a neutrino created as type $`\nu _1`$ oscillates to type $`\nu _2`$ is:

$$P(\nu _1\nu _2)=\mathrm{sin}^22\alpha \mathrm{sin}^2\left(\frac{1.27\mathrm{\Delta }m^2L}{E_\nu }\right),$$
(1)
where $`\mathrm{\Delta }m^2`$ is the mass squared difference between the mass eigenstates in $`\mathrm{eV}^2`$, $`\alpha `$ is the mixing angle, $`E_\nu `$ is the incoming neutrino energy in $`\mathrm{GeV}`$, and $`L`$ is the distance between the points of creation and detection in km.

Data from the Super-Kamiokande atmospheric neutrino experiment  have been interpreted as evidence for $`\nu _\mu \nu _\tau `$ oscillations with $`\mathrm{sin}^22\alpha >0.88`$ and $`1.6\times 10^3<\mathrm{\Delta }m^2<4\times 10^3`$ $`\mathrm{eV}^2`$. The LSND experiment has reported  a signal consistent with $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations with $`\mathrm{sin}^22\alpha 10^2`$ and $`\mathrm{\Delta }m^2\stackrel{>}{}1`$ $`\mathrm{eV}^2`$ . The solar neutrino experiments, and most recently SNO  have reported evidence for oscillations of $`\nu _e(\nu _\mu ,\nu _\tau )`$ with $`\mathrm{\Delta }m^2<10^3\mathrm{eV}^2`$. Within a three-generation mixing scenario and under the assumption that the $`\mathrm{\Delta }m^2`$ values for $`\nu `$ and $`\overline{\nu }`$ are the same, it is not possible to simultaneously accommodate the Super-Kamiokande, LSND, and SNO results. Therefore, experimental searches for oscillations with both $`\nu `$ and $`\overline{\nu }`$ beams are of interest. In this letter, we report on a search for oscillations in both the $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ channels using a new sign-selected neutrino beam.

High-purity $`\nu `$ and $`\overline{\nu }`$ beams are provided by the new Sign-Selected Quadrupole Train (SSQT) beamline at the Fermilab Tevatron during the 1996-1997 fixed target run. Hadrons are produced when the $`800`$ GeV primary proton beam interacts in a BeO target located 1436 m upstream of the neutrino detector. Sign-selected secondary particles of specified charge (mean energy of about 250 GeV) are directed in a 221 m beamline towards a 320 m decay region, while oppositely charged (and neutral) mesons are stopped in beam dumps. Two-body decays of the focused pions yield $`\nu _\mu `$ ($`\overline{\nu }_\mu `$) with a mean energy of $``$75 GeV. Two-body decays of the focused kaons yield $`\nu _\mu `$ ($`\overline{\nu }_\mu `$) with a mean energy of $``$200 GeV. Muons are stopped in a 915 m steel/earth shield.

The energy and spatial distributions of $`\nu _\mu `$ ($`\overline{\nu }_\mu `$) CC events in the detector provide a determination of the flux of pions and kaons in the decay channel (used in the determination of the predicted $`\nu _e`$ and $`\overline{\nu }_e`$ fluxes). For $`\nu _\mu `$ running mode, the predicted energy spectra for $`\nu _\mu `$, $`\overline{\nu }_\mu `$, and ($`\nu _e`$+$`\overline{\nu }_e`$) CC events are shown in Figure 1(a). The corresponding spectra for $`\overline{\nu }_\mu `$ running mode are shown in Figure 1(b). The $`\nu _\mu `$ ($`\overline{\nu }_\mu `$)beam contains a 1.7% (1.6%) $`\nu _e`$’s($`\overline{\nu }_e`$’s) $`93\%`$ and $`70\%`$ of which are produced from $`K^\pm \pi ^0e^\pm \stackrel{_{()}}{\nu _e}`$, for $`\nu `$ and $`\overline{\nu }`$ modes, respectively. The proton beam is incident on the production target at an angle such that forward neutral kaons do not point at the detector. This greatly reduces the electron neutrino flux from neutral kaon decays (which is more difficult to model). The error in the predicted electron neutrino flux is reduced from $`4.1\%`$ (in CCFR ) to $`1.4\%`$ (NuTeV).

The NuTeV detector  is an upgrade of the CCFR detector . It consists of an 18 m long, 690 ton total absorption target-calorimeter with a mean density of $`4.2\mathrm{g}/\mathrm{cm}^3`$. Muon energy is measured by a 10 m long iron toroidal spectrometer. The target consists of 168 steel plates, each $`3\mathrm{m}\times 3\mathrm{m}\times 5.15\mathrm{cm}`$, instrumented with liquid scintillation counters placed every two steel plates and drift chambers every four plates. The separation between consecutive scintillation counters corresponds to six radiation lengths. The energy resolution of the target calorimeter is $`\mathrm{\Delta }E_h/E_h0.85/\sqrt{E_h}`$(GeV), and $`\mathrm{\Delta }E_e/E_e0.50/\sqrt{E_e}`$(GeV) for hadrons and electrons, respectively. The muon momentum resolution is $`\mathrm{\Delta }p_\mu /p_\mu =0.11`$. The NuTeV detector is calibrated continuously every accelerator cycle (once a minute) with beams of electrons, muons, and hadrons during the slow spill part of the cycle.

While the neutrinos arrived in gates a few $`msec`$ wide, the calibration beam arrived in a different gate 20 $`sec`$ long, followed by an off-spill cosmic ray gate for background measurement. These continuous test beam calibrations yield a reduction in the hadron energy scale error from $`1\%`$ (in CCFR ) to $`0.43\%`$ (in NuTeV ).

The event sample used in this analysis is similar to that used in the recent precise NuTeV measurement of the electroweak mixing angle  with additional fiducial cuts, and $`E_{\mathrm{cal}}`$ $`>`$ 30 GeV. The data sample consists of 1.5 x 106 $`\nu `$ events and 0.35 x 106 $`\overline{\nu }`$ events with a mean visible energy in the calorimeter ($`E_{\mathrm{cal}}`$) of 74 GeV and 56 GeV, respectively. The observed neutrino events are separated into CC and NC candidates. Both CC and NC interactions initiate a cascade of hadrons in the target that is registered in both the scintillation counters and drift chambers. Muon neutrino CC events are distinguished by the presence of a final state muon, which typically penetrates well beyond the hadronic shower and deposits energy in a large number of consecutive scintillation counters. NC events usually have no final state muon and deposit energy over a range of counters typical of a hadronic shower (about ten counters $``$ 1 m of steel). For each event, the length ($`L`$) is defined as the number of scintillation counters between the interaction vertex and the last counter consistent with at least single muon energy deposition. A pure sample of $`\nu _\mu N\mu ^{}X`$ $`\nu _\mu `$ charged current events is obtained from a $`\mathrm{`}`$long’ sample with $`L29`$ for $`\nu `$ running mode ($`L28`$ for $`\overline{\nu }`$). The “short’ event sample consists of events with $`L28`$ for $`\nu `$ running mode ($`L27`$ for $`\overline{\nu }`$). Events with a $`\mathrm{`}`$short’ length are primarily NC induced and originate from:

1. $`\nu _{\mu ,e}N\nu _{\mu ,e}X`$; $`\nu _{\mu ,e}`$ NC events; ($`65\%`$);
2. $`\nu _\mu N\mu ^{}X`$; $`\nu _\mu `$ short CC events with muons which range out or exit the side of the calorimeter ($`25\%`$ for $`\nu `$, $`15\%`$ for $`\overline{\nu }`$);
3. $`\nu _eNeX`$ $`\nu _e`$ CC events ($`10\%`$ for $`\nu `$, $`15\%`$ for $`\overline{\nu }`$);
4. $`\mu N\mu X`$; steep cosmic ray interactions ($`2\%`$ for $`\nu `$ and $`9\%`$ for $`\overline{\nu }`$).

The electron produced in a $`\nu _e`$ CC event (source 3) deposits energy in a few counters immediately downstream of the interaction vertex; this changes the longitudinal energy deposition profile of the shower. The energy profile is characterized by the ratio of the sum of the energy deposited in the first three scintillation counters to the total visible energy in the calorimeter $`E_{\mathrm{cal}}`$:

$$\eta _3\frac{E_1+E_2+E_3}{E_{\mathrm{cal}}},$$
(2)
where $`E_i`$ is the energy deposited in the $`i^{th}`$ scintillation counter downstream of the interaction vertex, and $`E_{\mathrm{cal}}`$ is the sum of the energy in the 20 scintillation counters downstream (plus 1 upstream) of the vertex. We similarly define $`\eta _2`$ to be the ratio of the sum of the energy deposited in the first two scintillation counters to $`E_{\mathrm{cal}}`$.

Although the electromagnetic shower component is typically much shorter than a hadron shower, NC and $`\nu _e`$ CC interactions cannot be separated on an event-by-event basis. However, the difference in the shower profiles can be used to perform a statistical extraction of the number of $`\nu _e`$ CC events in the ”$`\mathrm{`}`$short’ sample using the $`\eta _3`$ distribution.

We first assume that (for the same final state hadron energy) hadron showers produced in NC and CC interactions are the same. Any difference in the $`\eta _3`$ (or $`\eta _2`$) distributions of $`\mathrm{`}`$long’ and $`\mathrm{`}`$short’ events is attributed to the presence of $`\nu _e`$ CC interactions in the $`\mathrm{`}`$short’ sample. To compare directly the $`\mathrm{`}`$long’ and $`\mathrm{`}`$short’ events, a muon track from the data is added to the $`\mathrm{`}`$short’ events to compensate for the absence of a muon in NC events. The fraction f of $`\nu _\mu `$ CC (source 2) events with a low energy muon contained in the $`\mathrm{`}`$short’ sample which now have two muon tracks is estimated from a detailed Monte Carlo of the experiment. A simulated sample of such events is obtained by choosing $`\mathrm{`}`$long’ events with the appropriate energy distribution from the data to which a second short muon track is added in software. The length of the short track and the angular distribution are obtained from a Monte Carlo of $`\nu _\mu `$ CC events.

A sample of $`\nu _e`$ CC interactions in our detector is simulated by adding a GEANT  generated electromagnetic shower of the appropriate electron energy to the calorimetry counter energies in events in the $`\mathrm{`}`$long’ data sample. There is good agreement between the GEANT simulation of electrons and the test beam data. The energy distribution of the electron neutrinos and the fractional energy transfer $`y`$ (in each event) are generated using a detailed Monte Carlo simulation of the experiment. The CC cross section model is tuned to agree with the measured CCFR differential cross sections  for $`\nu _\mu `$ events. Since the hadron showers in the $`\mathrm{`}`$long’ sample already have a muon track, this sample ( $`\nu _eCC+\mu `$) can be compared directly with the $`\mathrm{`}`$short’ and $`\mathrm{`}`$long’ samples.

The $`\mathrm{`}`$long’ and $`\mathrm{`}`$short’ $`\eta _3`$ distributions are further corrected by subtracting the contamination due to cosmic ray events. The cosmic ray component (source 4), which is only important for very low energy bins, is well measured using the off-spill cosmic ray data. Additionally, the $`\eta _3`$ distribution of short $`\nu _\mu `$ CC events (source 2), normalized to the predicted fraction f, is subtracted from the $`\mathrm{`}`$short’ event sample.

To extract the number of $`\nu _e`$ CC events in each $`E_{\mathrm{cal}}`$ bin, we fit the corrected shape of the observed $`\eta _3`$ distribution for the $`\mathrm{`}`$short’ sample to a combination of $`\nu _\mu `$ CC and $`\nu _e`$ CC distributions with appropriate muon additions:

$$[\mathrm{short}+\mu ]=\alpha [\nu _\mu \mathrm{CC}]+\beta [\nu _\mathrm{e}\mathrm{CC}+\mu ].$$
(3)
At this point we improve  over the previous CCFR analysis  by correcting for additional effects. First, the hadron shower in CC events also includes the contribution of photons that are radiated from the muon during the CC scattering process. These photons are not present in hadron showers of NC events. A correction is applied by using the PYTHIA Monte Carlo to generate the spectrum of photons radiated by the muon in CC events. The parameters in PYTHIA that govern the emission of photons are tuned to yield agreement with the radiative corrections formalism of deRujula . Second, the procedures of adding a muon in software to the $`\mathrm{`}`$short’ sample (to model a $`\mathrm{`}`$long’ event) and of modeling electron neutrino events by adding GEANT-generated electromagnetic showers to hadron showers from $`\mathrm{`}`$long’ events are corrected for imperfect modeling using a full LEPTO/GEANT/GHEISHA simulation of the experiment and analysis procedures. In the GEANT simulation of neutrino events, the Lund Model is used to generate the initial particle composition of hadron showers. The entire experimental procedure is simulated with the beam Monte Carlo ($`\nu _e`$,$`\overline{\nu }_e`$) flux as input. Modeling corrections are extracted for each $`E_{\mathrm{cal}}`$ bin from the small difference between the extracted ($`\nu _e`$,$`\overline{\nu }_e`$) flux using simulated data and the flux extracted using perfect modeling in the Monte Carlo.

The absolute flux of $`\nu _e`$’s is taken as the average of the results from analyses done using the $`\eta _3`$ and $`\eta _2`$ variables (the statistical error from the $`\eta _3`$ analysis is used) . An additional systematic error (2.3% in neutrino mode and 0.6% in antineutrino mode) is included to account for this difference. This systematic error for antineutrinos is smaller because the final state positron carries a much larger fraction of the energy in $`\overline{\nu }_e`$ CC events.

For the oscillation search we measure the absolute flux of $`\nu _e`$’s at the detector and compare it to the flux predicted by a detailed beamline simulation . In order to extract limits on oscillations, the data are fitted by forming a $`\chi ^2`$ which incorporates the Monte Carlo generated effect of oscillations and terms with coefficients accounting for systematic uncertainties. A best fit $`\mathrm{sin}^22\alpha `$ is determined for each $`\mathrm{\Delta }m^2`$ by minimizing the $`\chi ^2`$ as a function of $`\mathrm{sin}^22\alpha `$ and these systematic coefficients. Figure 2 show the ratios of the measured rate of ($`\nu _e`$,$`\overline{\nu }_e`$) CC events to the Monte Carlo predictions minus 1. The inner errors are the statistical errors. In order to show the magnitude of the systematic errors, the outer errors include all systematic errors added in quadrature (in the analysis all correlations are taken into account). The shaded area in the figure corresponds to the systematic error band from uncertainties in the predicted electron flux (primarily from the error in the $`K^\pm `$ branching ratio) and uncertainties in the measured flux at the detector (primarily from the $`\eta _3`$,$`\eta _2`$ difference, and the $`2\%`$ error in the electron energy scale). The curves correspond to the predictions for oscillations with $`\mathrm{sin}^22\alpha =0.01`$, and $`\mathrm{\Delta }m^2`$ of $`100`$ and $`1000`$ $`\mathrm{eV}^2`$.

At all $`\mathrm{\Delta }m^2`$, the data are consistent with no observed $`\nu _\mu \nu _e`$ oscillations (i.e. the best fit values of $`\mathrm{sin}^22\alpha `$ are within one standard deviation of zero). The frequentist approach  is used to set a 90% confidence upper limit for each $`\mathrm{\Delta }m^2`$. The limit in $`\mathrm{sin}^22\alpha `$ corresponds to a shift of 1.64 units in $`\chi ^2`$ from the minimum (including all systematic errors). The 90% confidence upper limit is shown in Fig. 3(a) for $`\nu _\mu \nu _e`$. Also shown are limits from BNL-E734  and BNL-E776 . For $`\mathrm{sin}^22\alpha =1`$, $`\mathrm{\Delta }m^2>2.4`$ $`\mathrm{eV}^2`$ is excluded, and for $`\mathrm{\Delta }m^21000`$ $`\mathrm{eV}^2`$, $`\mathrm{sin}^22\alpha >1.6\times 10^3`$ (the best fit is at $`1000`$ $`\mathrm{eV}^2`$ is $`\mathrm{sin}^22\alpha =(0.4\pm 0.9\times 10^3`$). In the large $`\mathrm{\Delta }m^2`$ region, NuTeV provides improved limits for $`\nu _\mu \nu _e`$ oscillations.

Similarly, the limit for $`\overline{\nu }_\mu \overline{\nu }_e`$ is shown Fig. 3(b). Also shown are the LSND  results and limits from KARMEN . For the case of $`\mathrm{sin}^22\alpha =1`$, $`\mathrm{\Delta }m^2>2.6`$ $`\mathrm{eV}^2`$ is excluded, and for $`\mathrm{\Delta }m^21000`$ $`\mathrm{eV}^2`$, $`\mathrm{sin}^22\alpha >1.1\times 10^3`$ (the best fit is at $`1000`$ $`\mathrm{eV}^2`$ is $`\mathrm{sin}^22\alpha =(0.3\pm 1.1\times 10^3`$). In the $`\overline{\nu }_\mu `$ mode, our results exclude the high $`\mathrm{\Delta }m^2`$ end of $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations parameters favored by the LSND experiment, without the need to assume that the oscillation parameters for $`\nu `$ and $`\overline{\nu }`$ are the same. These are the most stringent experimental limits  for $`\nu _\mu (\overline{\nu }_\mu )\nu _e(\overline{\nu }_e)`$ oscillations in the large $`\mathrm{\Delta }m^2`$ region.

This work was supported by the U.S. Department of Energy, the National Science Foundation, and the Alfred P. Sloan foundation.


text:
# 1 Introduction

## 1 Introduction

Large detectors of high energy ($`>1\mathrm{GeV}`$) neutrinos, like MACRO in the Gran Sasso Laboratory or SuperKamiokande  have been operated in the last decade of the past century, but no extrasolar astrophysical source has been detected. The occasion of this paper is related to the decommissioning of the MACRO detector at Gran Sasso, that is the largest neutrino telescope operated before the year 2k. The proposal for the detector was accepted in November 1984 . Physics runs started in March 1989, as soon as the first supermodule of the detector has been built  The full detector  has been run from April, 1994 to December 31, 2000 when it has been decommissioned. A sample of a 1000 up-ward going muons has allowed the measurement of the flux of geophysical neutrinos collected up to September, 1999 , producing also interesting result . The purpose of this paper is to compare the upper limits obtained by MACRO with current models of neutrino emission from astrophysical objects. First in the following §2 we will discuss the physics of neutrino detection in detectors of up-ward going muons. In particular in §2.1 we discuss the extrapolation of the differential neutrino cross section to the range of energies of interest for astrophysical neutrino detection, in §2.2 the conversion of neutrinos into muons in the rock and finally in §2.3 we give the rationale for a new statistical method for assessing a model independent upper limit to the differential flux of neutrinos inferred from the measured up-ward going muon flux. In §3 we present some results obtained from the MACRO survey of neutrino induced up-ward going muons. In §3.1 we present the constraint to the neutrino emission from the Crab nebula, in §3.2 that for a possible neutralino decay in the Galactic center, in §3.3 the constraint to diffuse neutrino background which can be inferred from the cumulative search of coincidences with blazars and in §3.4 we discuss the same topic for the GRB deriving an upper limit to the diffuse flux of neutrinos correlated with GRB . Finally in §4 we discuss the prospect for detection of neutrinos in the next future by the planned detectors.

## 2 Neutrino detection

### 2.1 Extrapolation of neutrino cross section

The detection process for muon neutrinos and anti neutrinos is the Charged Current (CC) scattering illustrated in Fig. 1

with double differential cross section

$`{\displaystyle \frac{d^2\sigma _{\nu (\overline{\nu })}}{dxdy}}={\displaystyle \frac{G_F^2(sm_N^2)}{2\pi }}{\displaystyle \frac{M_W^4}{(Q^2+M_W^2)^2}}`$ (1)
$`\times \{{\displaystyle \frac{y^2}{2}}\mathrm{\hspace{0.17em}2}xF_1^{\nu (\overline{\nu })}(x,Q^2)+(1y)F_2^{\nu (\overline{\nu })}(x,Q^2)+\pm (y{\displaystyle \frac{y^2}{2}})xF_3^{\nu (\overline{\nu })}(x,Q^2)\}`$
where $`y=1E_\mu /E_\nu `$ is the fraction of neutrino energy lost in the Lab system, $`x=Q^2/2m_N(E_\nu E_\mu )`$ the fraction of parton momentum carried by the struck parton, and $`Q^2=xy(sm_N^2)`$ the square of the transferred four momentum. The last term is positive for $`\nu `$’s and negative for $`\overline{\nu }`$’s. The three functions $`F_1,F_2`$ and $`F_3`$ are the structure functions (SF’s) of the nucleon. In the framework of the quark parton model (QPM) the parton are identified with quarks, and the SF’s can be expressed in terms of the probability density functions (PDF) of each quark types . The PDF’s are the functions $`q(x,Q^2)`$ where $`q=u,d,s,c,b,t`$ and the relative anti-quarks, which assign the probability that a quark carries a momentum fraction $`x`$. In particular we have for an isoscalar target

$`F_2^{\nu N(\overline{\nu }N)}`$ $`=`$ $`x\left\{u+\overline{u}+d+d\overline{d}+2s+2c+2b+2t\right\}`$
$`xF_3^{\nu N(\overline{\nu }N)}`$ $`=`$ $`x\left\{u\overline{u}+dd\overline{d}\pm 2s2c\pm 2b2t\right\}`$
At very high energies the SF’s $`F_2`$ and $`2xF_1`$, are equal, but if the transverse momentum of the $`W`$ is not negligible, they are related by the equation

$$2xF_1(x,Q^2)=F_2(x,Q^2)\frac{1+4m_N^2x^2/Q^2}{1+R(x,Q^2)}$$
(2)
where $`R(x,Q^2)=\sigma _L/\sigma _T`$ is the ratio of the longitudinal to transversal cross section of the $`W`$.

The structure functions $`F_2`$ and $`xF_3`$ and the function $`R`$ have been measured experimentally up to $`E_\nu 500\mathrm{GeV}`$ by the CCFR/NuTeV collaboration . To calculate the neutrino cross sections at higher energies, as we are interested in this paper, the SF’s should be extrapolated. The propagator term in Eq. (1) $`M_W^4/(Q^2+M_W^2)^20`$ for $`Q^2M_W^2`$, therefore the cross section at any energy will vanish for $`Q^2\mathrm{max}[xy(sm_N^2),M_W^2]`$. However the double differential cross section of Eq (1) must be integrated to obtain the differential cross section $`d\sigma /dy`$ for a given $`y`$ in the interval $`Q^2/y(sm_N^2)x1`$. If the minimum detectable energy for the muon in the detector is $`E_{th}`$ the lower integration limit is $`x_{min}Q^2/2m_NE_\nu `$. In Fig. 2 we have reported the area of the plane $`\mathrm{log}_{10}(1/x),Q^2`$ corresponding to $`E_\nu =10^6\mathrm{GeV}`$ where Eq. (1) must be integrated, compared with the area where the SF’s have been measured.

In the same plot of Fig. 2, we have reported the lines $`W^2=Q^2(1/x1)=4m_h^2`$ where $`m_h`$ are the masses of the heavy quarks $`c,b`$ and $`t`$. Above these lines the production of the heavy quarks will give an important contribution to the cross section. In fact it is well established that above $`100\mathrm{GeV}`$ the contribution of the quark production to the total cross section is of the order of 20% . From the plot of Fig. 2 we are forced to conclude that the experimental knowledge of the physics underlying the extrapolation of the cross section is, up to now, rather incomplete. In the frame work of quantum cromodynamics (QCD), the gauge field theory which describes the strong interactions of colored quarks and gluons, the cross sections can be calculated perturbatively in ascending order of the strong coupling constant $`\alpha _s`$. In practice this corresponds to the solution of a set of partial differential equations $`F_i(x,Q^2)/\mathrm{ln}(Q^2)`$ (the DGLAP equations) at fixed $`x`$, starting from the knowledge of the SF’s $`F_i(x,Q^2)`$ with $`i=1,2,3`$ at a given lower value $`Q_0^2`$ (see e.g. Sterman G. et al. 1999 and references therein). This method has been checked, using the measured differential cross sections. For example we report in Fig. 3 the fit to the measured differential cross section obtained using PDF’s calculated by the CTEQ5 collaboration , with evolution equations resulting from NLO analysis from a starting scale $`Q_0^2=1\mathrm{GeV}^2`$. It is evident from this figure that the agreement is reasonable but not perfect .

Nevertheless the extrapolation at higher energies of the neutrino cross section is still subjected to a non negligible uncertainty because, as is shown in Fig. 2, the domain of integration at higher energies extends to smaller and smaller values of $`x`$. Therefore, even if the evolution of SF’s at fixed $`x`$ to higher $`Q^2`$ can be calculated with the DGLAP equations, still we lack of information about the SF’s at low x values. Various estimates, based on different physical assumption , as shown in Fig. 4 disagree by factors up to 2-4 for $`E_\nu 10,000\mathrm{GeV}`$. The experimental point has been obtained by the H1 Collaboration .

The higher values of the cross section is given by extrapolation of the SF’s for $`x<10^5`$ based on the Regge theory .

In the same Fig. 4 we have reported the cross section of the reaction $`e^{}p\nu X`$ by the H1 experiment at HERA (Ahmed et al. 1994) which seems to be in better agreement with the lower cross section estimate. In the rest of this paper we will take as a baseline extrapolation of the cross section the more conservative one, shown as a solid line in Fig. 4, based on the CTEQ5 PDF’s.

### 2.2 Evaluation of conversion probability

The probability $`P_{\mu \nu }(E_\nu )`$ of detecting a muon with energy $`E_{th}`$ in the detector, originated by a neutrino with energy $`E_\nu `$, is given, in the continuous slowing down approximation, by the integral

$$P_{\mu \nu }(E_\nu )=\frac{1}{m_N}_0^{1E_{th}/E_\nu }\left[a+b(1y)E_\nu \right]^1_0^1\frac{d^2\sigma _\nu }{dxdy}𝑑y𝑑x$$
(3)
where $`m_N`$ is the mass of the nucleon, $`a`$ and $`b`$ are the coefficient of the muon energy loss $`dE_\mu /dz=a(E_\mu )+b(E_\mu )E_\mu `$ .

### 2.3 Model independent upper limits

The differential up-ward going muon flux in the detector will be $`d\mathrm{\Phi }_\mu /dE_\nu =P_{\mu \nu }d\mathrm{\Phi }_\nu /dE_\nu `$ where $`d\mathrm{\Phi }_\nu /dE_\nu `$ is the spectral distribution of the neutrino flux. However what is measured in an up-ward going muons neutrino telescope, like MACRO, is only the integrated muon flux

$$\mathrm{\Phi }_\mu ^{obs}(>E_{th})=_{E_{th}}^{\mathrm{}}P_{\mu \nu }\frac{d\mathrm{\Phi }_\nu }{dE_\nu }𝑑E_\nu $$
(4)
The conversion of the muon flux into a neutrino flux is subjected to the ambiguities arising from the lack of knowledge of the shape of the neutrino spectrum itself. The usual approach followed in the literature  is to give the upper limit to the neutrino fluxes assuming a spectrum $`d\mathrm{\Phi }_\nu /dE_\nu E_\nu ^2`$. In this case the integrated neutrino flux can be estimated using an average conversion probability

$$\mathrm{\Phi }_\nu (>E_{th})=\frac{\mathrm{\Phi }_\mu ^{obs}(>E_{th})}{\frac{1}{E_{th}}_{E_{th}}^{\mathrm{}}E_\nu ^2P_{\mu \nu }𝑑E_\nu }$$
(5)
We propose here to apply the maximum entropy principle  to infer the unknown neutrino differential flux, with the less prior assumption on the spectral shape. In fact, even if the maximum entropy principle has been questioned for the lack of solid theoretical foundations, it gives very significant results in several fields of applications of spectral deconvolution . According to this principle we assume that the probability distribution function of the energy of the neutrino that could have produced a muon in the MACRO detector will be uniform. In this case we have the ansatz $`P_{\mu \nu }d\mathrm{\Phi }_\nu ^{ME}/dE_\nu =\mathrm{const}`$. Therefore, normalizing to the integrated neutrino flux of Eq. (5), we have derived the estimate of the upper limit to the neutrino flux which will be shown in the following.

## 3 Results of the MACRO survey

### 3.1 Neutrino emission from the Crab nebula

Gamma ray spectrum of the Crab Nebula extends up to 25 TeV. Above $`500\mathrm{GeV}`$ the spectrum is a simple power law with $`\alpha 2.5`$. The total $`\gamma `$-ray luminosity is above $`500\mathrm{GeV}`$ $`L_\gamma 1.5\times 10^{26}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. The detection of neutrino emission from the Crab could support the possibility that this high energy emission is due to hadronic interactions. Bednarek & Protheroe (1997) have analyzed the consequences of acceleration of heavy nuclei in the pulsar magnetosphere as a possible mechanism of energetic radiation from the Crab Nebula. The MACRO limit for the Crab Nebula is

$$\mathrm{\Phi }_\mu ^{Crab}(>1.2\mathrm{GeV})2.5\times 10^{14}\mathrm{cm}^2\mathrm{s}^1(90\%\mathrm{C}.\mathrm{L}.)$$
(6)
In Fig. 7 we compare the limit on the neutrino flux, with two extreme models. The present emission is determined by the trapped high energy proton’s density in the nebula. Therefore the present emission is mainly due to acceleration in the early times after supernova explosion. Magnetic dipole losses determine the pulsar period at present, but the initial period is determined largely by gravitational losses and it has been probably shorter than 10 ms. Hence the two models consider two initial periods: 5 ms, and 10 ms.

### 3.2 Neutralino annihilation in the Galactic Center

Intense emission from the Galactic center is predicted if cold dark matter is present there, as in current models of the dark galactic halo . In particular it could be expected, according to Gondolo & Silk (1999) that a massive black hole at the galactic center could redistribute the WIMP’s into a cusp. The effect of this redistribution would be a strong enhancements of the self-annihilation rate of neutralinos. In Fig. 8 we report the predicted up-ward going muon flux in the two cases.

After the complete run MACRO  the 90% C.L. upper limits on up-ward going muons from the Galactic Center is $`\mathrm{\Phi }_\mu 0.34\times 10^{14}\mathrm{cm}^2\mathrm{s}^1`$. This new limit is reported in Fig. 8. It is clear from this figure that the predicted enhancement has not been observed. One possible explanation of this negative result has been given from a subsequent more careful evaluation, in a recent paper, of the dynamics of the cusp formation . The MACRO negative observation supports the possibility that the accreting matter could not spiral fast enough by dynamical friction only. Thus within a Hubble time only a mild enhancement could have taken place.

### 3.3 Diffuse emission from Blazars

The second EGRET catalog of high-energy $`\gamma `$-ray sources  contains 40 high confidence identification of AGN and all appear to be blazars. This is the main reason for which several authors  have proposed this type of object as potential powerful sources of high energy neutrinos. The proton blazar model  both protons and electrons are accelerated and protons interact with synchrotron radiation produced by electrons photoproducing pions that decay into $`\gamma `$-rays and neutrinos.

The up-ward going muons distribution, shown in Fig. 6 has been searched for angular coincidences in a 3 cone with each of the 181 blazars, listed in a recent catalogue , with a declination $`\delta 40^{}`$. No significant excess has been found for any individual source, with a cumulative average upper limit to the muon flux $`\mathrm{\Phi }_\mu (1.2\mathrm{GeV})5.44\times 10^{16}\mathrm{cm}^2\mathrm{s}^1`$. A bias in favour of northern declinations is clearly present in the sample. In fact 178 out of the 233 objects listed in the catalog, have $`\delta >0^{}`$ and there are no known BL Lacs with $`\delta <57^{}`$. We can estimate, very roughly, that the MACRO survey has covered the declination band $`57^{}\delta 40^{}`$, corresponding to $`3.2\pi \mathrm{sr}`$ of sky. From this figure we can estimate the U.L. to the diffuse up-ward going muons from blazars

$$\frac{d\mathrm{\Phi }_\mu ^{BLac}}{d\mathrm{\Omega }}0.98\times 10^{14}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1(90\%\mathrm{C}.\mathrm{L}.)$$
(7)
Applying the maximum entropy method we obtain the upper limit to the diffuse neutrino flux reported in Fig. 9. In the same plot we have reported the model independent upper bound that can be derived if the UHE cosmic rays are also accelerated in the same source . The higher prediction is obtained, presuming a blazar’s luminosity cosmological evolution following that of QSO’s, which may be described as $`f_{QSO}=(1+z)^3`$ at redshift $`z<1.9`$ and constant above. The lower prediction is for no evolution of the blazar’s density.

### 3.4 Search for neutrinos from GRB’s

Several mechanisms for production of intense HE and UHE neutrino and $`\gamma `$ burst associated with the main pulse, in the $`0.1\mathrm{MeV}`$ have been proposed in the literature . It is somewhat surprising however that a phenomenon that was discovered in the keV-MeV region could be the herald of big activity in the multi-TeraVolt region.

We have searched the catalogues BATSE 3B & 4B  containing 2527 gamma ray bursts from 21 Apr. 1991 to 5 Oct. 1999 . They overlap in time with 1085 upward-going muons collected by MACRO during this period. In Fig. 10 it is shown a plot of angular distance among the GRB and up-ward going muons in MACRO as a function of arrival time difference.

We observe that if neutrino are massive, as suggested by observation of neutrino oscillations, the neutrino burst from cosmological GRB will be delayed respect to the observation of the $`\gamma `$ rays observed by the GRB itself. The time delay will be

$$\delta t=\frac{1\beta }{\beta }(t_zt_0)$$
(8)
where $`\beta =\sqrt{E_\nu ^2m_\nu ^2c^4}/E_\nu `$. For $`E_\nu m_\nu c^2`$ we can put $`(1\beta )/\beta =1/2(E_\nu /m_\nu c^2)^2+𝒪(E_\nu /m_\nu c^2)^4`$.

In Fig. 11 we show the expected delay of the massive neutrino signal, for GRB in the range $`z=0.52`$ as a function of the Lorentz factor of the neutrino. In the same figure we have reported the normalized Lorentz factor distribution (for a source spectrum $`dN_\nu /dE_\nu E_\nu ^2`$) of neutrinos that produce upward going muons in the MACRO detector. From this figure we see that a maximum delay of $`1000\mathrm{s}`$ could be expected for very massive neutrinos ($`m_\nu =10\mathrm{eV}`$).

Therefore the time window to be searched for should to be expanded up to $`1000`$ s after the detection of the GRB for the muon to be detected in MACRO. However according to the current ideas on the production of the GRB explosion, it should not be possible to have the neutrino emission long before the GRB. Therefore we can assume that the coincidences between muons detected in the time window preceding the GRB should be very likely accidental. From the lower panel of Fig. 10 we estimate that the accidental coincidences among up-going muons and GRB with $`\mathrm{\Delta }\theta 20^{}`$ will be $`0.4`$ in the overlap period. From the upper panel of the same figure we could candidate two events as possible true coincidences, one event after 39.4 s from 4B950922 at an angular distance of 17.6 and another very horizontal event in coincidence with the 4B940527 inside 280 s at 14.9. However the cumulative Poissonian probability that those two events are accidental is $`7.9\times 10^3`$. We assume then that those events are expected accidental, and calculate that the 90% confidence interval for the unobserved true coincident upward-going muons  will be $`5.41`$. The corresponding upper limit to the muon fluence from the average GRB is $`\mathrm{\Phi }_\mu (1.2\mathrm{GeV})0.79\times 10^9\mathrm{cm}^2(90\%\mathrm{C}.\mathrm{L}.)`$ for neutrino masses $`m_\nu 10\mathrm{eV}`$. From this cumulative limit we can derive an upper limit to the diffuse neutrino background from GRB’s, taking into account the fact that the MACRO live time in the period from 21 Apr. 1991 to 5 Oct. 1999 has been $`4.62\mathrm{y}`$, that neutrinos could have been detected for GRB’s with declination $`\delta 40^{}`$ (corresponding to $`3.53\pi \mathrm{sr}`$ of sky) and that the exposure factor for BATSE in this declination range  is $`55\%`$ we can convert the upper limit for the average burst to a diffuse background

$$\frac{d\mathrm{\Phi }_\mu ^{GRB}}{d\mathrm{\Omega }}2.25\times 10^{15}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1(90\%\mathrm{C}.\mathrm{L}.)$$
(9)
We obtain also in this case the upper limit to the diffuse neutrino emission from GRB’s reported in Fig. 12. In the same figure we have reported the prediction of diffuse neutrino background produced by a cosmological distribution of GRB’s . Also we have reported the flux expected from the intriguing possibility that before the jet could emerge from the dense stellar progenitor of the GRB, the accelerated protons could interact with the thermal X-ray photon leading to neutrinos with energy $`5\mathrm{TeV}`$ . In this model the intensity of diffuse neutrino background depends from the amount of collimation of the jet $`\mathrm{\Omega }_{jet}/4\pi `$, because the $`1000`$ GRB’s observed per year, are only a small fraction of the actual explosions.

## 4 Conclusions

In Fig. 13 we show a gallery of past, present and future HE neutrino telescopes, compared with a massive detector with 1 km3 of volume. From this figure we have a direct impression of the smallness of the detectors which were operated in the past century, compared with the ones under construction. It is also worth noticing that the preliminary results obtained by Amanda, from April to October, 1997 which are reported in the upper panel of Fig. 6, include 1097 up-ward going muons. A number comparable to the number of events detected by MACRO in $`4.6`$ live years. This shows directly that the acceptance of the AMANDA-II detector is already more then a factor 10 greater then the acceptance of MACRO.

In the previous §3 we have practically shown how far was MACRO from effectively testing theoretical production models. From Fig 7 it can be crudely estimated that an acceptance $`10^4\mathrm{MACRO}1\mathrm{km}^2`$ will be needed to challenge possibly the predictions for galactic neutrino sources. The situation could be slightly better for the statistical detection of emission from selected classes of objects, like for example the association of up-ward going muons with blazars, where a possible detection seems to be attainable with an acceptance $`0.1\mathrm{km}^2`$.


text:
# Searching for the statistically equilibrated systems formed in heavy ion collisions

## I Introduction

It is known from more than 15 years that violent heavy ion collisions lead to an advanced disassembly of the compound system known under the name of nuclear fragmentation. The asymptotically resulted fragments are the only experimental link with the primordial process. Good agreements between various observables related to these fragments and results of various models assuming uniform population of the system phase space Gross1 ; Bondorf ; Randrup  led to the conclusion that a statistically equilibrated nuclear source is at the origin of the fragmentation process. The source size, its excitation energy and its volume are thus quantities which can only be indirectly evaluated by comparisons between experimental data and statistical multifragmentation models predictions. The comparison process is complicated by various quantities polluting the pure statistical signals such as preequilibrium particle emission, collective radial expansion, Coulomb propagation of the break-up fragments, secondary particle emissions. A complete theoretical reproduction of experimental fragmentation event topology (never achieved so far) would help to answer questions fundamental for determining the thermodynamics taking place in such systems. For example, the identification of the experimental freeze-out volume would allow to locate the process in the phase diagram and, therefore, would help solving a problem which made the subject of a strong debate over the last decade: is there any phase transition taking place in nuclear matter?

The aim of the present paper is to investigate whether an unambiguous determination of the experimental statistically equilibrated source can be achieved by means of the microcanonical model from Ref. noi1 ; temp3 . To this aim, the model is refined by including the conservation of the angular momentum and completed with a Coulomb propagation stage of the primary fragments with the possibility of superimposing radial flow. A complete set of experimental fragmentation data concerning the reactions Xe+Sn at 32 MeV/u and Gd+U at 36 MeV/u recently published in Ref. indra1 ; indra2  is used for comparison with the model results. The paper is structured as follows: Section II gives a brief presentation of the employed microcanonical model and describes its new refinements and improvements. In Section III a comparison between the model prediction and the above mentioned experimental data is presented. The influence of various model parameters on different observables is discussed in Section IV. Conclusions are drawn in Section V.

## II The model: brief review and improvements

The microcanonical multifragmentation model used in the present study was proposed in Ref. noi1 ; temp3 . The modifications from the previous version of the model concern the primary break-up stage, momentum generation, inclusion of collective radial expansion, Coulomb propagation, and the secondary emission stage.
i) Break-up stage
This stage concerns the disassembly of a statistically equilibrated nuclear source $`(A,Z,E,V)`$ (the mass number, the atomic number, the excitation energy and the freeze-out volume respectively). The basic assumption of the model is equal probability between all configurations $`C:\{A_i,Z_i,ϵ_i,𝐫_i,𝐩_i,i=1,\mathrm{},N\}`$ (the mass number, the atomic number, the excitation energy, the position and the momentum of each fragment $`i`$ of the configuration $`C`$, composed of $`N`$ fragments). Fragments are assumed to be spherical, are not allowed to overlap each other and are placed into a spherical recipient of volume $`V`$. In the previous version of the model noi1  the system was considered subject to the standard microcanonical constraints: $`_iA_i=A`$, $`_iZ_i=Z`$, $`_i𝐩_i=𝐏`$ (=0 in the c.m. frame), $`E`$ \- constant. An extra constraint will be considered herein: the conservation of the angular momentum, $`_i𝐫_i\times 𝐩_i=𝐋`$. The integration over fragment momenta can be analytically performed subject to the above mentioned constraints:

$`{\displaystyle \underset{i=1}{\overset{N}{}}\mathrm{d}𝐩_i\delta \left(HE\right)\delta \left(\underset{i}{}𝐩_i𝐏\right)\delta \left(\underset{i}{}𝐫_i\times 𝐩_i𝐋\right)}`$ (1)
$`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{\Gamma }\left(\frac{3}{2}(N2)\right)}}\left({\displaystyle \frac{_im_i}{_im_i}}\right)^{3/2}{\displaystyle \frac{1}{\sqrt{det𝐈}}}\left[2\pi \left(K{\displaystyle \frac{P^2}{2M}}{\displaystyle \frac{1}{2}}𝐋^T𝐈^1𝐋\right)\right]^{\frac{3}{2}(N2)1}`$
Here $`M=_im_i`$, $`𝐈`$ is the inertial tensor of the system: $`I_{\alpha \beta }=_im_i\left(r_i^2\delta _{\alpha \beta }r_i^\alpha r_i^\beta \right)`$ with $`\alpha ,\beta =1,2,3`$, $`𝐋^T𝐈^1𝐋=_{\alpha \beta =1}^3L_\alpha I_{\alpha \beta }^1L_\beta `$, $`H=_ip_i^2/(2m_i)+_{i<j}V_{ij}+_iϵ_i_iB_i`$ and $`K=E_{i<j}V_{ij}_iϵ_i+_iB_i`$ ($`V_{ij}`$ stands for the Coulomb interaction between fragments $`i`$ and $`j`$). Of course, the model corresponds to the c.m. frame (i.e. $`𝐫_i`$ with $`i=1,\mathrm{},N`$ are positions in the c.m. frame) where $`𝐏=0`$ 111For each particular fragmentation event the origin of the system is considered in its center of mass and not in the center of the spherical “recipient” frame. Therefore, at each step of the Metropolis type simulation one has to consider the transformation $`𝐫_i=𝐫_0^i𝐫_{\mathrm{c}.\mathrm{m}.}`$, $`i=1,\mathrm{},N`$. Here $`𝐫_0^i`$, $`i=1,\mathrm{},N`$ represent the positions of the fragments subject to the origin of the spherical recipient, generated as in noi1 .. If one also imposes $`𝐋=0`$ (which is the hypothesis employed in the present calculations) then the last factor in eq. (1) becomes $`(2\pi K)^{3/2N2}`$. The only modifications brought by the inclusion of the angular momentum conservation constraint to the statistical weight of a configuration $`C^{}:\{A_i,Z_i,ϵ_i,𝐫_i,i=1,\mathrm{},N\}`$ ($`W_C^{}`$) from Ref. noi1  are the inclusion of the factor $`1/\sqrt{det𝐈}`$ and the replacement of $`N`$ with $`N1`$. These new weights can be employed in a Metropolis-type simulation which allows the determination of the average value of any system observable in the very same manner as in Ref. noi1 .

Fragments with $`A4`$ are considered without excitation degrees of freedom except for the $`\alpha `$ particle for which the few levels larger than 20 MeV with $`\mathrm{\Gamma }2.01`$ MeV have been considered. These fragments are weighted in $`W_C^{}`$ with their energy levels degeneracies. Larger fragments carry internal excitation. For them, the following level density formula is included in the statistical weight of a configuration $`C^{}`$ (see Ref. temp3 ):

$$\rho (ϵ)=\frac{\sqrt{\pi }}{12a^{1/4}ϵ^{5/4}}\mathrm{exp}(2\sqrt{aϵ})\mathrm{exp}(ϵ/\tau )$$
(2)
with $`a=0.114A+0.098A^{2/3}`$ MeV-1 iljinov  and $`\tau =9`$ MeV. The factor $`\mathrm{exp}(ϵ/\tau )`$ is introduced to account for the dramatic decrease of the excited levels lifetime at high excitation energies Randrup .
ii) Generation of the primary decay fragment momenta
While integration over the fragments’ momenta has been carried out in order to simplify the Metropolis simulation, to produce events similar to the experimental ones one has to generate momenta for each given configuration $`C^{}`$. Each of these events is characterized by a kinetic energy $`K`$. Therefore, we have to deal with the following computational task: generate uniformly momenta for a system composed of $`N`$ particles such that $`_ip_i^2/2m_i=K`$, $`_i𝐩_i=0`$ and $`_i𝐫_i\times 𝐩_i=𝐋`$. This problem was solved in an elegant way in Ref. randrup . There it is proved that the following generation gives the right sampling of the system momenta, in agreement with the above mentioned constraints: Pick a preliminary set of $`N`$ particle momenta from a arbitrary canonical distribution, then eliminate the overall translational and rotational motion by making the transformation $`𝐩_i𝐩_im_i(𝐏^{}/M+𝝎^{\mathbf{}}\times 𝐫_i)`$ (where $`𝐏^{}=_i𝐩_i`$ and $`𝝎^{\mathbf{}}=_i𝐫_i\times 𝐩_i𝐈^1`$), then spin the system such as to match the desired angular momentum $`𝐋`$: $`𝐩_i𝐩_i+m_i𝝎\times 𝐫_i)`$ with $`𝝎=𝐋^T𝐈^1`$ and finally, renormalize the momenta such as to match the available energy $`K\frac{1}{2}𝝎𝐋`$.
iii) Radial flow
After generating the fragment momenta corresponding to the primary decay, inclusion of nonequilibrium effects such as collective radial expansion (flow) can be easily superimposed. Here we use the following parameterization for the flow velocity of fragment $`i`$: $`𝐯_f^i=𝐯_0(r_i/R)^\alpha `$, with $`𝐯_0=v_0(𝐫_i/r_i)`$ and $`\alpha `$ a real number defining the flow profile. After including the radial flow the momentum of the $`i`$th fragment reads: $`𝐩_i𝐩_i+m_i𝐯_f^i`$.
iv) Coulomb propagation
After break-up, hot primary fragments are supposed to suffer an expansive propagation under their mutual Coulomb interaction. This stage can be easily simulated by integrating the corresponding Newtonian equations of motion. Integration has been carried out up to 500 fm/c, when Coulomb interaction between fragments can be neglected.
v) Secondary decays
Since primary fragments carry internal excitation, a secondary decay stage was introduced in Ref. temp3 . Depending on the fragment excitation, secondary break-up processes or particle evaporation were considered. Since the above classification is rather arbitrary, here we resume to the second process, treated using the standard Weisskopf evaporation scheme. As in temp3 , the range of the evaporated particles is considered up to $`A=16`$. Evaporation events are simulated using standard Monte Carlo temp3 .

For each primary break-up event (i), the steps ii) - v) are performed. The resulted fragmentation events can be readily compared with experimental ones (after removing the particles coming from the preequilibrium stage).

## III Comparison of the model predictions with recent experimental data

In this section comparison between the model results and the recent experimental data concerning the reactions Xe+Sn at 32 MeV/u and Gd+U at 36 MeV/u are presented. For the above reactions, various experimental fragment size and kinematic observables have been published recently indra1 ; indra2 .

The aim of the study is to extract information about the physical phenomenon taking place by trying to fit the entire set of experimental data using the microcanonical multifragmentation model. It is known vfree  that various freeze-out assumptions may induce important differences in the statistical models’ results. For drastically reducing this uncertainty in the present study we employ two opposite freeze-out scenarios:
1) The standard working hypothesis of the model: Fragments are idealized as hard spheres placed into a spherical freeze-out recipient; fragments are not allowed to overlap each other or the recipient wall. The generation of the fragments’ positions is described in Ref. noi1 .
2) The hardcore interaction is switched off: Integration over the fragments’ positions may be approximately carried out by assuming as in Ref. fai-ran  that each fragment (from a configuration composed of $`N`$ fragments) is blocking the volume $`V_0/N`$ ($`V_0`$ is the volume of the nuclear system at normal density) for the rest of the fragments (as in fai-ran ) and for itself as well. Coulomb energy is approximated by the Wigner Seitz formula, being therefore independent of the fragments’ positions 222For a $`N`$ fragment partition the Wigner Seitz Coulomb interaction energy writes: $`V_{WS}=\frac{3}{5}\left(\frac{Z^2}{R}_{i=1}^N\frac{Z_i^2}{R_C^i}\right)e^2`$ where $`R`$ is the radius of the freeze-out recipient and $`R_C^i`$ is the radius of the Wigner Seitz cell corresponding to the fragment $`i`$, $`R_C^i=r_0(nA_i)^{1/3}`$ where $`r_0=1.2`$ fm, and $`n=V/V_0`$.. The integration over the fragments’ positions writes:

$$\underset{i=1}{\overset{N}{}}\mathrm{d}𝐫_i\frac{1}{\sqrt{det𝐈}}\underset{j<i}{}\theta _{ij}V_{free}\underset{i=1}{\overset{N}{}}\mathrm{d}𝐫_i\frac{1}{\sqrt{det𝐈}}$$
(3)
with:

$$V_{free}=\underset{i=1}{\overset{N}{}}\left(Vi\frac{V_0}{N}\right).$$
(4)
The factor $`_{j<i}\theta _{ij}`$ from eq. (3) is just a formal expression of the fragment blocking constraint formulated before. With this approximation fragment positions have to be generated into the spherical freeze-out volume, without any hardcore constraint only for evaluating $`det𝐈`$ and the statistical weight of a configuration will get (as in Ref. vfree ) the extra factor $`V_{free}`$. Having the weights of each configuration $`C^{}`$ the corresponding Metropolis type simulation is straight forward noi1 .

Using the above defined working hypotheses, comparisons between the model’s results and the corresponding experimental data have been performed. Both fragment size distributions and kinematic distributions have been considered. Very good agreement between calculations and experimental data can be observed for both reactions, in both considered working hypotheses for all the considered observables.

The comparisons between theoretical and experimental fragment size distributions corresponding to both Xe+Sn at 32 MeV/u and Gd+U at 36 MeV/u reactions are as follows: Charge multiplicities are given in Fig. 1. Multiplicities of largest charge ($`Z_{max1}`$), second largest charge ($`Z_{max2}`$) and third largest charge ($`Z_{max3}`$) from one fragmentation event are presented in Fig. 2.

In Fig. 3 intermediate mass fragment (IMF) distributions and $`Z_{bound3}`$ distributions are given (IMF fragments are considered to have $`Z5`$ as in indra2 ; $`Z_{bound3}`$ is defined as the sum over all fragments’ atomic numbers from one fragmentation event which are greater or equal to 3; $`Z_{bound5}`$ is similarly defined but the $`Z`$ limit is fixed to 5). $`Z_{bound5}`$ multiplicities are compared with experimental data in Fig. 4 only for the reaction Gd+U at 36 MeV/u. It can be noticed that for the $`Z_{bound}`$ distributions the calculated width of the distributions are slightly smaller than the experimental ones. This suggests a (small) fluctuation of the experimental equilibrated source.

This will become more clear in the next section. An excellent agreement for all considered observables related to fragment size was, therefore, obtained.

Comparisons between theoretical and experimental average fragment kinetic energy function of fragment charge are given in Fig. 5, left column for the the two considered reactions. Very good agreements are to be noticed. A deeper insight to the kinematics of the process can be obtained by analyzing the reduced fragments velocity correlations gross . Fragments with $`5Z20`$ as in indra-bologna  have been used for constructing the correlations. This observable was proven to reflect subtle topological features of the fragmentation events. The beginning of the correlation function is often called “Coulomb hole” and reflects the strength of the Coulomb repulsion between fragments (obviously related to the size of the freeze-out volume) but, as shown in the next section, is also related to the amount of radial expansion. Good agreement between calculations and experimental data can be observed here as well (see Fig. 5). While in the case of the Gd+U reaction the experimental velocity correlations are well reproduced in both freeze-out hypotheses, in the case of the first reaction a better agreement with experimental data is obtained for the hypothesis (2) (“without hard-core”).

In fact, this result is natural and suggests that deviations from the spherical fragments with hard-core interaction idealization (deformed fragments, surface diffusivity, etc.) are present in real multifragmentation. The deviation effect is stronger in the case of the Xe+Sn since the corresponding equilibrated system is smaller and, therefore the fragment partitions are more symmetric (as shown in vfree , the free volume is smaller for more symmetric partitions and larger for more asymmetric ones). The simultaneous description of both average fragment kinetic energies as a function of the fragment charge and of the reduced velocity correlations is remarkable and reflects a good description of the event topology. Here it is worth noticing that in Ref. indra-bologna  the velocity correlations for both reactions could not be reproduced via the SMM model.

The values of the model input parameters for which the comparisons with the experimental data have been performed are listed in Table 1. A remarkable fact is that the fitting parameters corresponding to the freeze-out hypotheses (1) and (2) are quite close ones to the others for each of the considered reactions. This tends to indicate a relatively small “region”, specific to each reaction in which the parameters of the real equilibrated nuclear source are situated. Quite small differences from one hypothesis to the other seem to be present in the case of the freeze-out volume. This result may look intriguing given the fact that freeze-out volume is directly related to both working hypotheses. In reality, the relatively large freeze-out volumes obtained are minimizing the influence of the hard-core interactions to the free-volume ($`V_{free}`$) (see e.g. vfree ), the values of of the freeze-out volume being mainly dictated by the Coulomb interaction (the remaining quantity entering the system’s density of states depending on the freeze-out volume). Unlike the free volume which in models can only be more or less arbitrary parametrized, Coulomb interaction can be accurately evaluated for each fragment partition (either by explicitly accounting for the interfragment interaction, or by employing the Wigner-Seitz approach). This contributes to the solidity of the present evaluations.

The differences in the excitation energies corresponding to the two hypotheses are due to the Coulomb energy variation from one case to another, caused by the slightly different values of the freeze-out volumes (see Table 1). For example, in the Gd+U case the energy variation (from (1) to (2)) is slightly larger than in the Xe+Sn case since the relative variation of the freeze out volume is larger as well. Freeze-out volume is always smaller in the working hypothesis (2). This is related to the different dependencies $`V_{free}(V)`$ resulting from the underlying hypotheses (more symmetric partitions and thus a more advanced fragmentation are specific to the without hardcore case and therefore both the freeze-out volume and excitation energies need to be smaller in the second hypothesis). The flow parameters were chosen such as to insure a good reproduction of the kinetic observables (see Fig. 5). The way in which each parameter is influencing the fitting observables is discussed in the next section.

## IV Influence of the model parameters on various fitting observables

While many fragment size and kinetic observables were simultaneously fitted in order to deduce the parameters of the equilibrated sources corresponding to the two reactions, the following question arises: are they unique? We address this question in the present section.

We start with the remark that the size of the equilibrated source appears to be dictated by the $`Z_{bound}`$ distribution. We illustrate this by representing in Fig. 6 the $`Z_{bound3}`$ distributions corresponding to three different sources: (215,90), (230,94) and (245,102).

Their excitation energies (5.4, 5.3, 5.3 MeV/u respectively) are chosen such as to insure a good fit of the charge distributions corresponding to the Xe+Sn at 32 MeV/u reaction; their freeze-out volume ($`V/V_0=9`$) is chosen as to provide a good description of the $`Z_{max}`$, $`Z_{max2}`$ and $`Z_{max3}`$ distributions. It can be easily observed that the best fit on the data is given by the (230,94) source. A smaller source deviates the maximum of the $`Z_{bound3}`$ distribution to the left while a larger one deviates it to the right (see Fig. 6).

The freeze-out volume of the equilibrated sources appears to be dictated by the distributions of the largest, second largest and third largest charge in one fragmentation event. Obviously, these distributions provide a measure of the degree of asymmetry of a given fragment partition. The degree of asymmetry of a partition is influenced by the freeze-out volume through both the Coulomb interaction and the free volume vfree . This fact is evidenced in Fig. 7 where multiplicities of $`Z_{max1}`$, $`Z_{max2}`$ and $`Z_{max3}`$ are represented for two different freeze out volumes in the freeze-out hypotheses (1) and (2) as follows. For the first hypothesis the source was considered (230,94) at two different volumes: $`V_0/5.5`$ and $`V_0/9`$ while for the second one the source was taken to be (220,92) at other two different volumes $`V_0/3`$ and $`V_0/8.5`$. The corresponding excitation energies are taken such as a good reproduction of the experimental charge distributions corresponding to the Xe+Sn reaction at 32 MeV/u to be achieved (the excitations energies are respectively: 6.4, 5.3, 6.1, 4.78 MeV/u). It can be noticed that deviations to the left of the maxima (indicating more asymmetric fragment partitions) of the $`Z_{max2}`$ and $`Z_{max3}`$ distributions are occurring when the freeze-out volume is decreased.

This behavior can be observed irrespective to the considered freeze-out hypothesis or to the chosen mass of the source. While, as demonstrated by Fig. 7 lower values of the freeze-out volume are leading to deviations of the $`Z_{max1}`$, $`Z_{max2}`$ and $`Z_{max3}`$ calculated distributions from the corresponding experimental data, larger values of the freeze-out volume will result in smaller excitation energies and larger flow energies leading to deviations from the experimental data of the kinematic observables represented in Fig. 5.

Given the previously described monotonical behavior of the “distance” from the calculated curves to the experimental ones, and also the very good reproduction of all the considered experimental observables for the parameters (given in Table I) corresponding to the case in which the $`Z_{bound3}`$ distribution is fitted best (see Fig. 6) it results that the corresponding sets of parameters provide the global minimum of the associated $`\chi ^2`$ function. Indeed, choosing a different value of $`A`$ would generate a deviation of the $`Z_{bound3}`$ from the experimental one and, obviously the global reproduction of the data would be worst.

Finally, the way in which the radial flow parameters $`E_{fl}`$ and $`\alpha `$ are influencing the experimental data fitting is presented in Fig. 8 for the nuclear source (230,94), $`V=V_0/9`$, $`E_{ex}=5.3`$ MeV/u and the first working freeze-out hypothesis. Three different cases are being considered: ($`\alpha =2`$, $`E_{fl}=1.4`$ MeV/u), ($`\alpha =1`$, $`E_{fl}=1.4`$ MeV/u) and ($`\alpha =2`$, $`E_{fl}=0.7`$ MeV/u). As shown in Fig. 8 the $`\alpha `$ parameter defining the flow profile appears to influence the location of the maximum of the average kinetic energy versus $`Z`$ ($`E_{c.m.}(Z)`$) and the value of the flow energy the height of the $`E_{c.m.}(Z)`$ dependency. Smaller effects are to be noticed in the reduced velocity correlation function and they concern mainly the width of the “Coulomb hole”. The effects are intuitive: smaller flow leads to a narrower Coulomb hole while larger ones lead to a larger hole. The parameters $`\alpha =2`$ and $`E_{fl}=1.4`$ MeV/u appear to give the best description of the kinetic experimental data corresponding to the Xe+Sn at 36 MeV/u reaction.

Therefore, it results that, within some (small) inherent uncertainties, the obtained parameters corresponding to the two working hypotheses ((1) and (2)) are uniquely determined.

## V Conclusions

In summary, the microcanonical multifragmentation model noi1  has been further refined and improved. Conservation of angular momentum has been included. Also, in order to provide events similar with the experimental ones, the model is completed with a stage of generation of momenta for the primary break-up fragments. Radial flow velocities can be superimposed here. Finally, an expansion stage of the hot primary fragments under their mutual Coulomb interaction is added. Using the new version of the model, the recent experimental data of the INDRA collaboration concerning the reactions Xe+Sn at 32 MeV/u and Gd+U at 36 MeV/u are analyzed. Both fragment size distribution observables and kinematic observables could be simultaneously fitted allowing thus the identification of the equilibrated sources parameters corresponding to the two reactions. Two different break-up scenarios were used for performing the analysis. The two types of calculations give close results for each of the considered reactions suggesting independence of the obtained results on the particular hypothesis. This is due to the fact that for relatively large values of the freeze-out volume (such as those obtained by fitting the data, 8-9 $`V_0`$) the influence of the hard-core is small and thus the results are rather independent on the extent to which these effects are included in a specific $`V_{free}`$ parametrization. In this case, the values of the freeze-out volume are dictated by the Coulomb interaction (accurately evaluated independently to any model assumptions) which makes the evaluation particularly robust. For each reaction a narrow freeze-out volume region was identified (8.5-9 $`V_0`$ for Xe+Sn and 8-8.5 $`V_0`$ for Gd+U). By analyzing the way in which variations of the model parameters are influencing the deviation between the calculated curves and the experimental ones, it is proved that the set of the parameters providing the best fit to the data is unique (within inherent uncertainties). It is shown that the size of the source is dictated by the $`Z_{bound}`$ distributions while its freeze-out radius can be deduced by fitting the experimental distributions of $`Z_{max1}`$, $`Z_{max2}`$ and $`Z_{max3}`$ and the kinematic observables. The way in which the flow energy and the parameter $`\alpha `$, defining the flow profile is influencing the fragment average kinetic energy versus $`Z`$ and also the reduced velocity correlation is also discussed. The above parameters are shown to be clearly identifiable as well. To our knowledge, this is the first complete reproduction of an entire set of experimental data (including both fragment size distribution observables and kinematic observables) by means of a microcanonical multifragmentation model. Moreover, it was shown for the first time that the parameters of the experimentally obtained equilibrated source can be uniquely determined by means of such a model.

###### Acknowledgements.

The authors acknowledge the support of the Alexander von Humboldt Foundation during this work.