text:
# Near-infrared [Fe II] emission from supernova remnants and the supernova rate of starburst galaxies

## 1 Introduction

Early near-infrared (IR) spectroscopic observations of galactic SNRs have shown that these objects present remarkably strong \[Fe II\] features (e.g. Graham, Wright & Longmore 1987). This has been proposed to be mainly due to the existence of an extended postshock region in which the ionization conditions are such that Fe+ can be efficiently excited by electron collisions (Mouri, Kawara & Taniguchi 2000). In contrast, the spatial extent of this zone, where hydrogen is partially ionized, is small in photoionized regions, resulting in little emission (e.g. Luhman, Engelbracht & Luhman 1998).

This peculiarity can be used as a diagnostic probe of supernova activity in distant starburst galaxies where the SNRs are generally unresolved and are thought to produce much of the integrated near-IR \[Fe II\] emission (e.g. Vanzi & Rieke 1997). Of particular interest is the fact that using near-IR \[Fe II\] lines to estimate supernova activity can potentially be superior to the usual radio technique, in that it more clearly distinguishes H II regions from SNRs.

From estimates of the total \[Fe II\] luminosity of a starburst galaxy ($`_{[\mathrm{Fe}\mathrm{II}]}`$) and assumptions regarding the typical values of the \[Fe II\] luminosity ($`L_{[\mathrm{Fe}\mathrm{II}]}`$) and \[Fe II\]-emitting lifetime ($`t_{[\mathrm{Fe}\mathrm{II}]}`$) of a single SNR, one can define the supernova rate, $`\eta `$, as:

$$\eta =\frac{_{[\mathrm{Fe}\mathrm{II}]}}{t_{[\mathrm{Fe}\mathrm{II}]}L_{[\mathrm{Fe}\mathrm{II}]}}\mathrm{yr}^1$$
(1)
This quantity can be used, for instance, to set constraints on the star formation history of starburst galaxies (e.g. Leitherer & Heckman 1995; Kotilainen et al. 1996).

Unfortunately, the supernova rate derived by this technique is considerably uncertain, as it scales with $`t_{[\mathrm{Fe}\mathrm{II}]}`$ and $`L_{[\mathrm{Fe}\mathrm{II}]}`$ which are poorly known. The \[Fe II\]-emitting lifetime of a SNR is generally taken to be of the order of 1–2 $`\times `$ 104 yr, but this estimate lies on weak grounds (e.g. van der Werf et al. 1993). On the other hand, the few \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosities quoted in the literature for SNRs in our Galaxy, the LMC or M33 span a wide range of values (from 0.3 to 720 L; Graham et al. 1987; Oliva, Moorwood & Danziger 1989, 1990; Lumsden & Puxley 1995 – hereafter LP). One issue of particular importance is to understand why SNRs in starburst galaxies exhibit dramatically higher \[Fe II\] luminosities (up to 1.6 $`\times `$ 105 L in M82; Greenhouse et al. 1997). One order of magnitude higher values are also observed in NGC 253 (Forbes et al. 1993). Before using near-IR \[Fe II\] features to estimate supernova activity in star-forming galaxies that presumably exhibit widely different properties (e.g. metallicity), one must first understand what parameters control the \[Fe II\] properties of individual SNRs. Many factors (e.g. evolutionary status of the SNR, density or metallicity of the local ISM) are susceptible to play a role in producing this observed luminosity range, although their relative importance has yet to be established.

## 2 A near-IR \[Fe II\] line-imaging survey of extragalactic SNRs

In an effort to address this issue, we present the results of a project in which narrowband images are used to secure \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosities for a large sample of extragalactic SNRs spanning a wide range of physical properties and evolutionary status. Previous efforts to determine the \[Fe II\] luminosities of individual SNRs have mainly concentrated on aperture photometry of galactic remnants. The shortcomings of such an approach are numerous: (i) those observations sample, in most cases, only a very small part of the SNR, rendering the total value of the \[Fe II\] luminosity considerably uncertain; (ii) because of the different apertures used, the published \[Fe II\] line luminosities do not constitute a homogeneous database based on which meaningful statistical studies can be conducted; (iii) the \[Fe II\] luminosities are entached of large uncertainties, as the distance to these SNRs and the extinction corrections are generally poorly known. On the other hand, because of the unavoidable limited size of the aperture, surveys in nearby galaxies carried out using long-slit spectroscopy exclude the largest objects and are therefore biased towards the youngest SNRs (see LP). Line-imaging observations of SNR populations in local group galaxies do not suffer from these limitations.

M33 is an obvious target in this respect because of its relatively face-on orientation ($`i`$ $``$ 55; Garcia-Gomez & Athanassoula 1991), well-known distance ($`D`$ $``$ 840 kpc; Freedman, Wilson & Madore 1991) and substantial abundance gradient that allows an investigation of metallicity effects (Smith et al. 1993). Furthermore, it has the largest catalogued extragalactic SNR population (Gordon et al. 1999). Of particular importance is the fact that a large body of observational data exists for SNRs in this galaxy; as will be shown below, this will allow us to relate the near-IR \[Fe II\] emission to observations in other parts of the electromagnetic spectrum. The latest catalogue of optically-selected SNR candidates in M33 comprises 98 objects (Gordon et al. 1998); many of which being detected at radio frequencies (Gordon et al. 1999, and references therein). Among this sample, 72 objects have been spectroscopically confirmed and therefore also have well-defined optical-line properties (Gordon et al. 1998, and references therein). In addition, a ROSAT survey of M33 has associated some X-ray sources with known SNRs (Long et al. 1996).

## 3 Observations and Reduction Procedure

Our observations were obtained at the Canada-France-Hawaii Telescope (CFHT) in 1997 and 1998, using the IR cameras MONICA and REDEYE, respectively. For both runs, the cameras were equipped with a 256 $`\times `$ 256 pixels HgCdTe NICMOS-3 detector array. At the f/8 Cassegrain focus of the CFHT, the image scales are 0.246 (MONICA) and 0.5 arcsec pixel-1 (REDEYE). For the assumed distance to M33 ($`D`$ $``$ 840 kpc), this translates into a physical scale of 1 and 2 pc pixel-1, respectively. In this configuration, the fields of view are about 63 (MONICA) and 128 arcsec (REDEYE). The interference filters used are a narrowband filter (FWHM $``$ 0.02 $`\mu `$m) centred on \[Fe II\] $`\lambda `$1.644 $`\mu `$m (a4D7/2 $``$ a4F9/2) and a wideband H filter. The \[Fe II\] filter is sufficiently wide that the broadening due to motions of the remnant gas, the rotation or the redshift of M33 do not shift lines to significantly different parts of the transmission curve. In total 42 objects drawn from the optically-selected catalogue of Gordon et al. (1998) have been observed. The vast majority has been spectroscopically confirmed (Blair & Kirshner 1985; Smith et al. 1993; Gordon et al. 1998). Particular attention has been paid to observe SNRs with widely different properties. The journal of observations is presented in Table 1.

For each SNR we have obtained 5 images: one centred on the object and 4 offset to the NO, NE, SO and SE by 22 and 45 arcsec for MONICA and REDEYE, respectively. These images were later combined in a mosaic centred on the SNR. Two type of domeflats have been obtained (dome lights “on” and “off”). The latter has been subtracted from the former in order to remove the contribution from the thermal background. These domeflats have been used to correct for changes in the quantum efficiency of the pixels in each object image. A mean sky frame has been created by applying a pixel-to-pixel median filtering procedure with proper masking of sources to the data frames themselves. No off-source images have been obtained. This procedure is justified in our case, as we are not attempting to detect uniform large-scale emission, but rather compact sources on the arcsecond scale (see Fig.1). After subtraction of this mean sky frame to the individual object frames and multiplication by a bad pixel mask in order to account for dead pixels, the object images were aligned and coadded. The astrometric calibration was carried out with reference to positions of stars in the USNO-A2.0 or GSC1.2 catalogues. The positions are accurate to within 3 arcsec.

## 4 Overview of the survey

Among our sample of 42 objects, only 7 have been firmly detected (we address in Section 6.1 the physical reasons for such a low detection rate). Fig. 1 shows the \[Fe II\] and H-band images of this subset. The \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosities are quoted in Table 2, along with the 3$`\sigma `$ upper limits for the remaining SNRs (col. ). The targets were observed at similar airmasses as the standard stars, eliminating the need to correct for atmospheric extinction. This is legitimate in the H band, where such corrections are small. Using the optical extinction measurements of Blair & Kirshner (1985) and the interstellar reddening law of Rieke & Lebofsky (1985), it is also found that reddening corrections are negligible. The contribution of Br12 $`\lambda `$1.641 $`\mu `$m and \[Si I\] $`\lambda `$1.645 $`\mu `$m to the total flux in the \[Fe II\] filter has been neglected (see Oliva et al. 1989). The non-detection of any significant flux in the H-band images also indicates that the contribution from continuum emission is negligible (see also LP). The flux calibration was achieved by using faint standard stars in the list of Hunt et al. (1998). The internal accuracy is typically 0.011 mag in the H band.

Some of our objects have been observed in near-IR spectroscopy by LP. We have adopted their values whenever appropriate (i.e. when the object was detected in their survey but not in ours or when they set a more stringent upper limit). Four objects have been detected in both surveys and can be used to assess the robustness of our luminosity calibration (this is in particular necessary when considering that some of our observations were acquired while some cirrus clouds were present; see Table 1). In all cases, both values are found to agree within the uncertainties. The internal consistency of our data can also be examined by comparing repeated observations of SNR28 taken under different photometric conditions, as detailed in Table 1. Once again, these quantities are found to be statistically undistinguishable. These checks confirm the reliability of our derived luminosities.

The contour maps of the detected SNRs are shown in Fig. 2. Although the limited spatial resolution of our images prevents us to examine this aspect in detail, a comparison with Hubble Space Telescope (HST) or ground-based optical images (Blair & Davidsen 1993; Gordon et al. 1998) reveals roughly similar \[Fe II\]- and optical-line morphologies. We note significant differences between the diameters determined from our \[Fe II\] images and from optical observations (Gordon et al. 1998). This might result from difficulties in subjectively assigning diameters for objects with complex morphologies (see Blair & Davidsen 1993).

## 5 Correlations between the \[Fe II\] emission and other SNR properties

In order to pin down the physical processes that are responsible for the \[Fe II\] emission, we sought for correlations between the \[Fe II\] luminosities and various quantities intrinsic to the SNRs (e.g. optical-line properties). These data were gathered from the literature and are quoted in Table 2.

A number of statistical methods have been developed to investigate relationships between astronomical data sets with censored data (e.g. with upper limits). We made use here of the generalized Kendall’s $`\tau `$ correlation technique (Isobe, Feigelson & Nelson 1986). Table 3 summarizes the results of this statistical analysis for all quantities under consideration. When a significant correlation is present (we consider this to be the case when the false alarm probability is less than 1 per cent), we give the results of a linear regression fit to the data for the detected objects (such regression techniques for censored data yield unreliable results in the presence of heavy censoring). We also show in Fig. 3 the dependence of $`L_{[\mathrm{Fe}\mathrm{II}]}`$ on the quantity in question. Survival analysis techniques yield more reliable results when the censored points are randomly distributed, a condition that is generally not fulfilled here (see Fig.3). In order to assess the robustness of our results, we performed the same analysis when only considering firm detections. This procedure broadly confirms the correlations previously found (albeit with a lower level of significance), and suggests that the results inferred by taking into account all the available information (i.e. detections and upper limits) can be regarded as robust. We discuss below the implications of this statistical analysis on our understanding of \[Fe II\] emission in SNRs.

### 5.1 Optical-line properties

