text:
# Development of Atmospheric Monitoring System at Akeno Observatory for the Telescope Array Project

## 1 Introduction

The air fluorescence technique for air shower observation has many advantages, (for example, direct measurement of shower longitudinal development, and stereo geometrical reconstruction with multiple eyes). These advantages will bring technical breakthroughs in EHE Cosmic Ray physics.

On the other hand, this fluorescence technique has a serious problem. In this technique the light yielded from air showers after the transmission is measured in atmosphere at a 10 $``$ 60 km distance. Atmosphere can be considered a part of the detector. Without any monitoring, the atmospheric conditions may cause large uncertainty in the experiment and result in serious systematic error. Therefore, the calibration of the light transmittance in the atmosphere is essential for the air fluorescence experiment.

Charged particles in the air shower excite air molecules and yield fluorescence light of between 330nm and 400nm. This fluorescence light is also scattered by molecules and aerosols in the air before it reaches the detector. The scattering process caused by molecule is called Rayleigh scattering. The attenuation length of this scattering process $`X_r`$ is 2974 $`g/cm^2`$ for a 400nm wave length. Amount of scattered photons with wave length of $`\lambda `$ can thus be calculated by following equation

$$\frac{dN_r}{dl}=\rho \frac{N_r}{X_r}(\frac{400nm}{\lambda })^4$$
(1)
where $`\rho `$ is the atmospheric density. In this equation, the number of scattered photons is proportional to $`\lambda ^4`$. This scattering process seriously effects the observation of fluorescence light in the air.

For example, when we consider the typical measurement of fluorescence light from an air shower at a distance of 30 km from detector, with a detector altitude of 1.5km and an elevation angle of 30 degrees, the fluorescence light has to pass through about 2200 $`g/cm^2`$ in the air and 70% of the photons will be scattered by this Rayleigh process.

There is another scattering process, Mie scattering. There are many kinds of scattering materials in the atmosphere, such as, water vapors, mists, clouds, dusts, blown small sands, artificial fumes, exhaust gas from cars, and smoke from forest fires. Furthermore, ’Kosa’ from the Gobi desert in China is spread over a wide area of the northern hemisphere. These atmospheric aerosol’s densities and their components change significantly with location and time. The scattering of photons by these aerosols is called Mie Scattering. In the Utah desert, the typical scale height of Mie scattering $`H_M`$ is 1.2 km and horizontal attenuation length $`L_M`$ is 20 km. The amount of scattering photons can be expressed by the following equation

$$\frac{dN_r}{dl}=\frac{N_r}{L_M}\times exp(\frac{h}{H_M})$$
(2)
Therefore, under the previously described conditions, approximately 27% of the photons contained in the fluorescence light from the air shower are scattered by the Mie process. When we take into account Rayleigh and Mie scattering, about 22% of the photons can reach our detector.

The Telescope Array (TA) will observe the air showers which flash far from the detector, for example, at more than 50km. Figure 1 shows the results of Monte-Carlo simulation for TA. According to this simulation study, the energy resolution of TA is better than 6% if atmospheric transmittance is accurately taken into account. However, if we apply the 20% shifted value in the attenuation length, the estimated cosmic ray energy has a systematic error of nearly 10%. To estimate the primary composition, TA will use $`X_{max}`$ which is the maximum position in the longitudinal shower development. The intrinsic resolution of $`X_{max}`$ is estimated to be 20g in the TA experiment using Monte Carlo simulation. If the assumed attenuation length differs by 20%, we estimate $`X_{max}`$ with a 10g systematic error. The correction of atmospheric transmittance affects not only energy determination but also air shower reconstruction. To realize this high resolution power of TA, we need to measure the attenuation length of the atmosphere with an accuracy of a few percent or better.

For the purpose of the atmospheric monitoring, TA and High Resolution Fly’s Eyes are developing various methods which use a laser or flashers. One of well-known technique for atmospheric monitoring is LIDAR which has been developed for environmental science. Using this technique, atmospheric conditions can be monitored remotely by measuring back-scattered light from a pulsed laser beam.

We are developing a steerable LIDAR system to study the atmospheric monitoring method in Akeno Observatory. In the present paper, we will describe the details of this system and discuss how to accurately measure transmittance in the atmosphere.

## 2 LIDAR equation

In general, the intensity of the back-scattered light detected by LIDAR is expressed by the following equation (LIDAR equation)

$$P(R)=\frac{P_0\kappa A(c\tau /2)\beta (R)T^2(R)Y(R)}{R^2}+P_b$$
(3)
where
$`P(R)`$ : intensity of the detected light
$`P_0`$ : laser intensity
$`P_b`$ : intensity of back ground photons
$`c`$ : light velocity
$`R`$ : distance from laser to target
$`\tau `$ : integration time
$`Y(R)`$ : geometrical efficiency of the beam track and receiver
$`A`$ : aperture of receiver
$`\kappa `$ : detector efficiency
$`\beta `$ : back-scatter coefficient
$`T(R)`$ : transmittance (=$`exp(_0^R\alpha 𝑑r)`$)
$`\alpha `$ : extinction coefficient (=1/attenuation length)
We define the following parameters for convenience

$$P_0\kappa A\frac{c\tau }{2}Y(R)C$$
(4)
$$X(R)R^2(P(R)P_b)$$
(5)
$$S(R)ln(X(R))$$
(6)
The LIDAR equation can be then written as follows:

$$X(R)=C\beta (R)T^2(R)=C\beta (R)exp(2_0^R\alpha 𝑑r)$$
(7)
If $`\alpha `$ is constant until $`R`$,

$$X(R)=C\beta exp(2R\alpha )$$
(8)
Then we obtain the differential LIDAR equation

$$\frac{dS(R)}{dR}=\frac{1}{\beta (R)}\frac{d\beta (R)}{dR}2\alpha (R)$$
(9)
There are two variable parameters in this LIDAR equation. One is the extinction coefficient $`\alpha `$ that represents the amount of scattered photons in a scattering volume. This value corresponds to the reciprocal of attenuation length. The other variable parameter is the back-scattering coefficient $`\beta `$ that represents the amount of scattered photons in the backward direction. This value is proportional to the density of scattering matter and its differential cross-section. Since the LIDAR equation has two variables, it can not be solved in general.

If the atmosphere is uniform, $`\beta `$ is constant and $`d\beta /dR=0`$. Then

$$\alpha =\frac{1}{2}\frac{dS}{dR}$$
(10)
Thus we can measure the $`\alpha `$ using this equation.

If we assume a relationship between $`\alpha `$ and $`\beta `$, the variable parameters can be reduced. For this purpose, an empirical relationship is proposed as follows:

$$\beta =const\alpha ^k$$
(11)
In this assumption, $`k`$ depends on the atmospheric condition and wave length of the laser. It is empirically estimated to be 0.6$``$1.0. Substituting this relationship into equation(9), we obtain

$$\frac{dS}{dR}=\frac{k}{\alpha }\frac{d\alpha }{dR}2\alpha $$
(12)
This is a Bernoulli-type equation.

Viezee adopted the boundary condition near the detector and obtained the following solution.

$$\alpha (R)=\frac{X(R)^{1/k}}{\frac{X(R_c)^{1/k}}{\alpha (R_c)}\frac{2}{k}_{R_c}^RX^{1/k}(r)𝑑r}$$
(13)
Since the second term of the denominator is negative, the denominator approaches zero with an increase of $`R`$. This solution often diverges, and is unstable.

Klett proposed a stable solution by adopting the boundary condition at the farthest point $`R_c`$, and obtained following solution

$$\alpha (R)=\frac{exp[(S(R)S(R_c))/k]}{\alpha (R_c)+\frac{2}{k}_R^{R_c}exp[(S(R)S(R_c))/k]𝑑r}=\frac{X(R)^{1/k}}{\frac{X(R_c)^{1/k}}{\alpha (R_c)}+\frac{2}{k}_R^{R_c}X^{1/k}(r)𝑑r}$$
(14)
In this solution, the integral term of the denominator becomes larger with R. Thus the contribution of uncertainty of the critical value $`\alpha (R_c)`$ becomes smaller. Therefore this solution converges on the correct value. Using Klett’s method, we can analyze $`\alpha `$, if we can assign certain values to the parameter ’$`k`$’ and the critical value $`\alpha (R_c)`$.

Furthermore, Fernald proposed a similar solution which considered atmospheric compositions. This solution assumes that $`k=1`$ and that scattering parameter $`S_a\alpha _M/\beta _M`$ is constant. Then $`\alpha `$ can be analyzed in successive steps using the Klett method. After the invention of these methods, analysis of the LIDAR signal based on elastic scattering made rapid progress.

## 3 Experiment

The steerable LIDAR system was developed at Akeno Observatory in Yamanashi Prefecture, Japan; its geographical coordinate is 900m a.s.l., $`138.5^{}`$N, $`35.78^{}`$E. This observatory utilizes are many cosmic ray detectors including AGASA and fluorescence detectors. Using these detectors, hybrid observation of cosmic rays are performed.

This LIDAR system is situated on the roof of the Akeno Main laboratory. All of the devices are housed in the astronomical dome shown in Figure 2. This astronomical dome consists of a cubic-shaped room of 1.8mW$`\times `$1.8mD$`\times `$1.0mH and a steerable dome. The dome window can open to 1m, and its direction must be rotated according to the movement of the laser. Two motors are equipped for opening and closing the dome and for dome rotation. Each motor is controlled by an electronic switch which operated by a local computer.

The steerable LIDAR system is illustrated in Figure 3. Inside the dome, laser and optical devices are installed on a optical table which is fixed on the concrete pad. This laser is a flash lamp pumped, Q-switched, and water cooled Nd:YAG type. Its specifications are shown in Table 1. The third harmonic beam with a wavelength of 355nm is used in this experiment. This wavelength is close to the major lines of air fluorescence light. The maximum beam intensity is 7 mJ and its maximum repetition rate is 10 Hz. The beam intensity, repetition rate, and shooting duration can be controlled by the computer through RS232C.

The laser beam is split into two directions after phase conversion by the circular polarizer. One has 10% intensity which is used for energy calibration, while the other 90% is transported to the shooting system. To stabilize the laser and to prevent condensation, these optical devices including the laser and the optical table are isolated from the outside atmosphere by a fireproof curtain and the temperature is kept constant by an air conditioner.

An alt-azimuth mount is adopted for the shooting system. Encoders are directly mounted on each axis. The resolution of these encoders is 5/1000 degrees. Each axis is driven by AC stepping motors with 0.0072 degree steps. The maximum speed of this shooting system is more than 10 degrees/sec.

An infrared camera is also mounted on the azimuthal axis. Using this infrared camera, clouds can be recognized as hot regions in the cool night sky. The typical cloud temperature is higher than a few degrees; however, the typical night sky temperature is lower than $`10`$ degrees. This infrared camera can measure the temperature between $`20`$ and 20 degrees and its field of view is 20 $`\times `$ 40 degrees. To increase the data rate, the video image from this camera is captured as a ppm format image file through a VIDEO capture board.

Two mirrors are mounted on the elevation axis. One is a shooting mirror for the laser beam. The other is a receiving mirror for back-scattered light. Maximum beam intensity shot to the sky is about 5 mJ. The receiving mirror is adjusted to parallel the laser beam. A parabolic mirror of 16 cm in diameter is adopted for this receiving system. Because the shooting and receiving mirrors are mounted on the same part, the laser beam and the receiver can be moved simultaneously. A one-inch PMT is located at the focus of 16cm diameter mirror. The field of view is about 1 degree. PMT gain is adjusted to about $`10^6`$. The signal from the PMT is acquired by a 100MHz band-width Digital Oscilloscope (TDS33014 made by Tektronix cooperation) which measures the time profile of the back-scattered light. The time range of the signal sweep is set to 10$`\mu `$sec/DIV and one sweep consists of 500 wards with a time resolution of 200nsec. Since the timing of the back-scattered light in the LIDAR system corresponds to the round-trip time to the target, 200nsec time resolution is equivalent to 30m in spatial resolution $`\mathrm{\Delta }R`$. When we set 1 bit to 0.4mV (it corresponds to 10mV/DIV with 256 bits digitizer), then 1 bit corresponds to about 2 photo-electrons. An average of 16 shots is calculated by this oscilloscope and recorded on the local computer.

