text:
# Sensitivity of an underwater acoustic array to ultra-high energy neutrinos

## I Introduction

The understanding of the highest energy cosmic rays represents one of the most challenging fields of modern physics. While abundant high quality data has greatly enhanced our knowledge of cosmic rays of energies up to several tens of GeV, the study of higher energies is limited by the low fluxes available. Yet the study of particles with energies in excess of $`10^{18}`$ eV (ultra-high energy cosmic rays, or UHECR) promises boundless opportunities of discovery. To date some 12 cosmic ray events have been observed with energy in excess of $`10^{20}`$ eV Watson (2000). While the acceleration mechanisms at these energies are not completely understood Olinto (2000), even more fundamental problems arise from the apparent inconsistency of data with the Greisen-Zatsepin-Kuzmin (GZK) cutoff. Such cutoff Greisen (1966); Zatsepin and Kuzmin (1966) is expected to limit the maximum energy of protons of cosmological origin somewhere below $`10^{20}`$ eV, because of the finite ($`50`$ Mpc) inelastic collision length of such particles in the cosmic microwave background radiation (CMBR). Indeed the inelastic scattering of UHE protons off CMBR should at the same time suppress the proton flux and breed neutrinos from decay-products of pions. While the present data seem to indicate that photons and neutrinos are not the main component of UHECR Watson (2000), reliable identification of the primary particle type and its energy are essential parameters for the study of this problem.

It was recently pointed out Halzen (2000) that a substantial UHE neutrino component may accompany UHE protons and nuclei, due to neutrino production in cosmic beam dumps. Neutrino production at ultra-high energies appears in fireball models of gamma-ray bursts Vietri (1998); Waxman (1995); Waxman and Bahcall (1997), active galactic nuclei Gaisser et al. (1995); Protheroe (1996a, b, 1998); Mannheim (1995); Stecker and Salamon (1996) and, to a lesser degree, in Galactic mechanisms Domokos et al. (1993) and, as mentioned, the GZK process Yoshida and Teshima (1993). Another hypothetical source is the decay of heavy objects predicted by some theories, known as the “top-down” models of UHECR production Olinto (2000); Sigl et al. (1997, 1999).

Weakly interacting neutrinos could, unlike UHE gamma rays and protons, reach us from distant and powerful sources, opening a deeper horizon for astrophysics, cosmology and, possibly, high-energy particle physics Gaisser et al. (1995); Gandhi et al. (1996). While atmospheric neutrinos represent an irreducible background for Earth-based detectors, such background is expected to be modest because of the extremely long decay length of pions at the energies of interest.

Three techniques have been used until now to detect UHE cosmic-rays, all involving the showering of the primary particle in the Earth’s atmosphere. The shower is then detected either by observing the fluorescence or Čerenkov light induced by the ionizing tracks in the air, or by directly detecting the charged particles in the shower tail with scintillation counters scattered on the ground O’Halloran et al. (1998); Weekes (1996); Watson (2000). While some of the these techniques provide the largest acceptances obtained in particle detectors, in general the study of UHECR is still hampered by the very low flux one has to be sensitive to. Typical fluxes are $`100`$ km-2y-1 above $`10^{18}`$ eV, $`1`$ km-2y-1 above $`10^{19}`$ eV, and $`1`$ km-2century-1 above $`10^{20}`$ eV Watson (2000). Dedicated neutrino telescopes using Čerenkov light under water and the Antarctic ice-cap are for the time being optimized for the TeV to EeV energy-region Halzen (2000).

While much has been learned from the above detectors, it is important to explore, in parallel, new methods that could either increase the flux sensitivity, and hence raise the energy threshold for detection, or help constraining the primary particle identification and its energy. Alternative methods being discussed Gaisser et al. (1995) include the detection of radio Čerenkov emission from the lunar soil Gorham et al. (1999) (sensitive only to neutrinos that can cross the moon and interact upon exiting the satellite from its near side) and active radar detection of showers in the atmosphere Gorham (2000) (which has different systematics for particle identification since it is more sensitive to horizontal showers than other techniques). As discussed in this paper, another possibility consists in the detection of showers by means of the acoustic energy released in the medium where they develop. While, at least in principle, this technique could use as a radiator either the soil of the moon or a large body of water on the earth, here we concentrate on this second case that is more practical and, as we will show, may be possible to test on a large scale with a very modest effort. It is important to realize that these two radiators would provide very different information, the first being sensitive to UHECR of any type, while the second being sensitive essentially only to neutrinos because of the filtering effect of the earth’s atmosphere. While the idea of taking advantage of the very high mechanical $`Q`$ for small oscillation amplitudes and ultra-low seismic noise characteristic of the moon was first mentioned in Learned (1989), acoustic detection of ionizing particles in water was proposed in Askarjan (1957) and then developed Learned (1979) in connection with the DUMAND project Bradner (1977); Parvulesku and Curtis (1980); Roberts (1992). Although the primary goal of DUMAND was optical Čerenkov detection of muon tracks in deep ocean water, the hydrophones, originally conceived to monitor photomultiplier positions, were proposed to be used for acoustic detection of neutrinos of $`>10^{16}`$ eV energy. These early ideas were nicely complemented by experimental data collected at accelerators Sulak et al. (1979). We will review the past work on the subject and present a model useful to study UHE neutrino detection in sea water. We will then use this model to simulate signals and develop an algorithm to optimally filter hydrophone data to extract the expected signals from statistical noise. While most the assumptions for the study are rather general, we will apply our results to the case of a large, high frequency test array that the US Navy operates off the coast of Florida.

We note here that UHE neutrino fluxes substantially larger than the UHECR flux quoted by the experiments utilizing air showers could have gone unobserved. The neutrino cross section due to neutral and charged current interactions in the UHE regime is calculated in Kwiecinski et al. (1999). Even for $`E_\nu 10^{20}`$ eV the probability of interaction in traversing entire atmosphere’s depth is only of the order of $`10^5`$. Many models Protheroe (1996a, b, 1998); Mannheim (1995); Stecker and Salamon (1996); Waxman and Bahcall (1997); Sigl et al. (1997, 1999), predict neutrino fluxes higher than the flux of the UHECR observed in the atmosphere. For example, fireball gamma-ray burst models predict a neutrino flux at energies $`E_\nu >10^{20}`$ eV as high as $``$1 km-2y-1, active galactic nucleus models $`0.1`$ km-2y-1, and top-down models $`10`$ km-2y-1, as quoted in Kwiecinski et al. (1999).

## II Acoustic signature of particles in water

At the energies of interest here, neutrinos interact in water by a deep inelastic scattering on quarks inside oxygen and hydrogen nuclei. The scattering produces a lepton and an hadronic shower, with similar energies shared among these two components. The structure of hadronic showers is the same for all three flavors of (anti)neutrino. The behavior of the leptons, however, is different. In the case of $`\nu _\mathrm{e}`$’s the lepton energy goes in an electromagnetic shower and is essentially detected together with the hadronic energy. For $`\nu _\mu `$’s calculations by Mitsui Mitsui (1992) show that the mean-free-path in water between catastrophic bremsstrahlung and direct pair production events with energy transfer greater than $`10^{19}`$ eV is large (several km) with respect to the vertical size of the detectors considered here. The residual ionization along the track is low, so that muons are virtually undetectable by acoustic methods Parvulesku and Curtis (1980). Tau neutrinos create $`\tau `$ leptons, whose mean free path turns out to be long enough to leave the detector volume for $`E_{\nu _\tau }>10^{17}`$ eV.

At $`E_\nu 10^{20}`$ eV in water a hadronic shower deposits 90% of its energy in a cylinder of some 20 cm radius and 20 m length. Sound is produced in water mainly by heating localized along the shower, resulting in volume expansion, as first suggested in Askarjan (1957). The shower development (and hence the energy deposition) occurs at the velocity of light and can be regarded as instantaneous for the purpose of acoustic phenomena. As already mentioned the sound generation has been confirmed experimentally using artificial particles Sulak et al. (1979). Although it is not completely clear to what extent these accelerator experiments, simulating $``$1020 eV energies using bunches of over 1011 protons of $``$200 MeV, can describe the details of the lower energy-density UHE neutrino interactions, their accuracy is probably sufficient for the present study.

We analyze the acoustic signal production following Learned (1979). Let the energy deposited per unit volume per unit time be given by a function $`E(𝐫,t)`$. The total neutrino energy is $`E_0=_VE(𝐫)d^3𝐫`$. The wave equation for the pressure pulse produced $`p`$ is:

$$^2\left(p+\frac{1}{\omega _0}\dot{p}\right)\frac{1}{c^2}\ddot{p}=\frac{\beta }{C_p}\frac{E}{t}$$
(1)
Where we use the parameters for sea water: $`c1500`$ m/s is the speed of sound, $`\beta 1.2\times 10^3`$ K-1 is the bulk coefficient of thermal expansion, $`C_p3.8\times 10^3`$ J kg-1K-1 is the specific heat at constant pressure, and $`\omega _02.5\times 10^{10}`$ s-1 is the characteristic attenuation frequency. $`\omega _0`$ is, strictly speaking, a function of frequency Fisher and Simmons (1977), as can be seen from the plot of the attenuation coefficient $`\alpha (\mathrm{dB}/\mathrm{km})=(10^4/\mathrm{ln}10\omega _0c)(2\pi f)^2`$ in Figure 1. For simplicity of calculations, we assume that $`\omega _0`$ is a constant in the frequency range $`f=10100`$ kHz characteristic of the signal.

We note here that the coefficient of thermal expansion $`\beta `$ depends upon the water temperature. It vanishes at $`3^{}`$C for typical sea water of 3.5% salinity Bolz and Tuve (1973). In the case of vanishing $`\beta `$ other mechanisms of energy coupling to acoustic modes have been proposed Hunter et al. (1981); Lyamshev (1992). Extreme temperature profiles as a function of depth are given in Figure 2 for the tropical waters of the site discussed below. As we can see, although $`\beta `$ decreases with depth, it does not reach zero, so we will restrict ourselves to the case of thermal emission mechanism.

The instantaneous nature of the heating mechanism can be expressed by $`E(𝐫^{},t)=E(𝐫^{})\delta (t)`$. The pressure wave can then be calculated at location $`𝐫`$ as a function of time $`t`$ as

$$p(𝐫,t)=_VE(𝐫^{})G(𝐫𝐫^{},t)d^3𝐫^{}$$