A significant correlation is found between all optical-line luminosities and $`L_{[\mathrm{Fe}\mathrm{II}]}`$ (Table 3). A good correlation between the \[Fe II\] and H$`\beta `$ surface brightnesses has already been reported for a sample of LMC and Galactic SNRs by Oliva et al. (1989). Although the optical and \[Fe II\] line-emitting regions are likely to differ in detail, this correlation suggests that the sites of \[Fe II\] emission are bright optical filaments whose emission is excited by radiative shocks. Support for this interpretation comes from observations of RCW 103 and N49 which show the \[Fe II\]- and optical-line morphologies to be virtually identical (Burton & Spyromilio 1993; Dickel et al. 1995). In both cases, the optical line-emitting regions are believed to trace the interaction of a radiative blast wave with a neighbouring molecular cloud (e.g. Banas et al. 1997; Oliva et al. 1999).

### 5.2 Density effects

A correlation is also found between the electronic density of the postshock gas (as derived from the optical \[S II\] doublet) and $`L_{[\mathrm{Fe}\mathrm{II}]}`$. Near-IR \[Fe II\] features in SNRs are collisionally excited via electron impact in the dense recombination zone behind the shock front. The critical electron number density of \[Fe II\] $`\lambda `$1.644 $`\mu `$m is (Blietz et al. 1994):

$$n_c(T_4)4.2\times 10^4T_4^{0.69}\mathrm{cm}^3$$
(2)
In this expression, $`T_4`$ is the electronic temperature in units of 104 K. For temperatures typical of \[Fe II\] line-emitting regions in SNRs ($`T_e`$ $``$ 6 500 K; Oliva et al. 1989), collisional de-excitation becomes important for: $`n_e`$ $``$ 3.1 $`\times `$ 104 cm-3. The postshock densities are of the order of 103 cm-3 for the \[Fe II\] brightest remnants (Fig.3$`b`$). However, the densities derived from the optical \[S II\] doublet have been generally found to be systematically lower (by about a factor of 5) than those determined from the near-IR \[Fe II\] lines, suggesting that the \[Fe II\] features are produced further downstream behind the shock (e.g. Oliva et al. 1989). The electronic densities prevailing in the \[Fe II\]-line emitting regions of the SNRs in our sample are thus less than 5 $`\times `$ 103 cm-3, but are still much below the critical density for collisional de-excitation. In a two-level approximation, the luminosity of the \[Fe II\] $`\lambda `$1.644 $`\mu `$m transition is given by (Blietz et al. 1994):