All of these systems are controlled by the local computer using a Linux operating system and operated by the remote computer via a network.

## 4 Observation

The Telescope Array is planed to build 10 stations with 40km separation in the Utah desert USA. The location of those stations will be spread over a very wide area, their operation should be done remotely via a network linked by microwave or optic fiber. This would constitute a severe constraint for our atmospheric monitoring system. For this reason, we developed a fully automated system at Akeno observatory. For example, if it begins to rain during the observation, the system is shut down automatically.

The intensity of the back-scattered light reaching the detector from a scattering point is proportional to the inverse of the square of the distance, the density of the atmosphere, and the intensity of the laser beam. For this reason, the amount of returned light from near the detector is huge and that from a distant point is faint. Therefore a very wide dynamic range, typically 4 orders of magnitude, is required for the light receiving detector. If we use a high intensity laser beam, the PMT is saturated by photons scattered near the detector. If it is weak, the number of measurements have to be increased in order to measure the signal from a distant place, thus lengthening the observation time.

To increase the dynamic range, we use low and high sensitivity channels. Each channel uses two settings on the digital oscilloscope as illustrated in Figure 4. For the high sensitivity channel, the laser beam is shot at 5 mJ into the sky and oscilloscope ranges are adjusted to 1mV/DIV and 10mV/DIV. In this case, the PMT is saturated by scattered photons near the detector and it takes more than 9 $`\mu `$sec to recover sensitivity. Data points with signal amplitudes of less than 2 mV are not used in the analysis in order to obtain a reasonable signal to noise ratio, and to reduce the uncertainty of the pedestal. This value corresponds to about 1 p.e., i.e., 5 bits for 1mV/DIV. For the low sensitivity channel, the laser intensity is reduced to about 1/12 and the oscilloscope ranges are adjusted to 10mV/DIV and 200mV/DIV. In this case, the PMT recovers by 8 $`\mu `$sec. Using this channel, we can measure the back-scattered light from a distance between 1.2 km and 3 km.

The signal profiles of 16 shots are averaged by the digital oscilloscope and then recorded on the local computer. This measurement is repeated 48 times in vertical and horizontal directions and 16 times for every 5 degrees between 85 and 5 degrees. In other words, 1024 shots are measured in the vertical and horizontal directions, and 256 shots are measured in other directions.

Examples of measurements are shown in Figure 5. In this figure, the results in the vertical and horizontal directions are shown. In the horizontal measurement, if the atmosphere depends on only altitude ( meaning that the atmosphere has a 1-D structure and is uniform), $`S(R)log(X(R))`$ is in proportion to R and extinction coefficient $`\alpha `$ is analyzed in terms of the slope using equation(10). If the atmosphere is not uniform in the horizontal direction, the slope is not constant, depending on the atmospheric conditions. We can see S(R) decline linearly to 10 km under good weather conditions as shown in Figure. There is a mountain at a distance of 10km and the light scattered by this mountain can be seen in this figure. The atmosphere is not uniform under bad weather and S(R) does not decline linearly. In measurements in the vertical direction, it can be clearly seen that clouds cause the magnitude of the signal strength to vary rapidly. We can measure the atmosphere up to an altitude of 12km under a clear sky.

## 5 Analysis

In this section, we discuss how to estimate the transmittance of the atmosphere using the results of the measurements described above.

One of the simplest cases is horizontal measurement. If the atmosphere is uniform, $`\alpha `$ can be analyzed using equation(10) at the level of the detector as mentioned previously. Figure 6 shows examples of this horizontal measurement. Data which is comparatively near the detector is used. It can be seen that the attenuation of the laser beam clearly depends on the weather. It can also be confirmed that the atmosphere is uniform in this region. Figure 7 shows the time profile of $`\alpha `$ in comparison with humidity which is measured simultaneously. Of course, the transmittance of the atmosphere depends not only on humidity; however, in this figure we can see that variation of $`\alpha `$ follows tiny variations of humidity. This method appears to be very convincing and accurate.

As described section 2, the LIDAR equation cannot be solved because it has two variable parameters. A solution which uses an assumption of the empirical relationship between $`\alpha `$ and $`\beta `$ is proposed by Klett. To estimate $`\alpha `$ using Klett’s method, we need to determine two parameters: one being the critical value at the highest altitude point, the other being parameter $`k`$ which represents the relationship between $`\alpha `$ and $`\beta `$. To determine these parameters, we use a simple assumptions as follows:

* In a measurement in the vertical direction, the smallest $`\alpha `$ corresponds to the value expected from Rayleigh scattering. If there are no clouds, this value can be estimated at its highest point.
* In a measurement in a different direction, a smaller $`\alpha `$ (meaning the value in a cloudless region) corresponds to the $`\alpha `$ which is measured at a direction of +5 degrees at the same altitude.
* $`\alpha `$ at points lower than 100m can be estimated by the $`\alpha `$ obtained in the horizontal measurement.

In the first assumption, we estimate the critical value based on Rayleigh scattering. This means that we assume the atmosphere is clear and the effect of Mie scattering can be ignored at high altitudes, which is about 10 km with the present system. Transmittance is analyzed in successive steps from the critical point to the lower point. For this reason, errors in the critical value affect the measurement of the lower points. This assumption is reasonable during good weather conditions. However, if the weather is not clear, this assumption causes significant systematic error.

At small elevation angles, it became difficult to measure the transmittance at high altitude. In the second assumption, the measurement at larger elevation angle is used for estimation of the critical value of the smaller one. In other words, we assume uniform atmosphere at the same altitude for the measurements of adjacent elevations. If atmospheric conditions changes rapidly or are localized, it is impossible to estimate the alpha because of this assumption.

In the Telescope Array experiment, we expect that no observations will be performed or that the data will be discarded under cloudy conditions or if the atmospheric conditions in the field of view change rapidly. Under these conditions we estimate the observation’s duty cycle to be 7$``$10% at the Utah site. For this reason, the assumptions made in this analysis should be acceptable. Based on these assumptions, we can determine the parameter ’$`k`$’ and the critical value.

Figure 8 shows the examples of estimated $`\alpha `$ as a function of altitude at each zenith angle. As we can see, the atmosphere is almost uniform if there are no clouds. In the lower left panel, a cloud can be seen at around an altitude of 6 km. This cloud obscures the signal from higher altitudes. Since the boundary condition at high altitudes cannot be determined, the absolute value of $`\alpha `$ cannot be analyzed correctly. Figure 9 shows the results obtained on a very clear night. The atmosphere of this night can be almost explained by pure molecular at high altitude. Also the atmosphere was very uniform as is illustrated. It can be clearly seen that the statistical error is very small at all points. For example, the transmittance between the detector and a 10 km distance in the vertical direction is estimated to be 73.27% with a statistical error of 0.05%. Figure 10 shows the result of measurements taken on the next night. The weather was hazy and the atmosphere was not uniform at lower heights. The transmittance between the detector and 10km in the vertical direction is 57.91% with a statistical error of 0.07% It takes approximately 20 minutes to measure one azimuthal direction.

## 6 Systematic Error

The systematic error in the present analysis based on the Klett’s method has been investigated using a simple Monte-Carlo simulation . In this simulation, the signal profile of back-scattered light is calculated using a simple atmospheric model. The number of Rayleigh-scattered photons follows molecular density only, while Mie scattering is calculated by assuming Mie parameters: scale height $`H_M`$, horizontal attenuation length $`L_M`$, and back-scatter coefficient $`\beta _M`$. An example of the atmospheric model with typical Mie parameters ($`L_M`$=20km, $`H_M`$=1.2km, $`\beta _M`$=0.05) is shown in Figure 11. In this calculation, the laser’s elevation angle is assumed to be 10 degrees.

This artificial data is analyzed using the same programs as above, employing real data analysis. In Figure 12, examples of the results are shown. The expected alpha which is calculated based on the input atmospheric model and the estimated one using Klett’s method are compared. In Figure 12.a the standard atmospheric model is assumed with typical Mie parameters as in the previous figure.

In Figure 12.b, Mie parameters are the same as the standard model except for $`\beta _M`$ which is transformed from 0.1 to 0.05 at 2 km height. The Klett’s assumption of a power-law relationship between $`\alpha `$ and $`\beta `$ does not take into account variation of the aerosol particle types. If the aerosol particle size is large, forward-scattered light is predominated in Mie scattering. If it is small, the phase function of Mie scattering becomes similar to that of Rayleigh scattering. In this panel the total cross section is the same as the standard model, but the phase function of the Mie scattering changes rapidly at 2 km height. The changing of aerosol types can cause significant systematic errors in this analysis as we can see.

In Figure 12.c, all Mie parameters are the same as those of the standard model. But the critical value which is determined based on Rayleigh scattering at the furthest point is assumed to have a 10% systematic error. This critical value is estimated under the assumption that all of the backscattering at the highest point is Rayleigh scattering. In this case, this critical value can be estimated based on molecular density which depends on temperature, pressure and the type of atmospheric model used. In this panel, we assume that the critical value is overestimated beyond the true value by 10%. A similar situation can arise when we ignore the effect of Mie scattering to obtain the critical value under an unclear atmosphere.

In Figure 12.d, all Mie parameters are the same as in the standard model. But multiple scattered light is considered under the assumption that 10% of the observed backscattering light is multiple scattered light. A more complete study of the effect of multiple scattering lies outside the scope of this paper. There is no effect to the estimated value under this simple model.

Figure 13 shows the systematic error of the total attenuation length between the laser and measurement point as a function of height. According to these figures, the systematic error caused by Klett’s method is less than 5% at all of the measurement points in the standard atmospheric mode. Most significant effect is caused by variation of the aerosol type. These systematic errors should be dependent on the laser’s elevation angle, the Mie parameters, and the type of atmospheric model used. Further investigation of the simulation study is described in .

The detector constants (for example: mirror reflectivity, PMT gain, beam intensity, and so on) do not contribute to the systematic errors in this analysis except for the linearity of the PMT. The effect of multiple scattering is ignored in this measurement. The most significant systematic error is caused by the uncertainty of the critical value. The altitude of 10 km, where we established the boundary condition, may not be high enough to ignore the effect of Mie scattering. This problem can be solved if we use a larger mirror or a higher intensity laser.

## 7 Summary

We have developed a steerable LIDAR system for atmospheric monitoring by the Telescope Array. The system consists of a 5 mJ pulse laser and 16 cm diameter mirror. Using this system, a technique for atmospheric monitoring was developed.

First, the extinction coefficient $`\alpha `$ at the level of the detector is measured. Then $`\alpha `$ is estimated at all directions using Klett’s method. The transmittance of the night sky to a distance of more than 10 km can be measured successfully. The statistical error in this analysis is less than 2% under a clear sky. It takes approximately 20 minutes to measure one azimuthal direction.

The systematic error which is caused by this method of analysis is estimated to be approximately 0.5% in a typical atmospheric model. The detector constants do not contribute to the systematic errors in this analysis except for the linearity of the PMT. The most significant systematic error is caused by uncertainty of the critical value and variation of the aerosol type. Some of these problems can be solved if we use a larger mirror or a higher intensity laser.

