text:
# GRB Afterglows from Anisotropic Jets

## 1 Introduction

Although the study of gamma-ray bursts (GRBs) has been revolutionized due to observations of their multiwavelength afterglows and particularly due to determinations of their redshifts in the past few years, the question whether the GRB emission is spherical or jet-like has been unsolved. This question has important implications on almost all aspects of the GRB phenomenon, e.g., the total energy that is released in an explosion, the event rate, the physical ejection mechanism and the afterglow decay rate (for a review see Sari 2000).

GRB 990123, the most energetic GRB well-studied, has an isotropic gamma-ray release of $`3.4\times 10^{54}`$ ergs, which corresponds to the rest-mass energy of $`1.9M_{}`$ (Kulkarni et al. 1999). From such an energetics, one deduced for the first time that the GRB emission may be highly collimated in this case with a half opening angle of $`\theta _m0.2`$ so that the intrinsic explosive energy could be reduced to $`10^{51}`$$`10^{52}`$ ergs, which is still consistent with the stellar death models, e.g., collapsars (MacFadyen, Woosley & Heger 2000). Another well-known feature of this burst is that its R-band afterglow light curve began to steepen at about day 2 after the burst (Kulkarni et al. 1999; Castro-Tirado et al. 1999; Fruchter et al. 1999). This steepening was first argued to be due to the possibility that a jet has evolved from the spherical-like phase to the sideways-expansion phase (Rhoads 1999; Sari, Piran & Halpern 1999) or that we have observed the edge of the jet (Panaitescu & Mészáros 1999; Mészáros & Rees 1999). Alternative explanation for this steepening was subsequently proposed to be that a shock expanding in a dense medium has undergone the transition from the relativistic phase to the non-relativistic phase (Dai & Lu 1999, 2000). Since then, the dynamics of expanding jets has been numerically studied by many authors (Panaitescu & Mészáros 1999; Moderski, Sikora & Bulik 2000; Huang et al. 2000a, b; Kumar & Panaitescu 2000; Gou et al. 2000; Wei & Lu 2000). Observationally, one marked break has also been observed to appear in the optical afterglow light curves of GRB 990510 (Harrison et al. 1999; Stanek et al. 1999), GRB 991216 (Halpern et al. 2000; Sagar et al. 2000a) and GRB 000301C (Rhoads & Fruchter 2000; Masetti et al. 2000; Jensen et al. 2000; Sagar et al. 2000b; Berger et al. 2000). Holland et al. (2000) recently collected, re-analyzed and explained all of the published photometry for GRB 990123 and GRB 990510. All these studies have assumed that both the energy per unit solid angle and the bulk Lorentz factor at any angle within jets are independent of the angle. Hereafter we refer to such jets as isotropic ones.

In fact, both the energy per unit solid angle ($`dE/d\mathrm{\Omega }`$) and the bulk Lorentz factor within a realistic jet may strongly depend on the angle $`\theta `$. For example, $`dE/d\mathrm{\Omega }\theta ^k`$ where $`k3`$ within the energetic jets expected in the collapsar model for GRBs of MacFadyen et al. (2000). We call such jets anisotropic ones. Salmonson (2000) argued that if GRBs arise from anisotropic jets, then the range of the viewing angle (viz., the angle between the line of sight and the jet axis) produces a range in the observed properties of GRBs, i.e. the lag-luminosity relationship. In this paper, we calculate the emission from anisotropic jets expanding both in the interstellar medium (ISM) and in the wind medium. In particular, we consider two kinds of wind: one ($`nr^{3/2}`$) possibly from a typical red supergiant star and another ($`nr^2`$) possibly from a Wolf-Rayet star. To our knowledge, this is the first numerical work about evolution of anisotropic jets and their emission. A preliminary analysis of afterglows from anisotropic jets was made by Mészáros, Rees & Wijers (1998), but our results are different from their analytical ones because we take into account the light travel effects related to different sub-jets within the jet. We also carry out a detailed modelling of the optical afterglow flux data of GRB 991208 given by Sagar et al. (2000a) and Castro-Tirado et al. (2000). Finally, we argue that this burst might arise from a highly anisotropic jet expanding in the wind from a red supergiant to explain the observed radio-to-optical-band afterglow data.

## 2 The Model

An energetic anisotropic jet is supposed to be ejected in an explosion. We assume that such a jet is adiabatic, i.e., the radiative energy is a negligible fraction of the total energy of the jet. This assumption is valid if the energy density of the electrons accelerated by a shock, produced by the interaction of the jet with its surrounding medium, is a small fraction $`ϵ_e`$ of the total energy density of the shocked medium or if most of the electrons are adiabatic, i.e., their radiative cooling timescale is larger than that of the jet expansion (Dai, Huang & Lu 1999; Böttcher & Dermer 2000).

### 2.1 Dynamics

The central question in the next study is how the Lorentz factor $`\gamma `$ at angle $`\theta `$ within the jet evolves with the observer time $`t`$ because all other quantities that appear in the observed flux are functions of $`\gamma `$ and of the jet radius and medium density. To obtain $`\gamma `$, we assume that the kinetic energy per unit solid angle within the jet has the following distribution:

$$\epsilon (\theta )\frac{dE}{d\mathrm{\Omega }}=\{\begin{array}{cc}\epsilon _0\hfill & \mathrm{if}\theta \theta _0\hfill \\ \epsilon _0\left(\frac{\theta }{\theta _0}\right)^k\hfill & \mathrm{if}\theta _0<\theta <\theta _m,\hfill \end{array}$$
(1)
where $`\theta _m`$ is the half opening angle of the jet. The ejected mass per unit solid angle, $`M_{\mathrm{ej}}(\theta )=M_0`$, is assumed to be independent of the angle (at least over some range of angles). A motivation for this assumption is that a different distribution of the ejected mass affects only evolution of the jet at very early times, but at later times the jet expands based on the Blandford-McKee’s (1976) self-similar solution during the relativistic phase and the Sedov-Taylor self-similar solution during the non-relativistic phase. Thus, the initial Lorentz factor at angle $`\theta `$ can be written as

$$\gamma _0(\theta )=\{\begin{array}{cc}\frac{\epsilon _0}{M_0c^2}+1\hfill & \mathrm{if}\theta \theta _0\hfill \\ \left(\frac{\epsilon _0}{M_0c^2}\right)\left(\frac{\theta }{\theta _0}\right)^k+1\hfill & \mathrm{if}\theta _0<\theta <\theta _m.\hfill \end{array}$$
(2)
The evolution of the Lorentz factor $`\gamma (\theta ,r)`$ can be calculated from the following equation

$$M(r)\gamma ^2+M_0\gamma =M(r)+M_0+\epsilon (\theta )/c^2,$$
(3)
where $`M(r)`$ is the mass of the swept-up medium per unit solid angle. Equation (3) expresses the conservation of energy and it applies to only the case without heating of the ejected material by the reverse shock. In this paper we neglect sideways expansion of the jet. Evidence for this consideration is that an observation of GRB 990123 on 7 February 2000 by HST (Fruchter et al. 2000), combined with the previous observations, indicates that the steepening of the R-band afterglow light curve of this burst during a period of $`378`$ days is still roughly consistent with the edge effect of an ultra-relativistic jet of a fixed opening angle expanding in a homogeneous medium proposed by Mészáros & Rees (1999). Further evidence is that, as argued by Jaunsen et al. (2000), the rapid initial decline, the sharp break in the optical light curve, and the spectral properties of the GRB 980519 afterglow are best interpreted as being due to an ultra-relativistic jet of a fixed opening angle expanding in a wind medium with $`nr^2`$.

The external medium density is assumed to be a power-law function of radius: $`n(r)=Ar^s`$, where $`A=n_{}\times 1\mathrm{cm}^3`$ for an $`s=0`$ (homogeneous) medium (e.g., ISM), $`A=3\times 10^{35}A_{}\mathrm{cm}^1`$ for an $`s=2`$ wind possibly from a Wolf-Rayet star as the GRB progenitor (Chevalier & Li 1999, 2000), and $`A=10^{30}\overline{A}_{}\mathrm{cm}^{3/2}`$ for an $`s=3/2`$ wind possibly from a typical red supergiant star (Fransson, Lundquist & Chevalier 1996). For convenience, $`x`$ is defined as the radius $`r`$ scaled to

$$r_0=\left[\frac{(3s)\epsilon }{m_pc^2A\gamma _0^2}\right]^{1/(3s)}=\{\begin{array}{ccc}2.7\times 10^{17}\mathrm{cm}\epsilon _{53}^{1/3}n_{}^{1/3}\gamma _{0,2}^{2/3}\hfill & \mathrm{if}s=0\hfill & \\ 2.2\times 10^{16}\mathrm{cm}\epsilon _{53}A_{}^1\gamma _{0,2}^2\hfill & \mathrm{if}s=2\hfill & \\ 4.6\times 10^{14}\mathrm{cm}\epsilon _{53}^{2/3}\overline{A}_{}^{2/3}\gamma _{0,2}^{4/3}\hfill & \mathrm{if}s=3/2,\hfill & \end{array}$$
(4)
where $`\epsilon _{53}=\epsilon (\theta )/10^{53}\mathrm{ergs}`$ and $`\gamma _{0,2}=\gamma _0(\theta )/10^2`$. Furthermore, we have

$$M(r)=\frac{m_pn(r)r^3}{3s},$$
(5)
where $`m_p`$ is the proton mass.

The solution of equation (3) is

$$\gamma (\theta ,r)=\frac{1}{2M}\left[\sqrt{M_0^2+4M(M+M_0+\epsilon /c^2)}M_0\right].$$
(6)
The dimensionless radius $`x`$ and the time $`t^{}`$ in the frame comoving with the jet evolve with the observer time $`t`$ based on

$$\frac{dx}{dt}=\frac{c\beta }{r_0(1\beta \mu )},$$
(7)
$$\frac{dt^{}}{dt}=\frac{1}{\gamma (1\beta \mu )},$$
(8)
where $`\beta =\sqrt{11/\gamma ^2}`$, $`\mu =\mathrm{cos}\theta `$, and the line of sight has been taken to be the jet axis. The solution of equations (7) and (8) combined with (6) gives the dynamics of a sub-jet at angle $`\theta `$ within the jet.

### 2.2 Synchrotron Radiation

We consider synchrotron radiation of the electrons accelerated by the shock. To calculate the spectrum and light curve, one needs to determine three crucial frequencies: the synchrotron self-absorption frequency ($`\nu _a^{}`$), the typical synchrotron frequency ($`\nu _m^{}`$) and the cooling frequency ($`\nu _c^{}`$). We assume a power law distribution of the electrons accelerated by the shock: $`dn_e/d\gamma _e\gamma _e^p`$ for $`\gamma _e\gamma _m`$, where $`p`$ is the index of the electron energy distribution and $`\gamma _m=[(p2)/(p1)](m_p/m_e)ϵ_e\gamma `$ is the minimum Lorentz factor. If $`ϵ_B`$ is assumed to be the ratio of the magnetic to thermal energy densities of the shocked medium, the magnetic field strength in the shocked medium can be approximated by $`B^{}=[32\pi ϵ_B\gamma (\gamma 1)nm_pc^2]^{1/2}`$.

Under these assumptions, the typical synchrotron frequency can be written: $`\nu _m^{}=\gamma _m^2eB^{}/(2\pi m_ec)`$. The cooling frequency, $`\nu _c^{}=\gamma _c^2eB^{}/(2\pi m_ec)`$, is clearly determined by the Lorentz factor $`\gamma _c`$, at which an electron is cooling on the dynamical expansion timescale ($`t^{}`$, measured in the comoving frame). According to Sari, Piran & Narayan (1998), $`\gamma _c=6\pi m_ec/(\sigma _TB^2t^{})`$ where $`\sigma _T`$ is the Thomson scattering cross-section. Scaled with the typical values of the involved parameters, these two frequencies, measured in the comoving frame, are derived,

$$\nu _m^{}=\{\begin{array}{ccc}3.7\times 10^9\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^2\left(\frac{ϵ_B}{0.01}\right)^{1/2}\left(\frac{p2}{p1}\right)^2n_{}^{1/2}\gamma ^{1/2}(\gamma 1)^{5/2}\hfill & \mathrm{if}s=0\hfill & \\ 9.5\times 10^{10}\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^2\left(\frac{ϵ_B}{0.01}\right)^{1/2}\left(\frac{p2}{p1}\right)^2\epsilon _{53}^1A_{}^{3/2}\gamma _{0,2}^2\gamma ^{1/2}(\gamma 1)^{5/2}x^1\hfill & \mathrm{if}s=2\hfill & \\ 3.7\times 10^{13}\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^2\left(\frac{ϵ_B}{0.01}\right)^{1/2}\left(\frac{p2}{p1}\right)^2\epsilon _{53}^{1/2}\overline{A}_{}\gamma _{0,2}\gamma ^{1/2}(\gamma 1)^{5/2}x^{3/4}\hfill & \mathrm{if}s=3/2;\hfill & \end{array}$$
(9)
$$\nu _c^{}=\{\begin{array}{ccc}2.8\times 10^{28}\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{3/2}n_{}^{3/2}\gamma ^{3/2}(\gamma 1)^{3/2}t^2\hfill & \mathrm{if}s=0\hfill & \\ 1.7\times 10^{24}\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{3/2}\epsilon _{53}^3A_{}^{9/2}\gamma _{0,2}^6\gamma ^{3/2}(\gamma 1)^{3/2}x^3t^2\hfill & \mathrm{if}s=2\hfill & \\ 2.8\times 10^{16}\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{3/2}\epsilon _{53}^{3/2}\overline{A}_{}^3\gamma _{0,2}^3\gamma ^{3/2}(\gamma 1)^{3/2}x^{9/4}t^2\hfill & \mathrm{if}s=3/2,\hfill & \end{array}$$
(10)
where $`t^{}`$ is in units of 1 s. The synchrotron self-absorption frequency is usually smaller than $`\nu _m^{}`$ and $`\nu _c^{}`$ a few hours after the burst. In this case, the optical depth due to synchrotron self-absorption at $`\nu _a^{}`$ can be approximated by

$$\tau _{ab}(\nu _a^{})\frac{5}{3s}\frac{enr}{B^{}\gamma _p^5}\left(\frac{\nu _a^{}}{\nu _p^{}}\right)^{5/3}1,$$
(11)
where $`\gamma _p=\mathrm{min}(\gamma _m,\gamma _c)`$ and $`\nu _p^{}=\mathrm{min}(\nu _m^{},\nu _c^{})`$ (see, e.g., Panaitescu & Kumar 2000). At early times of the afterglow, all the electrons accelerated behind the shock cool in a timescale smaller than the dynamical expansion timescale of the jet at angle $`\theta `$, i.e., they are in the fast cooling regime, implying $`\nu _c^{}<\nu _m^{}`$, in which case we find the self-absorption frequency,

$$\nu _a^{}=\{\begin{array}{ccc}0.15\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{6/5}\epsilon _{53}^{1/5}n_{}^{8/5}\gamma _{0,2}^{2/5}\gamma ^{6/5}(\gamma 1)^{6/5}x^{3/5}t^{}\hfill & \mathrm{if}s=0\hfill & \\ 7.5\times 10^3\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{6/5}\epsilon _{53}^3A_{}^{24/5}\gamma _{0,2}^6\gamma ^{6/5}(\gamma 1)^{6/5}x^3t^{}\hfill & \mathrm{if}s=2\hfill & \\ 1.3\times 10^{12}\mathrm{Hz}\left(\frac{ϵ_B}{0.01}\right)^{6/5}\epsilon _{53}^{7/5}\overline{A}_{}^{16/5}\gamma _{0,2}^{14/5}\gamma ^{6/5}(\gamma 1)^{6/5}x^{21/10}t^{}\hfill & \mathrm{if}s=3/2.\hfill & \end{array}$$
(12)
At later times, $`\nu _m^{}>\nu _c^{}`$, in which case the self-absorption frequency is given by

$$\nu _a^{}=\{\begin{array}{ccc}4.2\times 10^8\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^1\left(\frac{ϵ_B}{0.01}\right)^{1/5}\left(\frac{p2}{p1}\right)^1\epsilon _{53}^{1/5}n_{}^{3/5}\gamma _{0,2}^{2/5}\gamma ^{1/5}(\gamma 1)^{4/5}x^{3/5}\hfill & \mathrm{if}s=0\hfill & \\ 3.1\times 10^{10}\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^1\left(\frac{ϵ_B}{0.01}\right)^{1/5}\left(\frac{p2}{p1}\right)^1\epsilon _{53}^1A_{}^{9/5}\gamma _{0,2}^2\gamma ^{1/5}(\gamma 1)^{4/5}x^1\hfill & \mathrm{if}s=2\hfill & \\ 3.5\times 10^{13}\mathrm{Hz}\left(\frac{ϵ_e}{0.1}\right)^1\left(\frac{ϵ_B}{0.01}\right)^{1/5}\left(\frac{p2}{p1}\right)^1\epsilon _{53}^{2/5}\overline{A}_{}^{6/5}\gamma _{0,2}^{4/5}\gamma ^{1/5}(\gamma 1)^{4/5}x^{3/2}\hfill & \mathrm{if}s=3/2.\hfill & \end{array}$$
(13)
The synchrotron peak specific luminosity of a ring ($`\theta \theta +d\theta `$) in the comoving frame is given by

$$dL_m^{}=\{\begin{array}{ccc}6.1\times 10^{29}\mathrm{ergs}\mathrm{s}^1\mathrm{Hz}^1\left(\frac{ϵ_B}{0.01}\right)^{1/2}\epsilon _{53}n_{}^{1/2}\gamma _{0,2}^2\gamma ^{1/2}(\gamma 1)^{1/2}x^3\mathrm{sin}\theta d\theta \hfill & \mathrm{if}s=0\hfill & \\ 1.6\times 10^{31}\mathrm{ergs}\mathrm{s}^1\mathrm{Hz}^1\left(\frac{ϵ_B}{0.01}\right)^{1/2}A_{}^{3/2}\gamma ^{1/2}(\gamma 1)^{1/2}\mathrm{sin}\theta d\theta \hfill & \mathrm{if}s=2\hfill & \\ 6.1\times 10^{33}\mathrm{ergs}\mathrm{s}^1\mathrm{Hz}^1\left(\frac{ϵ_B}{0.01}\right)^{1/2}\epsilon _{53}^{1/2}\overline{A}_{}\gamma _{0,2}^1\gamma ^{1/2}(\gamma 1)^{1/2}x^{3/4}\mathrm{sin}\theta d\theta \hfill & \mathrm{if}s=3/2.\hfill & \end{array}$$
(14)
In the case of $`\nu _c^{}<\nu _m^{}`$, the spectrum is (Sari et al. 1998)

$$dL_\nu ^{}^{}=\{\begin{array}{cccc}dL_m^{}(\nu _a^{}/\nu _c^{})^{1/3}(\nu ^{}/\nu _a^{})^2\hfill & \mathrm{if}\nu ^{}<\nu _a^{}\hfill & & \\ dL_m^{}(\nu ^{}/\nu _c^{})^{1/3}\hfill & \mathrm{if}\nu _a^{}<\nu ^{}<\nu _c^{}\hfill & & \\ dL_m^{}(\nu ^{}/\nu _c^{})^{1/2}\hfill & \mathrm{if}\nu _c^{}<\nu ^{}<\nu _m^{}\hfill & & \\ dL_m^{}(\nu _m^{}/\nu _c^{})^{1/2}(\nu ^{}/\nu _m^{})^{p/2}\hfill & \mathrm{if}\nu ^{}>\nu _m^{}.\hfill & & \end{array}$$
(15)
If $`\nu _c^{}>\nu _m^{}`$, the spectrum becomes (Sari et al. 1998)

$$dL_\nu ^{}^{}=\{\begin{array}{cccc}dL_m^{}(\nu _a^{}/\nu _c^{})^{1/3}(\nu ^{}/\nu _a^{})^2\hfill & \mathrm{if}\nu ^{}<\nu _a^{}\hfill & & \\ dL_m^{}(\nu ^{}/\nu _m^{})^{1/3}\hfill & \mathrm{if}\nu _a^{}<\nu ^{}<\nu _m^{}\hfill & & \\ dL_m^{}(\nu ^{}/\nu _m^{})^{(p1)/2}\hfill & \mathrm{if}\nu _m^{}<\nu ^{}<\nu _c^{}\hfill & & \\ dL_m^{}(\nu _c^{}/\nu _m^{})^{(p1)/2}(\nu ^{}/\nu _c^{})^{p/2}\hfill & \mathrm{if}\nu ^{}>\nu _c^{}.\hfill & & \end{array}$$
(16)
The observed total flux density at observed frequency $`\nu `$ is given by

$$F_\nu =_0^{\theta _m}\frac{dL_{\nu \gamma (1\beta \mu )}^{}}{4\pi D_L^2\gamma ^3(1\beta \mu )^3},$$
(17)
where $`D_L`$ is the luminosity distance to the source. For a flat universe with $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, the luminosity distance to the source $`D_L=2c/H_0(1+z\sqrt{1+z})=1.8\times 10^{28}\mathrm{cm}[(1+z)/2]^{1/2}[(\sqrt{1+z}1)/(\sqrt{2}1)]`$.

## 3 Numerical Results

Integrating equations (7), (8), and (17) numerically, we can easily obtain an afterglow light curve. Figures 1-3 exhibit different afterglow light curves at R-band ($`\nu =4.4\times 10^{14}`$ Hz) for different $`k`$ in the homogeneous ISM ($`n_{}=1`$), in the wind medium with $`s=2`$ ($`A_{}=1`$), and in the wind medium with $`s=3/2`$ ($`\overline{A}_{}=0.01`$), respectively. We choose the remaining parameters: $`\epsilon _0=10^{53}\mathrm{ergs}\mathrm{sr}^1`$, $`M_0=5\times 10^4M_{}\mathrm{sr}^1`$, $`\theta _0=M_0c^2/(\epsilon _0+M_0c^2)10^2`$ rad, $`\theta _m=0.3`$ rad, $`p=2.5`$, $`ϵ_e=0.1`$, $`ϵ_B=0.01`$, and $`z=1`$. In each figure, the lines from top to bottom correspond to $`k=0`$, $`0.5`$, 1, 2 and 3 respectively. From these figures, we can see the following features:

1. At very early times, the afterglow from a jet expanding in the ISM brightens rapidly but the afterglow flux for the wind case is approximately a constant, no matter whether the jet is isotropic or anisotropic. This may provide a way of distinguishing between the environments as well as GRB progenitors in the near future because rapid, accurate locations of HETE-2 will be able to allow follow-up observations on very early-time X-ray and optical afterglows.
2. For an isotropic jet ($`k=0`$) expanding in the ISM, there is a sharp break in the light curve at later times. The time at which this break occurs is close to the analytical result,

$$t_{\mathrm{jet}}2\times 10^6(\epsilon _{53}/n_{})^{1/3}(\theta _m/0.3)^{8/3}\mathrm{s},$$
(18)
at which the Lorentz factor of the jet is equal to the inverse of its half opening angle (Mészáros & Rees 1999). But for an isotropic jet expanding in the wind, the break of the light curve is weaker and smoother.
3. In each type of medium, the break of the light curve becomes weaker and smoother as $`k`$ increases. When $`k=2`$ or 3, the break seems to disappear but the afterglow decays rapidly. This result shows that the emission from expanding, highly anisotropic jets may provide an explanation for some rapidly fading afteglows whose light curves have no break (e.g., GRB 991208; see the next section). We note that in the preliminary analysis of Mészáros et al. (1998), an afterglow decays as one single power-law for any $`k`$. The reason for this difference will be discussed in the final section.
4. At very late times, the slopes of the light curves are approximately equal. We also find that the slopes become steeper with the electron distribution index $`p`$.

In Figures 4-6 we present the observed radiation spectra computed for a sequence of time for a jet with $`k=3`$ in three kinds of medium for the same parameters as in Figures 1-3. It can be seen that even for such an anisotropic jet considered here, the observed radiation spectra can still be divided into four parts, which can be well described from low-energy to high-energy bands by power-law functions $`F_\nu \nu ^\beta `$, with $`\beta =2`$, $`1/3`$, $`(p1)/2=0.75`$ and $`p/2=1.25`$ at $`t10^4`$ s, respectively. Two adjacent spectrum portions are still joined very sharply. These results may be due to that the observed radiation spectra are mainly contributed by the $`\theta <\theta _0`$ part of the jet.

## 4 The Afterglow of GRB 991208

GRB 991208 was a long burst with a duration of $`60`$ s and a fluence of $`10^4\mathrm{erg}\mathrm{cm}^2`$ ($`>25`$ keV) (Hurley et al. 1999). Its redshift was measured as $`z=0.7055\pm 0.0005`$ (Djorgovski et al. 1999), and thus its isotropic energy in $`\gamma `$-rays is estimated: $`E_{\mathrm{iso}}1.3\times 10^{53}`$ ergs. Because its afterglows at optical, millimeter and radio wavelengths were detected and unusually bright (Sagar et al. 2000a; Castro-Tirado et al. 2000; Galama et al. 2000), a detailed study of this burst may provide new clues regarding the origin of the GRB phenomenon.

Sagar et al. (2000a) gave the observed spectral and temporal indexes of the optical afterglow $`\beta _{\mathrm{ob}}=0.75\pm 0.03`$ and $`\alpha _{\mathrm{ob}}=2.2\pm 0.1`$ while Galama et al. (2000) presented well-sampled spectra and light curves between radio and millimeter wavelengths for a two-week period, and obtained, for the first time, the evolution of the synchrotron self-absorption frequency $`\nu _at^{0.15\pm 0.12}`$, the peak frequency $`\nu _mt^{1.7\pm 0.4}`$, and the peak flux density $`F_mt^{0.47\pm 0.11}`$. The existence of one single power law decay for the R-band afterglow implies that the afterglow might be produced by one of the following models: (i) an isotropic jet with sideways expansion, (ii) a relativistic fireball expanding in the ISM, and (iii) a relativistic fireball expanding in the wind. Galama et al. (2000) argued that model (i) can explain the observations, but models (ii) and (iii) are ruled out from $`\alpha _{\mathrm{ob}}`$ and $`\beta _{\mathrm{ob}}`$.

If $`p=2.5`$, inferred from $`\beta _{\mathrm{ob}}=(p1)/2=0.75\pm 0.03`$, Figures 7 and 8 provide two good fits to the observed R-band afterglow ligh curve based on our jet model with $`k=3`$ in cases of media (ISM and wind with $`s=2`$), respectively. This shows that GRB 991208 might arise from a highly anisotropic jet. However, we note that these two cases may not provide further fits to the observed radiation spectra because the evolution of the theoretical synchrotron self-absorption frequency, peak frequency and peak flux density is inconsistent with the observations by Galama et al. (2000).

Now we assume a generic wind case: $`nr^s`$. In this case, we have derived the synchrotron self-absorption frequency $`\nu _at^{3s/[5(4s)]}`$ if $`\nu _a<\nu _m`$, and the peak frequency $`\nu _mt^{3/2}`$, and the peak flux density $`F_mt^{s/(82s)}`$ in the ultra-relativistic phase (Dai & Lu 1998). If $`s=3/2`$, we easily find $`\nu _at^{0.36}`$, $`\nu _mt^{1.5}`$, and $`F_mt^{0.3}`$. These scaling laws are in approximate accord with those obtained from the observations. This preliminary result prompts us to re-consider a highly anisotropic jet expanding in an $`s=3/2`$ wind to fit the afterglow data of GRB 991208. In Figure 9 we first show a satisfactory fit to the observed R-band light curve. It is interesting to note that such a wind has been inferred in the circumstellar medium of SN 1993J (whose progenitor is a red supergiant) by Fransson et al. (1996), and is possibly caused by a variation of the mass-loss rate from the progenitor or by a non-spherical geometry. In Figure 10 we further present a fit to the radio-to-optical-band afterglow spectrum observed on 1999 December 15.5 UT for the same parameters as in Figure 9. Therefore, we can conclude that GRB 991208 might come from a highly anisotropic jet expanding in the wind environment ($`s=3/2`$) from a red supergiant.

## 5 Discussion and Conclusions

The model presented in section 2 clearly includes the following assumptions. First, the initial opening angle of a highly anisotropic jet (e.g., $`k=3`$), $`\theta _m`$, was taken to be $`0.3`$, and the observer’s angle ($`\theta _{\mathrm{obs}}`$) between the axis line of the jet and the line of sight to be zero. Because $`dL_m\epsilon (\theta )`$ in the ultra-relativistic limit, the emission from the jet mainly arises from the shock-accelerated electrons moving along the axis line. Thus, the flux density of the emission is weakly dependent of $`\theta _m`$, but strongly depends on $`\theta _{\mathrm{obs}}`$. For example, the flux density at $`\theta _{\mathrm{obs}}=0.3`$ is about four orders of magnitude smaller than that of the emission at $`\theta _{\mathrm{obs}}=0`$. GRB 980425 was likely such an off-axis burst surrounded by an $`s=2`$ wind medium (Nakamura 1999) because its gamma-ray luminosity of GRB 980425 was also approximately four orders of magnitude smaller than that of GRB 970228/980326 and its X-ray afterglow slowly declined. Second, we considered jets of a fixed opening angle and neglected the effect of sideways expansion based on two observational facts. As shown analytically by Rhoads (1999) and Sari, Piran & Halpern (1999), this effect can lead to a sharp break in an afterglow light curve. However, numerical studies by many authors (e.g, Panaitescu & Mészáros 1999; Moderski, Sikora & Bulik 2000; Huang et al. 2000a, b; Kumar & Panaitescu 2000; Gou et al. 2000; Wei & Lu 2000) show that the actual break, when two effects such as the equal-time surface and detailed dynamics of the jet are considered, is much weaker and smoother than that predicted analytically. Finally, we used the adiabatic solution for kinematics of the jet. This argument is based on the radiative efficiency of the jet, $`f=ϵ_e[t_{\mathrm{syn}}^1/(t_{\mathrm{syn}}^1+t_{\mathrm{exp}}^1)]`$, defined by Dai et al. (1999). At very early times, the cooling timescale due to synchrotron radiation, $`t_{\mathrm{syn}}^{}`$, may be much less than the expansion timescale, $`t_{\mathrm{exp}}^{}`$, and $`f0.1(ϵ_e/0.1)`$; but at late times, $`t_{\mathrm{syn}}^{}t_{\mathrm{exp}}^{}`$ and thus $`f0.1(ϵ_e/0.1)(t_{\mathrm{exp}}^{}/t_{\mathrm{syn}}^{})1`$. Therefore, it can be seen that the energy losses due to synchrotron radiation are insignificant during the whole evolution stage of the jet.

Our numerical results show that for an isotropic jet ($`k=0`$) expanding in the ISM, there is a sharp break in the light curve at late times, but for an isotropic jet expanding in the wind, the break of the light curve is weaker and smoother. This result is easily understood in terms of the analytical model of Dai & Lu (1998) and Chevalier & Li (1999, 2000): at early times of the jet evolution, $`\gamma >\theta _m^1`$, in which stage the jet is spherical-like and thus the decay index of the afterglow $`\alpha _1=[s/(82s)+3(p1)/4]`$ in the slow cooling regime; but at later times, $`\gamma <\theta _m^1`$, the decay index becomes $`\alpha _2=[(6s)/(82s)+3(p1)/4]`$ due to the edge effect. These indexes imply the steepening of the light curve $`\mathrm{\Delta }\alpha \alpha _1\alpha _2=(3s)/(4s)=0.5`$ for $`s=2`$, 0.6 for $`s=3/2`$, or 0.75 for $`s=0`$. Therefore, the break in the light curve of the afterglow from a jet expanding in a wind is weaker and smoother than that for a homogeneous medium case. Our results also show that for an anisotropic jet, one break seems to appear in the light curve at small $`k`$ but becomes weaker and smoother as $`k`$ increases. This conclusion is different from that presented by Mészáros et al. (1998), who found that an afterglow decays as one single power-law for any $`k`$. The reason for this difference is that we took into account the light travel effects related to different sub-jets within the jet but Mészáros et al. (1998) didn’t. It is easily understood that an anisotropic jet with a small $`k`$ can be treated as a quasi-isotropic jet with a half opening angle $`\theta _m`$, and thus the temporal decay of the flux density of radiation from such a jet will begin to expedite at $`t=t_{\mathrm{jet}}`$, which is given by equation (18) for the ISM case. However, for a highly anisotropic jet (e.g., $`k=2`$ or 3), the contribution of radiation from $`\theta >\theta _0`$ is in fact much smaller than that from $`\theta <\theta _0`$ due to the fact that the energy per unit solid angle ($`dE/d\mathrm{\Omega }`$) within this jet decreases rapidly with increasing $`\theta `$, and thus, the light curve of the resulting radiation should begin to steepen at a very early time, $`t230(\epsilon _{53}/n_{})^{1/3}(\theta _0/10^2)^{8/3}\mathrm{s}`$, which is also estimated from equation (18) in the ISM case but for $`\theta _m`$ being replaced by $`\theta _0`$. Therefore, one cannot see a broken light curve at late times for such a highly anisotropic jet.

In summary, the energy per unit solid angle within a realistic jet may be a power-law function of the angle ($`\theta ^k`$). Such anisotropic jets are expected in the collapsar model of MacFadyen et al. (2000). We numerically calculated light curves and spectra of the emission from anisotropic jets expanding either in the interstellar medium (ISM) or in the wind medium. We also took into account two kinds of wind: one ($`nr^{3/2}`$) possibly from a typical red supergiant and another ($`nr^2`$) possibly from a Wolf-Rayet star. Two of the main results of this work are that (i) at very early times, the afterglow from a jet expanding in the ISM brightens rapidly but the afterglow flux for the wind case is approximately a constant, no matter whether the jet is isotropic or anisotropic. Based on this conclusion, future observations led by HETE-2 on very early-time optical afterglows may be able to distinguish between the environments as well as GRB progenitors. (ii) In each type of medium, one break appears in the late-time afterglow light curve for small $`k`$, but becomes weaker and smoother as $`k`$ increases. When $`k2`$, the break seems to disappear but the afterglow decays rapidly. Therefore, one expects that the emission from expanding, highly anisotropic jets provides an explanation for some rapidly fading afteglows whose light curves have no break. We also presented good fits to the optical afterglow flux data of GRB 991208. Finally, to interpret the observed radio-to-optical-band afterglow data of this burst, we argued that it might arise from a highly anisotropic jet expanding in the wind ($`nr^{3/2}`$) from a red supergiant.

We would like to thank the referee for several valuable comments that allowed us to improve the manuscript, and D. M. Wei for reading carefully the manuscript. This work was supported by the National Natural Science Foundation of China (grant 19825109) and the National 973 Project.


format: text

text:
# A New Interpretation of Compensate Effect

## 1 Introduction

Recently Zhao et al. developed an interpretation of Hawking effect as a compensate effect of the coordinate scale transformation . In this understanding the temperature appears to be the $`x^1`$ component of a pure gauge. It is of importance that the zero-temperature space-time must be singled out. Otherwise an uncertainty in this scheme will arises. It is confident that all is determined whenever the metric of space-time is given. Therefore we look a scheme that the compensate effect is determined by the space-time itself. In section two a brief review of Zhao’s idea is given via a simple example. In section three the temperature is presented to be $`x^1`$ component of a pure gauge on the horizon in generalized tortoise coordinates. It is proved for a general case and unnecessary to point out a zero-temperature space-time.

## 2 Zhao’s Compensate Effect

It is known that via the technique of conformal flat the line element in $`x^0,x^1`$ subspace is conformal to Minkowski space near the event horizon. Zhao et al. analyze the conformal factor and present the compensate effect. We now use Schwarzschild metric in advanced Eddington-Finkelstein coordinate to go over the skeleton.

$$ds^2=(1\frac{2m}{r})dv^2+2dvdr+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2).$$
(1)
In tortoise coordinates

$$r_{}=r+\frac{1}{2\kappa }ln(\frac{r}{2m}1),$$
(2)
the line element can be rewritten as

$$ds^2=\mathrm{\Omega }_2^2ds_2^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2),$$
(3)
in which

$$\mathrm{\Omega }_2^2=1\frac{2m}{r}=\frac{e^{2\kappa (r_{}r)}}{2\kappa r}$$
(4)
$$ds_2^2=dv^2+2dvdr_{}.$$
(5)
Introducing the generalized null Kruskal coordinates

$$U=\frac{1}{\kappa }e^{\kappa u},V=\frac{1}{\kappa }e^{\kappa v},$$
(6)
where the retarded Eddington-Finkelstein coordinate is defined as $`u=v2r_{}`$. Its relation to the kruskal coordinates is

$$T=\frac{1}{2}(V+U),R=\frac{1}{2}(VU).$$
(7)
The line element Eq.(3)can be written as

$$ds^2=\mathrm{\Omega }_1^2ds_1^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2),$$
(8)
in which

$$\mathrm{\Omega }_1^2=\frac{e^{2\kappa r}}{2\kappa r}$$
(9)
$$ds_2^2=dV^2+2dVdR.$$
(10)
Since the line element keeps invariant, from Eqs.(3)(8), it is clear

$$ds_2^2=\frac{ds_1^2}{\mathrm{\Omega }^2},\mathrm{\Omega }^2=\frac{\mathrm{\Omega }_2^2}{\mathrm{\Omega }_1^2}.$$
(11)
Eq.(11)can be regarded as a conformal isometry.

Now consider a infinitesimal 2-dimensional vector $`B_\mu `$ whose proper length is

$$L=|ds|=\sqrt{B_\mu B^\mu }=\sqrt{g_{\mu \nu }B^\mu B^\nu },\mu ,\nu =0,1$$
(12)
whose coordinate lengths are defined respectively

$$l_1=|ds_1|,l_2=|ds_2|.$$
(13)
Obviously

$$L=\mathrm{\Omega }_1l_1=\mathrm{\Omega }_2l_2.$$
(14)
Since $`L`$ is invariant under parallel transfer in two dimensional subspace ,

$$\delta lnl_1=dln\mathrm{\Omega }_1,\delta lnl_2=dln\mathrm{\Omega }_2$$
(15)
is obtained. On the other hand, the contractions of the affine connection of 2-dimensional subspace have the relations

$$\mathrm{\Gamma }_{\alpha \mu }^\alpha dx^\mu =d(ln\sqrt{g_1})=2d(ln\mathrm{\Omega }_1)$$
(16)
$$\stackrel{~}{\mathrm{\Gamma }}_{\alpha \mu }^\alpha dx^\mu =d(ln\sqrt{g_2})=2d(ln\mathrm{\Omega }_2)$$
(17)
where $`g_1`$ and $`g_2`$ are the metric determinants. Define 1-form $`A_1`$ and $`A_2`$

$$A_1=\delta lnl_1,A_2=\delta lnl_2,$$
(18)
which rapidly results in

$$A_1=\frac{1}{2}\mathrm{\Gamma }_{\alpha \mu }^\alpha ,A_2=\frac{1}{2}\stackrel{~}{\mathrm{\Gamma }}_{\alpha \mu }^\alpha .$$
(19)
It can be find that $`A_1`$ and $`A_2`$ are respectively the relative change rates of the coordinate lengths $`l_1`$ and $`l_2`$. The coordinate length has the same scale as the coordinate time in the two dimensional subspace, so both $`A_1`$ and $`A_2`$ reflect the relative change rates of the coordinate times. Since Eq.(25) the $`A`$ transform as a connection in the scale transformation

$$A_2=A_1+dln\mathrm{\Omega }.$$
(20)
The field strength 2-form

$$F=0$$
(21)
is immediately obtained since connection $`A`$ is an exact-form and $`dd=0`$. It is easily obtained the property of pure gauge $`A`$

$$\kappa =\frac{ln\mathrm{\Omega }}{r_{}}.$$
(22)
It is the $`x^1`$ component of connection that can be regarded as Hawking temperature.

## 3 Another Interpretation of Conformal Factor

It is wonderful that the Hawking effect can be regarded as the compensate effect of the coordinates. There is also disadvantage: the zero-temperature space-time must be pointed out manually. Otherwise this scheme can not determine the final result uniquely since the method itself does not point out which space has the temperature $`T=\frac{\kappa }{2\pi K_B}`$. Here we propose another interpretation. The advantages of the new interpretation is that the $`\kappa `$ is determined by the space-time itself and the zero-temperature space-time is unnecessarily to be pointed out. We directly give a general proof in advanced Eddington-Finkelstein coordinates.

Consider the most general case that $`\xi =\xi (x^0,x^2,x^3)`$ in Eddington coordinates. Tortoise coordinates transformation is

$$x_{}^1=x^1+\frac{1}{2\kappa }ln(x^1\xi )$$
(23)
with other components invariant.

$$dx_{}^1=(1+\frac{1}{ϵ})dx^1\frac{\xi _\nu ^{}}{ϵ}dx^\nu ,$$
(24)
in which $`ϵ=2\kappa (x^1\xi )`$ and $`\xi _\nu ^{}=\frac{\xi }{x^\nu }`$.

The metric is then obtained in terms of tortoise coordinates

$$ds^2=(\frac{ϵg_{11}\xi _0^{}}{(1+ϵ)^2}+\frac{ϵg_{10}}{1+ϵ})[\frac{\frac{g_{11}}{(1+ϵ)^2}\xi _0^{}\xi _0^{}+\frac{2g_{10}}{1+ϵ}\xi _0^{}+g_{00}}{\frac{ϵg_{11}}{(1+ϵ)^2}\xi _0^{}+g_{10}\frac{ϵ}{1+ϵ}}dx^0dx^0+2dx^0dx_{}^1]+(others).$$
(25)
The technique of conformally flat require around the horizon the coefficient of $`dx^0dx^0`$ in $`[]`$ in Eq.(25) being $`1`$

$$\frac{\frac{g_{11}\xi _0^{}\xi _0^{}+2g_{10}\xi _0^{}+g_{00}}{ϵ}+2g_{00}+2g_{10}\xi _0^{}+ϵg_{00}}{g_{10}+g_{11}\xi _0^{}+ϵg_{10}}=1.$$
(26)
So the line element in the two dimensional subspace $`x^0,x_{}^1`$ is conformal to Minkowski space

$$ds^2=\mathrm{\Omega }^2(dx^0dx^0+2dx^0dx_{}^1),$$
(27)
where the conformal factor $`\mathrm{\Omega }`$ is

$$\mathrm{\Omega }^2=\frac{ϵg_{11}\xi _0^{}}{(1+ϵ)^2}+\frac{ϵg_{10}}{1+ϵ}.$$
(28)
Putting Eqs.(26)(28), the technique of conformally flatness at the event horizon obtains

$$ds^2=\mathrm{\Omega }^2|_{x^1\xi }(dx^0dx^0+2dx^0dx_{}^1)+(otherterms),$$
(29)
in which when $`x^1`$ approach the horizon $`\xi `$ only keeps the first order of $`ϵ`$, $`\mathrm{\Omega }^2`$ equals to

$$\mathrm{\Omega }^2|_{x^1\xi }=ϵ(g_{10}+g_{11}\xi _0^{}).$$
(30)
Finally the relation of $`\kappa `$ and when $`x^1\xi `$ the $`x^1`$ component of connection $`A`$ is obtained as

$$\frac{ln\mathrm{\Omega }}{x_{}^1}=\frac{1}{2}[\frac{lnϵ}{x_{}^1}+\frac{ln(g_{10}+g_{11}\xi _0^{})}{x_{}^1}]=\frac{1}{2}\frac{lnϵ}{x_{}^1}=\kappa ,$$
(31)
where the condition is supposed that the components $`g_{10}`$ and $`g_{11}`$ is not divergent on the event horizon $`\xi `$, i.e. their derivatives with tortoise coordinate $`x_{}^1`$ vanish when $`x^1\xi `$. In fact this hypothesis is not strong. For many models of black holes in advanced Eddington-Finkelstein coordinates, $`g_{11}=0`$ and $`g_{10}`$ is a constant.

It is emphasized that although the Eq.(31) always holds, generally speaking, for the most general horizon $`\xi =\xi (x^0,x^2,x^3)`$ the $`\kappa `$ determined by this way is not the one determined by the method of Damour-Ruffini-Zhao


text:
# The shape of pulsar radio beams

## 1 Introduction

Pulsars are generally believed to be rotating neutron stars in which the observed pulses result from one or more emission beams sweeping across the Earth as the star rotates. Observations of radio polarisation \[Radhakrishnan & Cooke 1969\] led to the magnetic-pole model in which the emission beam was centred on the magnetic axis of a predominantly dipole magnetic field. On the assumption that the emission is directed radially, the observed pulse profile reflects the variations in emission intensity along a line of constant rotational latitude. Although more than 1000 pulsars are now known \[Taylor, Manchester & Lyne 1993, Lyne et al. 2000\] the shape of pulsar radio beams remains controversial.

There are two main areas of uncertainty. One concerns the outline shape of the radio beam. Originally assumed to be circular, some investigations (e.g. Narayan & Vivekanand 1983) argued for an elliptical beam extend in the latitude direction with a large axial ratio. In contrast, Biggs (1990) suggested that the beam was slightly compressed in the latitudinal direction. Most other investigations (e.g. Lyne & Manchester 1988; Björnsson 1998; Gil & Han 1996) have concluded that the emission beam is essentially circular, and we will assume this in the present investigation.

The second area of uncertainty concerns the form of the beam pattern. Early observations showed that there are two or more pulse components in many pulsars (e.g. Lyne, Smith & Graham 1971; Manchester 1971). Ruderman & Sutherland (1975) presented a detailed model for pulsar radio emission in which the beam had the shape of a hollow cone, more intense around the periphery, corresponding to the last open field lines emanating from the polar cap region. The radio emission was attributed to curvature radiation by positrons moving these field lines. Such radiation is linearly polarised in the plane of curvature of the magnetic field and hence the model naturally explained the smooth sweep of polarisation position angle seen in many pulsars, as well as the common occurrence of double-peaked pulse profiles.

Backer (1976) first discussed the more-or-less central pulse component seen in many pulsars and suggested that it resulted from an axial beam. Rankin (1983) made the distinction between this central or ‘core’ emission and the outer or ‘conal’ emission, and showed that these two components of the profile had rather different properties. At least in short-period pulsars, core emission has a steeper spectrum than conal emission \[Rankin 1983, Lyne & Manchester 1988\] and the fluctuation properties of the two types of emission are very different \[Rankin 1986\]. This model was subsequently extended to have two or more coaxial conal emission zones, either to account for multiple-component profiles \[Gil & Krawczyk 1997, Qiao & Lin 1998\] or for the appearance of components at different apparent radii from the conal axis \[Rankin 1993a, Kramer et al. 1994, Mitra & Deshpande 1999\]. These latter analyses depend on a knowledge of the angle between the beam and rotation axes ($`\alpha `$), usually computed from the observed width of a ‘core’ component. Unfortunately, the determination of this angle is very model dependent and this results in large uncertainties in the derived radii.

Based on a large sample of pulse shape and polarisation data, Lyne & Manchester (1988) found that pulse components in multiple-component profiles were not symmetrically located within the beam boundary. Furthermore, they were generally of very different intensities and, in some cases, missing altogether. Their results were consistent with components being randomly located within the beam boundary, leading to the ‘patchy-beam’ model. Gould (1994) and Gould, Lyne & Smith (2000) confirmed these conclusions with an even larger data set.

Manchester (1995) suggested an interpretation of these results in which the observed pulse profile is the product of a ‘window function’, which is a function of pulse period and radio frequency but common to all pulsars, and a ‘source function’ which determines the strength of the emission at a given point within the beam and is different for every pulsar. For example, the source function may be determined by the density or energy of the plasma beams along different field lines, whereas the window function is determined by the properties of the emission mechanism.

In this paper, we determine the two-dimensional shape of the window function at frequencies around 1 GHz by averaging the observed pulse profiles of multi-component pulsars having reliable polarisation data. This analysis depends only on the observed width of the pulse profile and the ‘impact parameter’, that is the minimum angle between the line of sight and the beam axis, usually given the symbol $`\beta `$. In contrast to the beam inclination angle $`\alpha `$, the impact parameter is normally well determined from the maximum rate of change of position angle at the profile centre \[Lyne & Manchester 1988\]. We assume that emission is seen across the entire polar cap and choose only pulsars with profiles where this appears to be the case. In Section 2 we describe the data set and our methods of determining the two-dimensional beam pattern. Results are presented and discussed in Section 3 and Section 4 summarises our conclusions.

## 2 The data set and analysis method

In order to reliably determine the two-dimensional beam shape we choose pulsars which have pulse profiles with two or more pulse components, good signal/noise ratio and good quality polarisation data.

For most pulsars with two or more pulse components, we see emission from right across the polar cap. The polarisation information allows a measurement of the normalised impact angle, $`\beta _n`$, that is, the impact parameter $`\beta `$ expressed as a fraction of the beam radius. This parameter is computed assuming $`\alpha =90\mathrm{°}`$, but it is not very sensitive to the actual value of $`\alpha `$ (Lyne & Manchester 1988). The polarisation data also allow a check on the assumption that we are seeing emission right to both edges of the beam. Only those pulsars for which the maximum rate of change of position angle is close to the pulse centre were included in the sample. Single component pulsars are not included in the sample as it is often difficult to reliably determine if they are predominantly core emission (low $`\beta _n`$) or from the outer edge of the beam (high $`\beta _n`$).

Millisecond pulsars were also excluded from the sample as they generally have very wide and complex profiles and polarisation variations which do not fit the simple rotating-vector model (e.g. Navarro et al. 1997, Stairs et al. 1999). It is therefore difficult or impossible to obtain reliable values of $`\beta _n`$ for these pulsars. Very short-period pulsars also commonly have wide or interpulse profiles and are excluded from the sample for similar reasons.

Pulse profiles in a given pulsar usually evolve with frequency, with some components becoming stronger and other becoming weaker or even disappearing at higher or lower frequencies. At low frequencies, the core component tends to dominate the profile, whereas at high frequencies the conal emission is generally dominant. To ensure relative uniformity in the data set, we restricted our attention to profiles at frequencies between 600 and 1600 MHz. The results we obtain therefore represent pulsar emission at frequencies around 1 GHz.

We checked the profiles available on the pulsar profile database 111See http://www.mpifr-bonn.mpg.de/div/pulsar/data/ of European Pulsar Network maintained at Max-Planck-Institut für Radioastronomie \[Lorimer et al. 1998\], which contains pulsar profiles from more than 50 papers, including the large datasets from Gould & Lyne (1998), Manchester, Han & Qiao (1998) and Hoensbroech & Xilouris (1997). Normalised impact parameters were obtained from Lyne & Manchester (1988), Rankin (1993), Gould (1994) and Manchester et al. (1998). The final data set of 87 pulsars which satisfy the above criteria is listed in Table 1.

To determine the average beamshape, we first defined the duration of each pulse profile by taking the points at which the signal rose above $`3\sigma `$, where $`\sigma `$ is the rms deviation of the off-pulse noise. Next, the mean intensity of the observed emission within the pulse was normalised to unity. The normalised impact parameter, $`\beta _n`$, gives the offset of the locus of line of sight across the polar cap in the latitude direction as a fraction of the beam radius. We assumed a circular beam of unity radius and a straight trajectory of the emission point across the polar cap. The beamshape is assumed to be symmetric about the equator, i.e., the sign of $`\beta _n`$ is ignored.

The profile intensity was interpolated on to an $`xy`$ array representing a semi-circular beam of unity radius along a line at $`y=\beta _n`$ between $`x`$ values of $`\pm \sqrt{1\beta _n^2}`$ using a five-point polynomial interpolation routine \[Press et al. 1992\]. To allow for uncertainties in the value of $`\beta _n`$, possible non-linearities in the beam trajectory across the polar cap, and to represent the finite width of a subpulse beam, the pulse profile is broadened in latitude $`(y)`$ assuming a guassian form, $`\mathrm{exp}[(\beta \beta _n)^2/0.1]`$, between $`\mathrm{\Delta }\beta =\pm 0.3`$. Examples of computed beam components for two pulsars are shown in Fig. 1.

## 3 Results and Discussion

Integrating all pulse profiles in the sample gives the beam shape shown in Fig. 2. This figure represents the probabilty of observing beam components at each location. However, it is strongly affected by the sample distribution of $`\beta _n`$. As shown in Fig. 3, this has a strong peak at $`\beta _n0.25`$, indicating that multi-component pulse profiles are most likely to be detected at about this impact parameter. The deficit at $`\beta _n0`$ is an observational selection effect which results from the smearing of rapid position angle changes near the profile centre due to finite subpulse beamwidths and instrumental broadening of the profile. Few multi-component profiles are observed at high $`\beta _n`$.

The beam distribution shown in Fig. 2 was normalised to correct for the non-uniform distribution of $`\beta _n`$ using the gaussian-smoothed distribution shown in Fig. 3, giving the final average beam shape shown in Fig. 4. Since there are only one or two pulsars in bins with $`\beta _n>0.8`$ (Fig. 3), there is large statistical uncertainty in the average profile shape for these bins. Therefore we do not plot this region of the normalised beam pattern.

The overall impression given by Fig. 4 is that of a ‘patchy’ beam. There are some systematic features, but in general the intensity is relatively uniform across the whole beam. The central or ‘core’ feature is significantly displaced toward later longitudes. The next most prominent property of the average beam after the core component is the presence of a few patches of slightly enhanced emission. Most of these are located within an irregular and rather broad annular maximum at a normalised radius of about 0.7 which is somewhat stronger on the leading side of the beam. Underlying this is more-or-less uniform emission over the whole of the emission beam. The decline at the beam edge is not very sharp.

Although there is no evidence for double or multiple cones in Fig. 4, our present data set is dominated by two- and three-component pulsars. A larger sample of pulsars with more than three components is needed to reliably distinguish multiple-cone models from a patchy beam models. However, our results suggest that, if multiple cones exist, they are at different radii relative to the beam radius in different pulsars. They also show that the conal emission is not confined to a single annular region at the beam boundary.

These results are consistent with the idea that components in pulse profiles are largely determined by the distribution of sources across the polar cap – the ‘source function’. The rather smooth distribution distribution of intensity in the mean beam shape suggests that these source regions are randomly distributed for different pulsars as suggested by Lyne & Manchester (1988) and Manchester (1995) – the ‘patchy-beam’ model. The number of identifiable source components across the profile is mostly limited by the finite width of subpulse beams. Except in a few cases with very high signal/noise ratio, e.g. PSR B0740$``$28 \[Kramer 1994\], it is generally not possible to identify more than four or five peaks or components, and often only one or two. Some very wide components identified by gaussian fitting, e.g. PSRs B2319+60, B2021+51 \[Kramer 1994\], are most probably regions of distributed emission producing overlapping subpulse beams.

We believe Fig. 4 is a good representation of the mean radio beam (for $`\beta _n<0.8`$) emitted at frequencies around 1 GHz by ‘typical’ pulsars, that is, the pulsars of medium or long period which dominate the sample of known radio pulsars. It represents the ‘window function’ in the model of Manchester (1995) which is determined by the effective gain or efficiency of the radio emission process.

## 4 Conclusions

We have computed the average radio beamshape at frequencies about 1 GHz of pulsars with medium to long periods. This beamshape has a peak near its centre and a mild, broad and rather irregular enhancement at a normalised beam radius of about 0.7, but is otherwise rather uniform. The decline at the beam edge is gradual.

These results suggest that the presence and location of profile components are determined by a ‘source function’ which varies randomly from pulsar to pulsar. The summing of these randomly distributed components results in a relatively uniform average beam profile which we interpret as the ‘window function’ representing the properties of the emission process common to all longer-period pulsars.

## Acknowledgments

We thank Prof. Qiao Guojun and Zhao Yongheng for helpful discussions and the referee for helpful comments. This work was initiated during JLH’s visit to the ATNF during 1997. JLH thanks the Su Shu Huang Astrophysics Research Foundation of CAS and the exchange program between CAS and CSIRO for support of visits in 1997 and 1999, respectively. His work in China is supported by the National Natural Science Foundation of China and the National Key Basic Research Science Foundation. The profile data were obtained from the pulsar database of the European Pulsar Network at the Max-Planck-Institut für Radioastronomie. The Australia Telescope is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO.


text:
# References

DEDUCTION OF PLANCK’S FORMULA FROM MULTIPHOTON STATES

Luis J. Boya111Permanent Address: Departmento de Física Teórica, Facultad de Ciencias. Universidad de Zaragoza. E-50009 Zaragoza, Spain. E-mail: luisjo@posta.unizar.es, Ian M. Duck222Physics Department, Rice University, Houston, TX. 77251-1892  and E. C. G. Sudarshan 333E-mail: sudarshan@utaphy.ph.utexas.edu

Center for Particle Physics

Department of Physics

The University of Texas at Austin

Austin, Texas 78712-1081, USA

Key Words: photons, radiation law, old quantum theory.

## Abstract

We obtain the black body radiation formula of Planck by considering independent contributions of multiphoton entities

1. Einstein introduced the light quantum hypothesis in 1905  by noticing that the entropy calculated from the Wien formula for the black body radiation

$$I_W(\nu ,T)=a\nu ^3\mathrm{exp}(b\nu /T),$$
(1)
which is the high-frequency, low-density limit of the complete Planck’s formula

$$I(\nu ,T)=\frac{8\pi h\nu ^3}{c^3}\frac{1}{e^{(h\nu /kT)}1},$$
(2)
this entropy had the same form as the entropy of a “gas” of independent particles, light quanta (later called photons). Viceversa, starting from a gas of independent light quanta one arrives only to Wien’s formula, as shown since 1910 by many people .

The question arises as how to get the complete radiation formula from the corpuscular concept of the photon. In this letter we present what we think is about the simplest derivation of Planck’s radiation formula, starting from the idea of “multiphoton states” , i.e. molecules of light compound of $`n`$ photons with zero binding energy.

We know of course that the most satisfactory deduction of Planck’s formula is obtained nowadays by incorporating the right B.E. statistics into the corpuscular picture, as first shown by S.N.Bose  and interpreted and extended by Einstein in 1924/25 . Our calculation makes use of well-defined multiphoton states as the natural way of taking in account the identical, indistinguishable character of the photons.

The idea of multiphoton “molecules” arose around 1910, after P. Debye presented  a simple, ondulatory method to obtain the radiation law; some historical comments are included below in this letter.

2. We start from the photon entity as particles of energy $`ϵ=h\nu `$ (they are truly massless particles of helicity $`\pm 1`$ with energy and momentum $`ϵ=h\nu =pc,p`$, but we want to think and write in the spirit of the 1910s).

We make now this additional simple hypothesis : light of definite direction and frequency $`\nu `$ presents itself in units (“molecules”) of $`0,1,2,\mathrm{},n,\mathrm{}`$photons, with energy $`nh\nu `$, that is, with zero binding energy, and they contribute independently to the energy density.

The rationale for this assumption is, of course, in the spirit of the atomistic or democritean point of view, that all the photons of same frequency are created equal (identical), and that in an occamian, economical description, we should abstain from considering situations in which the photons are individually distinct (indistinguishability); that is, there is a “state” with two photons different, for the counting, from the naive juxtaposition of two distinguishable photons; same thing for the three photon states, etc. So the multiphoton molecules should be thought of as entities only to the effect of counting, and not as claiming for the existence of some force which really binds the photons together.

In this way one gets the correct radiation formula at once. We are unaware of any complete proof of this statement, so we now show the deduction, leaving some historical remarks for later.

For the $`n`$-photon “ molecules” of energy $`nh\nu ,(n=0,1,2,\mathrm{})`$ the probability is of course given by the Boltzmann factor

$$p(n)=c\mathrm{exp}(nh\nu /kT)$$
(3)
$$\underset{n=0}{\overset{\mathrm{}}{}}p(n)=1\text{gives}c=1x,x\mathrm{exp}(h\nu /kT).$$
(4)
The multiphoton states still have frequency $`\nu `$, so the common phase space factor is $`2`$(polarization)$`\times 4\pi `$(solid angle)$`\times \nu ^2`$ (radial 3D factor); hence the density of energy per inverse wavelength is

$`I(\nu ,T)`$ $`=`$ $`8\pi (\nu ^2/c^3)\mathrm{\Sigma }_0^{\mathrm{}}nh\nu (1x)\mathrm{exp}(nh\nu /kT)`$ (5)
$`=`$ $`{\displaystyle \frac{8\pi \nu ^3}{c^3}}{\displaystyle \frac{hx}{1x}}`$
$`=`$ $`{\displaystyle \frac{8\pi \nu ^3h}{c^3}}{\displaystyle \frac{1}{\mathrm{exp}(h\nu /kT)1}}`$
i.e., the correct complete radiation formula.

The two traditional limits of the radiation formula do have a particle interpretation, as already remarked by Wolfke . E.g. at low temperature only the lightest “molecule” is excited, that is, the one-photon mode, of energy $`h\nu `$, that leads to the Wien’s limit law (1), which was Einstein’s starting point .

On the other hand, for high densities or temperatures all the $`n`$-photon states contribute together and sort of coalesce, giving a scale-invariant power law distribution $`I\nu ^2`$, i.e. the classical electromagnetic (Rayleigh-Jeans)limit formula.

3. If we now start from (5) and expand the denominator, we get

$$I(\nu ,T)=8\pi (\nu ^2/c^3)(h\nu )\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{exp}(nh\nu /kT)$$
(6)
already obtained by several people  in the 1910s; superficially, (6) looks like an addditive contribution from each $`n`$-photon state, and it was this analogy which prompted Ehrenfest , Wolfke , de Broglie  and others  to advance the “$`n`$-photon molecule” concept; notice however that the literal interpretation of (6) as a cooperative multiphoton formula falls short of the correct derivation in (5) on three accounts:

1. There is no contribution from the zero-photon state
2. The $`n`$photon state contributes with energy $`h\nu `$, instead of the correct $`nh\nu `$
3. Probabilities $`\mathrm{exp}(nh\nu /kT)`$ do not add up to one.

Still, it is remarkable that the expansions (5) and (6) are really identical!

4. Why our simple derivation was not, to the best of our knowledge, obtained and exploited many decades ago, is for the historians to discuss; perhaps light quanta correlations due to indiscernibility was too novel at the time. Recall, however, that the Gibbs paradox was well known at the time, with the proposed solution by Gibbs himself . Notice that this treatment of photons is totally correct even from the modern point of view; in particular the zero-point energy does not show up. Also we are really doing second quantization. Namely “first quantization” just amounts to counting frequency (or wavenumber) modes, because $`\nu =pc/h`$; the “second” quantization amounts to the existence of these photon molecules as separate entities. We finish by remarking that zero-energy bound states appear frequently in modern guise as Bogomol’nyi-Prasad-Sommerfield (BPS) states, and moreover uncharged PBS states, like the photon, are massless.


text:
# Neutrinos from Gamma Ray Bursts

## Introduction

Gamma Ray Bursts (GRBs) have for decades persisted as a new and unexplained phenomenon in astrophysics, prompting numerous speculations about their nature and origin Fishman\_97 . After the discovery of X-ray and optical counterparts, the cosmological origin of GRBs is now beyond doubt, and the physics scenario employs the cataclysmic release of $`E10^{54}`$ ergs, one solar rest mass equivalent, under near vacuum conditions. Quasi-thermal neutrinos at $`20`$ MeV energies are the expected prime carriers during the first stages of energy release Kumar\_99 , and the difficulties in converting this energy into a high-entropy fireball during a neutron star merger have been pointed out Ruffert\_98 . The difficulties could be overcome by assuming a disk phase with Poynting flux extraction during coalescenceMeszaros\_97 , or a hypernova scenario Paczynski\_98 , and possibly by considering the interaction between the neutrino magnetic moment (which arises in any model explaining finite neutrino masses) and the magnetic field gradients in the collapsing object. Owing to the large distances of GRBs and the strong stellar background at the nuclear binding energy scale, there is little hope for detecting the initial neutrino emission component. There could be a very powerful neutrino emission component in the GeV range due to internal dissipation driven by proton-neutron diffusion Derishev\_99 ; Meszaros\_00 , and the background at these energies is significantly lower.

At high energies, neutrino detection becomes easier, since the weak interaction cross section, the muon range, and the angular resolution all increase with energy. As the consequence of a low mass-load (high entropy) in GRBs, the exploding matter can attain ultra-relativistic velocities, typically corresponding to Lorentz factors $`\mathrm{\Gamma }300`$, and can form relativistic shocks which are efficient in accelerating particles to even higher energies Meszaros\_92 . This is known to occur in the relativistic jets expelled from Active Galactic Nuclei (AGN) Rees\_67 , which transport some $`E10^{60}`$ ergs over their lifetime and have typical Lorentz factors of $`\mathrm{\Gamma }10`$. As a matter of fact, GRB and AGN jets involve almost identical concepts in their theoretical modeling.

It is thus not surprising that after the hypothesis of the association of GRBs and ultrahigh-energy cosmic rays (UHECRs) by Milgrom and Usov Milgrom\_95  an effort was launched to parallel the advanced theory of accelerated protons in AGN Protheroe\_83 ; Biermann\_87 ; Sikora\_87 ; Mannheim\_89 ; Stecker\_91 ; Mannheim\_93 ; Rachen\_93 ; Mannheim\_95 ; Mannheim\_98 ; Rachen\_98  for GRBs Waxman\_95 ; Vietri\_95 ; Waxman\_96 ; Waxman\_97 ; Vietri\_98 ; Totani\_98 . If the proton acceleration mechanism operates at the gyro time scale or faster, both AGN and GRBs represent possible classes of sources for UHECRs. The difference is that the AGN produce roughly two orders of magnitude more cumulative flux in the observed gamma ray band, and thus are likely much stronger sources of cosmic rays and neutrinos, too. Moreover, spatial coincidences between GRBs and UHECRs and the equality of local GRB and UHECR emissivities appear as a poor heuristic basis for the claim of a common origin in view of the fact that GRBs are typically too distant to allow for cosmic rays to reach Earth unimpeded by energy losses in interactions with the microwave background.

While this can be regarded as an argument disfavoring the hypothesis MPR\_00 ; Stecker\_00 ; Pugliese\_00 , there is one more aspect which deserves a closer look. This aspect is related with the shape of the multi-wavelength spectra of relativistic jet sources. The canonical spectral shape consists of two components, one peaking at lower energies, and one peaking at higher energies. The first component has a flux density that is independent on the photon energy $`E`$ (in the radio/mm range in AGN and in the X/$`\gamma `$-ray range in GRBs) which steepens by the factor $`E^1`$ at higher energies. This part of the spectrum is successfully modeled as the synchrotron radiation from accelerated electrons. If protons are accelerated in the sources, too, then there must be a second peak due to the hadronically-induced radiation at high energies in the canonical spectrum. Such a component has been found in the spectra of AGN jets in the MeV-to-TeV range, albeit it could also be interpreted as being due to inverse-Compton scattering of low-energy photons by the synchrotron-emitting electrons Mannheim\_98 . Naturalness would then provoke the assumption that GRB spectra have a corresponding peak at ultrahigh energies which has escaped detection so far.

This is consistent with the fact that due to the typical high redshifts of GRBs, such an emission component would be attenuated by pair production in traversing the intergalactic radiation fields Stecker\_98 ; Kneiske\_00 . We could be fooled by beaming statistics, if the dominant isotropic gamma ray output of GRBs went into a diffuse channel, e.g. due to cascading in the intergalactic medium wiping out directional memory and temporal correlation Totani\_98 .

The issue highlights the importance of neutrino observations. Neutrinos suffer practically no energy losses (except for adiabatic losses owing to the expansion of the Universe) and can therefore be used to measure the flux of cosmic rays escaping from GRBs. Prompt UHE protons escaping from relativistic outflows suffer from synchrotron and adiabatic losses Rachen\_98 , and therefore the neutron-escape fraction of cosmic rays from internal shocks should be the dominant one. Protons from external shocks would have a too steep spectrum to be important at UHE energies for reasonable energetics Kirk\_00 . The neutron-fraction of the escaping cosmic ray flux can be inferred from measurements of neutrinos in the 100 TeV range. These neutrinos would be correlated with the electromagnetic bursts. Cosmic rays also give rise to an uncorrelated secondary neutrino flux during propagation through the microwave background due to photo-production of pions.

The paper first provides a general bound on the neutrino emission from GRB-like sources, constrained by the observed diffuse gamma-ray and cosmic ray proton fluxes. Secondly, more specific assumptions regarding the spectral shape of the emitted protons are employed to find out the expected flux of high-energy neutrinos following the Milgrom & Usov hypothesis. Finally, this would have consequences for the properties of high energy gamma-ray emission from GRB which are briefly discussed.

## Upper limit for GRB-like photo-hadronic sources of neutrinos

Recently, Mannheim, Protheroe, and Rachen MPR\_00  (MPR) have derived an upper limit for the neutrino flux from extragalactic sources, in the framework of photo-production models in which the neutrinos originate from pions produced in hadronically-induced interactions of accelerated protons with low-energy photons. With every neutrino, there is a minimum number of photons and neutrons associated, and the flux of those is not allowed to exceed observed omnidirectional fluxes of cosmic ray protons and gamma rays. In their Fig. 3, MPR show this limit for sources with conditions typical for AGN (steady sources with flux density spectral index $`\alpha =1`$). Here, the focus is on the short-lived GRBs in which the synchrotron target photons have a flat spectrum ($`\alpha =0`$). The relation between cosmic ray and neutrino fluxes at ultrahigh-energies (above $`10^7`$ GeV proton energy in the comoving frame) is then given by

$$Q_{\nu _\mu }(E)=416Q_{\mathrm{cr}}(25E)$$
(1)
i.e. the neutrino yield can be up to factor of 5 larger than in the case of AGN Muecke\_99 .

The experimental upper limit on a possible extragalactic proton contribution to the local flux of cosmic rays is consistent with an extrapolation of the observed, presumably protonic, spectrum at $`10^{10}`$ GeV down to $`10^6`$ GeV using a spectral index of 2.75. Owing to the steepness of this spectrum, saturating the upper limit implies an increasing emissivity of the putative GRB-like sources with decreasing energy. The electromagnetic and neutrino luminosities of GRB-like sources are the same. Below PeV energies, their electromagnetic output begins to reach the level of the observed extragalactic gamma ray background and the upper limit levels off. Generally, since gamma-rays from pion decays and from pair production are reprocessed by electromagnetic cascades, only the bolometric (integrated) fluxes can be compared.

Equation (1) still ignores the effect of cosmic ray attenuation, which is important for sources with strong redshift evolution. A strong redshift dependence $`(1+z)^{3.4}`$ of the cumulative source emissivity is assumed for the GRBs. This corresponds to the observed dependence for star formation in galaxies and AGN Boyle\_98 . Above $`10^{10}`$ GeV, the upper limit spectrum steepens rapidly due to the inescapable energy losses of pions and muons in compact, variable sources Rachen\_98 . Below PeV energies, the increase is only by a factor of 2, since the constraint from the diffuse gamma-ray background imposed by $`L_\gamma L_\nu `$ (where $`L_i`$ denotes the bolometric luminosity emitted in species $`i`$) becomes the stronger one, and AGN-like conditions have $`L_\gamma =2L_\nu `$ (due to the contribution of gamma-rays from Bethe-Heitler pairs). For realistic model assumptions, the turnover at $`10^{10}`$ GeV would occur at even lower energies (see next section). Figure 1 also shows previously predicted fluxes for GRBs based on a simplified treatment of the photo-production kinematics Waxman\_97  and for AGN Mannheim\_95  for comparison.

## Generic GRB neutrino flux

For the calculation of an upper limit for GRB-like sources, explicit model assumptions required for a generic GRB model have not been made. A crucial assumption is that of the slope of the energy distribution of accelerated protons.

The external shock where the ejecta impact the external medium is highly relativistic and should thus produce particle spectra $`E^{2.2}`$ Kirk\_00 . If protons accelerated at the external shock were to explain the UHE cosmic rays, the power in the low-energy part of the population would have to be unrealistically high, and the GRBs radiatively extremely inefficient.

Internal GRB shocks are expected to be only mildly relativistic owing to the velocity differences of superseding shells being smaller than those with the external medium, and therefore a proton spectrum $`E^2`$ seems possible. Such an energy distribution delivers a large fraction of the energy in the proton population at UHE energies. The shape of the neutrino spectrum follows $`E^1`$ below $`E_\mathrm{b}=100(\mathrm{\Gamma }/300)^2`$ TeV and $`E^2`$ above, until the spectrum steepens above 10 PeV due to energy losses of the muons in the expanding ejecta. The actual spectrum will be slightly more curved (see the lower thick solid curve in Fig. 1), since the photo-hadronic interactions are dominated by the $`\mathrm{\Delta }`$-resonance in the 100 TeV range, whereas by multiple pion production above PeV lab frame energies. A more accurate treatment using the SOPHIA code Muecke\_99  must await future work.

Since the protons must suffer synchrotron and adiabatic losses upon escape Rachen\_98 , the neutrons produced in photon collisions dominate the escaping fraction of cosmic ray baryons at ultrahigh energies. The generic GRB neutrino flux following from the assumption that GRBs produce the UHE cosmic rays above $`10^9`$ GeV (exempting the super-Greisen events), is shown in Fig. 1. The corresponding integral gamma-ray flux is given by $`F_\gamma =\frac{3}{2}(2\times 10^7\times \mathrm{ln}[100])`$ GeV cm-2 s-1 sr-1 which corresponds to 10% of the extragalactic gamma ray background between 3 MeV and 30 GeV ($`1.5\times 10^5`$ in the same units). Comparing with the extragalactic gamma ray background between 1 GeV and 30 GeV ($`4.2\times 10^6`$), the number is 33%. The factor $`\mathrm{ln}[100]`$ takes into account the bandwidth of the neutrino spectrum (the limit for AGN-like sources shown in Fig. 1 takes into account a factor $`\mathrm{exp}[1]`$ as the bolometric correction for trial spectra of the form $`E^1\mathrm{exp}[E/E_{\mathrm{max}}]`$).

Note, however, that the secondaries produced during propagation of UHE cosmic rays in interactions with the microwave background add another (albeit uncorrelated, and therefore truly diffuse) neutrino component at EeV energies (similar to the dotted curve in Fig. 1) and also add an electromagnetic component emerging in the observed gamma-ray band. This component has an energy flux level comparable to the minimum of the upper limit for AGN-like sources at $`10^9`$ GeV, since photo-production occurs close to the threshold energy in interactions with photons from the microwave background. The spectrum is rather sharply peaked at UHE energies and has a bolometric correction of $`\mathrm{exp}[1]`$. The corresponding gamma-ray flux is $`2\times 10^7`$ GeV cm-2 s-1 sr-1, so that the combined GRB contribution to the extragalactic gamma ray background would be $`1.6\times 10^6`$ in the same units (11% of the 3 MeV-30 GeV background or 38% of the 1 GeV-30 GeV background).

Constraints on the total isotropic energy per burst, their radiative efficiency, and the burst rate have not been imposed here. It is clear, however, that if GRBs shall explain the UHE cosmic rays, the current physical interpretation of GRBs is severely challenged by the dramatically increased energy demands.

## Implications for high-energy gamma-ray tails

The powerful high-energy gamma-ray emission implied by the Milgrom & Usov hypothesis raises the possibility for falsification by gamma-ray observations. However, the current generation of ground-based telescopes are sensitive only above 300 GeV corresponding to a pair-attenuation length of $`z=0.3`$. Only a small number of bursts is expected to occur at redshifts less than 0.3. From an analysis of the 2nd BATSE catalogue data, the expected rate of 100 per yearMannheim\_96  is too low for the small field of view of imaging air Cherenkov telescopes to make a chance detection. The large-field-of-view detectors MILAGRITO, Tibet, and HEGRA-AIROBICC have not found strong evidence for GRBs, but their thresholds are in the TeV domain where less than one GRB per year is expected to occur. It is interesting to note, however, that the groups have all reported marginal evidence for single GRBs which is statistically conceivable Atkins\_00 ; Padilla\_98 ; Amenomori\_98 . The implied flux is consistent with the Milgrom & Usov hypothesis Totani\_99 . This high-energy tail would have an extremely hard spectrum rising above 10 GeV, since EGRET has not discovered evidence for an upturn. Clearly, this is a challenging possibility to be tested by AGILE, GLAST, or MAGIC.

## Summary and conclusions

Non-thermal power law emission and afterglows over a broad range of wavelengths indicate the acceleration of charged particles to relativistic energies in GRBs. If GRBs are the sources of UHE cosmic rays, they must produce at least 10% of the extragalactic gamma ray background between 3 MeV and 30 GeV and a diffuse flux of neutrinos in the 100 TeV to 10 PeV regime which is a factor of 20 below the AMANDA-$`\nu `$2000 limit on isotropic neutrinos. The time stamp associated with GRBs increases the sensitivity further compared with the sensitivity for an isotropic, uncorrelated background. Low-redshift GRBs would also have powerful high-energy tails in their gamma ray spectra which are consistent with the marginal detections reported by the MILAGRO, HEGRA, and Tibet Collaborations. The next years will bring rapid progress to the understanding of GRBs and their high-energy emissions, with improved measurements of GRB neutrinos, redshifts, and high-energy gamma-ray tails using AMANDA-ICECUBE/ANTARES/NESTOR, HETE-2 follow-up programs, and MAGIC/AGILE/GLAST.

###### Acknowledgements.

I am deeply indebted to Jörg Rachen and Raymond Protheroe for their invaluable contributions to this research. Superior support through a Heisenberg-Fellowship granted by the DFG is gratefully acknowledged.


text:
# Optimization with Extremal Dynamics

## Abstract

We explore a new general-purpose heuristic for finding high-quality solutions to hard optimization problems. The method, called extremal optimization, is inspired by self-organized criticality, a concept introduced to describe emergent complexity in physical systems. Extremal optimization successively replaces extremely undesirable variables of a single sub-optimal solution with new, random ones. Large fluctuations ensue, that efficiently explore many local optima. With only one adjustable parameter, the heuristic’s performance has proven competitive with more elaborate methods, especially near phase transitions which are believed to coincide with the hardest instances. We use extremal optimization to elucidate the phase transition in the 3-coloring problem, and we provide independent confirmation of previously reported extrapolations for the ground-state energy of $`\pm J`$ spin glasses in $`d=3`$ and $`4`$. PACS number(s): 02.60.Pn, 05.65.+b, 75.10.Nr, 64.60.Cn.

Many natural systems have, without any centralized organizing facility, developed into complex structures that optimize their use of resources in sophisticated ways Bakbook . Biological evolution has formed efficient and strongly interdependent networks in which resources rarely go to waste. Even the morphology of inanimate landscapes exhibits patterns that seem to serve a purpose, such as the efficient drainage of water Somfai ; rivers .

Natural systems that exhibit self-organizing qualities often possess a common feature: a large number of strongly coupled entities with similar properties. Hence, at some coarse level they permit a statistical description. An external resource (sunlight in the case of evolution) drives the system which then takes its direction purely by chance. Like descending water breaking through the weakest of all barriers in its wake, biological species are coupled in a global comparative process that persistently washes away the least fit. In this process, unlikely but highly adapted structures surface inadvertently. Optimal adaptation thus emerges naturally, from the dynamics, simply through a selection against the extremely “bad”. In fact, this process prevents the inflexibility inevitable in a controlled breeding of the “good”.

Various models relying on extremal processes have been proposed to explain the phenomenon of self-organization PMB . In particular, the Bak-Sneppen model of biological evolution is based on this principle BS ; BoPa1 . Assuming an unspecified interdependency between species, it produces salient nontrivial features of pale-ontological data such as broadly distributed lifetimes of species, large extinction events, and punctuated equilibrium.

In the Bak-Sneppen model, species are located on the sites of a lattice, and have an associated “fitness” value between 0 and 1. At each time step, the one species with the smallest value (poorest degree of adaptation) is selected for a random update, having its fitness replaced by a new value drawn randomly from a flat distribution on the interval $`[0,1]`$. But the change in fitness of one species impacts the fitness of interrelated species. Therefore, all of the species at neighboring lattice sites have their fitness replaced with new random numbers as well. After a sufficient number of steps, the system reaches a highly correlated state known as self-organized criticality (SOC) BTW . In that state, almost all species have reached a fitness above a certain threshold. These species, however, possess punctuated equilibrium: only one’s weakened neighbor can undermine one’s own fitness. This coevolutionary activity gives rise to chain reactions called “avalanches”, large fluctuations that rearrange major parts of the system, potentially making any configuration accessible.

Although coevolution does not have optimization as its exclusive goal, it serves as a powerful paradigm. We have used it as motivation for a new approach to optimization BoPe1 . The heuristic we have introduced, called extremal optimization (EO), follows the spirit of the Bak-Sneppen model, merely updating those variables having the “worst” values in a solution and replacing them by random values without ever explicitly improving them.

Previously, several heuristics inspired by natural processes have been proposed, notably simulated annealing Science  and genetic algorithms Holland . These have aroused considerable interest among physicists, in the development of physically motivated heuristics Houdayer ; Dittes ; Moebius , and in their applications to physical problems Pal ; Hartmann\_d3 ; Hartmann\_d4 . EO adds a distinctly new approach to the approximation of hard optimization problems by utilizing large fluctuations inherent to systems driven far from equilibrium. In this Letter we demonstrate EO’s generality, the simplicity of its implementation, and its results, for a physical problem as well as for a classic combinatorial optimization problem. A tutorial introduction is given in Ref. CISE .

The physical problem we consider is a spin glass MPV . It consists of a $`d`$-dimensional hyper-cubic lattice of length $`L`$, with a spin variable $`x_i\{1,1\}`$ at each site $`i`$, $`1in`$ $`(=L^d)`$. A spin is connected to each of its nearest neighbors $`j`$ via a bond variable $`J_{ij}\{1,1\}`$, assigned at random. The configuration space $`\mathrm{\Omega }`$ consist of all configurations $`S=(x_1,\mathrm{},x_n)\mathrm{\Omega }`$ where $`|\mathrm{\Omega }|=2^n`$.

We wish to minimize the cost function, or Hamiltonian

$`C(S)=H(𝐱)={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}J_{ij}x_ix_j}.`$ (1)
Due to frustration MPV , ground state configurations $`S_{\mathrm{min}}`$ are hard to find, and it has been shown that for $`d>2`$ the problem is among the hardest optimization problems Barahona .

To find low-energy configurations, EO assigns a fitness to each spin variable $`x_i`$

$`\lambda _i=x_i\left({\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}J_{ij}x_j\right),`$ (2)
so that

$`C(S)={\displaystyle \underset{i=1}{\overset{n}{}}}\lambda _i.`$ (3)
Unlike in genetic algorithms, the fitness does not characterize an entire configuration, but simply a single variable of the configuration. In similarity to the Bak-Sneppen model, EO then proceeds to search $`\mathrm{\Omega }`$ by sequentially changing variables with “bad” fitness at each update step. The simplest “neighborhood” $`N(S)`$ for an update consists of all configurations $`S^{}N(S)`$ that could be reached from $`S`$ through the flip of a single spin. After each update, the fitnesses of the changed variable and of all its neighbors are reevaluated according to Eq. (2).

The basic EO algorithm proceeds as follows:

Initialize configuration $`S`$ at will; set $`S_{\mathrm{best}}:=S`$.    For the “current” configuration $`S`$,    evaluate $`\lambda _i`$ for each variable $`x_i`$,    find $`j`$ satisfying $`\lambda _j\lambda _i`$ for all $`i`$, i.e., $`x_j`$ has the “worst fitness”,    choose $`S^{}N(S)`$ so that $`x_j`$ must change,    accept $`S:=S^{}`$ unconditionally,    if $`C(S)<C(S_{\mathrm{best}})`$ then set $`S_{\mathrm{best}}:=S`$.      Repeat at step (2) as long as desired.    Return $`S_{\mathrm{best}}`$ and $`C(S_{\mathrm{best}})`$.

Initial tests have shown that this basic algorithm is quite competitive for optimization problems where EO can choose randomly among many $`S^{}N(S)`$ satisfying step (2c), such as graph partitioning BoPe1 . But, in cases such as the single spin-flip neighborhood above, focusing on only the worst fitness \[step (2b)\] leads to a deterministic process, leaving no choice in step (2c). To avoid such “dead ends” and to improve results BoPe1 , we introduce a single parameter into the algorithm. We rank all the variables $`x_i`$ according to fitness $`\lambda _i`$, i.e., find a permutation $`\mathrm{\Pi }`$ of the variable labels $`i`$ with

$`\lambda _{\mathrm{\Pi }(1)}\lambda _{\mathrm{\Pi }(2)}\mathrm{}\lambda _{\mathrm{\Pi }(n)}.`$ (4)
The worst variable $`x_j`$ \[step (2b)\] is of rank 1, $`j=\mathrm{\Pi }(1)`$, and the best variable is of rank $`n`$. Now, consider a probability distribution over the ranks $`k`$,

$`P_kk^\tau ,1kn,`$ (5)
for a given value of the parameter $`\tau `$. At each update, select a rank $`k`$ according to $`P_k`$. Then, modify step (2c) so that the variable $`x_i`$ with $`i=\mathrm{\Pi }(k)`$ changes its state.

For $`\tau =0`$, this “$`\tau `$-EO” algorithm is simply a random walk through $`\mathrm{\Omega }`$. Conversely, for $`\tau \mathrm{}`$, the process approaches a deterministic local search, only updating the lowest-ranked variable, and is bound to reach a dead end (see Fig. 1). However, for finite values of $`\tau `$ the choice of a scale-free distribution for $`P_k`$ in Eq. (5) ensures that no rank gets excluded from further evolution, while maintaining a clear bias against variables with bad fitness.

For the spin glass, we obtained our best solutions for $`\tau 1^+`$, as shown in Fig. 1. Generally, over many optimization problems, the preferred value seems to scale slowly as $`\tau 11/\mathrm{ln}n`$ for increasing $`n`$. Experiments with random graphs support this expectation BoPe2 , although the dependence on $`n`$ can be practically negligible. For instance, at fixed $`\tau =1.4`$, $`\tau `$-EO already reproduced many testbed results for the partitioning of graphs of sizes $`n=1000`$ to $`n10^5`$ BoPe1 .

We have run the $`\tau `$-EO algorithm with $`\tau =1.15`$ on a large number of realizations of the $`J_{ij}`$, for $`n=L^d`$ with $`L=5,6,7,8,9,10,12`$ in $`d=3`$, and with $`L=3,4,5,6,7`$ in $`d=4`$. To reduce variances, we fixed $`|J_{ij}|1`$. For each instance, we have run EO with 5 restarts from random initial conditions, retaining only the lowest energy state obtained, and then averaging over instances. Inspection of the convergence results for the genetic algorithms in Refs. Pal ; Houdayer  suggest a runtime scaling at least as $`n^3`$$`n^4`$ for consistent performance. Indeed, using $`n^4/100`$ updates enables EO to reproduce its lowest energy states on about 80% to 95% of the restarts, for each $`n`$. Our results are listed in Table 1. A fit of our data with $`e_d(n)=e_d(\mathrm{})+A/n`$ for $`n\mathrm{}`$ predicts $`e_3(\mathrm{})=1.7865(3)`$ for $`d=3`$ and $`e_4(\mathrm{})=2.093(1)`$ for $`d=4`$. Both values are consistent with the findings of Refs. Pal ; Hartmann\_d3 ; Hartmann\_d4 , providing independent confirmation of those results with far less parameter tuning.

In the future, we intend to use EO to explore the properties of states near the ground state. As the results on 3-coloring in Ref. BGIP  suggest, EO’s continued fluctuations through near-optimal configurations even after first reaching near-optimal states may provide an efficient means of exploring a configuration space widely.

To demonstrate the versatility of EO (see also Ref. PPSN ), we now turn to a popular combinatorial optimization problem, graph coloring G+J . Specifically, we consider random graph 3-coloring AI ; Culberson . A random graph is generated by connecting any pair of its $`n`$ vertices by an edge, with probability $`p`$ Bollobas . In $`K`$-coloring, given $`K`$ different colors used to label the vertices of a graph, we need to find a coloring that minimizes the number of “monochromatic” edges connecting vertices of identical color. We implement EO by defining $`\lambda _i`$ for each vertex to be $`1/2`$ times the number of monochromatic edges attached to it. Then, Eq. (3) represents exactly the cost function, counting the number of monochromatic edges present. As a simple neighborhood definition, at each update we merely change the color of one “bad” vertex selected according to $`\tau `$-EO \[step (2c)\].

As a special challenge, we have used $`\tau `$-EO for 3-coloring to investigate the phase transition that occurs under variation of the average vertex degree $`\alpha =pn`$ Cheeseman ; Culberson . Random graphs with small $`\alpha `$ can almost always be colored at zero cost, while graphs with large $`\alpha `$ are typically very homogeneous with a high but easily approximated cost. Located at some point $`\alpha _{\mathrm{crit}}`$ between these extremes, there is a sharp phase transition to non-zero cost solutions. Such a critical point appears in many combinatorial optimization problems, and has been conjectured to harbor those instances that are the hardest to solve computationally Cheeseman . Previously we have shown EOperc  that $`\tau `$-EO significantly outperforms simulated annealing near the phase transition of the bipartitioning problem of random graphs.

Using EO we can estimate the value of $`\alpha _{\mathrm{crit}}`$ for 3-coloring. To this end, we have averaged the cost EO obtains as a function of the vertex degree $`\alpha `$. We generated $`10000`$, $`5000`$, $`1300`$, and $`650`$ instances for $`n=32`$, 64, 128, and 256, respectively, for values of $`3.6\alpha 6`$. Since $`n`$ is relatively small and the runs were chosen to be very long ($`100n^2`$ updates), we found optimal performance at $`\tau =2.7`$. Such excessively long runs were used as part of a study to find all minimal-cost solutions for each instance, in order to determine their overlap (or “backbone”). Elsewhere BGIP  we show that this backbone appears to undergo a first-order phase transition as conjectured in Ref. Monasson .

Finite size scaling with the ansatz

$`C(\alpha ,n)f\left[\left(\alpha \alpha _{\mathrm{crit}}\right)n^{1/\nu }\right]`$ (6)
applied to the results depicted in Fig. 2 predicts $`\alpha _{\mathrm{crit}}4.72(1)`$ and $`\nu 1.53(5)`$. These are the most precise estimates to date. The numerical value of $`\nu `$ suggests that in fact $`\nu =3/2`$, just as for the percolation transition of random graphs Bollobas  and the transition in 3-satisfiability KS .

In conclusion, we have presented a new optimization method, called extremal optimization due to its derivation from extremally driven statistical systems. At each update step, the algorithm assigns fitnesses to variables $`x_i`$, and then generates moves by randomly updating an “unfit” variable. EO gives no consideration to the move’s outcome. Large fluctuations in the cost function can accumulate over many updates; only the bias against poor fitnesses guides EO back towards improved solutions.

A drawback to EO is that a general definition of fitness for individual variables may prove ambiguous or even impossible. Also, each variable $`x_i`$ could have a large number of states to choose from, as in $`K`$-coloring with large $`K`$; random updates would then be more likely to remove than to create well-adapted variables (this is notably the case for the traveling salesman problem BoPe1 ). And in highly connected systems (e. g. for $`\alpha 1`$ in $`K`$-coloring), EO may be slowed down considerably by continual fitness recalculations \[step (2a)\].

However, extremal optimization is readily applicable to problems whose cost can be decomposed into contributions from individual degrees of freedom. It is easily implemented and, using very few parameters, it can prove highly competitive. We have shown this on the spin glass Hamiltonian, obtaining $`d=3`$ and $`4`$ ground-state energies that are consistent with the best known results. We have also used EO to explore the phase transition in random graph 3-coloring. Its results enable us to provide, by way of finite size scaling, the first sound estimates of critical values for this problem.

###### Acknowledgements.

We would like to thank the participants of the 1999 Telluride Workshop on Energy Landscapes for valuable input. Special thanks to Paolo Sibani, Jesper Dall, Sigismund Kobe, and Gabriel Istrate. This work was supported in part by the URC at Emory University and by an LDRD grant from Los Alamos National Laboratory.