$$L_{[\mathrm{Fe}\mathrm{II}]}=\pi ^2\frac{T_4^{0.94}e^{1.57/T_4}𝒩_{\mathrm{Fe}^+}}{1+n_c(T_4)/n_e}\left(\frac{d}{\mathrm{cm}}\right)^2\mathrm{ergs}\mathrm{s}^1$$
(3)
Here $`𝒩_{\mathrm{Fe}^+}`$ is the column density of Fe+ in units of 1016 cm-2, $`n_c`$ is the critical density given by equation (2) and $`d`$ is the diameter of the SNR. In view of the similar size of the detected SNRs (see Fig.2), we will ignore in the following any dependence of $`L_{[\mathrm{Fe}\mathrm{II}]}`$ on the latter quantity. While $`L_{[\mathrm{Fe}\mathrm{II}]}`$ is only weakly dependent on the temperature (it varies within a factor of only 1.5 over 5 decades in density between 6 000 and 20 000 K), it can be expressed to a reasonable degree of accuracy as a linear function of $`n_e`$ for densities up to 104 cm-3. The relationship found in our data between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`n_e`$ ($`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ $`n_e^{1.244\pm 0.188}`$) is close to what might be expected from theoritical considerations. We will return to this point in the following.

The postshock densities derived for the detected SNRs are 2 to 4 orders of magnitude higher than canonical values for the ISM ($`n_0`$ $``$ 0.1–10 cm-3). We have previously argued that the \[Fe II\] brightest SNRs are in the radiative expansion phase. Large compression factors (i.e. ratio of postshock to preshock density) can be achieved for strong radiative shocks, with values in the range 10–50 being typical (e.g. Vancura et al. 1992). While we cannot rule out the possibility that the high postshock densities inferred for the detected SNRs are due to very large compression factors, we rather favour the idea that they are evolving in an ambient medium of high density. The tendency for the objects with high radio fluxes to be strong \[Fe II\] emitters supports this picture (Table 3). Since thermal X-ray emission from the shock-heated gas (which should dominate X-ray excitation in our sample) scales with ambient density (e.g. Magnier et al. 1997), X-ray data can also be used to assess the relevance of this interpretation. However, we find no evidence in our data for a relationship between the X-ray and the \[Fe II\] luminosities. A better knowledge of the X-ray properties of SNRs in M33 might shed some new light on this issue.

### 5.3 Metallicity effects

Abundances derived from modelling the optical-line properties of SNRs reflect to a large extent the chemical composition of the swept-up interstellar material (e.g. Russell & Dopita 1990), and can be characterized by the “metallicity index”, $`A`$, defined as: $`\mathrm{log}(A)`$ = 0.8 + 1/3\[$`\mathrm{log}(N/H)`$ \+ $`\mathrm{log}(O/H)`$ \+ $`\mathrm{log}(S/H)`$\] (Dopita et al. 1984). We find a correlation between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`A`$ which might indicate that \[Fe II\] emission is enhanced in regions characterized by high ISM chemical abundances. The column density of Fe+, $`𝒩_{\mathrm{Fe}^+}`$, in equation (3) can be rewritten as:

$$𝒩_{\mathrm{Fe}^+}=f_1\delta X_{\mathrm{Fe}}^{}𝒩_\mathrm{H}\mathrm{cm}^2$$
(4)
where $`f_1`$ is the ionization fraction of Fe+, $`\delta `$ is the gas phase iron abundance, $`X_{\mathrm{Fe}}^{}`$ is the solar abundance of iron relative to hydrogen and $`𝒩_\mathrm{H}`$ is the column density of neutral and singly ionized hydrogen. It can be seen that the correlation found in our data between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`A`$ is consistent with theoritical predictions.

One major difficulty with this interpretation lies, however, in the existence of a positive correlation between $`n_e`$ and $`A`$ (see Table 3 and Fig.4). While there is no evidence in our data for a negative correlation between $`n_e`$ and galactocentric distance, this might suggest that both the cosmic abundances (see Smith et al. 1993) and the density of the ISM increase when progressing towards the inner parts of M33. In view of this unexpected result, it is not clear whether the correlation between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`A`$ is fundamental or is mainly a density effect. It is conceivable, however, that a spread in the metallicity of the ambient ISM among our sample of SNRs contributes to the large scatter observed in the $`L_{[\mathrm{Fe}\mathrm{II}]}`$-$`n_e`$ relation. Variations in the depletion factor might also be relevant in this respect, as both theory and observations suggest that the fraction of iron-bearing grains returned to the gas phase by shocks may significantly vary from one SNR to another (Jones, Tielens & Hollenbach 1996; Oliva et al. 1999). The positive correlation between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and the iron cosmic metal abundance hinted at by equation (4) should be kept in mind when comparing the \[Fe II\] properties of SNR populations in host galaxies with widely different metallicities.

## 6 Discussion

### 6.1 Near-IR \[Fe II\] emission from SNRs

In the light of the results presented above, we conclude that strong near-IR \[Fe II\] emission in SNRs is subordinated to the existence of a radiative blast wave propagating through an ambient material of high density (other physical processes than shocks might account for \[Fe II\] emission in very young SNRs; Graham et al. 1990). This picture is supported by the shock models of Hollenbach & McKee (1989). The \[Fe II\] brightest SNRs might be at the onset of the radiative phase, where line luminosities are believed to reach their maximum (Falle 1981; Cioffi & McKee 1988). We note that the relatively small size of our detected SNRs (Fig.2) is not in conflict with the fact that they have reached such an advanced evolutionary stage; their confined nature may cause them to evolve rapidly. The low detection rate of our survey may result from: (i) the finite duration of strong \[Fe II\] emission in SNRs (see Section 6.4) and (ii) the fact that optically-selected samples of SNRs are biased in favour of objects evolving in a teneous medium, possibly the warm component of the ISM (e.g. Pannuti et al. 2000).

The following empirical relation has been found in our data between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`n_e`$ (see Table 3): $`L_{[\mathrm{Fe}\mathrm{II}]}`$ (L) $``$ 0.257 $`n_e^{1.244\pm 0.188}`$ (cm-3). Since this exponent of 1.244 $`\pm `$ 0.188 is very close to the value of unity expected on theoritical grounds (see Section 5.2), we shall consider in the following that:

$$L_{[\mathrm{Fe}\mathrm{II}]}=(1.1\pm 0.3)\left(\frac{n_e}{\mathrm{cm}^3}\right)\mathrm{L}_{}$$
(5)
The factor 1.1 is simply the weighted average of $`L_{[\mathrm{Fe}\mathrm{II}]}`$/$`n_e`$. The equation above can be used to estimate the \[Fe II\] $`\lambda `$1.644 $`\mu `$m line luminosity of radiative SNRs from the sole knowledge of the electronic density in the postshock region (since this relation refers to the density in the \[S II\]-emitting region, care should be taken when using other density diagnostics). We caution, however, that this relation has been derived on the basis of few detections that only sample a limited density regime. A better knowledge of the \[Fe II\] properties of SNRs in host galaxies with widely different ISM properties is needed to fully assess the robustness of this relationship, especially with respect to metallicity effects.

### 6.2 Comparison with Galactic and LMC SNRs

Is this picture consistent with \[Fe II\] observations of SNRs in our Galaxy and in the LMC? These observations are summarized in Table 4, along with the typical postshock densities associated with the SNRs considered. Great care has been taken to ensure that the density estimates (as determined from the \[S II\] doublet) relate to the same regions observed in \[Fe II\]. One major difficulty in comparing the \[Fe II\] luminosities of SNRs in M33 with their counterparts in our Galaxy or in the LMC lies in the extended nature of the latter class of objects. Most \[Fe II\] luminosities quoted in Table 4 are derived from aperture photometry and, in most cases, only sample a very small fraction of the SNR line-emitting regions. We performed an aperture correction to estimate the total \[Fe II\] luminosities by considering the ratio of the total optical size of the remnant by the area sustained by the aperture used.111We note that these corrections are likely to be grossly overestimated since the \[Fe II\] measurements are generally performed on the brightest optical knots where we expect the highest level of \[Fe II\] emission (see Section 5.1). For instance, we derive from the observations of Rho et al. (2001) a total dereddened \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosity for IC 443 of about 55 L. This value is much smaller than our estimate of 1 200 L derived from the measurements of Graham et al. (1987) through a 35 arcsec $`\times `$ 35 arcsec aperture. We show in Fig. 5 the resulting \[Fe II\] luminosities of SNRs in M33, our Galaxy and the LMC, as a function of $`n_e`$. A reasonable agreement between the expected and observed \[Fe II\] luminosities of IC 443, the Cygnus Loop, N49, RCW 103 and N63A is found considering the crudeness of the corrections applied. Interestingly, the 3 SNRs (Crab, Kepler and N103B) that fail by more than 2 orders of magnitude to follow the relationship between $`L_{[\mathrm{Fe}\mathrm{II}]}`$ and $`n_e`$ are the objects in our sample of Galactic or LMC SNRs most unlikely to undergo radiative expansion. The Crab nebula and Kepler’s remnant are both very young and therefore are still in the free (or Sedov-Taylor) expansion phase. The same is true for N103B which is likely to be still undergoing adiabatic expansion (Dickel & Milne 1995). The weakness of the \[Fe II\] emission in N103B and in Kepler’s remnant illustrates the important point that the existence of radiative shocks is a necessary condition for strong \[Fe II\] emission, but not a high-density ambient ISM. As can be seen in Fig.3b, the \[Fe II\] luminosities of SNRs with similar densities may differ by as much as 2 orders of magnitude. This is best explained by evolutionary effects, with the low \[Fe II\] emitters being dominated by adiabatic shocks.

### 6.3 Implications for the \[Fe II\] emission of SNRs in starburst galaxies

The fact that the level of \[Fe II\] emission scales with density may help to understand the differences that exist between the near-IR \[Fe II\] properties of SNR populations in galaxies with moderate star formation activity (such as M33) and starburst galaxies. While the present study demonstrates that the brightest remnants in M33 have \[Fe II\] luminosities of the order of 700 L, SNRs in the starburst galaxy M82 have luminosities up to 1.6 $`\times `$ 105 L (Greenhouse et al. 1997). Much higher values are also observed in NGC 253 (Forbes et al. 1993). We compare below the \[Fe II\] luminosities expected within the framework of the results presented above with the observed values.

The electronic density in the \[Fe II\]-emitting regions of NGC 253 has been directly measured from density-sensitive near-IR \[Fe II\] line ratios and varies from 5 $`\times `$ 103 cm-3 (Engelbracht et al. 1998) to 1 $`\times `$ 104 cm-3 (Simpson et al. 1996). The \[Fe II\] emission in NGC 253 is believed to mainly arise from SNRs in the disk of the galaxy, and not from superwind- or AGN-related phenomena (Forbes et al. 1993). It is therefore reasonable to assume that these densities are typical of the postshock regions of SNRs in the central regions of this galaxy. Keeping in mind that our empirical relation (equation 5) relates to the density derived from the \[S II\] doublet and that the density in the \[Fe II\]-emitting region is about 5 times higher (e.g. Oliva et al. 1989), we estimate for these SNRs: $`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 1.1–2.2 $`\times `$ 103 L. This is in good agreement with the observed value: $`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 2.6 $`\times `$ 103 L (Forbes et al. 1993).

The ISM electronic densities prevailing in the central regions of M82 show a large spread: $`n_0`$ = 50–500 cm-3 (Förster-Schreiber et al. 2001, and references therein). We shall consider in the following: $`n_0`$ = 200 cm-3. This translates into: $`n_e`$ $``$ 6 $`\times `$ 103 cm-3, if one considers a canonical value for the compression factor of the SNR blast wave of 30. This estimate is consistent with the upper limit of 3.2 $`\times `$ 104 cm-3 derived by Lester et al. (1990) from \[Fe II\] lines. According to equation (5), \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosities approaching 104 L can easily be reached for these SNRs. However, it appears that density effects alone have difficulties in explaining the high luminosities of the \[Fe II\] brightest sources in M82 (up to 1.6 $`\times `$ 105 L; Greenhouse et al. 1997). Two effects might account for this discrepancy. First, the SNR might expand in the interclump medium of molecular clouds, with the result that we severely underestimate $`n_e`$ (Chevalier & Fransson 2001). Second, in regard of the high space density of SNRs in the compact star-forming regions of this galaxy (Huang et al. 1994), it is far from clear that the strong \[Fe II\] sources mapped by Greenhouse et al. (1997) can be identified with individual SNRs. The radius at which a SNR enters the radiative phase scales with the initial blast energy, $`E_{51}`$, and upstream hydrogen number density, $`n_\mathrm{H}`$, as: $`R_{\mathrm{rad}}`$ $``$ $`E_{51}^{5/17}`$$`n_\mathrm{H}^{7/17}`$ (Blondin et al. 1998). Regardless of the details of the original explosion, one therefore expects the bright \[Fe II\] SNRs in starburst galaxies to be less spatially extended (by a factor 2–5) than the detected SNRs in M33 (with $`\overline{d}`$ $``$ 21 pc). It is thus plausible that the exceptionally strong \[Fe II\] sources in M82 (with sizes up to 40 pc) are made up of several individual SNRs with diameters in the range 4–10 pc.

We conclude that the dichotomy between the \[Fe II\] properties of SNRs in quiescent and in starburst galaxies may be interpreted as a result of the different ISM densities prevailing in both type of galaxies. Although clear metallicity effects have not been uncovered from our data, we also note that a substantial enrichment in metals is expected in the hostile environments of starburst regions (because of high chemical abundances and/or processing via grain destruction of the large reservoir of dust grains) and might also potentially contribute to strong \[Fe II\] emission.

### 6.4 The \[Fe II\]-emitting lifetime of SNRs

To estimate the typical \[Fe II\]-emitting lifetime of a SNR, we consider the integrated \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosity corrected for extinction in the central 2.4 arcsec $`\times `$ 12 arcsec region of NGC 253: $`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 1.3 $`\times `$ 105 L (Engelbracht et al. 1998). A distance of 2.5 Mpc has been adopted (De Vaucouleurs 1978). We have also assumed a homogeneous mixture of gas and dust with an extinction towards the \[Fe II\]-emitting regions of $`A_V`$ = 8.4 mag (Engelbracht et al. 1998) and the interstellar extinction law of Rieke & Lebofsky (1985). The global radio supernova rate of NGC 253 is uncertain, but estimates range from 0.1 to 0.2 yr-1 (Ulvestad 2000, and references therein). A value $`\eta `$ $``$ 0.15 yr-1 is adopted in the following. To translate this into a normalized value for the central 2.4 arcsec $`\times `$ 12 arcsec region, we consider the ratio of the area sustained by this aperture by the total radio extent of the galaxy, which gives: 9 $`\times `$ 10-3 yr-1. The electron density in the \[Fe II\]-emitting regions lies in the range: 5–10 $`\times `$ 103 cm-3 (Simpson et al. 1996; Engelbracht et al. 1998). Using equation (5) by converting to densities in the \[S II\]-emitting regions, and substituting the result in equation (1) yields: $`t_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 0.7–1.4 $`\times `$ 104 yr.

Similar calculations can be applied to M82. The integrated \[Fe II\] $`\lambda `$1.644 $`\mu `$m luminosity corrected for extinction in the central 16 arcsec $`\times `$ 10 arcsec region is: $`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 2.7 $`\times `$ 106 L (Förster-Schreiber et al. 2001). A distance of 3.3 Mpc (Freedman & Madore 1988) and the similar extinction model as above with $`A_V`$ = 36 mag has been adopted (Förster-Schreiber et al. 2001). By considering this relatively small area in M82, we minimize the contribution of large-scale, superwind-induced emission (Greenhouse et al. 1997). Radio observations yield a global supernova rate: $`\eta `$ = 0.11 $`\pm `$ 0.05 yr-1, which translates into a normalized value for the central 16 arcsec $`\times `$ 10 arcsec region of about 0.016 yr-1 (Huang et al. 1994). Additionally, we assume the typical \[Fe II\] luminosity for SNRs in M82 determined previously: $`L_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 104 L. Substituting in equation (1) yields: $`t_{[\mathrm{Fe}\mathrm{II}]}`$ $``$ 1.7 $`\times `$ 104 yr.

In regard of the large uncertainties involved, the values of $`t_{[\mathrm{Fe}\mathrm{II}]}`$ for SNRs in NGC 253 and M82 appear to be in good agreement. The \[Fe II\]-emitting lifetime derived for NGC 253 is more robust, as the density of the \[Fe II\]-emitting gas has been directly measured in this galaxy. A value of 1 $`\times `$ 104 yr will thus be adopted in the following. This estimate is broadly consistent with the duration of strong optical emission in SNRs, as determined by evolutionary models (e.g. Franco et al. 1994).

### 6.5 The supernova rate of starburst galaxies

By using equation (5) and the \[Fe II\]-emitting lifetime derived above, the supernova rate given by equation (1) can be written as a function of the postshock electron density as:

$$\eta \frac{_{[\mathrm{Fe}\mathrm{II}]}/\mathrm{L}_{}}{1.1\times 10^4}\left(\frac{n_e}{\mathrm{cm}^3}\right)^1\mathrm{yr}^1$$
(6)
In the absence of density diagnostics for the postshock regions, this expression can be conveniently expressed as a function of the (more easily measurable) ISM number density, $`n_0`$, by an appropriate choice of the typical compression factor of the SNR blast wave.

In regard of the various assumptions made here and the fact that metallicity issues still remain to be fully addressed, the supernova rate we derive is still uncertain. In spite of these limitations, we note that our expression is consistent with previous estimates (van der Werf et al. 1993; Vanzi & Rieke 1997) if canonical ISM densities for starburst galaxies are considered ($`n_0`$ $``$ 100 cm-3). The relation proposed by van der Werf et al. (1993) has been found to give a good agreement in a sample of starburst galaxies between supernova rates derived from near-IR \[Fe II\] lines and from their star-forming properties (Calzetti 1997). This provides some support to the validity of equation (6) to provide robust estimates of the supernova rates in galaxies with widely different ISM properties and for which the contribution of superwinds- and AGN-related phenomena to the total \[Fe II\] output can be neglected. Because only radiative SNRs significantly contribute to $`_{[\mathrm{Fe}\mathrm{II}]}`$, however, these values are only lower limits.

## 7 Conclusions

The main conclusions of this study are the following:

* We have carried out the first near-IR \[Fe II\] line-imaging survey of extragalactic SNRs, presenting evidence for large intrinsic differences in their \[Fe II\] properties.
* A significant correlation is found between the optical- and \[Fe II\]-line luminosities, suggesting that the detected SNRs are dominated by radiative shocks.
* We suggest that the SNRs with strong \[Fe II\] emission are evolving in a dense ISM. Density effects are most likely the cause of the very strong level of \[Fe II\] emission from SNRs in starburst galaxies.
* The typical \[Fe II\]-emitting lifetime of a SNR in starburst galaxies is found to be of the order of 104 yr.
* We provide a new empirical expression for the supernova rate of starburst galaxies, as derived from their integrated near-IR \[Fe II\] luminosity.

## Acknowledgments

We wish to thank an anonymous referee for useful comments. We acknowledge the TACs of the CFHT for their generous time allocation. The Canada-France-Hawaii Telescope is operated by the National Research Council of Canada, the Centre Nationale de la Recherche Scientifique of France, and the University of Hawaii. The authors also wish to thank the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Fonds pour la Formation de Chercheurs et l’Aide à la Recherche (FCAR) of Québec for financial support. T. M. would like to thank the organizers of the conference “The interstellar medium in M31 and M33” (where the first results of this project were presented) for this enjoyable meeting, as well as the Wilhelm und Else Heraeus-Stiftung fundation for financial support.


format: text

text:
# Afterglow Light Curve Modulated by a Highly Magnetized Millisecond Pulsar

## 1 Introduction

Gamma-ray bursts (GRBs) are widely accepted to be produced when fast-moving, relativistic shells ejected from a central source in a relatively short period collide with slowly moving, yet relativistic shells that were ejected at an earlier time (Rees & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 1994; Paczy$`\stackrel{´}{\mathrm{n}}`$ski & Xu 1994; Kobayashi et al. 1997; Daigne & Mochkovitch 1998, 2000). In connection with the so-called internal shock model, the external shock model also prevails as a possible origin of the GRB afterglows. In the external shock model the relativistic material is assumed to be decelerated via interactions with the surrounding medium. A shock wave results in heating the ambient matter to relativistic temperatures, and emitting photons in longer wave lengths (Rees & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 1992; $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ & Rees 1993; Paczy$`\stackrel{´}{\mathrm{n}}`$ski & Rhoads 1993; Sari et al. 1996; Vietri 1997). GRBs and their afterglows seem to result from the dissipation of bulk energy in the relativistic outflows, which are in the form of a narrow beam rather than a spherical shell.

Even though the origin of the observed GRBs are still unknown, from the observations of several GRB afterglows the evidence of beamed GRBs has been accumulated (Sari et al. 1999; Halpern et al. 1999; Rhoads 1999). There are several works on models for the geometry of GRBs (e.g., Chang & Yi 2001 and references therein) and their environments (e.g., Scalo & Wheeler 2001 and references therein). Much of the current research on GRBs is aimed at determining the nature and the origin of the central engine (Duncan & Thompson 1992; Narayan et al. 1992; Woosley 1993; Katz 1994; Usov 1992, 1994a; Shaviv & Dar 1995; $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ & Rees 1997b; Yi & Blackman 1997; Blackman & Yi 1998; Paczy$`\stackrel{´}{\mathrm{n}}`$ski 1998; MacFadyen & Woosley 1999; Portegies Zwart et al. 1999; Li 2000; Wheeler et al. 2000; Zhang & Fryer 2001). Although the simple cosmological fireball afterglow model is in a good agreement with the observed light curves of a power law decay (e.g., Wijers et al. 1997; Waxman 1997a,b), the afterglow of GRB970228 observed by HST, for instance, shows a deviation from the simple power law behavior (Fruchter et al. 1997). Several works have been done to further investigate more subtle effects that can change the afterglow characteristics (Katz & Piran 1997; $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ & Rees 1997a, 1999; Rhoads 1999; Berger et al. 2000; Dai & Lu 2000; Kumar & Panaitescu 2000a, b; Panaitescu & Kumar 2000; Dai & Lu 2001).

Most fireball models assume that the energy injection into the fireball occurs in a short period of time compared with the lifetime of the afterglows (Rees & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 1998; Kumar & Piran 2000; Sari & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 2000). However, in some types of central engines, such as a fast rotating new-born pulsar with the strong magnetic field (magnetar), a significant energy input into the fireball may in principle continue for a significantly longer timescale, and accordingly the temporal decay of the afterglow will be slower. Hence, it is worthwhile to investigate a continuously fed fireball in more details as a probe of the central engine of GRBs. We here consider the central engine that emits both an initial impulsive energy input $`E_{\mathrm{imp}}`$ and a continuous power. In fact, recently there were such attempts to provide an explanation for the deviation of the afterglow light curve from the simple power law (Dai & Lu 1998a, b; Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 2001a, b). Even though previous studies have considered continuous injection from a highly magnetized millisecond pulsar to interpret the deviation of afterglow light curves of some GRBs from the power law behavior, one of the most important aspects of pulsar emissions is ignored, that is, beaming. Beaming and continuous powering are the clear consequences of the pulsar origin of the GRB afterglows. In this Letter, this feature is included and a change of the beaming due to a dynamical evolution of the new born pulsar is explicitly taken into account. We find that the deviations suggested by previous studies indeed occur, but that the shape of the light curve is significantly modified due to a beaming effect. If such observational evidence for beaming is found, the corresponding pulsar origin for GRBs would imply that activities resulting in too massive objects to be a neutron star might be ruled out as a central engine of GRBs. Therefore, the question of whether or not a signature of beaming can be found in the afterglows thus extremely important.

## 2 Afterglows Continuously Powered by Pulsar

The total luminosity emitted from a young millisecond pulsar (MSP) has two important terms: an electromagnetic (EM) radiation term and a gravitational wave (GW) radiation term. Given that the spin-down is mainly due to electromagnetic dipolar radiation and to gravitational wave radiation, the spin-down law is given by

$`I\mathrm{\Omega }\dot{\mathrm{\Omega }}={\displaystyle \frac{B_p^2R^6\mathrm{\Omega }^4\mathrm{sin}^2\alpha }{6c^3}}+{\displaystyle \frac{32Gϵ^2I^2\mathrm{\Omega }^6}{5c^5}},`$ (1)
where $`\mathrm{\Omega }`$ and $`\dot{\mathrm{\Omega }}`$ are the angular frequency and its time derivative, $`B_p`$ is the dipolar field strength at the magnetic poles, $`\alpha `$ is the angle between spin axis and the magnetic field axis, $`I`$ is the moment of inertia, $`R`$ is the pulsar radius, and $`ϵ`$ is the ellipticity of the compact star (Shapiro & Teukolsky 1983). When an MSP rotates with spin frequency larger than a certain critical frequency, $`\mathrm{\Omega }_{\mathrm{crit}}`$, nonaxisymmetric secular instability drives the pulsar into nonaxisymmetric configuration with nonzero quadrupole moment. This is the case where the GW luminosity dominates. In this case, the spin-down of the MSP occurs on a timescale $`\tau _{\mathrm{gw}}=I_{}\mathrm{\Omega }_{}^2/2L_{\mathrm{gw}}3\times 10^3ϵ^2I_{45}\mathrm{\Omega }_4^4`$ s, where $`I_{}=I_{45}10^{45}\mathrm{g}\mathrm{cm}^2`$ is the MSP moment of inertia, $`\mathrm{\Omega }_{}=\mathrm{\Omega }_410^4\mathrm{s}^1`$ is the MSP spin frequency, and we have made use of the gravitational wave energy-loss rate, $`L_{\mathrm{gw}}=32Gϵ^2I_{}^2\mathrm{\Omega }_{}^6/5c^5`$. When an MSP rotates with spin frequency smaller than the certain critical frequency, the spin-down due to the electromagnetic dipole radiation gives a timescale $`\tau _{\mathrm{em}}I_{}\mathrm{\Omega }_{}^2/2L_{\mathrm{em}}3\times 10^2I_{45}B_{15}^2R_6^6\mathrm{\Omega }_4^2\mathrm{s}`$, where $`R_{}=R_610^6`$ cm is the MSP radius, $`B_{}=B_{15}10^{15}`$ G is the MSP dipole magnetic field, and we have assumed the electromagnetic dipole energy-loss rate, $`L_{\mathrm{em}}=2\mu _{}^2\mathrm{\Omega }_{}^4/3c^3`$, where $`\mu _{}=B_{}R_{}^3`$ is the electromagnetic dipole moment. The decay solution $`\mathrm{\Omega }(T)`$ as a function of time $`T`$ includes both EM and GW emissions, but the energy input into the fireball is due to the EM dipolar emission $`L_{\mathrm{em}}`$ alone. Provided that the spin-down at a time is dominated by one or the other loss term, one can obtain approximate solutions. When GW radiation losses dominate the spin-down, we have $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{gw}}(1+T/\tau _{\mathrm{gw}})^{1/4}`$, or approximately $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{gw}}`$ for $`T\tau _{\mathrm{gw}}`$, and $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{gw}}(T/\tau _{\mathrm{gw}})^{1/4}`$ for $`T\tau _{\mathrm{gw}}`$, where $`\mathrm{\Omega }_{0,\mathrm{gw}}`$ being the initial spin frequency. When EM dipolar radiation losses dominate the spin-down, we have $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{em}}(1+T/\tau _{\mathrm{em}})^{1/2}`$, or approximately $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{em}}`$ for $`T\tau _{\mathrm{em}}`$, and $`\mathrm{\Omega }=\mathrm{\Omega }_{0,\mathrm{em}}(T/\tau _{\mathrm{em}})^{1/2}`$ for $`T\tau _{\mathrm{em}}`$, where $`\mathrm{\Omega }_{0,\mathrm{em}}`$ being the initial spin frequency. Assuming $`\mathrm{\Omega }_{\mathrm{crit}}>\mathrm{\Omega }_{\mathrm{dynamo}}`$, where $`\mathrm{\Omega }_{\mathrm{dynamo}}`$ is the spin required for dynamo action, there are two types of MSPs: a supercritical strong field rotator (SPS) with $`\mathrm{\Omega }>\mathrm{\Omega }_{\mathrm{crit}}>\mathrm{\Omega }_{\mathrm{dynamo}}`$, and a subcritical strong field rotator (SBS) with $`\mathrm{\Omega }_{\mathrm{crit}}>\mathrm{\Omega }>\mathrm{\Omega }_{\mathrm{dynamo}}`$. These two classes may give the bimodal distribution of short and long pulsar-induced GRBs (Yi & Blackman 1998). For an SPS pulsar, GW spin-down is important in initial period and promptly below $`\mathrm{\Omega }_{\mathrm{crit}}`$. As soon as $`\mathrm{\Omega }<\mathrm{\Omega }_{\mathrm{crit}}`$ the spin-down becomes dominated by the EM losses. On the other hand, for an SBS pulsar, the spin-down is always dominated by EM regime.