Based on these observational results, we are considering future plans regarding atmospheric monitoring for the Telescope Array. The signal to noise ratio is in proportion to beam intensity, mirror diameter, square root of observation time, and transmittance. $`5mJ\times 16cm\times \sqrt{20min}=358[mJcmmin^{1/2}]`$ is required to measure to a distance of 10km under a clear sky (in this case, the transmittance is about 60% in the vertical direction). To measure more than 50 km, we need $`358\times 5^2/0.6^5=115000[mJcmmin^{1/2}]`$. This means that a 270 mJ laser is necessary for a 3 m diameter mirror if we want to measure one azimuthal direction in 2 minutes. This value can be reduced by a narrow band filter, small field of view, and high altitude location since the main component of the noise is night sky background light. For a telescope with a 1 degree field of view and F=1, the observed night sky background is about 5 MHz in the ultraviolet region.

For the operation of the Telescope Array detector, independent and redundant measurements of the atmosphere should be done for the sake of verification and cross-checking. In fact, several measurements being performed by the Telescope Array and the HiRes groups . For example, These include the measurement of the large angle scattered light of a polarized laser beam and radio-controlled xenon flashers using the HiRes detector, monitoring the transmittance of the air around the detector using the Nepherometer, measurement of the zenith angle dependence of star light and so on. By combining these measurements and the LIDAR method, more reliable measurements can be performed in the future.

This work is supported by grants in aid #12304012 and #11691117 for scientific research from the JSPS(Japan Society for the Promotion of Science). The authors thank Dr. Lawrence R. Wiencke, Dr. Michael E. Roberts, and Prof. John A.J.Matthews for their thoughtful discussions.


format: text

text:
# Creation of a dipolar superfluid in optical lattices

## Abstract

We show that by loading a Bose-Einstein condensate (BEC) of two different atomic species into an optical lattice, it is possible to achieve a Mott-insulator phase with exactly one atom of each species per lattice site. A subsequent photo-association leads to the formation of one heteronuclear molecule with a large electric dipole moment, at each lattice site. The melting of such dipolar Mott-insulator creates a dipolar superfluid, and eventually a dipolar molecular BEC.

The physics of strongly correlated systems in ultra-cold trapped quantum gases has attracted a growing interest recently. The experimental progress on the trapping and cooling techniques, and the control of the inter-atomic potentials by means of Feshbach resonances allow to analyze situations beyond the validity of the mean-field approach. In particular, the possibility to achieve 1D bosonic gases opens the pathways toward a gas of impenetrable bosons . Also, a rapidly rotating trapped bosonic gas should exhibit effects similar to the fractional quantum-Hall effect . The field of strongly correlated atomic systems concerns also ultra-cold Bose gases with large scattering lengths , and the Mott-insulator (MI) to superfluid (SF) phase transition  in cold bosonic gases in optical lattices . The latter has been recently observed in the remarkable experiment .

The physics of ultra-cold molecular gases is also a very active research area. Cooling and trapping of molecules can be achieved by means of buffer-gas techniques , by employing the dipolar moments of polar molecules , and by means of photo-association of already existing ultra-cold atomic gases . Jaksch et al.  have recently proposed that a bosonic atomic gas, which is in a MI phase with exactly two atoms per site in an optical lattice, can be used to obtain via photo-association a superfluid gas of homonuclear molecules, and eventually a molecular BEC.

The influence of dipole-dipole forces on the properties of ultra-cold gases has also drawn a considerable attention. It has been shown that these forces significantly modify the ground state and excitations of trapped dipolar BEC’s . In addition, since dipole-dipole interactions can be quite strong relative to the short-range (contact) interactions, dipolar particles are considered to be promising candidates for the implementation of fast and robust quantum-computing schemes . Several sources of cold dipolar bosons have been proposed, including atoms with large magnetic moments , and externally induced electric dipoles . Perhaps the most promising perspective in this sense is provided by ultra-cold gases of polar molecules with large electric dipole moments, which could eventually dominate the physics of the system .

This Letter is devoted to the analysis of the generation of ultra-cold polar molecules in an optical lattice. We consider the loading of the lattice by a superfluid of two different atomic species, and analyze the transition into a MI phase with exactly one atom of each species per site. This transition is followed by the formation of polar dimers by photo-association, or adiabatic passage. Finally, the quantum melting of the molecular Mott phase into a polar molecular superfluid takes place.

We consider a miscible gas of two trapped atomic bosonic species at very low temperature $`T`$. Although our calculations could in principle be employed for arbitrary atoms, we have assumed in the calculations below the particular mixture of $`{}_{}{}^{41}\mathrm{K}`$ and $`{}_{}{}^{87}\mathrm{Rb}`$, whose simultaneous Bose-Einstein condensation has been recently realized . We additionally consider that the binary gas is confined in a 2D optical lattice, similarly to the 3D one-component case recently observed in Ref. . The transverse dimension is considered as harmonically confined. The physics of the system is governed by the two-species generalization of the Bose-Hubbard (BH) model, described by the Hamiltonian:

$`H`$ $`=`$ $`{\displaystyle \underset{<i,j>}{}}[J_aa_i^{}a_j+J_bb_i^{}b_j]+U_{ab}{\displaystyle \underset{i}{}}n_{ai}n_{bi}`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}[U_{0a}n_{ai}(n_{ai}1)+U_{0b}n_{bi}(n_{bi}1)],`$ (2)
where $`a_i`$, $`b_i`$ are the annihilation operators of $`{}_{}{}^{41}\mathrm{K}`$ and $`{}_{}{}^{87}\mathrm{Rb}`$ atoms at the lattice site $`i`$, which occupy a state described by the Wannier functions $`w_a(𝐫𝐫_i)`$, $`w_b(𝐫𝐫_i)`$ of the lowest energy band, localized on this site. This implies that energies involved in the system are small compared to the excitation energies to the second band. We denote the position of the local minimum of the optical potential as $`𝐫_i`$, and the number operator of $`\mathrm{K}`$ ($`\mathrm{Rb}`$) atoms at the site $`i`$ as $`n_{ai}=a_i^{}a_i`$ ($`n_{bi}=b_i^{}b_i`$). In Eq. (1) only the nearest-neighbor tunneling is considered, which is described by the parameter $`J_{a,b}=w_{a,b}^{}(𝐫𝐫_i)[\frac{\mathrm{}^2}{2m}^2+V_l(𝐫)]w_{a,b}(𝐫𝐫_j)\mathrm{d}^3r`$, where $`j`$ and $`i`$ are the indices of neighboring sites, and $`V_l(𝐫)=_{\xi =x,y}V_\xi ^0\mathrm{cos}^2(k_\xi \xi )`$ is the optical lattice potential (with wavevector $`𝐤`$) which we consider to be the same for both species. The optical potential for $`\mathrm{K}`$ and $`\mathrm{Rb}`$ is almost the same for the detunings corresponding to the wavelength $`1064\mathrm{n}\mathrm{m}`$ of a Nd:Yag laser . The same species interact via a short range pseudo-potential, which is described by the coefficients $`U_{0a},U_{0b}`$: $`U_{0a,b}=4\pi \mathrm{}^2a_{\mathrm{K},\mathrm{Rb}}|w_{a,b}(𝐫𝐫_i)|^4/m_{\mathrm{K},\mathrm{Rb}}\mathrm{d}^3r`$, where $`a_\mathrm{K}`$, $`m_\mathrm{K}`$ ($`a_{\mathrm{Rb}}`$, $`m_{\mathrm{Rb}}`$) are the scattering length and atomic mass of potassium (rubidium). The interactions between different species are characterized by the coefficient $`U_{ab}=2\pi \mathrm{}^2a_{\mathrm{KRb}}|w_a(𝐫𝐫_i)|^2|w_b(𝐫𝐫_i)|^2/m_{ab}\mathrm{d}^3r`$, where $`a_{\mathrm{KRb}}`$ is the $`\mathrm{K}\mathrm{Rb}`$ scattering length, and $`m_{ab}`$ is the reduced mass for $`\mathrm{K}`$ and $`\mathrm{Rb}`$ atoms. The scattering length $`a_{\mathrm{KRb}}`$ has been estimated recently to be of the order of $`8.6\mathrm{nm}`$ , that is the repulsive interactions between $`\mathrm{K}`$ and $`\mathrm{Rb}`$ are stronger than those between $`\mathrm{K}`$-$`\mathrm{K}`$ and $`\mathrm{Rb}`$-$`\mathrm{Rb}`$, which are characterized by $`a_\mathrm{K}=3.17\mathrm{nm}`$ and $`a_{\mathrm{Rb}}=5.24\mathrm{nm}`$, respectively. As a result, it is energetically favorable for $`\mathrm{K}`$ and $`\mathrm{Rb}`$ particles to remain immiscible, i.e. to stay at different lattice sites, rather than in the same ones. This situation can be however modified by means of Feshbach resonances, which can strongly modify the value of $`a_{\mathrm{KRb}}`$. We have thus assumed that the value of $`a_{\mathrm{KRb}}`$ can be reduced to $`0.5\sqrt{a_\mathrm{K}a_{\mathrm{Rb}}}`$. In that case the two components acquire a miscible phase.

In the 2D calculations presented below, we have considered square lattices up to $`80\times 80`$ sites, with periodic boundary conditions. In the first stage, we demonstrate the possibility of performing a transition from a state with two miscible superfluids to a Mott insulator phase with exactly one $`\mathrm{Rb}`$ and one $`\mathrm{K}`$ atom per site. To this aim we first obtain the ground state of the system using a variational Gutzwiller ansatz .

We approximate the ground-state wave function by $`|\mathrm{\Psi }_{MF}=_i|\varphi _i`$, where the product is over all lattice sites. The functions $`|\varphi _i`$ for each site are expressed in the Fock basis, $`|\varphi _i=_{n,m=0}^{\mathrm{}}f_{n,m}^{(i)}|n,m_i`$, where $`n`$, $`m`$ indicate the occupation number for the $`\mathrm{K}`$ and $`\mathrm{Rb}`$ atoms, respectively. For numerical reasons we have assumed that maximally $`3`$ atoms of each species can be present at any lattice site, and checked the validity of this assumption self-consistently. The coefficients $`\{f_{n,m}^{(i)}\}`$ are found by minimization of $`\mathrm{\Psi }_{MF}|H_i[\mu _an_{ai}+\mu _bn_{bi}]|\mathrm{\Psi }_{MF}`$, where the chemical potentials $`\mu _a,\mu _b`$ were chosen so that the mean numbers of $`\mathrm{K}`$, $`\mathrm{Rb}`$ atoms were equal to the number of lattice sites.

The ground state for weak lattice potentials corresponds to the coexistence of superfluids for both atomic species, and it was used as the initial condition for subsequent time-dependent calculations. This state is characterized by a non-vanishing value of the superfluid order parameters $`|a_i|`$, $`|b_i|`$ (Fig. 1). To drive the system from superfluid to MI phase we slowly increase the lattice potential. This process can be described using a dynamical Gutzwiller method, as discussed in Ref. . To this aim one allows the coefficients $`\{f_{n,m}^{(i)}\}`$ to be time dependent and to follow the dynamics derived from the time dependent variational principle. The resulting equation becomes:

$`i\dot{f}_{n,m}^{(i)}`$ $`=`$ $`\left[{\displaystyle \frac{U_{0a}}{2}}n(n1)+{\displaystyle \frac{U_{0b}}{2}}m(m1)\right]f_{n,m}^{(i)}`$ (3)
$`+`$ $`J_a\left[\mathrm{\Phi }_i^a\sqrt{n+1}f_{n+1,m}^{(i)}+\mathrm{\Phi }_i^a\sqrt{n}f_{n1,m}^{(i)}\right]`$ (4)
$`+`$ $`J_b\left[\mathrm{\Phi }_i^b\sqrt{m+1}f_{n,m+1}^{(i)}+\mathrm{\Phi }_i^b\sqrt{m}f_{n,m1}^{(i)}\right]`$ (5)
$`+`$ $`U_{ab}nmf_{n,m}^{(i)},`$ (6)
where

$$\mathrm{\Phi }_i^a=\underset{}{}\mathrm{\Psi }_{MF}|a_j|\mathrm{\Psi }_{MF},\mathrm{\Phi }_i^b=\underset{}{}\mathrm{\Psi }_{MF}|b_j|\mathrm{\Psi }_{MF}.$$

When the lattice potential increases, the system follows the changes quasi-adiabatically, and the two species enter sequentially the Mott phase with one atom per site. As shown in Fig. 1 this occurs first for the heavier $`\mathrm{Rb}`$ atoms at $`t0.75\mathrm{s}`$, and then at $`t1.5\mathrm{s}`$ for $`\mathrm{K}`$ atoms. The ratio $`U_{0a}/J_a`$ ($`U_{0b}/J_b`$) at the transition point equals approximately $`23.5`$ ($`26`$), which interestingly compares very well with the expected value, $`23.2`$, for a 2D single-component gas. Therefore the transition point is rather unaffected by the interactions between different species.

At the end of the process each site contains one $`\mathrm{Rb}`$ and one $`\mathrm{K}`$ atom, with a relative atom number fluctuations less than $`3\%`$. This, however, requires a sufficiently slow modification of the lattice potential, typically within a time scale of 2 seconds for a lattice wavelength $`\lambda =1064\mathrm{n}\mathrm{m}`$. Note that the time scale is provided by the inverse of the recoil energy of the lattice, and therefore it scales as $`\lambda ^2`$. We note also at this point that a transition 3 times faster leads to about two times larger number fluctuation in the final state. We want to stress that during the whole time evolution the mean number of particles was constant, and equal to the number of lattice sites, which gives us the most complete analogy, within the Gutzwiller approach, to the experimental situation.

Having obtained a Mott phase with one Rb and one K atom per site, a heteronuclear molecule can be generated at each lattice site. One possible method would be two-color Raman photoassociation to make ground state $`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ or $`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ dimers. Since heteronuclear species do not have gerade/ungerade symmetry, levels of either symmetry could be formed via intermediate excited states of $`\mathrm{\Omega }=0^+`$ or 1 symmetry, where $`\mathrm{\Omega }`$ labels projection of electronic plus spin angular momentum on the molecular axis. Since the $`\mathrm{\Omega }`$ states excited at long range typically have mixed singlet and triplet character, photoassociation to ground $`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ levels is even possible for collision of two doubly spin-polarized atoms, which collide via the ground $`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ potential. The $`1/R^6`$ dependence of the excited molecular potential for the case of heteronuclear molecules, in contrast to the $`1/R^3`$ dependence for the homonuclear case, leads to excitation at relatively short internuclear separation $`R`$. This difference may lead to reduced Franck-Condon factors for the first exciation step but enhanced factors for the second de-excitation step . In any case, these can be worked out once the potential energy curves of the KRb dimer are better known. A second possible method of making ground state molecules is available for heteronuclear dimers, namely, use of a pulse of microwave radiation to directly associate two colliding atoms in a one-color transition from a scattering state to a bound dimer state. However, the dipole matrix elements within the $`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ state are small, so this will need fairly high microwave power. Either method will result in an excited vibrational level not far below the dissociation limit. It should be possible to use a succession of Raman pulses to transfer the population from such a level to the ground $`v=0`$ vibrational level of the $`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ or $`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ state, in a manner similar to that described by Ref. ; microwave/infrared pulses may also be feasible for this heteronuclear species. The $`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ $`v=0`$ level has a small dipole moment of $`0.0046`$ atomic units (1 atomic unit is $`8.478\times 10^{30}`$ Cm), whereas it is much larger, $`0.30`$ atomic units, for the $`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ $`v=0`$ level . In the latter case, special care should be taken to prevent the heating due to black-body radiation, which could be avoided by an appropriate shielding.