where $`G(𝐫,t)`$ is the pressure pulse generated by a point source $`E=\delta (𝐫)\delta (t)`$, taking attenuation into account:

$$G(𝐫,t)=\frac{\beta }{4\pi C_p}\frac{(tr/c)}{r\sqrt{2\pi }\tau ^3}e^{(tr/c)^2/(2\tau ^2)}$$
(2)
where $`\tau =\sqrt{r/(\omega _0c)}`$.

The region of the energy deposition is elongated in the direction of the initial velocity of the neutrino, which constitutes the axis of the cascade. The acoustic emission is coherent in the plane perpendicular to this axis. Thus, the radiation diagram has a pancake shape, perpendicular to the shower axis, as shown in Figure 3. This effect was also observed in the 30 cm long and 4.5 cm wide cylinder energy deposition by particles in one of the accelerator experiments described in Sulak et al. (1979).

We calculate the time dependence of the pressure at a distance $`r=1`$ km from the origin, for different directions of observation. As shown in Figure 3, the angle $`\theta `$ is calculated in the forward direction from the perpendicular to the shower axis and the origin is at the starting point of the shower. Since the maximum energy deposition does not occur at the beginning of the shower, the maximum acoustic pulse results at some angle $`\theta >0`$.

To calculate the effect of a $`10^{20}`$ eV hadronic shower, we use a model based on the results of one-dimensional Monte Carlo simulations in ice Alvarez-Muñiz and Zas (1998). This model includes the interactions of $`\pi ^0`$ and other short-lived resonances, as well as the Landau-Pomeranchuk-Migdal (LPM) effect Landau and Pomeranchuk (1952, 1953); Migdal (1956, 1957), that are important at UHE. The transverse structure is modeled after Sokolsky (1989). The results of acoustic emission by a hadronic shower are presented in Figure 4.

The pressure pulse has a typical bipolar shape and is shown, at the maximum of the radiation pattern, in Figure 4a. The peak pressure and energy fluence for other angles of observation are shown in Figure 4b. The signal at other observation angles has a smaller amplitude and is stretched in time, still preserving the bipolar shape.

To calculate the case of an electromagnetic shower we use the Monte Carlo model LPMSHOWER Alvarez-Muñiz and Zas (1997) that includes the LPM effect. In this case, with all the energy in the electromagnetic channel, the LPM effect dominates the shower shape, resulting in very large fluctuations, and a significant elongation for initial energies $`>`$1 EeV. For $`E_e=10^{20}`$ eV the shower is $``$300 m long. The non-uniform energy deposition is the result of individual sub-showers starting at the locations of large energy loss. These sub-showers also create peaks in the acoustical radiation pattern. The properties of a typical shower are presented in Figure 5. The multi-peak structure is clearly visible in the amplitude and fluence diagrams as function of the angle. While the pressure pulse shape, taken at the direction of maximum emission, has, also in this case, a simple bipolar shape, the maximum occurs at a large angle respect to the previous case, because of the longer shower profile.

The signal from electron (anti)neutrinos is the superposition of these two cases in a proportion corresponding to the way the energy of the primary is shared between the hadronic and electromagnetic shower. It was found Gandhi et al. (1996) that at $`E_{\nu (\overline{\nu })}=10^{20}`$ eV the hadronic component accounts for $``$20% of the total energy. For an UHE primary, the axes of the hadronic and electromagnetic showers are practically parallel, the angle between them being of the order of $`10^6`$ radians Mitsui (1992).

We note here that the peak pressures predicted at 1000 m from $`10^{20}`$ eV showers are well within the sensitivity of good quality hydrophones. Indeed typical sensitivities for the frequency band of interest are $`10^3`$ Pa Sulak et al. (1979). However it is rather clear that the two factors limiting the power of this technique will be the ambient noise and the characteristic emission pattern described above. Such pattern substantially limits the solid angle accessible to each sensor and hence, together with the noise level, it will dictate the maximum tolerable spacing between detection sites.

## III The AUTEC array as a UHE neutrino detector

While the installation of hydrophones to complement underwater Čerenkov arrays has been discussed by several groups Domogatsky (2000); Butkevich et al. (1998), it is interesting to consider whether a very large existing array could be used in parasitic mode for UHE neutrino detection. Although the scientific community has recently started discussing the use of long-range, early warning military arrays for a variety of oceanographic purposes Richelson (1998), the low bandwidth of these systems, essentially designed to detect low-frequency ship noise that can propagate over very large distances, make them quite unsuitable for the type of signals discussed here. Another type of arrays exists, generally designed to track ships and weapons equipped with special high frequency “pingers” during limited-range naval exercises. Such arrays have a relatively high-density of hydrophones with typical bandwidths in the tens of kHz. The array considered here, at the Sites 3 and 4 of the Atlantic Undersea Test and Evaluation Center of the US Navy (AUTEC) Atlantic Undersea Test and Evaluation Center (2000), covers an approximate area of 250 km2 ($`5\times 15`$ nautical miles2) with depths between 1400 and 1600 m, as shown in Figure 6.

AUTEC is located in the “Tongue of the Ocean”, a $`50\times 200`$ km2 tract of deep sea, bounded to the west by Andros Island and to the south and east by large areas of very shallow banks in the Bahamas. This peculiar geographical configuration, with shipping access only from the north through the narrow Providence Channel, provides quiet conditions, because of the low boat traffic and sluggish currents.

Individual hydrophones at Sites 3 and 4 are mounted at the ends of 4.5 m long booms extending above the ocean floor. A total of 52 sensors cover the two sites, arranged on a triangular lattice with 2.5 km sides. Analog signals from each hydrophone are preamplified and brought to shore where digitizers and processors are located. The frequency response of the hydrophone and analog chain is flat to within $`\pm 5`$ dB in the range from 1 to 50 kHz, while the sampling rate of the digitizers is about 100 kHz. Accurate GPS time stamping is provided in the data stream. In its normal operation the system is capable to gather highly accurate 3-dimensional in-water tracking data. While a denser sensor spacing would of course be desirable, the shower-to-hydrophone distance of 1 km we use though this work is close to the worst case scenario of a neutrino interacting half-way between two sensors.

In Figure 7 we show the random noise levels at AUTEC for different wind conditions, along with the cumulative probability for such conditions to occur on range. The approximate frequency spectrum of the expected neutrino signals is also given for reference. The noise spectra in the Figure are due to the waves at the ocean surface induced by wind Wenz (1962); Knudsen et al. (1948); Urick (1975) and, above $`>`$10 kHz, to thermal noise Wenz (1962). In addition coherent noise from human activities and natural phenomena should be considered. A-priori we expect two types of artificial noise: the first, due to ship screws, is mainly confined to low frequency and hence easy to reject, while the second, due to the range “pingers”, is concentrated at a few well known frequencies and hence easy to filter-out. In general exercises are performed on range roughly 50% of the time, so that very little man-made noise is expected for a substantial fraction of every day (typically outside of working hours). More serious is probably the high frequency noise produced by marine mammals using sonar to localize their prey and snapping shrimps Epifanio (1997). The severity of these backgrounds depends on the season and can only be quantitatively understood by analyzing a substantial amount of data from the array. Here we limit ourselves to a detailed analysis of the random noise, formulating a filtering algorithm that can be implemented on the digital signal processors (DSP) that analyze on-line the time-series from each hydrophone. This system of data selection can be used essentially as a trigger signal to log an interval of the data stream for the entire array (or maybe only some subset of sensors near to the one producing the trigger). Further data reduction, possibly involving the correlation of signals from different sensors, can then be done off-line.

Scattering of sound off the ocean surface and bottom may alter somewhat the simple picture given above. For instance we expect that each hydrophone will generally record two pulses: the direct one and the one scattered from the ocean floor immediately around the sensor. The time delay will be of the order of $`h/c3`$ ms, where $`h=4.5`$ m is the height of the hydrophones above the bottom. This effect can be used to estimate the angle of incidence of the acoustic wave and hence of the shower axis. In addition scattering phenomena can make the event detectable by more than one hydrophone. The attenuation and re-emission patterns from scattering at the sea-bed and surface (where bubbles play an important role) can be numerically estimated Beckmann and Spizzichino (1963); Crowther (1980). While large attenuations are to be expected so that it appears unlikely that these phenomena can be used to trigger an event that for geometrical reasons was not directly visible, it is probable that scattered signals can be detected below threshold by sensors in the vicinity of the triggering hydrophone. While off-line study of these correlations would be particularly important in confirming the event and in measuring position, orientation and energy of the shower, in the rest of this paper we concentrate on the triggering function that has to rely on single sensors.

## IV The sensitivity of the detector

In order to analyze the problem of signal detection in some detail we assume a discrete data sample set $`x_k`$ ($`k=0,\mathrm{},M1`$) with a sampling rate of $`f_s=100`$ kHz, as provided by the AUTEC digitizers. We consider the detection from a single site, as it is natural given the small probability that the radiation pattern discussed in the previous section would intercept more than one hydrophone. Sub-threshold use of signals from neighboring hydrophones may be possible in off-line analysis but it is not discussed here. The data are the sum of the stationary Gaussian noise of a given power spectrum and a signal of a given shape, but unknown amplitude, starting at position $`k_0`$ in the data set, as illustrated in the simulated time series of length $`M=10000`$ in Figure 8. In this Figure a signal from a $`2\times 10^{19}`$ eV hadronic shower at 1 km with the optimal angular orientation to illuminate a hydrophone has been superimposed at a random time (chosen in this case to start at sample number $`k_0=5000`$) over the statistical noise spectrum relative to 13 knots wind from Figure 7. These conditions (or better) occur at the AUTEC site 60% of the time. The sampling frequency is high enough to provide negligible distortions.

We extract the signal using the algorithm described in Appendix, which is based on the digital filter of transfer function $`H_l`$ calculated on the basis of the signal shape and the Gaussian noise spectrum given in equation (8). The variable $`Y_k`$ resulting from the application of the filter

$$Y_k=\underset{l=k}{\overset{k+N1}{}}H_lx_{k+l}$$

has a Gaussian distribution and is used to assert the presence of the signal, using a threshold $`Y_{\mathrm{th}}`$ that can be chosen for a certain detection efficiency and false-alarm rate at a given signal amplitude $`A`$ ($`Y_k>Y_{\mathrm{th}}`$ indicating the presence of a signal). In practice the value of $`N`$ can be made small enough for efficient calculation on the commercial DSP processors used at AUTEC.