The differential energy conservation relation for the self-similar blast wave can be written as $`dE/dt=_0(t/t_0)^q^{}\kappa ^{}(E/t)`$, where $`E`$ and $`t`$ are the energy and time measured in the fixed frame and $`q^{}`$ and $`\kappa ^{}`$ are constants (Cohen & Piran 1999). The first term denotes the continuous luminosity injection, and the second term takes into account radiative energy losses in the blast wave. For $`t>t_0`$, the bulk Lorentz factor of the fireball scales with time as $`\mathrm{\Gamma }^2t^m`$, with $`m`$ and $`\kappa ^{}`$ related by $`\kappa ^{}=m3`$ (Cohen et al. 1998). If $`m=3`$, it corresponds to the adiabatic case (Blandford & McKee 1976). In the observer frame, the time $`T`$ is related to the fixed frame $`t`$ by $`dT=(1\beta )dtdt/2\mathrm{\Gamma }^2`$, and $`T=_0^t(2\mathrm{\Gamma }^2)^1𝑑tt/[2(m+1)\mathrm{\Gamma }^2]`$ when $`tt_0`$. The differential energy conservation relation in the observer frame is now given by $`dE/dT=L_0(T/T_0)^q\kappa (E/T)`$, and can be integrated as

$`E={\displaystyle \frac{L_0}{\kappa +q+1}}\left({\displaystyle \frac{T}{T_0}}\right)^qT+E_{\mathrm{imp}}\left({\displaystyle \frac{T}{T_0}}\right)^\kappa ,T>T_0,`$ (2)
where $`T_0=t_0/[2(m^{}+1)\mathrm{\Gamma }^2]`$, $`m^{}`$ being the self-similar index for $`t<t_0`$, $`L_0=2\mathrm{\Gamma }^2_0`$, $`q=(q^{}m)/(m+1)`$, and $`\kappa =\kappa ^{}/(m+1)`$. Setting $`T=T_0`$, the total energy at the beginning of the self-similar expansion is the sum of two terms, $`E_0=L_0T_0/(\kappa +q+1)+E_{\mathrm{imp}}`$. The first term is the accumulated energy from the continuous injection before the self-similar solution begins. The second term $`E_{\mathrm{imp}}`$ is the energy injected impulsively by the initial blast.

The total energy of the blast wave given by Eq. 2 may be dominated either by the continuous injection term ($`T^{q+1}`$) or by the initial impulsive term ($`T^\kappa `$), subject both to the relative values of the two indices and to the values of $`L_0`$ and $`E_{\mathrm{imp}}`$ (Dai & Lu 1998a, b; Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 2001a). One may classify three regimes according to the relative values of the two indices as discussed by Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ (2001a). We are interested in the case where $`q>1\kappa `$, since otherwise a pulsar signature is no longer observable even if there is a pulsar in the central engine. If $`q>1\kappa `$, the first term in Eq. 2 will eventually dominate over the second term after a critical $`T_c`$. The injection-dominated regime begins at a critical time $`T_c`$ defined by equating the injection and energy-loss terms in Eq. 2,