After the creation of the K-Rb dimers, the system becomes single-component. However, since the molecules have a large dipole moment and are considered to be oriented by an external electric field ($``$ $`100\mathrm{V}/\mathrm{cm}`$), we have to include the dipole-dipole interactions in our calculations. The corresponding Hamiltonian reads :

$`H`$ $`=`$ $`J{\displaystyle \underset{<i,j>}{}}b_i^{}b_j+{\displaystyle \frac{1}{2}}U_0{\displaystyle \underset{i}{}}n_i(n_i1)`$ (7)
$`+`$ $`{\displaystyle \frac{1}{2}}U_{\sigma _1}{\displaystyle \underset{<i,j>}{}}n_in_j+{\displaystyle \frac{1}{2}}U_{\sigma _2}{\displaystyle \underset{<<i,j>>}{}}n_in_j+\mathrm{},`$ (8)
where $`b_i`$ is the annihilation operator of a heteronuclear molecule at the lattice site $`i`$. The number operator $`n_i`$, the Wannier basis $`w(𝐫𝐫_i)`$ and the tunneling coefficient $`J`$ for molecules are defined similarly as in Eq. (1). However, the interaction part significantly differs from that of Eq. (1). The coefficients $`U_\sigma `$ are defined as follows: $`U_\sigma =|w(𝐫𝐫_i)|^2V_{int}(𝐫𝐫^{})|w(𝐫^{}𝐫_j)|^2\mathrm{d}^3r\mathrm{d}^3r^{}`$, where $`\sigma =|𝐫_i𝐫_j||𝐤|/4\pi `$ is the dimensionless intersite distance. In particular, $`U_0`$ determines the on-site interactions, $`U_{\sigma _1}`$ the nearest-neighbor interactions, $`U_{\sigma _2}`$ the interactions between the next-nearest neighbors, etc. Consequently, the respective summations in Eq. (7) must be carried out over appropriate pairs of sites which are marked by $`<\mathrm{}>`$ for the nearest neighbors, $`<<\mathrm{}>>`$ for the next-nearest neighbors, etc. In the 2D calculations presented below, we have taken into account interactions with up to $`4`$ neighbors ($`\sigma _1=1`$, $`\sigma _2=\sqrt{2}`$, $`\sigma _3=2`$, $`\sigma _4=\sqrt{5}`$), since the effects of interactions of a longer range are negligible. In the case of polarized dipoles the interaction potential is

$$V_{int}=d^2\frac{13\mathrm{cos}^2\theta }{|𝐫𝐫^{}|^3}+\frac{4\pi \mathrm{}^2a}{M}\delta (𝐫𝐫^{}),$$
(9)
where the first part of Eq. (9) provides the dipole-dipole interaction characterized by the dipole moment $`d`$ and the angle $`\theta `$ between the dipole direction and the vector $`𝐫𝐫^{}`$. The second part of Eq. (9) represents the short-range interactions given by the $`s`$-wave scattering length $`a`$ and the molecular mass $`M`$. We use the dipole moment as $`0.3`$ a.u. , and assume the $`s`$-wave scattering length of K-Rb molecule as $`a=\sqrt{a_\mathrm{K}a_{\mathrm{Rb}}}`$.

We have checked that the MI phase of these heteronuclear molecules is indeed the ground state of the system for a mean lattice filling factor equal to 1. We employ this state as the initial condition for our dynamical calculations of the melting of the Mott phase towards a molecular superfluid. It is, however, interesting to note, that by itself the MI phase of dipolar particles offers promising perspectives as a quantum computation device . The melting dynamics is analyzed by employing again a dynamical Gutzwiller ansatz, but with the Hamiltonian (7). During the dynamics we consider a first stage in which we reduce the lattice trapping potential up to $`4`$ in recoil energy units of $`\mathrm{K}\mathrm{Rb}`$ molecules. This first stage allows us to reduce the value of $`U_0/J`$, but the system remains still in the MI state mostly due to the large contribution of the repulsive dipolar on-site interactions to the $`U_0`$ coefficient. This process can be performed in a relatively fast time scale compared to the total duration of the melting process, since the MI gap guarantees the adiabaticity at this stage.

Later on, we have observed that the dynamics can be better controlled by modifying the transverse confinement. As shown in Ref. , the properties of the dipolar gas strongly depend on the aspect ratio of the on-site wavefunction. In particular, by reducing the transverse confinement, the on-site wavefunction becomes more elongated in transverse direction. As a consequence, $`U_0`$ decreases , and the system enters into the superfluid phase. Fig. 2 shows this process. Due to the non-completely adiabatic evolution the gas does not exactly follow the expected stationary result, although very clearly a superfluid molecular phase is accomplished. The ratio $`U_0/J23.2`$ at which the transition from MI to superfluid equals surprisingly to the value expected for a 2D single-component gas without a permanent dipole moment. As shown in Fig. 2 this dynamical transition into the superfluid phase occurs for realistic parameters in a time scale of $`1`$s. This transition can be realized, $`50`$ times faster, although the dynamics of the superfluid order parameter is in such case not as smooth as in Fig. 2.

In this Letter we have analyzed the formation of a superfluid of polar molecules. We have considered the loading of a BEC of 87Rb and 41K atomic species into an optical lattice. The modification of the lattice potential produces the dynamical transition into a MI phase with only one atom of each species per site. The atoms can then be photo-associated on-site. Once the polar molecules are created, a molecular MI is formed. By reducing the lattice potential and modifying the transverse confinement, such MI evolves into a molecular superfluid within a feasible time scale. Eventually, a molecular BEC is formed. The mechanism studied in this Letter, provides a perspective towards a dipolar BEC.

We acknowledge support from the Alexander von Humboldt Stiftung, the Deutsche Forschungsgemeinschaft, EU RTN Network ”Cold Quantum Gases”, the Polish KBN (grant no 2 P03B 124 22), ESF PESC Program BEC2000+ and the US Office of Naval Research. We thank G. Birkl, G. V. Shlyapnikov and E. Tiesinga for discussions.


text:
# 1 Vector diagram of the lie algebra 𝐶₂=𝑠𝑝(2,𝐑) (according to Cartan’s notation). It contains 10 vectors including the two null vectors in the center.

Heisenberg quantization for the

systems of identical particles and

the Pauli exclusion principle in

noncommutative spaces

S. A. Alavi

High Energy Physics Division, Department of Physics, University of Helsinki
and
Helsinki Institute of Physics, FIN-00014 Helsinki, Finland.
On leave of absence from : Department of Physics, Ferdowsi University of Mashhad, Mashhad, P. O. Box 1436, Iran
E-mail: $`Ali.Alavi\mathrm{@}Helsinki.fi`$
$`s_{}alialavi\mathrm{@}hotmail.com`$

1 Abstract. We study the Heisenberg quantization for the systems of identical particles in noncommtative spaces. We get fermions and bosons as a special cases of our argument, in the same way as commutative case and therefore we conclude that the Pauli exclusion principle is also valid in noncommutative spaces.

1 Introduction.

Recently there have been notable studies on the formulation and possible experimental consequences of extensions of the standard (usual) quantum mechanics in the noncommutative spaces \[1-13\]. Many physical problems have been studied in the framework of the noncommutative quantum mechanics (NCQM), see e.g. \[1-4\].
NCQM is formulated in the same way as the standard quantum mechanics SQM (quantum mechanics in commutative spaces), that is in terms of the same dynamical variables represented by operators in a Hilbert space and a state vector that evolves according to the Schroedinger equation :

$$i\frac{d}{dt}|\psi >=H_{nc}|\psi >.$$
(1)
we have taken in to account $`\mathrm{}=1`$ and $`H_{nc}H_\theta `$ denotes the Hamiltonian for a given system in the noncommutative space. In the literatures two approaches have been considered for constructing the NCQM :
a) $`H_\theta =H`$, so that the only difference between SQM and NCQM is the presence of a nonzero $`\theta `$ in the commutator of the position operators :