In Figure 9 we plot the variable $`Y_k`$ for the time series in Figure 8. We see that $`Y_k`$ at the time sample $`k=k_0=5000`$ has a value substantially higher that elsewhere. Using a threshold $`Y_{\mathrm{th}}=Y_{k_0}`$, we find the probability of false alarm to be $`5\times 10^{16}`$ at each data point. On the other hand the naive technique of trying to find the signal using an amplitude threshold in time domain would give a false alarm probability of 0.14 at each sample point, with a threshold set at amplitude of the signal. The effectiveness of the method, also evident by simple inspection of Figures 8 and 9, is based on the difference between the signal and noise spectra in Figure 7. This technique is equivalent to the application of a matched filter Bradner (1977); Cappellini et al. (1978). A value of $`N=7`$ was used in the calculation, giving $``$15 floating-point operations at each data point, which corresponds to a processor speed requirement of $`1.5`$ Mflop/s. The duration and amplitude of the signal can be estimated by maximizing a likelihood function.

We now use the probabilities of a signal miss and a false alarms found in the Appendix to quantitatively analyze the power of the array as a UHE neutrino detector. Here we stress the fact that our analysis does not take into account possible interference from coherent sources, as already mentioned above.

We define the detector efficiency as the fraction of neutrinos that are detected respect to those interacting in the volume of water $`=\mathrm{\Lambda }/\mathrm{\Lambda }_0`$. We can write the interaction rate $`\mathrm{\Lambda }_0`$ as

$$\mathrm{\Lambda }_0=F(E)N_A\rho _w\sigma (E)V_w$$

where $`F(E)`$ is the UHE neutrino flux at energy $`E`$, $`\sigma (E)`$ is the cross-section of interaction with a nucleon, $`N_A\rho _w`$ is the number of nucleons in unit volume of water, and $`V_w=ah`$ is the volume of the detector ($`a`$ and $`h`$ being the area and average depth of the array). The detected rate is limited by the effective volume of the radiation pattern $`V_E`$ in which the amplitude of the signal (relative to produced by a shower of energy $`10^{20}`$ eV at a distance of 1 km) is above the detection threshold $`A`$. We can write $`V_E`$ as:

$`\left({\displaystyle \frac{V_E}{1\mathrm{km}^3}}\right)`$ $`=`$ $`{\displaystyle \frac{2\pi \mathrm{\Delta }\theta }{3}}\left({\displaystyle \frac{R_E}{1\mathrm{km}}}\right)^3`$ (3)
$`=`$ $`{\displaystyle \frac{2\pi \mathrm{\Delta }\theta }{3}}\left[\left({\displaystyle \frac{E}{10^{20}\mathrm{eV}}}\right){\displaystyle \frac{1}{A}}\right]^3,`$
where $`\mathrm{\Delta }\theta `$ is the effective angular extent of the radiation lobe that can be found by integrating the curves in Figures 4b and 5b, and $`R_E`$ is its radial extent. We obtain $`\mathrm{\Delta }\theta 0.5^{}`$ ($`\mathrm{\Delta }\theta 1.5^{}`$) for showers induced by $`\nu _\mu `$, $`\overline{\nu }_\mu `$, $`\nu _\tau `$, $`\overline{\nu }_\tau `$ ($`\nu _e`$, $`\overline{\nu }_e`$). The difference comes about from the electron neutrinos having the electromagnetic part of the shower, while other flavors have only hadronic part. We note that the expression for $`V_E`$ includes a signal amplitude inversely proportional to the distance $`R_0`$, as it should in our regime of radiation Learned (1979).

The rate of detected neutrinos is then found by integrating over volume:

$$\mathrm{\Lambda }=N_hF(E)N_A\rho _w\sigma (E)_0^{V_{E,\mathrm{max}}}p_{\mathrm{detect}}(A,Y_{\mathrm{th}})𝑑V_E$$
(4)
where $`N_h`$ is the number of hydrophones and $`p_{\mathrm{detect}}`$ is the detection probability given in (11), with $`A`$ found from $`V_E`$ using (3). The upper limit of the integration $`V_{E,\mathrm{max}}`$ is determined from the condition that the acoustic radiation pattern is limited in size by the size of the detector, i.e. $`R_E<R_{\mathrm{max}}`$, where $`R_{\mathrm{max}}100`$ km. Finally, we get

$$=\frac{N_h}{V_w}_0^{V_{E,\mathrm{max}}}p_{\mathrm{detect}}(A,Y_{\mathrm{th}})𝑑V_E$$

We can now use this expression to set the value of $`Y_{\mathrm{th}}`$ for any given efficiency $``$ and then calculate the resulting false alarm rate as $`\lambda =N_hf_sp_{\mathrm{false}}`$, where $`p_{\mathrm{false}}`$ is determined from $`Y_{\mathrm{th}}`$ using (12). We plot the calculated false alarm rates for different efficiencies in Figure 10a. The solid lines correspond to $`\nu _e`$ or $`\overline{\nu }_e`$ that produce an electromagnetic and an hadronic showers, while the dashed lines are for other neutrino flavors that produce a narrower sound lobe from the shorter hadronic shower. As expected electromagnetic showers can afford a higher threshold that results in a lower false alarm rate.

## V Conclusions

We have reviewed the possibility of searching for ultra-high-energy neutrinos in cosmic rays using the acoustic emission from the electromagnetic and hadronic showers produced by the neutrino interactions in sea water. We have analyzed the expected statistical noise and devised an optimal algorithm to filter it out of the data stream in real-time. In a simulation we have applied this technique to trigger the data-acquisition of a very large, high-frequency multi-hydrophone array of the US Navy and found that one could trigger on events at or above $`10^{20}`$ eV with tolerable false alarm rates. Our algorithms are based on signals from individual sensors, as appropriate for a trigger system, and are optimized to run on the array’s digital signal processors. This would make a data taking campaign rather straightforward. Further off-line analysis would be needed to study the additional information that can be obtained by multi-hydrophone correlations. Although it is clear that an array with higher sensor density would allow lower energy thresholds and better redundancy, a test on the existent array analyzed here will provide information on coherent noise and give a definite assessment on the power of this technique.

###### Acknowledgements.

We would like to thank D. Bryant, J. Cecil, N. DiMarzio, T. Kelly-Bissonnette, D. Moretti (US Navy) for the many discussions on the neutrino detection capabilities of AUTEC. We are also indebted to J. Vandenbroucke (Stanford) for the help in programming the AUTEC signal processors. One of us (G.G.) is indebted to D. Kapolka (US Navy) for early guidance in understanding the high frequency capabilities of different arrays. We thank P. Gorham (JPL) for a critical reading of an early version of the manuscript. This work was supported, in part, by a Terman Fellowship from the Stanford University. Partial support was also provided by the Office of Naval Research under Grant No. N00014-93-1-0054.

## Appendix A The detection algorithm

In order to analyze the problem of signal detection in some detail we assume a discretized data sample from a single hydrophone. Since the amplitude of the signal can vary with respect to the distance to the source, we only fix the shape of the signal, without specifying its amplitude in advance. Let us assume that the signal shape after discretization with sampling frequency $`f_s`$, is represented as an array of N real numbers, $`F_n`$, where $`n=0,\mathrm{},N1`$, so that the signal is $`AF_n`$, where $`A`$ is the amplitude.

The discrete data sample $`x_k`$, $`k=0,\mathrm{},M1`$, which is discretized with sampling frequency $`f_s`$, is the sum of the stationary Gaussian noise $`w_k`$ of a given power spectrum and a signal embedded at a random position $`k_0=0,\mathrm{},MN`$ with an unknown amplitude $`A`$:

$$x_k=\{\begin{array}{cc}w_k,\hfill & k<k_0\mathrm{or}kk_0+N\hfill \\ w_k+AF_{kk_0},\hfill & k_0k<k_0+N.\hfill \end{array}$$

We will search for the signal position $`k_0`$ and the amplitude $`A`$ which maximizes the conditional probability of a signal to be present within a given data set. This is a regression problem that can be solved by finding the maximum in the appropriate likelihood function Hogg and Craig (1995). Let us take a subsample of the original data $`X_kx_{k_0+k}`$, where $`k=0,\mathrm{},N1`$, and analogously define a noise subsample $`W_kw_{k_0+k}`$. The likelihood of the signal presence at position $`k_0`$ is given by the Bayes formula

$`\mathrm{Prob}\left\{F,A[A,A+dA]\right|X\}`$ (5)
$`={\displaystyle \frac{\mathrm{Prob}\left\{X\right|F,A\}\mathrm{Prob}\left\{F\right\}f(A)dA}{\mathrm{Prob}\left\{X\right\}}}`$
The ratio to the probability of the absence of signal (which is found in a similar manner) is

$`{\displaystyle \frac{\mathrm{Prob}\left\{F,A[A,A+dA]\right|X\}}{\mathrm{Prob}\left\{\mathrm{no}\mathrm{signal}\right|X\}}}`$
$`={\displaystyle \frac{\mathrm{Prob}\left\{X\right|F,A\}\mathrm{Prob}\left\{F\right\}f(A)dA}{\mathrm{Prob}\left\{X\right|\mathrm{no}\mathrm{signal}\}\mathrm{Prob}\left\{\mathrm{no}\mathrm{signal}\right\}}}`$
We notice that some of these conditional probabilities are related to the distribution of the Gaussian noise:

$`\mathrm{Prob}\left\{X\right|F,A\}`$ $`=`$ $`f_W(XAF)dX`$
$`\mathrm{Prob}\left\{X\right|\mathrm{no}\mathrm{signal}\}`$ $`=`$ $`f_W(X)dX`$
where $`dX=_kdX_k`$. The probability distribution $`f_W(W)`$ of the Gaussian noise, defined so that $`f_W(W)dW`$ is the probability that the noise values $`W_k`$ are in the interval $`dW=_{k=0}^{N1}dW_k`$, is given by

$$f_W(W)=C\mathrm{exp}\left(\frac{1}{2}W^TK^1W\right)$$

where $`K`$ is the noise covariance matrix defined as $`K_{jk}=W_jW_k`$ that, for a stationary noise, has the property of being a function of the difference $`(kj)`$ only. If this portion of the time series is long enough so that the noise values separated by $`>N/2`$ points are not correlated, we can assume circular stationarity $`K_{jk}=R_{(kj)\mathrm{mod}N}`$, where $`R`$ is an array of length $`N`$.

The calculations are easiest in a basis in which the covariance matrix is diagonal. For a circularly stationary signal, this is the basis of Fourier components

$$\stackrel{~}{W}_k=\underset{n=0}{\overset{N1}{}}W_n\mathrm{exp}\left(\frac{2\pi i}{N}kn\right),k=0,\mathrm{},N1$$
(6)
The distribution of the noise Fourier components is defined so that $`f_{\stackrel{~}{W}}(\stackrel{~}{W})d\stackrel{~}{W}`$ is the probability of $`\stackrel{~}{W}`$ to be in the interval $`d\stackrel{~}{W}=d\stackrel{~}{W}_0d\stackrel{~}{W}_{N/2}_{k=1}^{N/21}d[\mathrm{Re}\stackrel{~}{W}_k]d[\mathrm{Im}\stackrel{~}{W}_k]`$. Hence the probability distribution is given by

$$f_{\stackrel{~}{W}}(\stackrel{~}{W})=\stackrel{~}{C}\mathrm{exp}\left(\frac{1}{2}\underset{k}{}|\stackrel{~}{W}_k|^2/S_k\right)$$

where $`S_k`$, $`k=0,\mathrm{},N1`$ is the noise spectrum defined as $`\stackrel{~}{W}_k^{}\stackrel{~}{W}_l=\delta _{kl}S_k`$. This is related to the covariance matrix as $`S=N\stackrel{~}{R}`$, with $`R`$ defined above.

Substituting $`f_{\stackrel{~}{W}}(\stackrel{~}{W})`$ we get the ratio of the likelihoods

$$\frac{\mathrm{Prob}\left\{F,A[A,A+dA]\right|X\}}{\mathrm{Prob}\left\{\mathrm{no}\mathrm{signal}\right|X\}}=\frac{\mathrm{Prob}\left\{F\right\}}{\mathrm{Prob}\left\{\mathrm{no}\mathrm{signal}\right\}}e^L^{}dA$$

where the argument of the exponential function is

$$L^{}=\frac{1}{2}\underset{k}{}\frac{1}{S_k}\left(|\stackrel{~}{X}_kA\stackrel{~}{F}_k|^2|\stackrel{~}{X}_k|^2\right)+\mathrm{ln}f(A)$$

Instead of maximizing the original probability, we maximize the log-likelihood $`L^{}`$.

Since $`A`$ is always positive and $`\mathrm{ln}f(A)\mathrm{const}`$ for $`A>0`$, we choose the optimal amplitude estimate $`\overline{A}`$ which maximizes $`L^{}`$:

$$\overline{A}=\mathrm{max}\{0,\frac{_k\mathrm{Re}(\stackrel{~}{X}_k\stackrel{~}{F}_k^{})/S_k}{_k|\stackrel{~}{F}_k|^2/S_k}\}$$

Substituting this value into $`L^{}`$, we form a random variable

$$L=\frac{1}{2}\underset{k}{}\frac{1}{S_k}\left(|\stackrel{~}{X}_k\overline{A}\stackrel{~}{F}_k|^2|\stackrel{~}{X}_k|^2\right)=L^{}(\overline{A})+\mathrm{const}.$$

Note that the maximization with respect to $`A`$ can be generalized to other parameters of the signal. Let $`F^{(m)}`$ be different signal shapes from a set given by an index $`m`$. Then to get the maximum likelihood of the signal shape we choose the value of $`m_0`$ such that $`L^{(m_0)}=\mathrm{max}\{L^{(m)}\}`$ (note that $`\mathrm{const}=\mathrm{ln}f(A)`$ is the same for all $`m`$).

It turns out that the equation for $`L`$ can be simplified in terms of the number of calculations required. Substituting into it the expression for $`\overline{A}`$, we derive

$$L=\{\begin{array}{cc}\frac{Y^2}{2\sigma _Y^2},\hfill & Y>0\hfill \\ 0,\hfill & Y0\hfill \end{array}$$

where $`Y=_k\mathrm{Re}(\stackrel{~}{X}_k\stackrel{~}{F}_k^{})/S_k`$ and $`\sigma _Y^2=_k|\stackrel{~}{F}_k|^2/S_k`$ can be shown to be the variance of $`Y`$. Substituting the Fourier transform definition (6) for $`\stackrel{~}{X}`$, we finally obtain

$$Y=\underset{l=0}{\overset{N1}{}}X_lH_l$$
(7)
where

$$H_l=\underset{k}{}(\stackrel{~}{F}_k/S_k)e^{2\pi ikl/N}.$$
(8)
The variable $`Y`$ is thus the result of applying a digital filter with response function $`H`$ or transfer function $`\stackrel{~}{H}_k=N\stackrel{~}{F}_k/S_k`$. The amount of calculation can be reduced even further, due to a fact that the quantity of $`H_k`$ that are significantly different from zero is $`N^{}<N`$. The number $`N^{}`$ is estimated to be $`\tau _{\mathrm{sig}}f_s`$ because $`H_k`$ has the time duration of the same order as the duration of the signal, $`\tau _{\mathrm{sig}}`$. This last fact is confirmed by a direct calculation. Thus, for typical values of $`f_s=10^5`$ Hz and $`\tau _{\mathrm{sig}}10^4`$ s, we have $`N^{}10`$. The number of floating-point operations (flops) required is then $`10^6`$ per second, which is feasible on modern DSP processors capable of more that 50 Mflops/s.

Signal detection can then be implemented at every hydrophone utilizing the fact that $`Y`$ has different distributions $`f_{Y,A}(Y)`$ in the cases presence ($`A0`$) and absence ($`A=0`$) of a signal. We chose a threshold value $`Y_{\mathrm{th}}`$ and decide that the signal is present (absent) when $`YY_{\mathrm{th}}`$ ($`Y<Y_{\mathrm{th}}`$). Since the distributions $`f_{Y,A}(Y)`$ for the cases of signal presence and absence in general overlap, there will be the possibility of false alarms and signal misses. The threshold value has then to be chosen in an appropriate manner.

We now estimate the probabilities of false alarms and misses for a certain threshold $`Y_{\mathrm{th}}`$ by estimating the probability distribution $`f_{Y,A}(Y)`$ of the variable $`Y`$ in the presence of a signal of amplitude $`A`$ (including the no-signal situation as a particular case $`A=0`$). Since $`Y`$ is a linear combination of Gaussian variables, the distribution sought is also Gaussian

$$f_{Y,A}(Y)=\frac{1}{\sqrt{2\pi }\sigma _Y}e^{(YA\sigma _Y^2)^2/(2\sigma _Y^2)}.$$
(9)
The probability of missing a signal is then

$`p_{\mathrm{miss}}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{Y_{\mathrm{th}}}}f_{Y,A}(Y)𝑑Y`$ (10)
$`=`$ $`1{\displaystyle \frac{1}{2}}\mathrm{erfc}\left({\displaystyle \frac{Y_{\mathrm{th}}A\sigma _Y^2}{\sqrt{2}\sigma _Y}}\right),`$
while the probability of the signal detection is

$`p_{\mathrm{detect}}=1p_{\mathrm{miss}}`$ $`=`$ $`{\displaystyle _{Y_{\mathrm{th}}}^{\mathrm{}}}f_{Y,A}(Y)𝑑Y`$ (11)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{erfc}\left({\displaystyle \frac{Y_{\mathrm{th}}A\sigma _Y^2}{\sqrt{2}\sigma _Y}}\right).`$
Finally, the false alarm probability is

$$p_{\mathrm{false}}=_{Y_{\mathrm{th}}}^{\mathrm{}}f_{Y,A=0}(Y)𝑑L=\frac{1}{2}\mathrm{erfc}\left(\frac{Y_{\mathrm{th}}}{\sqrt{2}\sigma _Y}\right)$$
(12)
In the expressions above $`\mathrm{erfc}\left(x\right)=1\mathrm{erf}\left(x\right)`$ is the complementary error function. The above expressions can be inverted to derive $`Y_{\mathrm{th}}`$ from $`p_{\mathrm{detect}}`$ and $`p_{\mathrm{false}}`$:

$$Y_{\mathrm{th}}=\{\begin{array}{c}A\sigma _Y^2+\sqrt{2}\sigma _Y\mathrm{erf}^1\left(12p_{\mathrm{detect}}\right)\hfill \\ \sqrt{2}\sigma _Y\mathrm{erf}^1\left(12p_{\mathrm{false}}\right)\hfill \end{array}$$
(13)
where $`\mathrm{erf}^1\left(x\right)`$ is the inverse error function.


format: text

text:
# References

## Acknowledgements

This work was supported by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF).

Correspondence and request for materials should be addressed to A.Z. (e-mail: anton.zeilinger@univie.ac.at)

Captions:

Figure 1: The wave front (top) and the intensity pattern (bottom) of the simplest Laguerre Gauss (LG$`{}_{}{}^{l}{}_{p}{}^{}`$) or doughnut mode. The index $`l`$ is referred to as the winding number and ($`p+1`$) is the number of radial nodes. Here we only consider cases of $`p=0`$. The customary Gaussian mode can be viewed as LG mode with $`l=0`$. The handedness of the helical wave fronts of the LG modes is linked to the sign of the index $`l`$ and can be chosen by convention. The azimuthal phase term $`e^{il\varphi }`$ of the LG modes results in helical wave fronts. The phase variation along a closed path around the beam center is $`2\pi l`$. Therefore in order to fulfill the wave equation the intensity has to vanish in the center of the beam.

Figure 2: Experimental setup for single-photon mode detection. After parametric down conversion each of the photons enters a mode detector consisting of a computer generated hologram and a mono-mode optical fiber. By diffraction at the hologram the incoming mode undergoes a mode transformation in a way that a LG mode can be transformed into a Gaussian mode. Since it has a smaller spatial extension than all LG modes, only the Gaussian mode can be coupled into the mono-mode fiber. Thus observation of a click projects the mode incident on the fiber coupler into the Gaussian mode. The hologram is a phase grating with $`\mathrm{\Delta }m`$ dislocations in the center blazed for first order diffraction. An incoming Gaussian laser beam passing through the dislocation of the hologram is diffracted by the grating and the n-th diffraction order becomes a LG mode with an index $`l=n\mathrm{\Delta }m`$ and vice versa. Intuitively speaking the phase dislocation exerts a “torque” onto the diffracted beam because of the difference of the local grating vectors in the upper and lower parts of the grating. This “torque” depends on the diffraction order n and on $`\mathrm{\Delta }m`$. Consequently the right and left diffraction orders gain different handedness. Reversing this process a photon with angular momentum $`\mathrm{\Delta }m\mathrm{}`$ before the grating can be detected by the mono-mode fiber detector placed in the first diffraction order. A photon with zero angular momentum (Gaussian mode) is detected by diffracting the beam at the border of the hologram faraway from the dislocation. All our measurements were performed in coincidence detection between the two down-converted photons.

Figure 3: Conservation of the orbital angular momentum.

Coincidence mode detections for photon 1 and photon 2 in 15 possible combinations of orthogonal states were performed. This was done for a pump beam having an orbital angular momentum of $`\mathrm{}`$, 0, and $`+\mathrm{}`$ per photon respectively. Coincidences was observed in all cases where the sum of the orbital angular momenta of the down converted photons were equal to the pump beams orbital angular momentum. The coincidence counts for each fixed value of the orbital angular momentum of photon 1 was normalized by the total number of coincidences varying the orbital angular momentum of photon 2.

Figure 4:

Experimental evidence of entanglement of photon states with phase singularities: The dislocation of the hologram in the beam of photon 1 is shifted out of the beam center step by step (top, middle, bottom). In these positions this hologram together with the mono-mode fiber detector projects the state of photon 1 into a coherent superposition of LG and Gaussian modes. The mode filter for photon 2 with the hologram taken out makes a scan of the second photon s mode in order to identify the location of its singularity with respect to the beam center. The coincidences show that the second photon is also detected in a superposition of the LG and the Gaussian mode. Classical correlation would yield a coincidence picture which is just a mixture of Gaussian and LG modes. In that case the intensity minimum would remain in the beam center but would become washed out. In the experiment a hologram with two dislocations in the first diffraction order was used. This results in a superposition of the l=0 and l=2 modes.


text:
# A high frequency radio study of G11.2-0.3, a historical supernova remnant with a flat spectrum core

## 1 Introduction

Only a small number of young supernova remnants (SNRs) has been identified in our Galaxy so far. One of them is the non–thermal source G11.2$``$0.3. Clark & Stephenson (clark (1977)) associated this SNR with the supernova explosion of A.D. 386 observed by Chinese astronomers during the Chin Dynasty. In the radio continuum G11.2$``$0.3 shows the typical structure of a shell–type SNR with a spectral index of $`\alpha 0.5`$ (S $`\nu ^\alpha `$). Green et al. (green (1988)) reported a clumpy shell and the lack of a sharp outer boundary based on high resolution VLA data. This indicates a SNR in an early evolutionary state. Green et al. (green (1988)) found no evidence for spectral index variations between 1.4 GHz and 5 GHz. However, Morsi & Reich (morsi (1987)) reported on a flat spectrum core based on early 32 GHz observations with the Effelsberg 100-m radio telescope.

H i–absorption measurements towards the SNR indicate a distance of about 5 kpc (Radhakrishnan et al., radha (1972); Green et al. green (1988) after reinterpretation of the absorption spectrum measured by Becker et al. becker1985 (1985)). At this distance G11.2$``$0.3 has a diameter of approximately 6 pc, which is comparable to other historical SNRs like Kepler and Cas A with $`\mathrm{D}4`$ pc and Tycho with $`\mathrm{D}6`$ pc.

X–ray observations of G11.2$``$0.3 show a soft shell–like structure which is most prominent in the southeast (Reynolds et al. reynolds (1994), see Fig. 1). In addition a hard spectrum core exists (Vasisht et al. vasisht (1996)), which was interpreted as the synchrotron nebula of an embedded pulsar. This pulsar was subsequently discovered by Torii et al. (torii1997 (1997)) through observations with the Advanced Satellite for Cosmology and Astrophysics (ASCA). It is a young pulsar with a period of 65 ms. The rotational energy loss rate was calculated to be $`\dot{\mathrm{E}}=8.8_{2.6}^{+9.9}\times 10^{36}`$ erg/s from follow–up X–ray observations (Torii et al. torii1999 (1999)). This is a rather low value for a young pulsar and indicates a low radio surface brightness for the surrounding synchrotron nebula (Kothes kothes (1998)).

In order to detail the properties of G11.2$``$0.3, in particular its central component and its polarization structure, we conducted a number of new observations at high radio frequencies.

## 2 Observations

The Effelsberg 100-m radio telescope was used to map fields around G11.2$``$0.3 at high radio frequencies. In 1996 observations at 4.85 GHz and 14.7 GHz were carried out as part of test measurements for newly installed receiver systems. In 1999 observations at 10.45 GHz and 32 GHz were added. At 32 GHz a new receiver system was used with improved sensitivity compared to Morsi & Reich’s (morsi (1987)) measurements at the same frequency. All receiver systems used were installed in the secondary focus of the telescope. The 4.85 GHz, 10.45 GHz, and 32 GHz observations were made by recording both circularly polarized components to obtain total intensity. IF correlation of the circular components results in linear polarization data. At 14.7 GHz just total intensity channels were available. The receiving systems consist of two feeds at 4.85 GHz, four feeds at 10.45 GHz and 14.7 GHz, and three feeds at 32 GHz. Observational parameters are listed in Table 1. All systems have highly stable cooled HEMT receivers at their front–ends and run in total power mode. Data reduction relies on application of the “software beam switching” technique to reject influences of weather effects (Morsi & Reich morsi2 (1986)).

All observations were made in the equatorial coordinate system (B1950). The scan direction was along the azimuth. The standard data reduction software package based on the NOD2 format has been applied. The multi–feed observations at 14.7 GHz and 32 GHz were restored using the algorithm described by Emerson et al. (emerson (1979)). At 4.85 GHz and 10.45 GHz data from different feeds were reduced separately and combined afterwards. Baseline improvement by unsharp masking (Sofue & Reich sofue (1979)) was applied to all observations. The “Plait” algorithm described by Emerson & Gräve (emerson88 (1988)) was used to combine the numerous coverages observed at 32 GHz at different parallactic angles by destriping the maps in the Fourier plane. This increases the signal–to–noise ratio of the final map substantially.

We made also use of an unpublished 23 GHz map with an angular resolution of 1$`\stackrel{}{.}`$2 observed in 1987 by W. Reich, P. Reich, and Y. Sofue with the Nobeyama 45-m radiotelescope. A cooled HEMT receiver was used and a map was made from two coverages at different parallactic angles. It was calibrated relative to NGC 7027 assuming 5.7 Jy. The reduction of the double–beam observations was described by P. Reich (pr (1987)).

## 3 Results

The total intensity maps at 10.45 GHz, 14.7 GHz, 23 GHz and 32 GHz are shown in Fig. 2. In Fig. 3 polarized intensities at 10.45 GHz and 32 GHz are displayed.

### 3.1 The radio emission from the shell

The overall structure of the radio continuum emission is most evident in the 32 GHz map (Fig. 2d), where the angular resolution is highest. The SNR shows an almost circular shell consisting of four prominent sectors, with the strongest towards the southeast. Three of those sectors show clear polarization with electric field vectors perpendicular to the expansion direction (Fig. 3). The southeastern sector is almost unpolarized.

### 3.2 The interior structure

Within the shell there are two slightly extended sources visible at 32 GHz (see Fig. 2d). On our lower frequency maps these sources appear unresolved and confused with the shell as a result of the lower angular resolution. Both inner sources seem to be polarized (see Fig. 3, Table 3) but a more prominent polarization feature is located slightly north of them. It has no distinct counterpart in total intensity and seems to be related to diffuse emission either from the shell or the inner structure of the SNR.

In order to measure the spectrum of the two compact sources we used the NVSS data at 1.4 GHz (Condon et al. condon (1998)) having an angular resolution of $`45\mathrm{}`$. We applied the “background filtering technique” (Sofue & Reich sofue (1979)) to separate the diffuse emission from the shell which is confused with the inner structure. We applied the method to the NVSS map and the 32 GHz map, convolved to $`45\mathrm{}`$, in exactly the same way to get comparable results. The resulting maps are displayed in Fig. 4.

The NVSS map shows two clearly separated sources, while at 32 GHz a bar-like structure with two emission peaks is visible. This results from an additional component showing up at 32 GHz, which is located between the two sources visible at 1.4 GHz. In the observations of Green et al. (green (1988)) there is a slightly extended source visible at the very center of the 4.76 GHz map which we believe has an inverse spectrum because it is not visible in their 1.46 GHz map. There are only the two extended diffuse structures which are also visible at 4.76 GHz.

### 3.3 The radio spectrum

The integrated spectrum of the SNR is plotted in Fig. 5. The flux density values obtained from the Effelsberg and the Nobeyama maps are listed in Table 2. Other values were taken from the references indicated. The fitted spectral index of $`\alpha =0.50`$ is a typical value of a mature supernova remnant.

However, G11.2$``$0.3 has a flat spectrum core component. Fig. 6 shows the radial variation of spectral index between 10.45 GHz and 32 GHz at 1$`\stackrel{}{.}`$15 angular resolution. The spectral index between these two frequencies calculated with the integrated flux densities from Table 2 is $`\alpha =0.45`$, just marginally flatter than the spectral fit in Fig. 5. The spectrum flattens from about $`\alpha =0.45`$ at the outer shell to about $`\alpha =0.3`$ in the inner part of the SNR. The $`\alpha `$ profile suggests a flatter radio spectrum in the interior of the remnant. The overall spectral index of $`\alpha =0.45`$ between 10.45 GHz and 32 GHz indicates that the flat component starts to flatten the overall spectrum of the source.

Integrating the flux density of the compact sources gives $`110\pm 20`$ mJy at 1.4 GHz and $`80\pm 15`$ mJy at 32 GHz. A spectral index of $`\alpha =0.10\pm 0.08`$ results. However, there are spectral variations. The central source is not visible at 1.4 GHz, thus indicating an inverted spectrum. This is supported by the observations of Green et al. (green (1988)). At 4.76 GHz a compact source with a peak flux density of about 2 mJy is visible which is absent in their 1.46 GHz map at the same angular resolution. The 4.885 GHz map of Becker et al. (becker1985 (1985)) at $`8\mathrm{}`$ angular resolution barely resolves that source, but a central enhancement of about 10 to 15 mJy/beam is seen. Of the outer sources, the spectrum of the eastern one seems slightly steeper than that of the western component. Clearly, sensitive higher angular resolution data are needed to determine the spectrum of the source components more precisely.

A separation of the shell emission and the central emission is difficult with the present data. Based on the pulsar data we estimate in Sect. 4.2 the flat spectrum radio emission to be about 1 Jy, which we assume to have a spectral index of $`\alpha =0.0`$. Subtracting the plerionic component from the integrated flux density values results in a spectral index of $`\alpha =0.57`$ for the shell emission. Fig. 5 shows the resulting total spectrum when adding both components. The differences to the spectral fit of the integrated flux densities are quite small for the present frequency range, but would become larger at higher frequencies.

### 3.4 The polarization properties

At 32 GHz a rotation measure (RM) of 2000 rad/m2 is required to rotate the polarization angle by about 9°. Therefore we can ignore Faraday rotation effects at 32 GHz. This means that the magnetic field responsible for the synchrotron emission from the shell of G11.2$``$0.3 has a radial orientation. The magnetic field structure is characteristic for a young SNR in which the ejected material, with its frozen–in radial magnetic field, dominates the expansion and the synchrotron emission.

The integrated linearly polarized flux density (PI, see Table 2) is about 2–3 %, which is quite low for a SNR. Nevertheless, at the three wavelengths where polarization was observed PI is strong enough to derive rotation measures (RM). At the common angular resolution of 2$`\stackrel{}{.}`$5 we calculated RMs for three components as marked in Fig. 7.

The southern and eastern component have low RMs, while the northern blob shows about $``$1530 rad/m2. The 4.75 GHz data suffers from ambiguities, but not the higher frequency data. Therefore a map of the lowest possible RM was calculated between 10.45 GHz and 32 GHz at 1$`\stackrel{}{.}`$15 angular resolution. We assumed the polarization angles at 32 GHz to be correct. At 10.45 GHz we checked the original angle $`\varphi `$, $`\varphi +\pi `$, and $`\varphi \pi `$. The angle which resulted in the lowest RM was taken. The RM map is displayed in Fig. 8. The PI map at 10.45 GHz shown in Fig. 9 is corrected for RM effects. The SNR clearly has a radial magnetic field except for the southeastern quadrant where the polarized intensity is very low. Otherwise the map is fully in agreement with the 32 GHz map, which is mostly unaffected by even high RMs. The reason for the excessive RM in the northern part of the SNR is unknown. Either a magnetic field enhancement or an enhanced electron density within this part of the remnant or in its immediate surroundings is required. The much lower RM observed towards the south and east of the SNR rules out a significant contribution from Faraday rotation within the interstellar medium along the line of sight between the SNR and us.

## 4 Discussion

A spectral index of $`\alpha =0.50`$ from the fit to the integrated flux densities of G11.2$``$0.3 indicates a spectrum flatter than all the other early phase SNRs. Their spectral indices as listed by Green (snrcat (2000)) are:

* Cas A: $`\alpha =0.77`$
* Kepler: $`\alpha =0.64`$
* Tycho: $`\alpha =0.61`$
* SN1006: $`\alpha =0.60`$

When taking into account the flat spectrum core the spectrum of the shell steepens to about $`\alpha 0.57`$. Since G11.2$``$0.3 is older than the four listed historical SNRs the spectral difference may reflect an evolutionary effect from free expansion to adiabatic expansion. It is still not fully understood what exactly happens during that transition phase and in particular under what circumstances it happens.

In the subsequent discussion we refer to six different components of the SNR. The shell of G11.2$``$0.3 is separated into the four prominent sectors as indicated in Fig. 10. Each sector has a width of $`90\mathrm{°}`$ starting with Sector I from $`0\mathrm{°}`$ (straight north) in clockwise direction. Within the SNR we refer to unresolved sources visible at 32 GHz, which are also indicated in Fig. 10. Properties of the individual components are listed in Table 3.

### 4.1 The structure of the shell

All the shell sectors look very similar in total intensity. As already mentioned by Green et al. (green (1988)), the structure of the shell appears rather clumpy without a sharp outer boundary. An exception is Sector III, which has a higher peak flux and is well defined by an outer boundary in contrast to the others (see Fig. 2). This indicates a later phase of evolution for this part of the SNR. This assumption is supported by the very low percentage polarization in Sector III, while all the other sectors show polarization between 5 and 7 % (Table 3).

In early stages of SNR evolution most of the expanding material is ejecta from the supernova explosion itself. A radial magnetic field is frozen into this material. In the free expansion phase the ejecta with the radial magnetic field dominate the hydrodynamics and the synchrotron emission. Since a radial magnetic field is observed in the shell of G11.2$``$0.3 we conclude that the SNR is in this phase of evolution. In the later development more ambient material is swept up by the expanding shock wave. This material is compressed and collected within a thin shell behind the shock wave, outside the ejecta. The ambient magnetic field component frozen into this shell is perpendicular to the expansion direction which means the expansion direction which means perpendicular to the magnetic field in the ejecta. The percentage polarization will decrease until the contribution of the swept up magnetic field to the synchrotron emission equals that of the ejecta. At this point we expect no linear polarization unless the observations are able to resolve these structures, which is not the case here. Sector III seems to be exactly at that stage while the other sectors are still in an earlier stage of evolution.

This advanced evolutionary phase of Sector III is also implied by the X–ray observations. The soft X–ray luminosity of a SNR depends on the amount of emitting mass which is of course increasing with time while the expanding shock wave sweeps up interstellar matter. The soft X–ray luminosity also depends on the emissivity of this material in the soft X–ray band which mainly depends on its temperature. In early phases the temperature exceeds $`10^8`$ K by a factor of five or even more. Such a high plasma temperature peaks in hard X–rays. During the free expansion phase the soft X–ray luminosity increases with decreasing temperature and an increase of emitting material. In the adiabatic expansion phase the temperature changes are rather small. Thus the soft X–ray luminosity increases just depending on the increasing mass in the expanding shell. In the free expansion phase the soft X–ray luminosity increases much faster than in later evolutionary phases of the SNR. This gives a very good indication of evolutionary differences within the remnant.

Sector III shows the strongest emission in the ROSAT HRI image (see Fig. 1). This implies that Sector III is the most developed part of the SNR. Apparently the ‘aging’ of the remnant seems to increase from Sector I to Sector III, which is likely due to different densities of the ambient medium. The required density increase points away from the Galactic Plane. This is either a local effect or the density peak at about 5 kpc distance is below $`b=0.3\mathrm{°}`$. However, no such density gradient is visible in the distribution of H i from the channel maps of the Effelsberg H i-survey by Braunsfurth & Rohlfs (brau (1984)) at about 9′ angular resolution. Therefore any density contrast has to be small. This assumption is supported by the nearly circular shape of the SNR, while a density variation would lead to different expansion velocities.

X–ray observations indicate that between 3 and 4 M of interstellar material has been swept–up so far (Vasisht et al. vasisht (1996); Reynolds et al. reynolds (1994)). A tangential magnetic field is frozen into this material, which interacts with the ejecta by forming a convective zone between both layers. In this zone material from both layers and their magnetic field lines are mixed by turbulence. Synchrotron emission is expected from all these zones. Since linear polarization with a radial magnetic field is observed, the mass of the ejecta should be considerably larger than the mass of the swept up material. Only a type II supernova explosion will eject such a large amount of material. The pulsar found within the remnant supports this conclusion. There is no indication for a stellar wind bubble around the SNR from the H i data (Braunsfurth & Rohlfs brau (1984)). This suggests an early B-type progenitor star with a mass of about 10–15 M. These stars have no high velocity winds or high mass loss rates. A supernova explosion of a B0 or a B1 star would fit all observations.

### 4.2 The interior structures

Torii et al. (torii1997 (1997)) detected an X–ray pulsar in the geometrical center of G11.2$``$0.3 within the hard spectrum X–ray synchrotron nebula previously discovered by Vasisht et al. (vasisht (1996)). The exact position of these components and therefore their physical relation to the central radio sources cannot be determined very precisely. The positional error is about $`\pm 1\mathrm{}`$, which makes it difficult to identify the synchrotron nebula of the pulsar with either of the three radio components within the remnant. However, it is most likely that the pulsar coincides with the central component visible in the ROSAT-HRI image (Reynolds et al. reynolds (1994)). In Fig. 1 a greyscale representation of this observation is shown together with white contours indicating the 32 GHz measurement in total power. Evidently the most likely radio counterpart of the hard X–ray synchrotron nebula is the central source (Source I) with an inverted spectrum. The flux density listed in Table 3 is, however, confused with a contribution from the western component showing up at 1.4 GHz (Fig. 4).

With the correlation between a pulsar’s rotational energy loss rate $`\dot{\mathrm{E}}`$, the radio surface brightness $`\mathrm{\Sigma }`$ and the diameter D of its synchrotron nebula found by Kothes (1998) we can estimate the flux density of this nebula in G11.2$``$0.3. With $`\dot{\mathrm{E}}=8.810^{36}`$ erg/s and a mean diameter of $`1\mathrm{}`$ at a distance of 5 kpc we get a flux density of about 1 Jy at 1 GHz. Assuming a spectral index of $`\alpha _{\mathrm{plerion}}=0.0`$ for this nebula we get $`\alpha _{\mathrm{shell}}=0.57`$ (see Fig. 5) for the shell.

The total flux density of the central component is much larger than the integrated flux density of the compact sources (Sect. 3.2.). Most of the flat spectrum emission is therefore diffuse. The compact polarization feature seen at 32 GHz (Fig. 3) is slightly offset to the north of the compact structures and therefore related to the diffuse emission. Its RM is about $``$800 \[rad/$`\mathrm{m}^2`$\], rather similar to that of the compact sources I and II. The central RMs follow the gradient of RM across G11.2$``$0.3 as seen in Fig. 8.

We note some similarities of G11.2$``$0.3 with the crab–like SNR G21.5$``$0.9, where a compact X–ray source (Slane et al. slane (2000)), which has no radio counterpart at 22 GHz, is seen at the center of axisymmetric radio structures suggesting some collimated outflow of particles (Fürst et al. fuerst (1988)). Certainly the angular resolution of the present data does not allow a decision as to whether this scenario also holds for G11.2$``$0.3. The possibility of confusion of the central emission with extragalactic background sources can be ruled out as the high resolution maps by Green et al. (green (1988)) and Becker et al. (becker1985 (1985)) resolve the eastern and western sources of the NVSS map into diffuse extended emission.

## 5 Conclusions

Our high frequency observations confirm the historical SNR G11.2$``$0.3 as a composite SNR. The structure of the centrally peaked flat spectrum component, which is powered by electron and magnetic field injection by the central pulsar, is partly resolved. At 32 GHz we have an indication of a central source with an inverted spectrum, which likely coincides with the pulsar. The two outer symmetrical components were previously resolved at high angular resolution and might result from a collimated outflow. However, this needs to be confirmed by more sensitive data. A strong gradient in RM across G11.2$``$0.3 is observed with highest absolute RM values in the northern part of the SNR. The SNR shows a radial magnetic field over most of the shell structure. The shell was found to be in the transition phase from free to adiabatic expansion, with differences between the four sectors of the shell. This suggests that the interacting interstellar material likely has different densities. Combining available X–ray and radio data we conclude that G11.2$``$0.3 is likely the remnant of a type II supernova explosion with an early type B progenitor star.

###### Acknowledgements.

RK acknowledges the support from Prof. Richard Wielebinski and the MPG through a PhD and a post-doctoral fellowship at the MPIfR, when most of the observations were made. This research made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. We wish to thank Ernst Fürst, Lloyd Higgs, and Tom Landecker for critical reading of the paper. We are grateful to the referee John Dickel for helpful comments.


text:
# Small x physics and the initial conditions in heavy ion collisions

## 1 Introduction

Why is small $`x`$ physics relevant to heavy ion collisions? The canonical picture of how matter is produced and evolves in a high energy nuclear collision was articulated by Bjorken . The “valence” components of the nuclear wavefunction (partons carrying a large fraction, per nucleon, of the nuclear momentum) interact weakly and populate primarily the fragmentation regions at large projectile and target rapidities. It is the “wee” virtual excitations of the nuclear ground state, gluons (mostly) and sea quarks, carrying a small fraction of the nuclear longitudinal momentum, that populate the central rapidity region, eventually producing the $`1000`$ hadrons that were measured in one unit of pseudorapidity at RHIC. If the nucleus is viewed as a superposition of Fock states, these wee partons belong to a state with a large number of gluons and sea quarks–each of these can therefore have only a small fraction of the nuclear momentum . The higher the energy of the nucleus, the larger the number of wee partons. In a nuclear collision, this Fock state is “destroyed” and the wee partons go on–shell creating the large energy densities in the relatively large space–time volumes associated with quark gluon plasma formation. Understanding the possible formation and dynamical evolution of this partonic state therefore requires understanding small $`x`$ physics at a deep level.

In the perturbative  QCD approach to small $`x`$ physics, the “cloud” of gluons is generated by Bremsstrahlung of a laddder of soft gluons off harder “valence” partons. The parton density grows rapidly with decreasing $`x`$ at a rate which, if extrapolated to $`x0`$, would violate the Froissart (unitarity) bound. This seems unlikely because at sufficiently small $`x`$ the Bremsstrahlung approach must break down. When the density of partons becomes large, parton recombination and screening effects become important, and may perturbatively “saturate” the rapid growth of parton distributions at small $`x`$ .

Since the density of partons is large, one can formulate the small $`x`$ problem as a classical effective field theory (EFT)  111The approach described here is by no means the only approach. For a review of the recent status of different approaches to small $`x`$ physics, see Ref. .. The small $`x`$ partons described by the EFT constitute a color glass condensate (CGC) . They are a glass because their properties are analogous to a spin glass. They form a condensate because they belong to a state with high occupation number of order $`1/\alpha _S`$.

The CGC contains a dimensionful scale, the saturation scale $`\mathrm{\Lambda }_s(x,A)`$. It is proportional to the gluon density per unit area and grows as $`\mathrm{\Lambda }_s^2A^{1/3}/x^\delta `$, where $`\delta 0.20.3`$ is the power that controls the rise in parton distributions in a nucleon. If $`\mathrm{\Lambda }_s\mathrm{\Lambda }_{QCD}`$, weak coupling methods are applicable and one can compute the parton distributions in the classical EFT. At RHIC, one can estimate $`\mathrm{\Lambda }_s=2`$ GeV, while at LHC, $`\mathrm{\Lambda }_s=4`$ GeV –weak coupling methods are therefore only marginally applicable at the former .

This talk is organized as follows. In section 2, we discuss the color field of a single nucleus. We will advertise the impressive recent theoretical developments in this approach. We will also discuss the striking phenomenological success of saturation models in describing the HERA data. In the following section, we will discuss the application of these ideas to nuclear scattering. Results from non–perturbative numerical simulations of the initial stages of nuclear collisions will be presented. In section 4, we will briefly confront saturation models against the recent RHIC data. We end in section 5 with an outlook describing relevant projects in progress.

## 2 The color field of a nucleus: QCD is a color glass condensate at high energies

The classical EFT approach to small $`x`$ physics has been discussed extensively in the literature . It contains a dimensionful scale $`\mathrm{\Lambda }_s`$ equal to the color charge squared per unit area of the valence sources. The EFT, formulated in the infinte momentum frame, in light cone gauge, is a (nearly) two-dimensional theory with a structure mathematically analogous to that of a disordered system of Ising spins in a random magnetic field . That’s where the color glass comes from 222Unlike the original model of Parisi and Sourlas, which fails in describing real spin glasses, this theory, on account of its peculiar quantum structure, may indeed be a glass. This remains to be proven though.. Further, since the occupation number of the gluonic fields is large–the field strength squared $`F_{\mu \nu }^21/\alpha _S`$–what we have is a color glass condensate (CGC) . A virtual photon (from an energetic leptonic probe) would therefore indeed “see” a colored glass on the time scales of its interaction.

The classical theory is solvable and classical correlation functions have been computed . The gluon number has the Weizsäcker–Williams distribution $`\mathrm{\Lambda }_s^2/k_t^2/\alpha _S`$ for $`k_t>>\mathrm{\Lambda }_s`$ and (due to an infinite resummation in $`\mathrm{\Lambda }_s/k_t`$) saturates as the distribution $`\mathrm{ln}(\mathrm{\Lambda }_s/k_t)/\alpha _S`$. This suggests that the typical momentum of the gluons is peaked around $`\mathrm{\Lambda }_s`$. Recently, Lam and Mahlon have done a lot of work elucidating the nature of the classical distributions when strict color neutrality is imposed .

When quantum corrections are included, a Wilsonian renormalization group (RG) picture emerges. The structure of the classical field remains intact while the scale $`\mathrm{\Lambda }_s`$ grows. This is because including hard gluons (harder than the $`x`$ scale of interest) in the source increases the typical size of color fluctuations. The scale $`\mathrm{\Lambda }_s(x)`$ is then given by the gluon density per unit area,

$$\mathrm{\Lambda }_s^2(x)=\frac{A^{1/3}}{x^\delta }.$$
(1)
Here $`\delta (\mathrm{\Lambda }_s)`$ denotes the power of the rise in the nucleon gluon distribution. This equation must be solved self–consistently for $`\mathrm{\Lambda }_s`$ at each $`x`$.

In particular, the weight function that gluonic correlation functions have to be averaged over satisfies a non–linear RG equation . There has been very important work recently in understanding and solving this RG equation . (The correlation functions obtained from this functional equation are identical to those obtained independently in different approaches by Balitsky and Kovchegov . For numerical solutions of these, see Ref. .) The equations have been solved in a mean field approximation–for large transverse momenta ($`k_t>>\mathrm{\Lambda }_s`$) the structure of the source is the simple Gaussian one of the McLerran–Venugopalan model, while at small transverse momenta, ($`k_t<<\mathrm{\Lambda }_s`$) it is entirely different. Much work remains to be done to uncover the complex structure of these equations–their phenomenological implications for heavy ion collisions and for deeply inelastic scattering  are tremendous.

Simple saturation models which incorporate the scale $`\mathrm{\Lambda }_s`$ have been applied to study the HERA data. With only three parameters, they fit the entire range of the HERA inclusive data (for $`x<0.01`$ and for all $`Q^2`$ upto $`450`$ GeV2). The model, with the same parameters, reproduces the HERA diffractive and (with some caveats) vector meson production data .

## 3 Classical approach to nuclear collisions: shattering the color glass condensate

In a nuclear collision, the color glass shatters producing a large number of gluons. The full quantum description of how this happens is a very difficult problem which has no solution yet. The standard pQCD factorization (“mini-jets”) approach to particle production is not of much help here because it is most unreliable at small $`x`$ and small $`Q^2`$. Inevitably, the better pQCD based models of nuclear scattering include saturation ideas in some form .

However, one can make quite a bit of headway in the classical EFT approach. Since the classical fields of the two nuclei are known, they provide the initial conditions for the evolution of the gluon fields produced in the nuclear collision. The small $`x`$ fields are described by the classical Yang-Mills equations

$$D_\mu F_{\mu \nu }=J_\nu $$
(2)
with the random sources on the two light cones: $`J_\nu =_{1,2}\delta _{\nu ,\pm }\delta (x_{})\rho _{1,2}(r_t).`$ The two signs correspond to two possible directions of motion along the beam axis $`z`$. As shown by Kovner, McLerran and Weigert (KMW) , low-$`x`$ fields in the central region of the collision obey sourceless Yang-Mills equations (this region is in the forward light cone of both nuclei) with the initial conditions in the $`A_\tau =0`$ gauge given by

$$A^i=A_1^i+A_2^i;A^\pm =\pm \frac{ig}{2}x^\pm [A_1^i,A_2^i].$$
(3)
Here the pure gauge fields $`A_{1,2}^i`$ are solutions of (2) for each of the two nuclei in the absence of the other nucleus.

In order to obtain the resulting gluon field configuration at late proper times, one needs to solve (2) with the initial condition (3). Since the latter depends on the random color source, averages over realizations of the source must be performed. KMW showed that in perturbation theory the gluon number distribution by transverse momentum (per unit rapidity) suffers from an infrared divergence. A reliable way to go beyond perturbation theory is to re-formulate the EFT on a lattice by discretizing the transverse plane. (Boost invariance in pseudorapidity is assumed, which reduces the theory to a $`2+1`$–dimensional theory.) The resulting lattice theory can then be solved numerically 333For alternative analytical approaches, see the work of Balitsky  and of Kovchegov .. We shall not dwell here on the details of the lattice formulation, which is described in detail in Ref. .

We have simplified the problem considerably since $`\mathrm{\Lambda }_s`$ and the linear size $`L`$ of the nucleus are the only parameters in the problem. We can write any dimensionful physical quantity $`q`$ as $`\mathrm{\Lambda }_s^df_q(\mathrm{\Lambda }_sL)`$, where $`d`$ is the dimension of $`q`$. All the non-trivial physical information is contained in the dimensionless function $`f_q(\mathrm{\Lambda }_sL)`$ . We can estimate the values of the product $`\mathrm{\Lambda }_sL`$ which correspond to key collider experiments. Assuming Au-Au collisions, we take $`L=11.6`$ fm (for a square nucleus!) and estimate the standard deviation $`\mathrm{\Lambda }_s`$ to be 2 GeV for RHIC and 4 GeV for LHC . Also, we have approximately $`g=2`$ for energies of interest. The rough estimate is then $`\mathrm{\Lambda }_sL120`$ for RHIC and $`\mathrm{\Lambda }_sL240`$ for LHC.

Our simulations were performed for the SU(2) gauge group–all results for physical quantities are based on extrapolations to SU(3). A straightforward computation is of the initial energy per unit transverse area per unit rapidity, deposited in the central region by the colliding nuclei . At late times, this is given by the non–perturbative formula

$$\frac{1}{\pi R^2}\frac{dE_t^i}{d\eta }|_{\eta =0}=\frac{1}{g^2}f_E\mathrm{\Lambda }_s^3,$$
(4)
where $`f_Ef_E(\mathrm{\Lambda }_sR)=0.21`$$`0.26`$. Using this formula, and assuming, in accordance with Ref. , the $`(N_c^21)/N_c`$ dependence of the energy on the number of colors $`N_c`$, we arrive at the values of 2700 GeV and of 25000 GeV for the initial transverse energy per unit rapidity at RHIC and at LHC, respectively . We can also estimate a formation time $`\tau _f`$ in our picture. For RHIC (LHC), we have $`\tau _f=0.3`$ fm ($`0.13`$ fm). This then enables one to estimate the initial energy density to be $`66.5`$ GeV/fm3 ($`1.3`$ TeV/fm3) at RHIC (LHC).

The number and distribution of produced gluons are of considerable interest as initial data for possible evolution of the gluon gas towards thermal equilibrium . Strictly speaking, a particle number is only well-defined in a free-field theory, and there is no unique extension of this notion to a general interacting case. For this reason, we use two different generalizations of the particle number to an interacting theory, each having the correct free-field limit. We verify that the two definitions agree in the weak-coupling regime corresponding to late proper times in the central region . Our first definition is straightforward. We impose the Coulomb gauge condition in the transverse plane: $`\stackrel{}{}_{}\stackrel{}{A}_{}=0`$ and and determine the momentum components of the resulting field configuration. Our second definition is based on the behavior of a free-field theory under relaxation. This “cooling” technique allows a gauge invariant determination of the total particle number. Unfortunately, it cannot be used to determine the momentum distributions leaving the Colomb gauge determination as the only one available.

We obtain

$$\frac{1}{\pi R^2}\frac{dN}{d\eta }|_{\eta =0}=\frac{1}{g^2}f_N\mathrm{\Lambda }_s^2,$$
(5)
where $`f_Nf_N(\mathrm{\Lambda }_sL)=0.13`$ ($`0.15`$) for RHIC (LHC) energies. The agreement between the cooling and Coulomb techniques at the larger values of $`\mathrm{\Lambda }_sL`$ relevant for RHIC and LHC is excellent. It is not as good at the smaller values: in general, the cooling number is more reliable. Naively extrapolating our results to SU(3), we find for $`Au`$-$`Au`$ central collisions at RHIC energies, $`dN/d\eta 950`$ for RHIC and $`dN/d\eta 4300`$ for LHC energies.

Our results for the number distribution, computed in Coulomb gauge, are as follows. We have verified that for large $`k_{}`$ our numerically obtained multiplicity agrees with the lattice analogue of the perturbative expression. At smaller $`k_{}`$, the distribution softens and converges to a constant value, unlike its perturbative counterpart. Notably, this qualitative change of the distribution occurs at $`k_{}\mathrm{\Lambda }_s`$. We tried to quantify this non-perturbative behavior by fitting the distribution to a variety of physically motivated functional forms. Surprisingly, we find that the shape of the distribution is closely reproduced by the two dimensional Bose-Einstein form $`n(k)=A/(\mathrm{exp}(\beta \omega _k)1)`$, with the inverse temperature $`\beta `$ of the order of 1 in units of $`\mathrm{\Lambda }_s`$, and with $`\omega _k`$ corresponding to a free massive dispersion relation, with the mass of the order of $`0.1\mathrm{\Lambda }_s`$. This is a remarkable result for a purely classical theory, whose meaning, beyond providing us with a convenient parametrization, is not yet clear.

We have studied the number distribution at different proper time slices. Our numerical results show clearly that the size of the mass gap decreases as the square root of the proper time–this is what one would expect for instance for a screening mass. Further, though the peak of the distribution softens, its slope remains unchanged–the inverse two dimensional “temperature” $`\beta `$ is a constant as a function of proper time. Thus the particle number is essentially unchanged–it increases slightly due to the “decay” of the mass gap.

Interestingly, all of these results (as well as the $`\mathrm{\Lambda }_s/k_t`$ form of the distribution for $`k_t<\mathrm{\Lambda }_s`$) are consistent with the theoretical analysis of Baier et al.  as communicated to us by Son . Baier et al. have studied the importance of $`23`$ processes in driving the system to equilibrium. The particular feature of their work, distinguishing it from previous ones, is the inclusion of the Landau–Pomeranchuk–Migdal (LPM) effect that results in the suppression of hard gluon radiation in QCD . In their kinetic theory analysis, for proper times $`\tau <1/\alpha _S^{3/2}`$, the occupation number $`f`$ is greater than unity. This is the overlap region between our classical simulations and kinetic theory and it is gratifying to see that they agree. This agreement further suggests that LPM interference effects are automatically included in the classical simulation.

What we have discussed above are the initial energy and number distributions of gluons produced in nuclear collisions. The subsequent evolution of the system (when the occupation number $`f`$ falls below unity) is beyond the scope of the classical approach–it does however provide the initial conditions for the evolution. The approach to equilibrium in the Mclerran–Venugopalan model was first investigated by Mueller . He initially considered only elastic scattering and derived a Landau–type equation for the single particle distributions. This equation has since been solved numerically and the approach to equilibrium studied quantitatively . One finds that elastic processes are very inefficient in driving the system to equilibrium (see also Ref. ). Baier et al. have argued that $`23`$ processes are not suppressed relative to $`22`$, and that they may actually be more efficient. For a more detailed discussion of thermalization, we refer the reader to their paper and to Son’s talk at this conference .

## 4 The CGC confronts RHIC data

Data from RHIC were presented at this conference and some of it has already been published. How does the CGC picture fare when confronted with the data? The predictions of various models for the RHIC multiplicity have been summarized nicely by Eskola in his QM2001 talk  (see also Ref. ). Our prediction in Ref.  does quite well. This however was fortuitous since from Eq. 5 it is clear that the multiplicity depends very sensitively on the saturation scale–the latter clearly has at least a 10% uncertainity.

A more sensitive test is the centrality dependence of the multiplicity. Kharzeev and Nardi plotted the charged particle multiplicity per participant pair at pseudorapidity $`\eta =0`$ versus the number of participants. They find that both the CGC and Glauber models give nearly identical predictions . These in turn agree well with the data .

A consequence of the CGC picture is that $`<p_t^2>dN/dy/\pi R^2`$, as is suggested by Eq. 5. This dependence has been verified recently for the E735 $`p\overline{p}`$ data at the Tevatron as well as the UA(1) data . It also works for the NA49 Heavy Ion CERN SPS data. First comparisions with the RHIC data are also promising .

In Ref. , we found that the momentum distributions of the produced gluons were nearly universal functions that scaled as $`ff(p_t/\mathrm{\Lambda }_s)`$. Whether the RHIC data reflects this can be easily checked because $`\mathrm{\Lambda }_s`$ should be different for different centralities. Juergen Schaffner–Bielich has checked that the RHIC data for charged particles satisfies this property by rescaling the $`p_t`$ scale by $`p_t\lambda p_t`$, where $`\lambda `$ was an apparently arbitrary scaling constant. What is truly remarkable however is that the scaling constant is not arbitrary–it corresponds closely to the ratio of the saturation scales at the two centralities. If this scaling constant is used instead to determine $`\mathrm{\Lambda }_s`$ at the different centralities and if then the number of charged particles per participant pairs at $`\eta =0`$ is plotted versus the number of participants (a la Kharzeev–Nardi), it agrees well with the data .

So far so good. However, the ratio of $`R_i=E_t/N_{\mathrm{gluon}}1.5\mathrm{\Lambda }_s`$ computed by us  is too large relative to the data . This is not necessarily a problem because $`R_i`$ is not a conserved quantity and depends on the complex dynamics of the space–time evolution. Our simulation only computes the initial value of this ratio. Both number changing processes a la Baier et al. and the expansion of the system should decrease this ratio. If the former is the dominant process, one might be throwing the baby out with the bath water because the $`\alpha _S`$ dependence of the multiplicity might change. A more mundane possibility is that the function $`f_E`$ which was computed in SU(2) might be significantly lower for SU(3). That’s one reason to press ahead with SU(3) simulations of the problem. As pointed out by Dumitru in his QM2001 talk, $`R_i`$ is a sensitive discriminant of the space–time dynamics in various models . Detailed space–time simulations are in progress .

Tests of the CGC scenario in the coming years include the energy dependence of $`\mathrm{\Lambda }_s`$, event by event fluctuations, and heavy quark production among others. We have recently computed the likelihood of the formation of CP–odd domains in the classical approach .

## 5 Outlook

The theoretical foundations of the classical approach to nuclear scattering clearly need to be on a firmer foundation. This is especially pressing in light of recent theoretical advances in small $`x`$ physics . There are several practical improvements that need to be done independently of the above. The extension to SU(3) and relaxing the strict boost invariance assumption are on the top of the agenda. A related interesting problem is: does the hot gluonic matter thermalize? At least asymptotically, this problem should have a solution. Recent attempts to address this problem are encouraging but much work (hopefully aided by new ideas) remains. Finally, as more and more RHIC data pour in over the years, we expect that they will provide a sensitive test of the ideas discussed here. We live in interesting times!