$`T_c=\mathrm{Max}\{1,\left[(\kappa +q+1){\displaystyle \frac{E_{\mathrm{imp}}}{L_0T_0}}\right]^{1/(\kappa +q+1)}\}.`$ (3)
If initially the impulsive term dominates, the critical time $`T_c`$ could be much longer than $`T_0`$, depending on the ratio of $`E_{\mathrm{imp}}`$ and $`L_0T_0`$. The case where $`T_cT_0`$ ensures that a self-similar solution has already formed when the continuous injection term dominates. Besides the condition of $`q>1\kappa `$, for detecting the signature of a pulsar the critical time $`T_c`$ should be greater than $`\tau _{\mathrm{em}}`$ or $`\tau _{\mathrm{gw}}`$ as discussed below.

To obtain the temporal decay index of the afterglow light curve for which the MSP is responsible we adopt the cylindrical geometry instead of the spherical geometry, which may accommodate elongated beaming configurations. The rotational axis of the MSP coincides with the $`z`$-axis of the geometry. For a fireball blastwave decelerated by a homogeneous external medium with particle number density $`n`$, the energy conservation equation at time $`t=r/c`$ is given by

$`E=2\pi R^2rnm_pc^2\mathrm{\Gamma }^2=2\pi R^2ctnm_pc^2\mathrm{\Gamma }^2,`$ $`t>t_0,`$ (4)
where $`r`$ is the distance of the blastwave in $`z`$ direction, $`R`$ is a radius of the beaming which can be characterized by that of the light cylinder, $`c/\mathrm{\Omega }`$, $`\mathrm{\Omega }`$ being given by the MSP evolution, and all other symbols have their usual meanings. One can work out the relationship between the temporal index $`\alpha `$ and the spectral index $`\beta `$, where $`F_\nu T^\alpha \nu ^\beta `$. To directly compare indices suggested by Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ (2001a), we derived the scaling laws for the slow cooling regime assuming that the reverse shock is mildly relativistic (Sari et al. 1998). We derive the scaling law for $`T\tau _{\mathrm{em}}`$ as an illustrating example. The energy E in Eq. 2 should have the same time dependence as the first term on the right-hand side of Eq. 1, giving $`T^{q+1}t\mathrm{\Gamma }^2`$, or $`m=q/(q+2)`$. Since in general, $`\mathrm{\Gamma }t^{m/2}r^{m/2}T^{m/2(m+1)}`$ and $`rT^{1/(m+1)}`$, the $`\mathrm{\Gamma }`$ dependence on $`T`$ is given by $`\mathrm{\Gamma }T^{q/4}`$. The relationship between the temporal index $`\alpha `$ and the spectral index $`\beta `$ can be obtained as below. For the forward shock, the synchrotron peak frequency $`\nu _m^f\mathrm{\Gamma }^4T^{2m/(m+1)}T^q`$, the peak flux $`F_{\nu _m}^fT^3\mathrm{\Gamma }^8T^{(3m)/(m+1)}T^{3+2q}`$. Thus, $`\alpha ^f=3+2qq\beta =(2m\beta ^f+3m)/(1+m)`$ since $`T^{\alpha ^f}\nu _m^{\beta ^f}T^{3+2q}T^{(3m)/(m+1)}`$. For the reverse shock, $`\nu _m^r=\nu _m^f/\mathrm{\Gamma }^2\mathrm{\Gamma }^2T^{m/(m+1)}T^{q/2}`$, the peak flux $`F_{\nu _m}^r=\mathrm{\Gamma }F_{\nu _m}^fT^3\mathrm{\Gamma }^9T^{(63m)/2(m+1)}T^{(12+9q)/4}`$. Thus, $`\alpha ^r=(2m\beta ^r+63m)/2(1+m)=(12+9q2q\beta )/4`$. Similarly, one can derive the temporal index for each regime. In Table 1 we have summarized temporal indices in various regimes for the forward shock and the reverse shock, respectively.

If the continuous injection term becomes dominant over the impulsive term after $`T_c`$, the afterglow light curves rises after $`T_c`$ and steepen again after some time, that is, about $`\tau _{\mathrm{em}}`$ or $`\tau _{\mathrm{gw}}`$. Therefore, there may be two types of afterglow patterns for continuous injections according to a appropriate combination. Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ (2001a) discussed conditions which allow to detect a signature for a pulsar and concluded that physical parameters are consistent with those of a magnetar in case of the afterglow features in case of GRB 000301c. In practice, however, $`\tau _{\mathrm{gw}}`$ is very short, can be even shorter than $`T_c`$ unless ambient matter density is very high, and is therefore unlikely to be observed. At around $`T_c`$ and $`\tau _{\mathrm{em}}`$, $`q`$ changes -1 to 0 and 0 to -2, respectively. These scaling laws represent a change from the standard adiabatic case to an EM-loss dominated regime as shown above (see also Zhang & $`\mathrm{M}\stackrel{´}{\mathrm{e}}\mathrm{sz}\stackrel{´}{\mathrm{a}}\mathrm{ros}`$ 2001a). For the forward shock, the temporal decay index changes around $`T_c`$ from $`3\beta ^f/2`$ to $`3`$ and returns to $`3\beta ^f3`$ after $`\tau _{\mathrm{em}}`$. For the reverse shock, the temporal decay index changes around $`T_c`$ from $`3\beta ^r/43/8`$ to $`3`$ and returns to $`3\beta ^r/215/4`$ after $`\tau _{\mathrm{em}}`$.

## 3 Discussions

A continuous energy injection signature in the GRB afterglow light curve may directly provide diagnostics about the nature of the injection as well as information on the GRB progenitor. Therefore, the question of whether or not a bump in the afterglow light curve is such an observational evidence for beaming from a pulsar is thus extremely important. In this sense, the correct geometry should be applied to the model calculation. We have discussed the case where a pulsar is continuously injecting energy cylindrically rather than spherically. We suggest that a possible explanation for the deviation of the afterglow light curve of GRB 970508 is due to a beaming from a central pulsar. The optical afterglow of GRB 970508 has been explained as evidence of gravitational lensing event (Loeb & Perna 1998; Dado et al. 2001). So far the gravitational lensing can only account for the fact that the GRB afterglow shows a rise and a fall. We show in this Letter that with a correctly assumed geometry of energy injection the afterglow of GRB 970508 can be explained by a strongly magnetized fast-rotating pulsar with the continuous energy injection in a beam. In this calculation we ignore effects of the lateral expansion of the jet which may occur in the later time of the jet evolution. The light curve should be modified and becomes that of the spherical geometry case if the sideways expansion occurs in a timescale comparable to those we discussed after a neutron star formed. For instance, if the jet is expanded much faster at earlier stage than the light cone evolution, the cylindrical geometry effect becomes less obvious.

If the bump is caused by such a pulsar indeed, it puts constraints on the GRB progenitor models. That is, the corresponding pulsar origin for GRBs would imply that activities resulting in too massive objects to be a neutron star, such as, neutron star mergers, black hole formation, should be ruled out as a central engine of GRBs. During the birth of the neutron star, an initial fireball may occur through electromagnetic processes (Usov 1992). This has led to models in which GRBs are powered by rapidly spinning compact objects with strong magnetic fields (Blackman & Yi 1998; Blackman et al. 1996; Usov 1994b; Yi & Blackman 1997, 1998).

###### Acknowledgements.

We are grateful to the referee, Robert Mochkovitch, for useful comments and suggestions. We are grateful to Ethan Vishniac for hospitality while visiting Johns Hopkins University where this work began.


text:
# Conservation Laws and Particle Production in Heavy Ion Collisions*footnote **footnote *Invited talk presented at International Symposium on Hadrons and Nuclei, Seoul, Korea, 2001.

## I Introduction

There are in general two approaches to describe integrated particle yields measured in ultrarelativistic heavy ion collisions: (i) the microscopic transport models  and (ii) the macroscopic statistical thermal models . In this article we will discuss the statistical approach and show that it provides a very satisfactory description of experimental data. We will emphasize the importance of the conservation laws in the particular strangeness conservation when modelling particle chemical freezeout conditions.

Within a statistical approach, the production of particles is commonly described using the grand canonical (GC) ensemble, where the charge conservation is controlled by the related chemical potential. In this description a net value of a given U(1) charge is conserved on the average. The (GC) approach can be only valid if the total number of particles carrying quantum number related with this symmetry is very large. In the opposite limit of a small particle multiplicities, conservation laws must be implemented exactly and locally, i.e., the canonical (C) ensemble for conservation laws must be used . The local conservation of quantum numbers in the canonical approach severely reduces the phase space available for particle productions. This treatment of charge conservation is of crucial importance in the description of particle multiplicities in proton induced processes , in $`e^+e^{}`$  as well as in central heavy ion collisions at low beam energies .

In this article we describe the exact strangeness conservation in the context of relativistic statistical thermodynamics. A kinetic theory for the time evolution of particle production and the approach to the grand canonical and the canonical equilibrium distribution will be also introduced. Finally the example of the applications of the statistical model in (C) ensemble is presented in the context of low energy central as well as in high energy peripheral heavy ion collisions.

## II Statistical model and particle multiplicity

The exact treatment of quantum numbers in statistical mechanics has been well established for some time now . It is in general obtained by projecting the partition function onto the desired values of the conserved charges by using group theoretical methods. For our purpose we shall only consider the conservation laws related to the abelian U(1) symmetry group. In particular, we concentrate on strangeness conservation.

The basic quantity in the statistical mechanics describing a thermal properties of a system is the partition function $`Z(T,V)`$. In the (GC) ensemble,

$$Z^{GC}(T,V,\mu _Q)Tr[e^{\beta (H\mu _QQ)}]$$
(1)
where $`Q`$ is the conserved charge, $`H`$ the hamiltonian of the system, $`\mu _Q`$ is the chemical potential which plays the role of the Lagrange multiplier which guarantees that the charge $`Q`$ is conserved on the average in the whole system. Finally $`\beta =1/T`$ is the inverse temperature.

In the (C) ensemble the charge $`Q`$ is conserved exactly. Thus there is no more chemical potential under the trace and instead we calculate the partition function summing only these states which are carrying exactly the quantum number $`Q`$, that is

$$Z_Q^C(T,V)Tr_Q[e^{\beta H}]$$
(2)
The canonical and the grand canonical partition functions are related through the following cluster decomposition,

$$Z^{GC}(T,V,\lambda )=\underset{Q=\mathrm{}}{\overset{+\mathrm{}}{}}\lambda ^QZ_Q^C(T,V)$$
(3)
where the fugacity $`\lambda \mathrm{exp}(\beta \mu _Q)`$ and the sum is taken over all possible values of the charge $`Q`$.

For the (GC) partition function, which is well behaving analytic function of the fugacity $`\lambda `$, the above relation can be inverted and the canonical partition function with a given value of the charge $`Q`$ reads,

$$Z_Q=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi e^{iQ\varphi }\stackrel{~}{Z}(T,V,\varphi )$$
(4)
where the generating function $`\stackrel{~}{Z}`$ is obtained from the grand canonical partition function replacing the fugacity parameter $`\lambda `$ by the factor $`e^{i\varphi }`$,

$$\stackrel{~}{Z}(T,V,\varphi )Z^{GC}(T,V,\lambda e^{i\varphi })$$
(5)
The form of the generating function $`\stackrel{~}{Z}`$ in the above equation is model dependent. Having in mind the application of the statistical description to particle production in heavy ion collisions we calculate $`\stackrel{~}{Z}`$ in the ideal gas approximation, however, including all particles and resonances . This is not an essential restriction, because, describing the freeze-out conditions we are dealing with a dilute system where the interactions should not influence particle production anymore. We neglect any medium effects on particle properties. In general, however, already in the low-density limit, the modifications of resonance width or particle dispersion relations could be of importance . For the sake of simplicity, we use classical statistics, i.e. we assume temperature and density regime so that all particles can be treated using Boltzmann distributions.