$$[\widehat{x}_i,\widehat{x}_j]=i\theta _{ij}[\widehat{x}_i,\widehat{p}_i]=i\delta _{ij}[\widehat{p}_i,\widehat{p}_j]=0.$$
(2)
b) By deriving the Hamiltonian from the Moyal analogue of the standard Schroedinger equation :

$$i\frac{}{t}\psi (x,t)=H(p=\frac{1}{i},x)\psi (x,t)H_\theta \psi (x,t),$$
(3)
where $`H(p,x)`$ is the same Hamiltonian as in the standard theory, and as we observe the $`\theta `$ \- dependence enters now through the star product . In , it has been shown that these two approaches lead to the same physical theory. For the Hamiltonian of the type :

$$H(\widehat{p},\widehat{x})=\frac{\widehat{p}^2}{2m}+V(\widehat{x}).$$
(4)
the modified Hamiltonian $`H_\theta `$ can be obtained by a shift in the argument of the potential  :

$$x_i=\widehat{x}_i+\frac{1}{2}\theta _{ij}\widehat{p}_j\widehat{p}_i=p_i,$$
(5)
which lead to

$$H_\theta =\frac{p^2}{2m}+V(x_i\frac{1}{2}\theta _{ij}p_j).$$
(6)
The variables $`x_i`$ and $`p_i`$ now, satisfy in the same commutation relations as the usual case :

$$[x_i,x_j]=[p_i,p_j]=0[x_i,p_j]=\delta _{ij}.$$
(7)
In this paper we study one of the most important problems in the quantum mechanics namely the Heisenberg quantization for the systems of identical particles and the Pauli exclusion principle in noncommutative spaces. It appears that the main difference is between one and two dimensions not between two and three dimensions and therefore we study only one and two dimensional spaces. Heisenberg quantization for the systems of identical particles in the framework of SQM has been studied in detailed in .The process of quantization in SQM consists of:
a) Identifying the classical observables and the Lie algebra which they form under poisson bracket.
b) Looking for the linear representations of this lie algebra such that classical observables $`A`$ and $`B`$ are represented by Hermitian operators $`A`$ and $`B`$ on some Hilbert space.
c)Replacing the Poisson bracket $`C=\{A,B\}`$ by the commutator :

$$\widehat{C}=\frac{1}{i}[\widehat{A},\widehat{B}]=\frac{1}{i}(\widehat{A}\widehat{B}\widehat{B}\widehat{A}).$$
(8)
The classical observables of a system of $`N`$ identical particles are the real valued functions on the $`N`$ particles phase space. Every observable must be symmetric as a function of the $`N`$ one-particle phase space variables, because if there would exist a nonsymmetric observable, it can be used to distinguish beween particles . Boson and fermion quantization are separate possibilities, because the symmetry of the quantum mechanical observables ensures that they do not mix symmetric and antisymmetric wave functions. In order to apply the Heisenberg quantization to a system of identical particles in NCQM, we follow the same way as SQM but we replace the commutator relations (7) by their noncommutative counterparts i.e. equ.(2).

2 Two particles in one dimensional spaces

In this section We apply the Heisenberg quantization to systems of two identical particles in commutative (usual) space. The fundemental relations are the cononical commutation relations :

$$[x,x]=[p,p]=0[x,p]=i.$$
(9)
To describe the two-particles systems, we introduce the relative coordinate and momentum :

$$x=x_{(1)}x_{(2)}p=\frac{1}{2}(p_{(1)}p_{(2)}),$$
(10)
where the subscript $`\mathrm{}=1,2`$ refers to the particles. If the two particles are not identical, the relative $`x`$ and $`p`$ of the two-particles system are observables and satisfy the same cononical relations as (9), but if particles are identical, then $`x`$ and $`p`$ are not observables any more, because they are antisymmetric under exchange of the particles indices. The observables can be :

$$A=\frac{1}{4}(p^2+x^2)B=\frac{1}{4}(p^2x^2)C=\frac{1}{4}(px+xp),$$
(11)
which are symmetric under the exchange of the particles indices, i.e. $`xx`$ and $`pp`$ . The observables $`A`$, $`B`$ and $`C`$ satisfy the $`sp(1,R)`$ algebra :

$$[A,B]=iC[A,C]=iB[B,C]=iA.$$
(12)
We can transform the above commutation relations in to the more familiar form by introducing the operators $`B_\pm `$ as follows :

$$B_\pm =B\pm iC.$$
(13)
Then we have :

$$[A,B_\pm ]=\pm B_\pm [B_+,B_{}]=2A.$$
(14)
We can define the irreducible representation of the algebra by introducing an orthonormal vector $`|\alpha _0,n>`$, $`n=0,1,2,\mathrm{}`$ as follows  :

$$\mathrm{\Gamma }|\alpha _0,n>=\alpha _0(\alpha _01)|\alpha _0,n>,$$
(15)
$$A|\alpha _0,n>=(\alpha _0+n)|\alpha _0,n>,$$
(16)
$$B_+|\alpha _0,n>=\sqrt{(n+1)(n+2\alpha _0)}|\alpha _0,n+1>,$$
(17)
$$B_{}|\alpha _0,n>=\sqrt{n(n1+2\alpha _0)}|\alpha _0,n1>,$$
(18)
where $`\mathrm{\Gamma }=A^2B^2C^2`$ is the Casimir operator and $`\alpha _0`$ is the minimum eigenvalue of $`A`$(an arbitrary constant which defines the representation). The cases $`\alpha _0=\frac{1}{4}`$ and $`\alpha _0=\frac{3}{4}`$ correspond to bosons and fermions respectively.

3 Two identical particles in two dimensional non-commtative spaces.

We first consider the problem in commutative spaces . We can describe the two particles systems by relative coordinates in the same way as in one dimensional case. We define complex quantities $`a_{j\pm }`$ as follows :

$$a_{j\pm }=\frac{1}{\sqrt{2}}(p_j\pm ix_j),j=1,2.$$
(19)
They satisfy in the following commutation relations :

$$[a_{j+},a_{k+}]=[a_j,a_k]=0[a_j,a_{k+}]=\delta _{jk}.$$
(20)
As we mentioned before, $`x`$ and $`p`$ are not observables. The generalized one-dimensional observables $`A`$, $`B`$ and $`C`$ are :

$$A_j=\frac{1}{4}(a_{j+}a_j+a_j,a_{j+}),B_{j\pm }=B_j\pm iC_j=\frac{1}{2}(a_{j\pm })^2.$$
(21)
In addition we have two-dimensional observables which are the real and imaginary parts of :

$$D_\pm =D_{re}\pm iD_{im}=a_{1\pm }a_{2\pm }E_\pm =E_{re}\pm iE_{im}=a_1a_{2\pm }.$$
(22)
There are two $`sp(1,R)`$ algebras $`A_1`$, $`B_{1\pm }`$ and $`A_2`$, $`B_{2\pm }`$  :

$$[A_j,B_{j\pm }]=\pm B_{j\pm }[B_j,B_{j+}]=2A_jj=1,2.$$
(23)
There are also two other algebras, one sp(1,R) algebra :

$$[A_1+A_2,D_\pm ]=\pm D_\pm [D_+,D_{}]=2(A_1+A_2).$$
(24)
and one $`su(2)`$ algebra :

$$[A_2A_1,E_\pm ]=\pm E_\pm [E_+,E_{}]=2(A_2A_1).$$
(25)
Now we study the same problem in two dimensional noncommutative spaces. The operators $`a_{j\pm }`$ are given by :

$$\widehat{a}_{j\pm }=\frac{1}{\sqrt{2}}(\widehat{p}_j\pm i\widehat{x}_j),$$
(26)
where $`\widehat{p}_{j\pm }`$ and $`\widehat{x}_j`$ satisfy in the Heisenberg commutation relations in noncommutative spaces i.e. equ.(2). We have the following commutation relations :

$$[\widehat{a}_{j+},\widehat{a}_{k+}]=[\widehat{a}_j,\widehat{a}_k]=\frac{i}{2}\theta _{jk},[\widehat{a}_j,\widehat{a}_{k+}]=\delta _{ij}+\frac{i}{2}\theta _{jk}.$$
(27)
The two $`sp(1,R)`$ algebras (23) are also valid in this case, but the commutation relations (24) and (25) are no longer valid and in addition $`A_1`$ and $`A_2`$ don’t commute with each other and can not have common eigenvectors. Now we use the variables $`x_i`$ and $`p_i`$ introduced in equ.(5), then we have :

$$a_{j\pm }=\frac{1}{\sqrt{2}}(p_j\pm ix_j),j=1,2.$$
(28)
where $`x_{\mathrm{}}`$ and $`p_{\mathrm{}}`$ satisfy in (7). One can show that $`a_{j\pm }`$ satisfy the same commutation relations as commutative case :

$$[a_{j+},a_{k+}]=[a_j,a_k]=0[a_j,a_{k+}]=\delta _{ij},$$
(29)
and therefore the operators $`A_j`$, $`B_{j\pm }`$, $`A_2\pm A_1`$, $`E_\pm `$ and $`D_\pm `$ satisfy in the same $`sp(1,R)`$ and $`su(2)`$ algebras as commutative space :

$$[A_j,B_{j\pm }]=\pm B_{j\pm }[B_j,B_{j+}]=2A_jj=1,2,$$
(30)
$$[A_1+A_2,D_\pm ]=\pm D_\pm [D_+,D_{}]=2(A_1+A_2),$$
(31)
$$[A_2A_1,E_\pm ]=\pm E_\pm [E_{},E_+]=2(A_2A_1).$$
(32)
We also have the following commutation relations :

$$[A_1,A_2]=0,$$
(33)
$$[E_{},D_+]=2B_{1+}[D_{},E_+]=2B_1,$$
(34)
$$[E_+,D_+]=2B_{2+}[D_{},E_{}]=2B_2,$$
(35)
$$[E_{},B_{2+}]=[E_+,B_{1+}]=D_+,$$
(36)
$$[E_{},B_1]=[E_+,B_2]=D_{},$$
(37)
$$[D_{},B_{1+}]=[D_+,B_2]=E_{},$$
(38)
$$[D_{},B_{2+}]=[D_+,B_1]=E_+.$$
(39)
We observe that by using the variables $`x_i`$ and $`p_i`$ defined by equ.(5), all of the commutation relations are the same as commutative case and the algebra governs physics of the system of identical particles and its representation remain unchanged.
The operators $`A_1`$ and $`A_2`$ together with the rasing and lowering operators $`B_{1\pm }`$, $`B_{2\pm }`$, $`D_\pm `$ and $`E_\pm `$ form a cartan basis for $`sp(2,R)`$ algebra.
The ”symplectic” algebra $`sp(2,R)`$ has 10 generators $`a_i`$, which are $`4\times 4`$ matrices and satisfy in the following relation :

$$\stackrel{~}{a}_iJ=Ja_i,$$
(40)
where $`\stackrel{~}{a}_i`$ is the transpose of the $`a_i`$, and the matrix $`J`$ is defined by :

$`J=\left[\begin{array}{cc}0& \mathrm{𝟏}\\ \mathrm{𝟏}& 0\end{array}\right]`$

where 1 is the $`2\times 2`$ unit matrix. The graphic representation of the root vectors of $`sp(2,R)`$ algebra can be drawn by introducing the angle $`\varphi `$ between the root vectors $`\alpha `$ and $`\beta `$ :

$$\mathrm{cos}(\varphi )=\frac{(\alpha ,\beta )}{\sqrt{(\alpha ,\alpha )(\beta ,\beta )}},$$
(41)
where $`(\alpha ,\beta )=\alpha ^i\beta _i`$ is the scalar product. Fig.1 shows the graphic representation of the root vectors.