In nucleus-nucleus collisions the absolute values of the baryon number, electric charge and strangeness are fixed by the initial conditions. Modelling particle production in statistical thermodynamics would in general require the canonical formulation of all these quantum numbers. From the previous analysis , however, it is clear, that in heavy ion collisions only strangeness should be treated exactly, whereas the conservation of baryon and electric charges can be described by the appropriate chemical potentials in the grand canonical ensemble.

Within the approximations described above and neglecting first the contributions from multi-strange baryons, the generating function in equation (4) has the following form,

$$\stackrel{~}{Z}(T,V,\mu _Q,\mu _B,\varphi )=\mathrm{exp}(N_{s=0}+N_{s=1}e^{i\varphi }+N_{s=1}e^{i\varphi })$$
(6)
where $`N_{s=0,\pm 1}`$ is defined as the sum over all particles and resonances having strangeness $`0,\pm 1`$,

$$N_{s=0,\pm 1}=\underset{k}{}Z_k^1$$
(7)
and $`Z_k^1`$ is the one-particle partition function defined as

$$Z_k^1\frac{Vg_k}{2\pi ^2}m_k^2TK_2(m_k/T)\mathrm{exp}(b_k\mu _B+q_k\mu _q)$$
(8)
with the mass $`m_k`$, spin-isospin degeneracy factor $`g_k`$, particle baryon number $`b_k`$ and electric charge $`q_k`$. The volume of the system is $`V`$ and the chemical potentials related with the electric charge and the baryon number are determined by $`\mu _q`$ and $`\mu _B`$ respectively.

With the particular form of the generating function equations (6,7,8) the $`\varphi `$-integration in equation (4) can be done analytically giving the canonical partition function for a gas with total strangeness $`S`$ :

$$Z_S(T,V,\mu _B,\mu _Q)=Z_0(T,V,\mu _B,\mu _Q)I_S(x)$$
(9)
where $`Z_0=\mathrm{exp}(N_{S=0})`$ is the partition function of all particles having zero strangeness and the argument of the Bessel function

$$x2\sqrt{S_1S_1}.$$
(10)
with $`S_{\pm 1}N_{s=\pm 1}`$. The parameter $`x`$ thus measures the total number of strange particles in thermal fireball.

The calculation of the particle density $`n_k`$ in the canonical formulation is straightforward. It amounts to the replacement

$$Z_k^1\lambda _kZ_k^1$$
(11)
of the corresponding one-particle partition function in equation (7) and taking the derivative of the canonical partition function equation (4) with respect to $`\lambda _k`$

$$n_k^C[\lambda _k\frac{}{\lambda _k}\mathrm{ln}Z_Q(\lambda _k)]_{\lambda _k=1}$$
(12)
As an example, we quote the result for the density of thermal kaons in the canonical formulation assuming that the total strangeness of the system $`S=0`$,

$$n_K^C=\frac{Z_K^1}{V}\frac{S_1}{\sqrt{S_1S_1}}\frac{I_1(x)}{I_0(x)}$$
(13)
Comparing the above formula with the result for thermal kaons density in the grand canonical ensemble, $`n_K^{GC}=(Z_K^1/V)\mathrm{exp}(\mu _s/T)`$, one can see that the canonical result can be obtained from the grand canonical one replacing the strangeness fugacity $`\lambda _S\mathrm{exp}(\mu _S/T)`$ in the following way:

$$n_K^C=n_K^{GC}\left(\lambda _S\frac{S_1}{\sqrt{S_1S_1}}\frac{I_1(x)}{I_0(x)}\right)$$
(14)
In the limit of large $`x`$ that is large volume and/or temperature the canonical and the grand canonical formulation are equivalent. For a small number of strange particles in a system, however, the differences are large. This can be seen in the most transparent way when comparing two limiting situations: the large and small $`x`$ limit of the above equation. In the limit $`x\mathrm{}`$ we have

$$\underset{x\mathrm{}}{lim}\frac{I_1(x)}{I_0(x)}1$$
(15)
and the kaon density is independent of the volume of the system as expected in the grand canonical ensemble. On the other hand in the limit of a small $`x`$ we have

$$\underset{x0}{lim}\frac{I_1(x)}{I_0(x)}\frac{x}{2}$$
(16)
and the particle density is linearly dependent on the volume. It is thus clear, that the major difference between the canonical and the grand canonical treatment of the conservation laws appears through different volume dependence of strange particle densities as well as strong suppression of thermal particle phase space. The relevant parameter, $`F_S`$, which measures the suppression of particle multiplicities from their grand canonical result is determined by the ratio of the Bessel functions

$$F_S\frac{I_1(x)}{I_0(x)}$$
(17)
with the argument $`x`$ defined in equation (10).

In Fig. 1 we show the canonical suppression factor $`F_S(x)`$ as a function of the argument $`x`$. To relate the initial volume of the system with the number of participant in A-A collisions one uses the approximate relation $`V1.9\pi A_{part}`$. The corresponding value of $`x`$ at SIS, AGS and SPS energy is calculated with the baryon chemical potential and temperature extracted from the measured particle multiplicity ratios . The results in Fig. 1 shows the importance of the canonical treatment of strangeness conservation at SIS energy. Here, the canonical suppression factor can be even larger than an order of magnitude.

For central Au-Au collisions at AGS or SPS energy this suppression is not relevant any more and the (GC)-formalism is adequate. In general, one expects, that the statistical interpretation of particle production in central heavy ion collisions requires the canonical treatment of strangeness conservation if the CMS collation energy $`\sqrt{s}<23`$GeV. This is mainly because at these energies the freeze-out temperature is still too low to maintain large argument expansion of the Bessel functions in equation (15). The canonical description of strangeness conservation can be, however, also of importance at the SPS energy were $`\sqrt{s}18`$GeV when one considers the peripheral heavy ion collisions. This is particularly true for multistrange particle production since the canonical suppression of the thermal particle phase-space increases with strangeness content of the particle .

### A Multistrange particle multiplicities

The extension of the canonical description to multi-strange particle multiplicities is straightforward. One needs first to extend the generating functional in Eq. 6 by including the contributions of multistrange baryons. In this case the canonical partition function constraint by the strangeness neutrality condition reads ,

$$Z_{S=0}^C=\frac{1}{2\pi }_\pi ^\pi 𝑑\varphi \mathrm{exp}(\underset{s=3}{\overset{3}{}}S_se^{is\varphi }),$$
(18)
where $`S_s_iZ_i^1`$ and the sum is taken over all particles and resonances carrying strangeness $`s`$. The one-particle partition function $`Z_i^1`$ is defined in Eq.(8).

With the particular form of the partition function given by Eq. (18) the density $`n_s`$ of particle $`i`$ with strangeness $`s`$ in volume $`V`$ is obtained by the replacement $`Z_i\lambda _iZ_i`$ in Eq. (18) and then taking an appropriate derivative :

$`n_{\pm s}`$ $`=`$ $`{\displaystyle \frac{N_{\pm s}}{V}}[{\displaystyle \frac{\lambda _i}{V}}{\displaystyle \frac{\mathrm{ln}Z_0}{\lambda _i}}]_{\lambda _i=1}Z_{\pm s}{\displaystyle \frac{(S_{\pm 1})^s}{(S_{+1}S_1)^{s/2}}}`$ (20)
$`\{I_s(x_1)I_0(x_2)+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}I_m(x_2)[I_{2m+|s|}(x_1)A^{m/2}+I_{2m|s|}(x_1)A^{m/2}]\}/Z_{S=0}`$
where

$$x_k2\sqrt{S_kS_k},A\frac{S_1^2S_2}{S_1^2S_2},$$
(21)
and the partition function

$$Z_{S=0}I_0(x_1)I_0(x_2)+\underset{m=1}{\overset{\mathrm{}}{}}I_{2m}(x_1)I_m(x_2)[A^{m/2}+A^{m/2}].$$
(22)
In the derivation of Eq. (20) we have neglected, after differentiation over particle fugacity, the term $`S_{\pm 3}`$. This approximation, however, due to small value of $`S_{\pm 3}`$ coefficients in comparison with $`S_{\pm 1}\mathrm{and}S_{\pm 2}`$ is quite satisfactory. The complete expression for the partition function without any approximation can be found in reference .

In the large system like Pb-Pb and for large collision energy, required to reach high $`T`$, the density $`n_s`$ of particle carrying strangeness $`s`$ is $`V`$ independent. In the opposite limit, however, this dependence is changed to $`n_sV^s`$, which can be verified from Eq. (20). Indeed for small $`x_i`$ we have approximately:

$$n_{\pm s}Z_{\pm s}\frac{(S_{\pm 1})^s}{(S_{+1}S_1)^{s/2}}\frac{I_s(x_1)}{I_0(x_1)}.$$
(23)
and when expanding the Bessel functions $`I_n(x)x^n`$ in Eq. (23) we see that $`n_sV^s`$.

From the above expression it is clear that strangeness suppression, which is measured by the ration $`I_s/I_0`$, is increasing with the strangeness content of the particle. Thus, there are two important ingredients of canonical modifications of multistrange particle density with respect to their (GC) value: (i) the density is volume dependent that is also centrality dependent and (ii) the thermal phase is suppressed and this suppression increases with strangeness content of the particle.

## III Time evolution and strangeness equilibration

In the last section we have formulated the statistical model for strange particle multiplicities $`<N_S>`$ assuming that the system is in thermal and chemical equilibrium. We have shown that dependently on the total number of strange particles there are two distinct equilibrium limits : if $`<N_s>`$ is small then we are in the canonical regime and in the opposite limit the canonical and grand canonical description coincides. In this section we consider the time evolution of the multiplicity of strange particles and formulate a kinetic master equations which distinguish between these two equilibrium limits. For a sake of illustration we consider a simple example of $`K^+K^{}`$ production in the environment of thermal pions in volume $`V`$ and temperature $`T`$ due to the following binary process, $`\pi ^+\pi ^{}K^+K^{}`$. We formulate for this example the kinetics for the time evolution of kaon multiplicities and their approach to chemical equilibrium .

In the standard formulation , the rate equation for this binary process is described by the following population equation:

$`{\displaystyle \frac{d<N_K>}{d\tau }}={\displaystyle \frac{G}{V}}<N_{\pi ^+}><N_\pi ^{}>{\displaystyle \frac{L}{V}}<N_{K^+}><N_K^{}>,`$ (24)
where $`G\sigma _Gv`$ and $`L\sigma _Lv`$ give the momentum-averaged cross sections for the gain $`\pi ^+\pi ^{}K^+K^{}`$ and the loss $`K^+K^{}\pi ^+\pi ^{}`$ process respectively. The value of $`<N_K>`$ represents the total number of produced kaons.

To include the possible correlations between the production of $`K^+`$ and $`K^{}`$ , let us define $`P_{i,j}`$ as the probability to find $`i`$ number of $`K^+`$ and $`j`$ number of $`K^{}`$ in an event. We also denote by $`P_i`$ as the probability to find $`i`$ number of $`K`$ pairs in an event. The average number of $`K`$ per event is defined as:

$`N_K={\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}iP_i.`$ (25)
We can now write the following general rate equation for the average kaon multiplicities:

$`{\displaystyle \frac{dN_K}{d\tau }}={\displaystyle \frac{G}{V}}N_{\pi ^+}N_\pi ^{}{\displaystyle \frac{L}{V}}{\displaystyle \underset{i,j}{}}ijP_{i,j}.`$ (26)
Due to the local conservation of quantum numbers, we have:

$`P_{i,j}`$ $`=`$ $`P_i\delta _{ij},`$ (27)
$`{\displaystyle \underset{i,j}{}}ijP_{i,j}`$ $`=`$ $`{\displaystyle \underset{i}{}}i^2P_iN^2=N^2+\delta N^2,`$ (28)
where $`\delta N^2`$ represents the event-by-event fluctuation of the number of $`K^+K^{}`$ pairs. Note that we always consider abundant $`\pi ^+`$ and $`\pi ^{}`$ so that we can neglect the number fluctuation of these particles and the change of their multiplicities due to the considered processes.