Like commutative case we can construct an arbitrary representation of the $`sp(2,𝐑)`$ by a common eigenvector of $`A_1`$ and $`A_2`$ :

$$A_1|a>=a_1|a>A_2|a>=a_2|a>,$$
(42)
where $`(a_1,a_2)`$ is a weight of the representation. All of the operators $`B_1,B_2,D_{}`$ and $`E_\pm `$ lower the eigenvalues of at least one of the $`A_1`$ or $`A_2`$ by either $`\frac{1}{2}`$ or $`1`$, and we can use this fact to obtain the eigenvectors and eigenvalues of $`A_1`$ and $`A_2`$, we will have  :

$$A_1|jkmn>=(a_{10}+j+\frac{mn}{2})|jkmn>,$$
(43)
$$A_2|jkmn>=(a_{20}+k+\frac{m+n}{2})|jkmn>,$$
(44)
where $`(a_{10},a_{20})`$ are the weights of the lowest leading weight vector $`|0>`$.
Like in one dimension we can get bosons and fermions as special cases of the Heisenberg quantization in two dimensions :
a) If $`a_{10}=a_{20}=\frac{1}{4}`$, the eigenfunction is symmetric in both $`x_1`$ and $`x_2`$ or antisymmetric in both and this corresponds to the bosons.
b) If $`a_{10}=\frac{3}{4}`$, $`a_{20}=\frac{1}{4}`$, the eigenfunction has opposite symmetry in $`x_1`$ and $`x_2`$ which corresponds to fermions, thus :

$`2(a_{10}a_{20})`$ = K = $`\{\begin{array}{cc}0\hfill & \text{for Bosons}\hfill \\ 1\hfill & \text{for Fermions}\hfill \end{array}`$

This means that a system of two identical fermions or bosons in two dimensional noncommutative spaces has a definite symmetry in each of the relative coordinates $`x_1`$ and $`x_2`$ (like commutative case) and therefore the Pauli exclusion principle is also valid in noncommutative spaces.

Acknowledgment.
I would like to thank Professor P. Prešnajder for his valuable comments. I am also very grateful to Professor Masud Chaichian for his warm hospitality during my visit to the University of Helsinki. I acknowledge Dr. S. M. Harun-or-Rashid for his help on preparation of the manuscript. This work is partialy supported by the Ministry of Scince, Research and Technology of Iran.

References.

1. M. Chaichian, M. M. Sheikh-Jabbari and A. Tureanu, Phys. Rev. Lett. 86, 2716 (2001).
2. M. Chaichian, A. Demichev, P. Pre$`\stackrel{˘}{s}`$najder, M. M. Sheikh-Jabbari, and A. Tureanu, Phys. Lett. B 527 (2002) 149.
3. A. Smailagic and E. Spallucci, Phys. Rev. D 65, 107701 (2002).
4. $`\ddot{O}`$. Dayi and A. Jellal, Phys. Lett. A 287, 349 (2001).
5. P-M Ho and H-C Kao, Phys. Rev. Lett 88, 151 602 (2001).
6. R. Banerjee, hep-th/0106280.
7. J. Gamboa, F, M$`\stackrel{´}{e}`$ndez, M. Loewe and J. C. Rojac, Mod. Phys. Lett. A 16, 2075 (2001).
8. V. P. Nair and A. P. Polychronakos, Phys. Lett. B 505, 267 (2001).
9. B. Muthukamar and R. Mitra, hep-th/0204149.
10. D. Kochan and M. Demetrian, hep-th/0102050.
11. L. Mezincesu, hep-th/0007046.
12. Y. Zunger, JHEP o1o4, 039 (2001).
13. O. Espinosa and P. Gaete, hep-th/0206066.
14. J. M. Leinaas and J. Myrheim, International Journal of Modern Physics A, Vol.8, No.21 (1993) 3649.
15. W. Heisenberg, Z. Phys. 38, 411 (1926).
16. P. A. M. Dirac, Proc. R. Soc. London A 112, 661 (1926).


text:
# Long-term Neutrino Afterglows from Gamma-Ray Bursts

## 1 Introduction

In the standard fireball shock model, GRBs and their afterglows are produced due to the dissipation of kinetic energy of relativistic ejecta with Lorentz factor of $`\gamma _01001000`$, so-called “fireball”, released from an unknown central engine (see Piran 1999; Paradijs, Kouveliotou, & Wijers 2000; Mészáros 2002 for detailed reviews). In this model, the shocks, resulting from the collisions between shells with different Lorentz factors in the ejecta (internal shocks) or between the ejecta and the surrounding medium (external shocks), accelerate the electrons responsible for prompt gamma-rays in GRBs and for X-ray, optical and radio emissions of afterglows. The same shocks should also accelerate protons. It is pointed out that the physical conditions in the shocks allow protons to be accelerated to $`ϵ_p10^{20}`$ eV (Waxman 1995a; Vietri 1995; Dermer & Humi 2001. See but Stecker 2000; Mannheim 2000; Scully & Stecker 2002). Furthermore, the spectrum and flux of ultra-high energy cosmic rays (UHECRs) are consistent with those expected from Fermi acceleration of protons in GRBs (Waxman 1995b), suggesting a common origin of GRBs and UHECRs. This picture will lead obviously to high energy neutrinos produced by $`\pi ^+`$ created in photo-pion interactions between the accelerated protons in the internal or external shocks and synchrotron photons from electrons accelerated by the same shocks (Waxman & Bahcall 1997, 2000; Halzen 1998; Vietri 1998ab; Dai & Lu 2001; see Waxman 2001 for a review). This $`p\gamma `$ reaction has a peak at the energy threshold of the photo-meson $`\mathrm{\Delta }`$-resonance, i.e. $`ϵ_\gamma ϵ_p=0.2\mathrm{GeV}^2\gamma ^2`$, and the resulting neutrinos own $`5\%`$ of the proton energy in the decay of charged pions and muons. This leads to neutrinos of $`10^{14}`$ eV for prompt gamma-rays from internal shocks (Waxman & Bahcall 1997; Halzen 1998), while a UV/optical flash from the external reverse shock implies higher energy neutrinos (e.g., $`10^{18}`$ eV) (Waxman & Bahcall 2000). Vietri (1998a, 1998b) studied the neutrino production in the case of external forward shocks, and pointed out that the neutrino energy may exceed $`10^{19}`$ eV. But Waxman & Bahcall (1999) argued that the contributions of external shocks to the neutrino flux is small.

In recent years, the standard afterglow model, which is significantly successful in interpreting observations, has been developed (Mészáros & Rees 1997; Sari, Piran & Narayan 1998; Wijers & Galama 1999). In this paper, we re-investigate the case of external shock neutrinos based on the present detailed afterglow model. We derive the maximum energy that protons could reach by Fermi acceleration and the efficiency of energy loss from accelerated protons to pions in the forward shock. By considering the cooling of secondary particles, $`\pi ^+`$ and $`\mu ^+`$, we derive the spectrum and the flux of the resulting neutrinos. We find that at late times, e.g., $`1`$ day after the GRB trigger, the forward shock may be able to produce an intense neutrino signal, even though the flux is smaller at earlier times, as stated by Waxman & Bahcall (1999).

We discuss in section 2 the standard afterglow model, including the dynamics and synchrotron radiation of the blast wave, and in section 3 the proton acceleration by the same shock. In section 4 we discuss the neutrino production rate from the shock-accelerated protons and the maximum energy that neutrinos may achieve. The neutrino spectrum, flux and detectability are derived in section 5. Conclusions and discussion are given in section 6.

## 2 Dynamics and emissions of afterglows

The fireball shock model suggests that after the prompt burst from internal shocks, the ejecta with total kinetic energy $`E`$ continues to expand and drives a relativistic blast wave into the surrounding medium. Once the shocked medium gets an energy similar to $`E`$, the blast wave will enter the Blandford-McKee (1976) self-similar evolution, where the radius $`r`$ and the Lorentz factor $`\gamma `$ evolve with observer time as (Wijers & Galama 1999)

$$r=4.07\times 10^{17}(E_{53}/n_0)^{1/4}t_d^{1/4}\mathrm{cm},$$
(1)
$$\gamma =8.87(E_{53}/n_0)^{1/8}t_d^{3/8},$$
(2)
where $`n`$ is the baryon number density of the medium, and $`t_d=t/1\mathrm{day}`$. Unless stated otherwise, the useful form $`Q=10^xQ_x`$ has been used in this paper.

Let $`\xi _e`$ and $`\xi _B`$ be the fractions of the post-shock energy density $`U=4\gamma ^2nm_pc^2`$ (in the comoving frame with the shocked material) that are carried by electrons and magnetic fields respectively. The characteristic electron Lorentz factor (in the comoving frame of the shocked material) is $`\gamma _m\xi _e\gamma m_p/m_e`$, where the post-shock energy per proton is $`\gamma m_pc^2`$. The strength of the magnetic field in the comoving frame is

$$B=0.34\xi _{B,2}^{1/2}E_{53}^{1/8}n_0^{3/8}t_d^{3/8}\mathrm{G}.$$
(3)
From the above equations, we can calculate the evolution of synchrotron spectrum from the blast wave. The characteristic energy of synchrotron photons from electrons with Lorentz factor $`\gamma _m`$ is

$$ϵ_m^{\mathrm{ob}}\mathrm{}\gamma \gamma _m^2\frac{eB}{m_ec}=0.09\xi _{e,1}^2\xi _{B,2}^{1/2}E_{53}^{1/2}t_d^{3/2}\mathrm{eV},$$
(4)
where and hereafter we denote the particle energy in the observer frame with superscript “ob”, and particle energy measured at the comoving frame without superscript. The characteristic energy of synchrotron photons from electrons with the cooling time, $`6\pi m_ec/\sigma _T\gamma _cB^2`$, which equals the blast wave expansion time ($`r/\gamma c`$), is given by

$$ϵ_c^{\mathrm{ob}}\mathrm{}\gamma \gamma _c^2\frac{eB}{m_ec}=3\xi _{B,2}^{3/2}n_0^1E_{53}^{1/2}t_d^{1/2}\mathrm{eV},$$
(5)
and the specific luminosity $`L_ϵ=dL/dϵ^{\mathrm{ob}}`$ at the peak is

$$L_{max}(2\pi \mathrm{})^1\gamma \frac{e^3B}{m_ec^2}N_e5\times 10^{58}(\xi _{B,2}n_0)^{1/2}E_{53}\mathrm{s}^1,$$
(6)
where $`N_e=\frac{4}{3}\pi r^3n`$ is the total number of the swept-up electrons. Eq. (6) implies that the peak specific luminosity is a time-independent constant. The electrons are accelerated by the shock to a power-law energy distribution, $`dN_e/d\gamma _e\gamma _e^p`$ for $`\gamma _e>\gamma _m`$, where $`p2.2`$ is fitted to the observations of afterglows. We take $`p=2`$ for simplicity here. After $`t1`$ hour, $`ϵ_c^{\mathrm{ob}}>ϵ_m^{\mathrm{ob}}`$. In this case, the specific luminosity $`L_ϵ`$ peaks at $`ϵ_m^{\mathrm{ob}}`$, $`L_{max}=L_m`$, and the synchrotron spectrum is a broken power law (Sari, Piran & Narayan 1998) as: $`L_ϵ=L_m(ϵ/ϵ_m)^{1/2}`$ for $`ϵ_m<ϵ<ϵ_c`$ and $`L_ϵϵ^1`$ for $`ϵ>ϵ_c`$. Note that in the above equations, the typical parameter values have been adopted as $`\xi _e=10^1\xi _{e,1}`$, $`\xi _B=10^2\xi _{B,2}`$, $`E=10^{53}E_{53}`$ erg and $`n=1n_0`$ cm-3, which are inferred from some afterglow observations (e.g. Wijers & Galama 1999; Granot, Piran & Sari 1999).

## 3 High-energy proton production

Following the analysis of Waxman (1995a) in the case of internal shocks, we investigate the physical condition in external forward shocks for high energy protons production. The typical Fermi acceleration time is $`t_a=fR_L/c`$, where the Larmor radius is given by $`R_L=ϵ_p/eB`$ and $`f`$ is of order unity (Hillas 1984). The protons are available for acceleration only within a comoving expansion time $`t_dr/\gamma c`$. This leads to the requirement that $`t_a<t_d`$. From Eqs. (1)-(3), we have

$$ϵ_p^{\mathrm{ob}}<4\times 10^{19}\xi _{B,2}^{1/2}E_{53}^{3/8}n_0^{1/8}t_d^{1/8}\mathrm{eV}.$$
(7)
If a magnetic field is nearly in equipartition (i.e., $`\xi _B1`$), the protons can be accelerated to well above $`10^{20}`$ eV. However, fits of afterglow observations give a typical value $`\xi _B=0.01`$, showing that the magnetic fields are too weak to accelerate protons up to $`10^{20}`$ eV in the blast waves. This upper limit is somewhat robust since its dependence on the medium density and the time ($`n^{1/8}t^{1/8}`$) is very weak.

The proton acceleration is affected by energy loss due to synchrotron radiation. The requirement that the proton synchrotron cooling time, $`t_{p,syn}=(6\pi m_p^4c^3/\sigma _Tm_e^2)ϵ_p^1B^2`$, should be larger than the acceleration time $`t_a`$ is

$$ϵ_p^{\mathrm{ob}}<3.2\times 10^{21}\xi _{B,2}^{1/4}E_{53}^{1/16}n_0^{5/16}t_d^{3/16}\mathrm{eV}.$$
(8)
The proton acceleration is also limited by energy loss due to interaction with afterglow photons. In the next section, we show that, provided a typical medium density $`n_0=1`$, the fraction of energy loss of a proton in the expansion time $`t_d`$ is usually not larger than one, which implies that the limit from energy loss due to photo-pion interaction is a less stringent constraint than that from $`t_a<t_d`$. Thus, the maximum energy of protons comes from Eq. (7), $`ϵ_{p,max}^{\mathrm{ob}}4\times 10^{19}`$ eV.

Gallant & Achterberg (1999) studied the shock-acceleration, also called first order Fermi-acceleration, in GRBs and argued that it may be impossible for protons to be shock-accelerated to above $`10^{18.5}`$ eV in an ultra-relativistic blast wave which is driven by the fireball ejecta into a typical interstellar medium. However we still expect protons to be accelerated to ultra-high energy in the blast waves by second order Fermi-acceleration, where the protons gain energy from the scattering at plasma waves (e.g. Dermer & Humi 2001).

## 4 Neutrino production

Denoting the photon number density in the comoving frame of the shocked medium by $`n_\gamma (ϵ)dϵ`$ and following Waxman & Bahcall (1997), we can write the fractional energy-loss rate of a proton with energy $`ϵ_p`$ due to pion production,

$`t_\pi ^1(ϵ_p)`$ $`{\displaystyle \frac{1}{ϵ_p}}{\displaystyle \frac{dϵ_p}{dt}}`$
$`=`$ $`{\displaystyle \frac{1}{2\gamma _p^2}}c{\displaystyle _{ϵ_0}^{\mathrm{}}}𝑑ϵ\sigma _\pi (ϵ)\xi (ϵ)ϵ{\displaystyle _{ϵ/2\gamma _p}^{\mathrm{}}}𝑑xx^2n(x),`$ (9)
where $`\gamma _p=ϵ_p/m_pc^2`$, $`\sigma _\pi (ϵ)`$ is the cross section for pion production for a photon with energy $`ϵ`$ in the proton rest frame, $`\xi (ϵ)`$ is the average fraction of energy lost to the pion, $`ϵ_0=0.15\mathrm{GeV}`$ is the threshold energy, and the proton number density is related to the observed specific luminosity by $`n(x)=L_ϵ(\gamma x)/(4\pi r^2c\gamma x)`$. Because of the $`\mathrm{\Delta }`$-resonance, the photo-meson production for $`ϵ_0/2ϵ_c\gamma _p<ϵ_0/2ϵ_m`$ is dominated by interaction with photons in the energy range $`ϵ_m<ϵϵ_c`$, where $`L_ϵϵ^{1/2}`$. For this photon spectrum, Eq. (9) leads to

$$t_\pi ^1(ϵ_p)\frac{2^{3/2}}{2.5}\frac{L_m}{4\pi r^2\gamma }\left(\frac{ϵ_{\mathrm{peak}}}{\gamma _pϵ_m}\right)^{1/2}\frac{\sigma _{\mathrm{peak}}\xi _{\mathrm{peak}}\mathrm{\Delta }ϵ}{ϵ_{\mathrm{peak}}}.$$
(10)
Here, $`\sigma _{\mathrm{peak}}5\times 10^{28}\mathrm{cm}^2`$ and $`\xi _{\mathrm{peak}}0.2`$ at the resonance $`ϵ=ϵ_{\mathrm{peak}}=0.3\mathrm{GeV}`$, and $`\mathrm{\Delta }ϵ0.2\mathrm{GeV}`$ is the peak width. The fraction of energy loss of a proton with observed energy $`ϵ_p^{\mathrm{ob}}`$ by pion production, $`f_\pi (ϵ_p^{\mathrm{ob}})`$, is estimated by $`t_\pi ^1(ϵ_p)`$ times the comoving expansion time of the shocked material ($`r/\gamma c`$). Thus, according to Eqs. (1), (2), (4) and (6), we obtain

$$f_\pi (ϵ_p^{\mathrm{ob}})0.07\xi _{e,1}\xi _{B,2}^{3/4}E_{53}^{5/8}n_0^{9/8}t_d^{1/8}(ϵ_{p,19}^{\mathrm{ob}})^{1/2}.$$
(11)
Eq. (11) is valid for protons in the energy range

$`6\times 10^{17}\xi _{B,2}^{3/2}E_{53}^{3/4}n_0^{3/4}t_d^{1/4}\mathrm{eV}<ϵ_p^{\mathrm{ob}}`$
$`<5\times 10^{19}\xi _{e,1}^2\xi _{B,2}^{1/2}E_{53}^{1/4}n_0^{1/4}t_d^{3/4}\mathrm{eV}.`$ (12)
At lower energy, protons interact mainly with photons in the energy range $`ϵ>ϵ_c`$, where $`L_ϵϵ^1`$, and $`f_\pi ϵ_p^{\mathrm{ob}}`$. Eqs. (11) and (12) imply that the fraction of energy loss, $`f_{\pi b}`$, at the “break” energy $`ϵ_{pb}^{\mathrm{ob}}=6\times 10^{17}\xi _{B,2}^{3/2}E_{53}^{3/4}n_0^{3/4}t_d^{1/4}\mathrm{eV}`$, is a time-independent constant,

$$f_{\pi b}0.02\xi _{e,1}\xi _{B,2}^{3/2}E_{53}n_0^{3/2}.$$
(13)
The maximum energy of the resultant neutrinos is analyzed as follows. This energy is first determined by the maximum energy that the accelerated protons achieve. From Eq. (7), the maximum energy of neutrinos is $`2\times 10^{18}`$ eV. This energy is also limited by the energy loss of secondary particles, i.e., pions and muons, because both pions and muons may suffer synchrotron and inverse-Compton (IC) losses before decay (Rachen & Mészáros 1998). For synchrotron losses of pions and muons to be negligible, their Lorentz factor (in the comoving frame of shocked medium) should not exceed the critical Lorentz factor $`\gamma _{}`$ given by $`\gamma _{}^2\tau _{}=\gamma _pt_{p,syn}(m_{}/m_p)^3`$ (Rachen & Mészáros 1998), where $``$ refers to either pions or muons, $`\tau _\pi =2.6\times 10^8`$ s and $`\tau _\mu =2.2\times 10^6`$ s are respective lifetimes in their rest frames, and $`m_{}/m_p0.1`$ for both pions and muons. Due to the longer lifetime compared to pions, muons suffer greater energy loss, resulting in lower upper limit of neutrino energy. Using Eq. (3), we have (in the observer frame) the critical energy for muons

$$ϵ_{\mu ,syn}^{\mathrm{ob}}=10^{20}\xi _{B,2}^{1/2}n_0^{1/2}\mathrm{eV},$$
(14)
which is independent of the total fireball energy $`E`$ and observer time $`t`$. For IC loss to be negligible, the critical Lorentz factor is given by $`\gamma _{}^2\tau _{}=\gamma _pt_{p,syn}(m_{}/m_p)^3(U_B/U_{ph})`$, where $`U_B=\xi _BU`$ and $`U_{ph}=\xi _e\eta U`$ are the energy densities of the magnetic field and afterglow photons, respectively, and $`\eta `$ is the fraction of the electron energy that is radiated away. Note that $`\eta =(ϵ_c/ϵ_m)^{(2p)/2}`$ for slow cooling electrons with $`ϵ_c>ϵ_m`$. We take $`p2.2`$. $`\eta `$ is not sensitive to other parameters as $`\eta 1`$. (In fact, $`\eta =1`$ shows the strongest IC loss and provides the strongest constraint). Thus, the critical Lorentz factor is $`\gamma _{}^{IC}\gamma _{}^{syn}(\xi _B/\xi _e)^{1/2}`$, or

$$ϵ_{\mu ,IC}^{\mathrm{ob}}=4\times 10^{19}\xi _{e,1}^{1/2}n_0^{1/2}\mathrm{eV},$$
(15)
which is independent of $`\xi _B`$. The IC loss is somewhat stronger than the synchrotron loss provided $`\xi _B/\xi _e<1`$. But, the energy limit from both losses exceeds the maximum proton energy from Eq. (7), so the neutrino production does not suffer the energy losses of secondary particles. The maximum neutrino energy is determined by the maximum proton energy, $`ϵ_{\nu ,max}^{\mathrm{ob}}2\times 10^{18}`$ eV.

## 5 Neutrino spectrum, flux and detectability

The photo-meson interactions include (1) production of $`\pi `$ mesons: $`p\gamma p+\pi ^0`$ and $`p\gamma n+\pi ^+`$, and (2) decay of $`\pi `$ mesons: $`\pi ^02\gamma `$ and $`\pi ^+\mu ^++\nu _\mu e^++\nu _e+\overline{\nu }_\mu +\nu _\mu `$. These processes produce neutrinos with energy $`5\%ϵ_p`$ (Waxman & Bahcall 1997). Eq. (12) implies that for a fixed time, the spectrum of neutrinos below $`ϵ_{\nu b}^{\mathrm{ob}}3\times 10^{16}\xi _{B,2}^{3/2}E_{53}^{3/4}n_0^{3/4}t_d^{1/4}\mathrm{eV}`$ is harder by one power of the energy than the proton spectrum, and by half a power of the energy at higher energy. Therefore, if the differential spectrum of accelerated protons is a power law form $`n(ϵ_p)ϵ_p^2`$ (Blandford & Eichler 1987), the differential neutrino spectrum is $`n(ϵ_\nu )ϵ_\nu ^1`$ below the break and $`n(ϵ_\nu )ϵ_\nu ^{3/2}`$ above the break. Unlike the prompt bursts of neutrinos from either internal or reverse shocks, which has a short duration similar to the GRBs, the neutrino emission from forward shocks last for a longer time due to the long-term forward shocks. Without the detailed knowledge of the evolution of the spectrum of accelerated protons in the forward shocks, it is far from calculating clearly the neutrino fluxes. However, as $`ϵ_{\nu b}^{\mathrm{ob}}(t^{1/4})`$ and $`ϵ_{p,max}^{\mathrm{ob}}(t^{1/8})`$ do not change much from one day to even weeks which is concerned about here, one can expect a time-integrated neutrino spectrum similar to the spectrum at $`t=1`$ day. For simplicity, the particle energies, hereafter in this section, refer to the value when $`t=1`$ day, e.g. $`ϵ_{\nu b}^{\mathrm{ob}}ϵ_{\nu b}^{\mathrm{ob}}(t_d=1)`$, $`ϵ_{pb}^{\mathrm{ob}}ϵ_{pb}^{\mathrm{ob}}(t_d=1)`$, etc.