Following Eqs.(26-28) the general rate equation for the average number of $`K^+K^{}`$ pairs can be written as:

$`{\displaystyle \frac{dN_K}{d\tau }}={\displaystyle \frac{G}{V}}N_{\pi ^+}N_\pi ^{}{\displaystyle \frac{L}{V}}N_K^2.`$ (29)
For abundant production of $`K^+K^{}`$ pairs where $`N_K1`$,

$`N_K^2N_K^2,`$ (30)
and Eq.(29) obviously reduces to the standard form:

$`{\displaystyle \frac{dN_K}{d\tau }}{\displaystyle \frac{G}{V}}N_{\pi ^+}N_\pi ^{}{\displaystyle \frac{L}{V}}N_K^2.`$ (31)
However, for rare production of $`K^+K^{}`$ pairs where $`N_K1`$, the rate equations (24) and (31) are no longer valid. We have instead

$`N_K^2N_K,`$ (32)
which reduces Eq.(29) to the following form :

$`{\displaystyle \frac{dN_K}{d\tau }}{\displaystyle \frac{G}{V}}N_{\pi ^+}N_\pi ^{}{\displaystyle \frac{L}{V}}N_K.`$ (33)
Thus, in the limit where $`N_K1`$, the absorption term depends on the pair number only linearly, instead of quadratically for the limit of $`N_K1`$. Thus, it is clear that the time evolutions and equilibrium values for kaon multiplicities are obviously different in the above limiting situations.

In the limit of large $`<N_K>`$ the equilibrium value for the number of $`K^+K^{}`$ pairs, which coincides with the multiplicity of $`K^+`$ and $`K^{}`$, is obtained from Eq.(29) as,

$`<N_K>_{\mathrm{eq}}^{\mathrm{GC}}={\displaystyle \frac{V}{2\pi ^2}}m_K^2TK_2(M_K/T)`$ (34)
thus, it is described by the (GC) result with vanishing chemical potential due to strangeness neutrality condition.

In the opposite limit where $`N_K1`$, the time evolution is described by Eq.(33), which has the following equilibrium solution:

$`N_{\mathrm{eq}}^\mathrm{C}=\left[{\displaystyle \frac{V}{2\pi ^2}}M_{K^+}^2TK_2(M_{K^+}/T)\right]\left[{\displaystyle \frac{V}{2\pi ^2}}M_K^{}^2TK_2(M_K^{}/T)\right].`$ (35)
The above equation demonstrates the locality of strangeness conservation. With each $`K^+`$ the $`K^{}`$ is produced in the same event in order to conserve strangeness exactly and locally. This is the result expected from the (C) formulation of the conservation laws as described in the previous section. We note that Eq.(35) is just the leading term in the expansion of the canonical result for multiplicities of particles which are carrying $`U(1)`$ charges. The general expression is given by Eq.(13) and Eq.(20).

Comparing Eq.(32) and Eq.(33), we first find that, for $`N_K1`$, the equilibrium value in the canonical formulation is far smaller than what is expected from the grand canonical result as

$`<N_K>_{\mathrm{eq}}^\mathrm{C}=[<N_K>_{\mathrm{eq}}^{\mathrm{GC}}]^2<N_K>_{\mathrm{eq}}^{\mathrm{GC}}.`$ (36)
This shows the importance of the canonical description of quantum number conservation when the multiplicity of particles carrying non-zero $`U(1)`$ charges is small. We also note that the volume dependence in the two cases differs. The particle density in the GC limit is independent of $`V`$ whereas in the opposite canonical limit the density scales linearly with $`V`$. Secondly, we note that the relaxation time for a canonical system is far shorter than what is expected from the grand canonical result . It is also clear from Eq.(28) that (C) and (GC) limits are essentially determined by the size of $`\delta N^2`$, the event-by-event fluctuation of the number of $`K^+K^{}`$ pairs. The grand canonical results correspond to small fluctuations, i.e., $`\delta N^2/N^2<<1`$; while the canonical description is necessary in the opposite limit.

## IV Model predictions versus experimental data

In heavy ion collisions the number of produced strange particles depends on the collision energy and the centrality of these collisions. At low collision energies like eg. in GSI/SIS the freezeout temperature is relatively low being of the order of $`5080`$MeV. Consequently the number of strange particles in the final state is very small. Following our discussion in the last sections it is clear that the statistical description require here the canonical approach. In Fig. 2 we show the experimental data on $`K^+`$ yield per participant $`A_{part}`$ as a function of $`A_{part}`$ measured in Au-Au collisions at $`E_{lab}1`$ A/GeV .

The volume parameter in the statistical model scale with the number of participant. We can thus directly compare the model with data by fixing thermal parameters, the temperature and the baryon chemical potential, such that to reproduce the measured particle multiplicity ratios . In Fig. 2 the results of the canonical model is shown by the full line. The results clearly indicates that strong almost quadratic dependence of kaon yield on the number of participant is well reproduced by the model. The quadratic dependence of particle multiplicity in the canonical regime is the basic property of the model. In Fig. 3 we calculate the ratio $`\mathrm{\Xi }/K^+`$ multiplicities for central $`NiNi`$ collisions in the temperature range which corresponds to the collision energy of $`E_{lab}2`$ A/GeV.

The yield of $`\mathrm{\Xi }`$ is seen to be substantially smaller then the yield of $`K^+`$. This result is not only related with the differences in particle masses but particularly it appears due to the canonical suppression. Since $`\mathrm{\Xi }`$ carry strangeness minus two it has to be produced together with two $`K^+`$ to locally and exactly neutralized strangeness. It is also interesting to note that the $`\mathrm{\Xi }/K^+`$ ratio is independent on the baryon chemical potential. This is because $`K^+`$ appears together with $`\mathrm{\Lambda }`$, thus it contains the same dependence on $`\mu _B`$ as multistrange baryons.

The importance of the canonical treatment of strangeness conservation is also seen in higher collision energies like at the SPS when considering centrality dependance of multistrange baryons. In peripheral collisions the yield of strange particles is small such that also here the canonical description should be applied. The canonical suppression of thermal particle phase-space increases with strangeness content of the particle. The exact conservation of strangeness requires that each particle carrying strangeness $`\overline{s}`$ has to appear e.g. with $`s`$ other particles of strangeness one to satisfy strangeness neutrality condition.

In Fig. 4 we calculate the multiplicity/participant of $`\mathrm{\Omega },\mathrm{\Xi },`$ and $`\mathrm{\Lambda }`$ relative to its value in a small system with only two participants . Thermal parameters were assumed here to be $`A_{part}`$ independent. Fig. 4 shows that the statistical model in (C) ensemble reproduces the basic features of WA97 data : the enhancement pattern and enhancement saturation for large $`A_{part}`$ indicating here that (GC) limit is reached. Fig. 4 also demonstrate different $`A_{part}`$ dependence of strange and multistrange baryons. For small $`A_{part}`$ this dependence is power like as describe by Eq.(22). The quantitative comparison of the model with the experimental data would require an additional assumption on the variation of $`\mu _B`$ with centrality to account for larger value of $`\overline{B}/B`$ ratios in p+A than in Pb+Pb collisions . The most recent results of NA57 , showing an abrupt change of the enhancement for $`\overline{\mathrm{\Xi }}`$ are, however, very unlikely to be reproducible in terms of the canonical approach.

## V Summary and conclusions

We have discussed the importance of the conservation laws in the application of the statistical model to the description of strangeness production in heavy ion collisions. We have presented the arguments that the more general treatment of strangeness conservation based on the canonical ensemble is required if one compares the model with experimental data for particle yields obtained in central A-A collisions at SIS energies or peripheral collisions at the SPS. In both situations the number of produced strange particles per event is still too small to use the asymptotic grand canonical ensemble. The time evolution of strangeness production and the approach to chemical equilibrium limit was discussed in the context of a kinetic approach. We have shown on few examples that the statistical model predictions are consistent with the experimental data. A more complete presentation of the model versus data can be found in .

## VI Acknowledgments

On of us (K.R.) acknowledge stimulating discussions with Su Houng Lee and the partial support of the Committee for Scientific Research (KBN-2P03B 03018). We also acknowledge stimulating discussions with P. Braun-Munzinger, B. Friman, V. Koch, Z. Lin, M. Stephanov, H. Satz and X.N. Wang.


text:
# First-principles Calculation of the Formation Energy in MgO–CaO Solid Solutions

## I Introduction

In the last decade significant progress was achieved in application of first principles methods to calculate properties of materials. The full-potential linearized-augmented-plane-wave (FLAPW) and ultra-soft pseudopotential methods allow one to map phonon potential surface, describe phase transitions, electromechanical properties of ferroelectrics  and magnetic properties of materials . However these methods are applicable to ordered systems with translational symmetry, whereas many technologically relevant materials do not meet this requirement. Such important properties as ultrahigh piezoelectric efficiency, colossal magnetoresistance, high oxygen permeation through fuel cell membranes occur only when corresponding oxides form solid solutions, most of which are disordered.

First-principles methods developed for disordered systems are mostly based on the Green’s function formalism such as the Korringa-Kohn-Rostoker (KKR) method  that allows one to take into account atomic disorder within the average $`t`$-matrix (ATA)  or coherent potential (CPA)  approximations. They have been successfully applied to metallic alloys . A few calculations were performed for copper oxides  with mixed covalence – ionic character of chemical bonding. In the present paper, we consider the simplest oxide solid solution, the system MgO-CaO. Magnesiowüstite (Mg,Fe)O is believed to be a major constituent of the Earth’s lower mantle, so that it would be an ideal system for study, but the FeO constituent leads to complex Mott insulator behavior. The study of CaO in MgO can be considered as a step in the direction of the understanding of solid solution in minerals, as well as important in understanding the behavior of the minor element Ca. Understanding of oxide solid solutions is also important from the perspective of ferroelectric solid solutions and high temperature superconductors.

Several approaches have been used previously to study the MgO-CaO system. The electronic structure and total energy were calculated for constitutive compounds of the system, MgO and CaO using FLAPW method . A potential model  was developed to compute the phase diagram of the MgO-CaO solid solution. However the approach was found to be inadequate to reproduce the properties with required accuracy. Authors of Ref.  applied a tight-binding method to calculate the formation energy of ordered MgO-CaO solutions. Tight-binding parameters were obtained from results of first principles pseudopotential calculations performed for MgO, CaO and MgCaO2 compounds and then used for computing the total energy of solutions with different compositions and crystal structures. Comparison with the pseudopotential results suggests that the tight-binding approach is an efficient tool to study ionic oxides. However those were still ordered systems.

In the present paper, we considered both ordered and fully disordered MgO-CaO solutions using the multiple scattering theory within the local density approximation. We calculated the electron structure, total energy, equilibrium lattice parameters and formation energy for ordered MgO, CaO and MgCaO2 compounds as well as for 4 disordered MgO-CaO solutions with different compositions. The obtained formation energy values were compared to results of Ref. .

## II Computational Method

The calculations were performed by means a computer code based on the local density approximation (LDA)  and multiple scattering theory. The code embodies the local self-consistent multiple scattering (LSMS) method  where a compound is divided into overlapping clusters – local interaction zones (LIZs) centered around atoms of different sorts. The multiple scattering problem is solved for each LIZ separately in the lattice site – angular momentum representation for the muffin-tin (MT) potential given the cluster Green’s functions, the local densities of electron states and valence charge density for central atom.

Our approach is self-consistent with respect to local charge densities and potentials. However it does not treat the disorder self-consistency. We use an actual single site scattering matrix $`t_l^i(E)`$ for the central atom of LIZ and average $`t`$-matrix: $`\stackrel{~}{t_l^i}=xt_l^a+(1x)t_l^b`$ for surrounding sites, if they belong to a sublattice with a substitutional disorder of the $`a`$\- and $`b`$-sort atoms with concentrations $`x`$ and $`1x`$. According to the Ref. , such an approach is more accurate than the regular average $`t`$-matrix approximation. In the considered solid solution the Mg and Ca atoms randomly occupy the metal sublattice sites. They have the same number of valence electrons. This also allows us to believe that our approximation is reasonable for the system. One more advantage of the approach is it can be explicitly implemented within the ordinary LSMS theory .

Solving the Poisson equation at each iteration, we use the actual charge density for the local contribution to the MT potential and the average charges for the contribution of the substitutionally disordered sublattice to the Madelung potential. The exchange and correlation parts of the potential is determined within LDA using the technique described in Ref. .

The total energy was calculated following procedure described in Ref.  and using an expression for multicomponent compound .

In the present calculations, each LIZ built around Mg, Ca and O sites contained 123 atoms. This provided good convergence over the LIZ size, since the difference in the total energies obtained for 123 and 93 atom LIZs was 3 and 5 mRy for MgO and CaO respectively. To estimate an error induced by MT potentials, the calculations were performed in two approximations – maximum MT radii proportional to the unit cell volume and fixed MT radii corresponding to a minimum volume considered for a given compound.

## III Results and Discussion

The density of electron states, charge density and total energy were calculated versus the unit cell volume for ordered compounds MgO, CaO and MgCaO2 and disordered solid solutions Mg1-xCaxO with $`x=0.2,0.5,0.7,0.8`$. For MgCaO2 the structure was chosen where cations were ordered by an alternate stacking of Mg and Ca planes along the  direction. Thus the cations form the Ll0 structure in Strukturbericht notation .

Based on the computed volume dependence of the total energy the equilibrium volumes were calculated for each considered compound. The results are shown in Fig. 1. The concentration dependence of unit cell volume is found to be close to linear one, though some wavy behavior is seen. The difference in results obtained by means of maximum and fixed MT radii approaches gives an order of magnitude for the error induced by muffin-tin approximation. For the MgO and CaO compounds the calculated equilibrium volumes are less than experimental ones that reflects the well known feature of LDA.

Using the calculated values of total energy we have determined the formation energy of the solutions determined as:

$$E_f=E(Mg_{1x}Ca_xO)[(1x)E(MgO)+xE(CaO)],$$

where E is the total energy of corresponding compound. The $`E_f`$ values compared to results of pseudopotential calculations  are shown in Fig. 2. Our result and result of Ref.  obtained for the ordered Ll0 phase of MgCaO2 are in a good quantitative agreement (0.285 eV and 0.282 eV respectively). The difference in formation energy of the ordered and disordered MgCaO2 is less than 10%. Such a small difference takes place since a decay of band states caused by random potentials mostly involves unoccupied states in these ionic compounds and is much lower than in case of materials with covalence bonding . The formation energy is found to be positive for all considered solid solutions suggesting that a phase separation is preferable for this system. This is in agreement with the experimental phase diagram , which shows a miscibility gap.

Charge transfer is one of the key mechanisms determining the total energy of complex oxides. Therefore, we have calculated effective charges on atoms in the considered solid solutions by integration of the valence electron density over Wigner-Seitz spheres. It is important at this stage to make a clear reasonable definition of the space belonging each nonequivalent atom that is determined by the ratio of the Wigner-Seitz radii $`r_O^{WS}/r_M^{WS}`$ (subscripts $`O`$ and $`M`$ denote oxygen and metal respectively). We suppose that in the present case it is convenient for interpretation purpose to keep this ratio independent of the metal composition. This ratio should also make some physical sense. Because the considered solutions are definitely ionic, the ratio should be related to the ionic radii. Thus we come to the Wigner-Seitz radii ratio defined as

$$r_O^{WS}/r_M^{WS}=2R_i^{OII}/[R_i^{MgII}+R_i^{CaII}]=1.24,$$

where $`R_i^{OII},R_i^{MgII}`$ and $`R_i^{CaII}`$ are ionic radii of O, Mg and Ca respectively. We have also determined the average cation charges in the solutions as $`Q_{av}=(1x)Q_{Mg}+xQ_{Ca}`$. The results are shown in Fig. 3. Since the Ca atom has the valence wave functions more extended than Mg, substitution of Mg with Ca leads to a noticeable increase in the electron charge on oxygen atoms (more than 0.2e, going from MgO to CaO). The concentration dependence of these charges as well as average cation charges are found to be linear.

In summary, we have calculated the electron structure, total energy, equilibrium lattice constants and formation energy for ordered and disordered solid solutions Mg1-xCaxO. A linear composition dependence has been found for lattice parameters and effective charges. The formation energy is positive for all considered materials that is in agreement with the experimental phase diagram. The results suggest that LSMS method and LDA can be an efficient tool to study properties of disordered ionic solutions.

Acknowledgements. This work was supported by Office of Naval Research grant N00014-97-1-0052. Computations were performed on the Cray SV1 at the Geophysical Laboratory, supported by NSF grant EAR-9975753 and the W. M. Keck Foundation.


text:
# Substructure of jets at HERA

## 1 Introduction

The internal structure of jets gives insight into the transition between partons produced in the hard scattering process and the experimentally observable jets of hadrons. The substructure of jets is expected to depend mainly on the type of primary parton, quark or gluon, and to a lesser extent on the particular hard scattering process. QCD predicts that at high transverse energy of the jet, $`E_T^{jet}`$, where fragmentation effects become negligible, the jet structure is driven by gluon emission from the primary parton and gluon jets are broader than quark jets due to the larger colour charge of the gluon. The substructure of jets has been studied in terms of the integrated jet shape , $`\mathrm{\Psi }(r)`$, and the differential jet shape , $`\rho (r)`$, where the energy flow inside a jet is considered, and the subjet multiplicity , $`n_{sbj}`$, where jet-like structures (subjets) within a given jet are studied. The data samples used in these analyses were collected with the ZEUS detector at HERA. During 1995-1997 (1999-2000) HERA operated with positrons of energy $`E_e=27.5`$ GeV colliding with protons of energy $`E_p=820`$ GeV ($`E_p=920`$ GeV).

## 2 Jet substructure in exclusive and charm-induced dijets in $`\gamma p`$

Jet shapes and subjet multiplicities in exclusive and charm-induced dijet samples in photoproduction (PHP) for the kinematic region $`Q^21`$ GeV2 and inelasticity $`0.2<y<0.85`$ have been measured. Dijets events were identified with the $`k_T`$-cluster algorithm  and selected with $`E_T^{jet_{1,2}}7,6`$ GeV and $`1<\eta ^{jet}<2`$. Charm quarks were tagged by identifying $`D^\pm `$ mesons through the $`K2\pi `$ decay mode using the $`\mathrm{\Delta }M`$ method . The jet with closest distance in azimuthal angle to the $`D^{}`$ meson was associated with the charm meson. The jet not associated with the $`D^{}`$, referred to as the “untagged-charm” jet, represents the unbiased, i.e. not influenced by the $`D^{}`$ selection, jet candidate whose internal properties are studied.

The measured $`n_{sbj}`$ as a function of the resolution scale, $`y_{cut}`$, for the exclusive dijet sample is shown in Fig. 1a. The predicted $`n_{sbj}`$ is larger for gluon-initiated jets than for quark-initiated jets and the measured data are located between the two curves, showing that the dijet sample is a mixture of quark and gluon jets. PYTHIA  gives a good description of the data. Fig. 1b shows $`n_{sbj}`$ as a function of $`y_{cut}`$ for the “untagged-charm” jet. The agreement between data and theory is very good and the predictions of charm-initiated jets are consistent with the measurements.

The dependence of the substructure of jets with $`\eta ^{jet}`$ has also been studied. In the exclusive dijet sample, $`n_{sbj}`$ increases with $`\eta ^{jet}`$ (Fig. 2a), which is consistent with the predicted increase in the fraction of gluon-initiated jets with $`\eta ^{jet}`$. In the “untagged-charm” jet sample, the results are consistent with a pure sample of quark jets for $`1<\eta ^{jet}<1.5`$. For the highest $`\eta ^{jet}`$ values, the data show a deviation from the prediction for quark-induced jets. Since the estimated gluon contamination to the charm-induced jet sample (mainly due to “charm excitation”) has its highest contribution in the forward region, the deviation could be explained by the increase of the gluon-jet fraction in the charm-enriched sample.

Since the dijet PHP sample consists of a mixture of quark and gluon jets, any measured observable $`𝒪`$ of the internal structure can be written as: $`𝒪_{dijet}=f_q𝒪_{quark}+f_g𝒪_{gluon}`$, where $`f_q`$, $`f_g=1f_q`$ are the fractions of quark and gluon jets. The measurements of the substructure of the charm-enriched sample at high transverse energies ($`E_T^{jet}>15`$ GeV) can be considered as measurements for a pure sample of quark jets ($`𝒪_{charm}=𝒪_{quark}`$). Taking the fractions $`f_q`$, $`f_g`$ from LO Monte Carlo, the substructure of gluon jets can be extracted. The fraction $`f_q`$ predicted for a dijet sample with $`E_T^{jet}>15`$ GeV and $`1<\eta ^{jet}<2`$ by PYTHIA ($`f_q=0.66`$) and HERWIG ($`f_q=0.69`$) are found to be similar. Fig. 2b shows the extracted gluon jet substructure, which is consistent with the QCD predictions.

## 3 Jet substructure in deep inelastic scattering

Quark initiated jets are expected to be predominant in charged current (CC) and neutral current (NC) deep inelastic scattering (DIS). In Fig. 3a, $`n_{sbj}`$ as a function of $`y_{cut}`$ in the “untagged-charm jet” PHP sample and an inclusive jet NC sample ($`Q^2>125`$ GeV2) are compared. The measurements are found to be similar and in very good agreement for resolution scales $`y_{cut}>0.01`$, where the charm mass effects are negligible. The charm-initiated jets are very similar to light quark jets. Fig. 3b shows the evolution of $`n_{sbj}`$ with $`E_T^{jet}`$ for a fixed $`y_{cut}=0.01`$ in NC and CC processes. The value of $`n_{sbj}`$ decreases as $`E_T^{jet}`$ increases. The agreement between both measurements indicates that the pattern of parton radiation within quark jets is to a large extent independent of the hard scattering process. In Fig. 3c $`n_{sbj}`$ for a fixed $`y_{cut}=0.01`$ as a function of $`\eta ^{jet}`$ in CC interactions ($`Q^2>200`$ GeV2) is shown. Both ARIADNE (CDM)  and LEPTO (MEPS)  give a good description of the data. The substructure in DIS processes shows no dependence with $`\eta ^{jet}`$.

In NC processes, the subjet multiplicity (Fig. 4a) as a function of $`E_T^{jet}`$ for a fixed $`y_{cut}=0.01`$ and the mean integrated jet shape (Fig. 4b) for a fixed radius of $`r=0.5`$ have been compared to NLO QCD predictions and $`\alpha _S`$ has been determined. The extracted value from $`n_{sbj}`$ for jets with $`E_T^{jet}>25`$ GeV is $`\alpha _S(M_Z)=0.1185\pm 0.0016(stat)_{0.0048}^{+0.0067}(exp)_{0.0071}^{+0.0089}(th)`$.
The extracted value from $`\mathrm{\Psi }(r)`$ for jets with $`E_T^{jet}>21`$ GeV is
$`\alpha _S(M_Z)=0.1179\pm 0.0017(stat)_{0.0065}^{+0.0054}(exp)_{0.0073}^{+0.0094}(th)`$.
Both values are in good agreement with the world average value .