We derive the neutrino fluxes as follows. Defining by $`\overline{f}_\pi n(ϵ_p)ϵ_pf_\pi 𝑑ϵ_p/n(ϵ_p)ϵ_p𝑑ϵ_p`$ the fraction of energy loss of all accelerated protons to neutrinos, and considering the neutrino spectrum given in the above paragraph, we obtain, after some calculations,

$$\overline{f}_\pi 2f_{\pi b}\frac{(ϵ_{p,max}^{\mathrm{ob}}/ϵ_{pb}^{\mathrm{ob}})^{1/2}}{\mathrm{ln}(ϵ_{p,max}^{\mathrm{ob}}/ϵ_{pl})},\mathrm{for}ϵ_{pl}ϵ_{pb}^{\mathrm{ob}}.$$
(16)
Here $`f_{\pi b}`$ is a constant for a certain GRB, therefore this value is mainly bound to $`ϵ_{pl}`$ which is determined by the lower threshold of a detector. In this equation, $`ϵ_{pb}^{\mathrm{ob}}`$ should be replaced by $`ϵ_{pl}`$ for $`ϵ_{pl}ϵ_{pb}^{\mathrm{ob}}`$. Taking $`ϵ_{pl}=ϵ_{pb}^{\mathrm{ob}}`$, we have $`\overline{f}_\pi 0.07`$.

The present day muon neutrino energy flux due to GRBs is approximately given by $`J0.25(c/4\pi )\overline{f}_\pi \dot{E}t_H`$, where $`t_H10`$ Gyr is the Hubble time, and $`\dot{E}`$ is the production rate of UHECRs in GRBs per unit volume. The factor $`0.25`$ here accounts for that about one half of the proton energy is lost to neutral pions which do not produce neutrinos, and about one half of the energy in charged pions is converted to $`\nu _\mu `$’s and $`\overline{\nu }_\mu `$’s. We take $`\dot{E}10^{44}\dot{E}_{44}\mathrm{erg}\mathrm{Mpc}^3\mathrm{yr}^1`$ (Waxman 1995b). Integrating the specific neutrino flux $`\mathrm{\Phi }_\nu `$ over the neutrino spectrum ($`\mathrm{\Phi }_\nu [ϵ_\nu ^{\mathrm{ob}}]^0`$ below $`ϵ_{\nu b}^{\mathrm{ob}}`$, and $`\mathrm{\Phi }_\nu [ϵ_\nu ^{\mathrm{ob}}]^{1/2}`$ above $`ϵ_{\nu b}^{\mathrm{ob}}`$), we have $`J=\mathrm{\Phi }_\nu 𝑑ϵ_\nu 2(ϵ_{\nu b}^{\mathrm{ob}}ϵ_{\nu ,max}^{\mathrm{ob}})^{1/2}\mathrm{\Phi }_{\nu b}`$. Combining these two equations about $`J`$ one obtains that the neutrino flux at $`ϵ_{\nu b}^{\mathrm{ob}}`$ is

$$\mathrm{\Phi }_{\nu b}2\times 10^{18}\frac{\overline{f}_\pi }{0.07}\dot{E}_{44}\left(\frac{ϵ_{\nu b}^{\mathrm{ob}}}{3\times 10^{16}\mathrm{eV}}\frac{ϵ_{\nu ,max}^{\mathrm{ob}}}{2\times 10^{18}\mathrm{eV}}\right)^{1/2}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1.$$
(17)
The resulting high-energy neutrinos may be observed by detecting the Cherenkov light emitted by upward moving muons produced by neutrino interactions below a detector on the surface of the Earth (Gaisser, Halzen & Stanev 1995; Gandhi et al. 1998). Planned 1 km3 detectors of high energy neutrinos include ICECUBE, ANTARES, NESTOR (Halzen 1999) and NuBE (Roy, Crawford & Trattner 1999). The probability that a neutrino could produce a high-energy muon in the detector is approximated by $`P_{\nu \mu }10^3(ϵ_\nu /10^{16}\mathrm{eV})^{1/2}`$. We obtain the observed neutrino event rate in a detector,

$`N_{\mathrm{events}}(>ϵ_{\nu b}^{\mathrm{ob}})`$ $`=2\pi {\displaystyle \frac{\mathrm{\Phi }_\nu }{ϵ_\nu ^{\mathrm{ob}}}P_{\nu \mu }𝑑ϵ_\nu ^{\mathrm{ob}}}2\pi \mathrm{\Phi }_{\nu b}P_{\nu \mu }(ϵ_{\nu b}^{\mathrm{ob}})\mathrm{ln}\left({\displaystyle \frac{ϵ_{\nu ,max}^{\mathrm{ob}}}{ϵ_{\nu b}^{\mathrm{ob}}}}\right)`$ (18)
$`0.04{\displaystyle \frac{\overline{f}_\pi }{0.07}}\dot{E}_{44}\left({\displaystyle \frac{ϵ_{\nu ,max}^{\mathrm{ob}}}{2\times 10^{18}\mathrm{eV}}}\right)^{1/2}\mathrm{km}^2\mathrm{yr}^1.`$
## 6 Conclusions and discussion

After a prompt $`\gamma `$-rays burst, the fireball ejecta continues to drive a relativistic blast wave into the surrounding medium. The shock wave accelerates the electrons to give rise to observed afterglow by synchrotron radiation. The protons in the shock are also expected to be Fermi-accelerated at the same time. However, the maximum proton energy may not well exceed $`10^{20}`$ eV based on the typical parameters of afterglows. The interaction of the protons with synchrotron photons leads to a significant fraction of proton energy loss to neutrino production at energy $`10^{16}`$$`10^{18}`$ eV. Though the synchrotron photon energy density in the blast wave is much lower than in the reverse shock, the efficiency of proton energy loss in the blast wave is still comparable to that in the reverse shock, because of the long expansion time of the blast wave available for energy loss of protons. Since both shocks carry comparable energies (Sari & Piran 1999), we assume the same production rate of high-energy cosmic rays, $`\dot{E}_{44}1`$ for blast waves as Waxman & Bahcall (2000) for reverse shocks. Thus, the expected detection rate of muon neutrinos $`0.04`$ km-2 yr-1 is also comparable to the reverse shock. These mean that besides the internal shock and the reverse shock, the blast wave is also one region suitable for intense neutrino production.

Waxman & Bahcall (1999) have estimated a much lower efficiency of the blast wave, i.e. $`f_\pi 10^4`$, using Eq. (4) of Waxman & Bahcall (1997). However, the parameter values they chose are only valid for the early phase when the external shock just forms. When a set of parameters, which are valid for afterglows one day after the burst, are chosen as follows: $`L_\gamma =10^{45}\mathrm{erg}/\mathrm{s},\mathrm{\Delta }t=10^5\mathrm{s},ϵ_{\gamma b}^{\mathrm{ob}}=0.1\mathrm{eV}\mathrm{and}\mathrm{\Gamma }=10`$, their equation gives an efficiency $`f_\pi 10^2`$ in general agreement with ours by orders of magnitude.

It may be that UHECRs can not be produced in an ultra-relativistic forward shock (Gallant & Achterberg 1999) during the transition of the shock to the self-similar expansion. However the hard X-rays from the forward shock will interact with the protons accelerated by the reverse shock, yielding a new component of neutrino emission besides those produced by interaction between protons and optical/UV photons from the reverse shock, which is proposed by Waxman & Bahcall (2000). At the moment of transition, especially for the photons at energy $`ϵ^{\mathrm{ob}}ϵ_m^{\mathrm{ob}}80`$ keV (cf. Eq. ), the specific luminosity from the forward shock dominates that from the reverse shock, which implies a higher neutrino flux at energy $`ϵ_\nu ^{\mathrm{ob}}10^{15}(\gamma _0/250)(ϵ_m^{\mathrm{ob}}/100\mathrm{k}\mathrm{e}\mathrm{V})^1`$ eV.

Different models for progenitors of GRBs include different environments: compact object mergers occur in the ISM and the explosions of massive stars in the preexisting stellar wind or denser medium. Dai & Lu (2001) have found that the neutrino production from reverse shocks in the wind case has a quite different spectrum and 10 times more detection rate compared with the ISM case, which provides a way to distinguish between GRB progenitor models. In the blast-wave case discussed in this paper, the proton energy loss rate shows a strong dependence on density $`n`$(Eq. ), we expect much higher fluxes for a denser medium. For $`n_4=1`$ typically for giant molecular clouds, following the analysis in this paper, the energy loss efficiency of a proton at one day after the the GRB trigger is given by $`f_\pi 8.6ϵ_{p,18}^{\mathrm{ob}}t_d^{1/4}`$ for $`ϵ_p^{\mathrm{ob}}<5\times 10^{18}t_d^{3/4}`$ eV, implying that protons of energy $`ϵ_p^{\mathrm{ob}}>10^{17}`$ eV lose all their energy to pion production and the neutrino spectrum for $`ϵ_\nu ^{\mathrm{ob}}>5\times 10^{15}`$ eV will trace the proton one. The maximum energy that protons can be accelerated to, which is now confined by the pion production, is $`4\times 10^{18}t_d^{3/4}`$ eV. This implies a maximum neutrino energy of $`10^{17}`$ eV.

Another consequence of the photo-meson interactions is a sub-TeV afterglow. The decay of neutral pions produced in the blast waves would lead to $`10^{18}`$ eV photons. Such high energy photons would suffer pair production in the afterglows because the optical depth for photons above $`1`$ TeV is important (Zhang & Mészáros 2001). So these ultra-high-energy photons would be degraded and leak out as a long-term sub-TeV afterglow, carrying away $`10\%`$ energy of protons, similar to the energy losses to neutrinos. For a GRB at typical redshift $`z=1`$, the intergalactic diffuse radiation field absorbs the photon above $`100`$ GeV (Mannheim, Hartmann & Funk 1996), allowing only a multi-GeV afterglow observed on the earth. An EGRET observation of GeV emission extended to 1.5 hours after the trigger of GRB 940217 (Hurley et al. 1994) has indicated that a high energy spectral component can extend to such high energy band and persist for a long period of time. However, the electron-IC emission in the blast wave can also contribute to this GeV observation. For the parameters, $`n_0\xi _{e,1}\xi _{B,2}1`$, adopted in this paper, the electron-IC emission is important and plausible to explain the observation (Zhang & Mészáros 2001). From the afterglow model presented in section 2, we can see that the peak energy in the IC emission spectrum, $`\gamma _c^2ϵ_c^{\mathrm{ob}}`$ (Sari & Esin 2001), is still marginally up to 100 GeV even one day after the GRB trigger, but we do not expect GeV IC emission a few days later, while the emission produced by neutral pion decay might still be possible. In principle, the GeV emission could also be due to the proton synchrotron emission, but it is not important for the parameters adopted here (Zhang & Mészáros 2001).

###### Acknowledgements.

The authors thank the referee for very good comments to improve this paper. This work was supported by the National Natural Science Foundation of China (grants 19825109 and 19973003) and the National 973 Project (NKBRSF G19990754).