text:
# A Systematic Analysis of Supernova Light in Gamma-Ray Burst Afterglows

## 1 Introduction

Observational and theoretical evidence suggest that the majority of GRB progenitors are stars at the endpoint in stellar evolution (e.g., Fryer, Woosley, & Hartmann 1999; Heger et al. 2003). Since the discovery of a near-by Type Ic supernova (SN 1998bw) in the error circle of the X-ray afterglow for GRB 980425 (Galama et al. 1998; Kulkarni et al. 1998), evidence is accumulating that core collapse supernovae are physically related to long-duration GRBs. The supernova picture is further supported by the observation that all GRB hosts are star-forming, and in some cases even star-bursting galaxies (e.g., Frail et al. 2002; Sokolov et al. 2001). Evidence for host extinction by cosmic dust in GRB afterglows (e.g., Castro-Tirado et al. 1999; Klose et al. 2000), and the discovery of an ensemble of optically ’dark bursts’ (for a recent discussion, see Fynbo et al. 2001; Klose et al. 2003; Lazzati, Covino, & Ghisellini 2002) also is consistent with the picture that GRB progenitors are young, massive stars (Groot et al. 1998; Paczyński 1998). Furthermore, for several GRB afterglows X-ray data suggest a period of nucleosynthesis preceding or accompanying the burst (e.g., Antonelli et al. 2000; Lazzati, Campana, & Ghisellini 1999; Mészáros & Rees 2001). The angular distribution of the afterglows with respect to their hosts also favors a physical relation of young, massive stars to GRBs (Bloom, Kulkarni, & Djorgovski 2002).

As a natural consequence of a physical relation between the explosion of massive stars and GRBs supernova light should contribute to the afterglow flux, and even dominate under favorable conditions. The most convincing example is GRB 030329 (Peterson & Price 2003) at $`z`$=0.1685 (Greiner et al. 2003a) with spectral confirmation of supernova light in its afterglow (Hjorth et al. 2003; Kawabata et al. 2003; Matheson et al. 2003; Stanek et al. 2003). In contrast to this unique spectroscopic evidence, several cases of photometric evidence for extra light in GRB afterglows have been reported, starting with GRB 980326 (Bloom et al. 1999). Various groups successfully fit SN 1998bw templates to explain these late-time bumps, the most convincing case being that of GRB 011121 (Bloom et al. 2002; Garnavich et al. 2003; Greiner et al. 2003b).

The goal of the present paper is to analyze the supernova bumps in GRB afterglow light curves using a systematic approach. While this was done also for several bursts by Dado, Dar, & de Rújula (2002a; and references therein) with respect to a verification of their cannonball model (Dar & de Rújula 2003), we tackle this issue in an independent and different way. First, from the numerical site, we have developed our own computational approach. This includes a numerical procedure to redshift SN 1998bw light curves (see § A.2) and to fit afterglow data within the context of the fireball model. Second, from the observational site, when necessary and possible we have performed late-time observations of some GRB host galaxies (§ 2). A considerable part of the data we have included in our study is based on observing runs where we have been involved. Additional data have been collected from the literature and checked for photometric consistency. Nearly two dozen afterglows had sufficient data quality, and a known redshift, to search for a late-time bump in the light curve (§ 3). Third, we concentrate our attention on a statistical analysis of the phenomenological parameters for this class of GRB-SNe (§ 4). In this respect our investigation goes beyond the approaches undertaken by others to explain late-time bumps in individual afterglow light curves.

## 2 Observations and data processing

Some of the GRB afterglows we analyzed had poorly sampled late-time data, which made it difficult to find or determine the parameters of a SN bump (GRBs 000418, 991208, and 010921). In order to perform late-time photometry of these GRB hosts, we carried out two observing runs at the Calar Alto 3.5-m telescope on March 13–14 and May 23–25, 2003. Observations were performed using the multi-purpose camera MOSCA, which uses a SITe 2k$`\times 4`$k CCD with a plate scale of 0.32 arcsec per pixel. The field of view is $`11\times \mathrm{\hspace{0.17em}11}`$ arcmin2. During the first observing run the seeing varied between 1.4 and 1.6 arcsec; in May the seeing was better than 0.8 arcsec. Data reduction followed standard procedures.

Most of the light curves we investigated have been followed in more than one photometric band. For each of the GRBs we chose the best-sampled light curve as a reference light curve for the fit in the other photometric bands. In the most cases this was the $`R`$ band light curve. In brief, our approach was as follows. First, we assumed that the afterglow slopes and break time are the same in all filters (Eq. A2), in reasonable agreement with observational data. Consequently, once we fit the reference light curve of an optical transient and deduced the corresponding afterglow parameters, we treated them as fixed parameters when fitting the light curves of the optical transient in other photometric bands. Thereby the fit procedure was based on a $`\chi ^2`$-minimization with a Levenberg - Marquardt iteration. Second, for the representation of the supernova component we used published $`UBVRI`$ data of the light curve of SN 1998bw (Galama et al. 1998) as a template, and redshifted them to the corresponding cosmological distance of the burster (§ A.2). These light curves are different from band to band. In addition, we allowed a variation of the SN luminosity with respect to SN 1998bw and a stretch in time (Eq. A1 in § A.1). Third, the host magnitude, which represents a constant component in the integral light of the optical transient, was usually considered as a free parameter. Only for GRB 011121 and 020405 we used host-subtracted magnitudes to fit the light curves.

Before performing a numerical fit, the observational data was corrected for Galactic extinction along the line of sight using the maps of Schlegel, Finkbeiner, & Davis (1998). This also holds for SN 1998bw, where we assumed $`E(BV)`$= 0.06 mag. We calculated the Galactic visual extinction according to $`A_V^{\mathrm{Gal}}`$ = 3.1 $`E(BV)`$, whereas the extinction in $`U`$ and $`B`$ were obtained via Rieke & Lebofsky (1985), and in $`R_c`$ and $`I_c`$ by means of the numerical functions compiled by Reichart (2001).

## 3 Results

Among the 36$`\pm 1`$ GRBs with detected optical afterglows up to the end of 2002111http://www.mpe.mpg.de/jcg/grbgen.html, 21 had sufficient data quality, and a known redshift, for a meaningful search for an underlying supernova component (Table 1). Among them in nine cases evidence for a late-time bump was found. The results are summarized in Table 2. A general inspection of this table makes clear that the burst ensemble with detected late-time bumps in their afterglows separates into a group with statistically significant evidence for a bump (990712, 991208, 011121, and 020405), mostly in at least two photometric bands, and a group with less significant bumps (970228, 980703, 000911, 010921, and 021211). Given the fact that evidence for these bumps has also been reported by other groups (with the only exception being GRB 010921), we feel confident that the results presented in Table 2 can be used for a first statistical approach to understand this type of GRB-SNe.

The most interesting result is that our numerical procedure found evidence for a late-time bump in all GRB afterglows with a measured redshift $`z\text{ }\stackrel{<}{}\text{ }`$0.7. We believe that the interpretation of these bumps as an underlying supernova component is the most natural and observationally most founded explanation. Among the higher redshifted bursts, we confirm the finding by Holland et al. (2001) of a possible bump in the afterglow of GRB 980703, the discovery by Lazzati et al. (2001) of a bump in the afterglow of GRB 000911 ($`z`$=1.06; Price et al. 2002b), and a bump in the afterglow of GRB 021211, which was discovered by Della Valle et al. (2003) and is also discussed by Dado, Dar, & de Rújula (2003b).

For five afterglows (GRB 970508, 991216, 000418, 010222, 020813) with $`0.7<z<1.5`$, we can only place upper limits on the luminosity of any underlying supernova component. These five upper limits have typical uncertainties of a factor of two. The only exception is GRB 020813 where the uncertainty is much larger, so that no firm conclusions can be drawn here. The remaining seven bursts in our sample (GRB 971214, 990123, 990510, 000301, 000926, 011211, 021004) have redshifts $`z>1.5`$ and, therefore, have not been investigated here, since this would have required a substantial extrapolation of the SN 1998bw template into the UV domain. Finally, the late-time bump clearly seen in the afterglow of GRB 980326 (Bloom et al. 1999) is not included in our study, because the redshift of the burster is not yet known.

As we have outlined in the previous section, we fit the SN component using the light curves of SN 1998bw as a template. Thereby, we allowed the luminosity and the light curve shape to be a free parameter. The former means a scaling of the luminosity of SN 1998bw by a factor $`k`$ (Eq. A1), whereby $`k`$ always refers to the corresponding wavelength region in the redshifted SN frame. Differences in the light curve shapes were modeled by means of a stretch factor $`s`$, which allows the supernova light curve to develop slower ($`s>1`$) or faster ($`s<1`$) than the one of SN 1998bw (Eq. A1). In this respect we follow Hjorth et al. (2003) to describe the light curve of GRB 030329/SN 2003dh. The advantage of such an approach is that we can use these two parameters to explore the entire ensemble of GRB-SNe in a statistical sense. Table 2 shows that for the bursts with the photometrically best sampled late-time bumps in their optical light curves (GRB 011121, 020405) the deduced parameters $`k`$ and $`s`$ are consistent with each other in different wavelength regions. A comparison of these parameters of the nine afterglows with late-time bumps, which are at different redshifts, seems to be a reasonable first approach in order to constrain the width of the luminosity distribution of GRB-SNe.

In Fig. 1 we display the deduced parameter $`k`$ (luminosity) for every individual GRB-SNe. We plot luminosity versus redshift just to look for a potential evolutionary effect (which is not apparent) and to separate the individual SNe from each other. While in Fig. 1a we allowed the stretch factor $`s`$ to be a free parameter during the numerical fit, Fig. 1b shows the results obtained when we fixed $`s=1`$. The reason for the latter was twofold. First, if $`s=1`$ we can constrain the luminosity of an underlying SN for those GRBs, where we do not detect a bump in the late-time light curve. This is not possible in a reasonable way, if we allow $`s`$ to be a free parameter. Second, sometimes the data base is too small in order to include also the stretch parameter in the fitting procedure, so that we have to fix $`s`$. Note that in Fig. 1 there are three small sets of bursts at redshifts 0.4, 0.7 and 1.0. Between them is no difference apparent neither in the luminosities of the GRB-SNe nor in the width of the luminosity distribution. What is apparent from a comparison of Fig. 1a,b is that introducing the stretch factor reduces the width of the luminosity distribution of the GRB-SNe and brings the luminosity of all SNe a little closer to those of SN 1998bw ($`k=1`$). The distribution of the deduced stretch factor itself is shown in Fig. 2. Although $`s`$ varies by a factor of two in both directions around $`s=1`$, within their individual 1$`\sigma `$ error bars most data are close to $`s=1`$. Finally, no correlation was found between the deduced SN luminosity (parameter $`k`$) and the stretch factor $`s`$ (Fig. 3).

## 4 Discussion

### 4.1 The luminosity distribution of the GRB-SNe

When plotting the parameter $`k`$ deduced for the $`R`$ band in the observer frame, the width of the luminosity distribution of the class of GRB-SNe (Fig. 4) is similar to what is observed for other classes of core collapse SNe (Richardson et al. 2002). The mean of $`k`$ in the $`R`$ band is 0.7 independent of whether or not we fix the stretch parameter at $`s=1`$, while our template SN 1998bw is at the bright end of the GRB-SN luminosity distribution. The latter conclusion is supported when we plot $`k`$ for the same wavelength region in the SN host frame, which is a better indicator of the spread in SN luminosities. Most of the GRB-SNe we explored have a data point for the photometric band centered around 395$`\pm `$10 nm in the corresponding host frame (Table 2). The distribution of $`k`$ now indicates again that SN 1998bw is among the most luminous GRB-SNe. It also indicates that in fact there is no peak around $`k`$=0.8 (Fig. 4), but we may so far have only sampled the bright part of the GRB-SN luminosity function. Some caution is of course required, given the partly large error bars of the $`k`$ factors, which are not shown in Figs. 4 and 5. While the conclusion that SN 1998bw was among the most luminous members of its class seems to be robust, the shape of the GRB-SNe luminosity function is still less well-determined. Extinction by interstellar dust in the host galaxies could in principle also affect these results, although only for GRB 010921 was a significant host extinction ( $`\stackrel{>}{}`$ 1 mag) reported (Price et al. 2003).

As we have outlined before, for redshifts $`z`$ $`\stackrel{<}{}`$ 0.7 all GRB afterglows show evidence for an underlying late-time bump. Within our context this means that we trace a complete set of GRB-SNe, i.e., not only the brightest members of this class. The width of this GRB-SNe luminosity distribution in the photometric band centered around 395$`\pm `$10 nm in the SN frame is $``$1 to 1.5 magnitudes. This wavelength region is roughly placed in the $`B`$ band, so that we can compare the corresponding luminosities with those of other Type Ib/c supernovae (Richardson et al. 2002), i.e., those class of SNe, which is believed to include the progenitors of GRBs. It turns out that the GRB-SNe do fit into a region between approximately $`M_B=19.5`$ and $`M_B=18`$ in figure 6 of Richardson et al., where no data on Type Ib/c SNe are known. If all GRB-SNe are indeed of type Ib/c this would favor the conclusion that the luminosity function of Type Ib/c SNe is rather described by a broad Gaussian than by a bimodal distribution (figures 6 and 7 in Richardson et al.).

The non-detection of a supernova bump in more than half of the investigated GRB afterglow light curves may be accounted for by several reasons, like a relatively bright host, or a faint supernova. In particular, finding a SN bump for high-$`z`$ bursts is an observational challenge. For $`z\text{ }\stackrel{>}{}\text{ }0.7`$ and $`k=1`$, this peak magnitude exceeds $`R_c=24`$, which poses a major challenge for 3-m class telescopes, given the usually very limited amount of target of opportunity time for such observations. It is thus no surprise that only for three of the GRBs above $`z=0.7`$ a supernova component was found (GRB 980703, GRB 000911, GRB 021211), even though we can not rule out that the SN ’drop-out’ is due to some evolutionary effect of the underlying burster population and their environment.

### 4.2 The SupraNova model

Vietri & Stella (2000) argued that GRBs are the result of delayed black hole formation, which implies that the core collapse and its subsequent supernova may significantly precede the burst. The delay could be of order months to years (Vietri & Stella 2000), or perhaps as short as hours (Woosley, Zhang, & Heger 2002). While constraining the latter possibility can not be accomplished with the data at hand, the longer time scales are easily constrained. For only two of the SN light curves the fit indeed improved if we allowed for a shift in time between the onset of the burst and the onset of the SN (GRB 990712, 011121). The offsets never exceeded 5 days, and were both negative and positive. However, the uncertainties in this parameter are large, due to the poorly sampled shape of the underlying supernova (e.g., Garnavich et al. 2003). The average shift is basically consistent with zero. Presumably, these shifts are due to an underlying correlation between luminosity and light curve shape, as observed in other types of supernovae (e.g., Candia et al. 2003; Stritzinger et al. 2002). This is just what the parameter $`s`$ takes into account. On the other hand, its clear that we have no information about this issue in those cases where we have not found evidence for a SN bump at late times. While this still leads open the possibility of two populations of bursters (collapsars & SupraNovae), we emphasize again that we find a late time bump in all afterglows with redshift $`z<0.7`$.

### 4.3 X-ray lines and supernova bumps

The identification of late-time bumps in afterglow light curves with SN light would benefit from observations in the X-ray band (for the cannonball model, e.g., Dado, Dar, & de Rújula 2003a). If the X-ray lines seen in some afterglows (e.g., Reeves et al. 2002) have their origin in the circumburster medium (e.g., Lazzati et al. 2001) and not in the exploding star (e.g., Mészáros & Rees 2001) this would be difficult to reconcile with the interpretation of a late-time bump with an underlying SN component. Unfortunately, the majority of bursts with high-resolution XMM-Newton or Chandra spectroscopic X-ray follow-up observations have no well-observed optical light curves. Among the afterglows with a detected optical late-time bump listed in Table 2 only for GRB 020405 such observations exist (Mirabal, Paerels, & Halpern 2003); no evidence for X-ray lines has been found there. BeppoSAX observed the afterglow of GRB 970228 (Frontera et al. 1998) with comparable low energy resolution, and no X-ray lines have been reported. Among those low-$`z`$ bursts with well-defined optical light curves and no evidence for a late-time bump in the data, only for GRB 970508 BeppoSAX X-ray follow-up observations have been published (Piro et al. 1999). Evidence for an iron line was found. Although one might add to the list of bursts with well-observed late-time light curves and additional spectral information in the X-ray band GRBs 990123, 000926, 010222, 011211, and 021004, the redshift of these bursts was  $`\stackrel{>}{}`$ 1.5, making it more or less hopeless to find an underlying SN component in the available data base (an upper limit for GRB 010222 is reported in Fig. 1b; Henden et al. 2004, in preparation). While it is very interesting that neither for GRB 020405 (Mirabal et al. 2003) nor for GRB 030329 (Tiengo et al. 2003) lines have been found in X-ray spectra of their afterglows, at present we cannot confirm a possible anti-correlation between the occurrence of X-ray lines and the appearance of SN light in GRB afterglows. This important issue remains to be investigated in the Swift era.

### 4.4 SN properties vs. afterglow parameters

Of particular interest is whether the burst and afterglow properties are to some degree related to the existence of an underlying SN component. For this reason we have investigated if the deduced SN luminosity is correlated with the corresponding energy release in the gamma-ray band (as given in Bloom, Frail, & Kulkarni 2003). No such correlation was found. We have also checked whether the afterglow parameters $`\alpha _1`$ (pre-break decay slope), $`\alpha _2`$ (post-break decay slope), and $`t_b`$ (break time; Eq. A2) from those GRBs with detected SN component are different from those without such a component. Again, no correlation was apparent, even though one should keep in mind that the data base is still very small. For illustration, in Fig. 6 we display the relation between the afterglow parameters $`\alpha _1`$ and $`\alpha _2`$ for all GRB afterglows we have investigated. While there is a ’forbidden region’ with $`\alpha _2\text{ }\stackrel{<}{}\text{ }\alpha _1+3/4`$ apparent in the data, with the border line representing the simplest jet model with no sideways expansion of the jet (Mészáros & Rees 1999), no bimodality in this distribution is visible. A tendency in our data that GRB afterglows with detected underlying SN component seem to prefer pre-break slopes $`\alpha _1>1`$ should not be overinterpreted, since the $`\alpha _1<1`$ sample includes several bursts with redshifts $`z>1.2`$, where the photometric detection of an underlying SN component is difficult. On the other hand, at least one selection effect does occur here. If a bright SN component is apparent in the data then the parameter $`\alpha _2`$ of the afterglow light curve is usually much more difficult to quantify, because the late-time evolution of the genuine afterglow is less well-defined. This problem is well seen in the afterglow light curve of GRB 011121 (e.g., Greiner et al. 2003) and GRB 030329 (e.g., Lipkin et al. 2004).

## 5 Summary

In an attempt to study the GRB-SN association, we have re-analyzed in a systematic way all GRB afterglow data published by the end of 2002. We have found that in nine cases evidence for extra light at late times is apparent in the optical afterglow light curves. In most cases this is seen in more than one photometric band. This extra light can be modeled well as supernova light peaking $`(1+z)15\mathrm{}20`$ days after a burst. Our main finding is that all GRB afterglows with redshift $`z\text{ }\stackrel{<}{}\text{ }0.7`$ showed evidence for extra light at later times. For larger redshifts the data base is usually not of sufficient quality, or the SN is simply too faint, in order to search for such a feature in the late-time afterglow light curve.

The cut off date of our sample (end of 2002) was chosen to ensure that all GRBs had published follow-up observations. Since that date, five new afterglows with redshifts have been established. All but one were at redshifts above 0.7, and again for the only nearby event (GRB 030329) a supernova component was established. This is consistent with the statistical inferences from the sample of earlier long-duration GRBs and leads us to conclude that the current world sample of GRB afterglow measurements provides strong statistical support for the link between (long-duration) GRBs and the final stages of massive star evolution. This conclusion basically agrees with earlier reports by Dado, Dar, & de Rújula (2002a; and references therein), and is essentially independent of the underlying GRB model. While so far only one event (GRB 030329) allowed a direct spectroscopic confirmation of this link, the larger photometric sample discussed here supports this idea by statistical means.

Based on our sample of nine GRB-supernovae we have performed a first statistical approach to get insight on the characteristic luminosities of this type of supernovae. We have found strong evidence that SN 1998bw is at the bright end of the GRB-SN luminosity distribution, with the latter matching well into what is known so far about the luminosities of the brightest members of other types of core collapse supernovae (Richardson et al. 2002). While GRB-SNe are not standard candles, their peak luminosities are comparable to those of Type Ia. In fact, within the context of the SN interpretation of the late-time bumps in afterglows, our results demonstrate once more that the first years in GRB research have already provided a first sample of high-$`z`$ core collapse SNe up to a redshift of 1. This sample might grow rapidly in the near future if indeed all long-duration GRBs tell us when and where in the universe a massive star explodes.

Some caution is of course required. First, there is some bias in the sample of bursts with detected optical afterglows. None of the bursts with a detected SN bump was classified as an X-ray rich burst or an X-ray flash, and for none of them were X-ray lines reported in the literature. In other words, it is still possible that bursts with SN bumps do not belong to these classes of events (but see Fynbo et al. 2004). Second, for most of the bursts discussed here evidence for a SN bump is based on a very small number of data points around the SN peak time (say, between 10 to 40 days after the burst), with the most critical cases being GRB 991208 and 010921. However, we see no reason why we should disregard these events.

In their discovery paper, Klebesadel, Strong, & Olson (1973) noted that a potential relation of GRBs to supernovae might still be an option to explain this new phenomenon. While the model they refer to (Colgate 1968) does not describe what is today believed to be the underlying GRB mechanism, historically it is nevertheless remarkable that the first paper ever about GRBs might have given the right hint on the underlying source population, followed by many years of trial and error.

S.K. and A.Z. acknowledge financial support by DFG grant KL 766/12-1 and from the German Academic Exchange Service (DAAD) under grant No. D/0103745. D.H.H. acknowledges support for this project under NSF grant INT-0128882. A.Z. acknowledges the receipt of a scholarship from the Friedrich-Schiller-University Jena, Germany. We thank Nicola Masetti and Eliana Palazzi, Bologna, for providing host-subtracted photometric data on the afterglow of GRB 020405. This work has profited from the GCN data base maintained by Scott Barthelmy at NASA and the GRB big table maintained by Jochen Greiner, Garching. We thank Kevin Lindsay (Clemson) for a careful review of the manuscript. We thank the referee for critical comments, which helped to improve the paper.

## Appendix A Numerical approach

### A.1 The light curve of the optical transient

We model the light curve of the optical transient (OT) following a GRB as a composite of afterglow (AG) light, supernova (SN) light, and constant light from the underlying host galaxy. The flux density, $`F_\nu `$, at a frequency $`\nu `$ is then given by

$$F_\nu ^{\mathrm{OT}}(t)=F_\nu ^{\mathrm{AG}}(t)+kF_\nu ^{\mathrm{SN}}(t/s)+F_\nu ^{\mathrm{host}}.$$
(A1)
Here, the parameter $`k`$ describes the observed brightness ratio (in the host frame) between the GRB-supernova, and the SN template (SN 1998bw) in the considered photometric band (in the observer frame). We allowed $`k`$ to be different in every photometric band, but within a band independent of frequency. The parameter $`s`$ is a stretch factor with respect to the used template. We have also explored the consequences of a shift in time between the onset of the burst and the onset of the supernova explosion, as implied by certain theoretical models (Vietri & Stella 2000). Then, in Eq. (A1) $`F_\nu ^{\mathrm{SN}}(t/s)`$ was replaced by $`F_\nu ^{\mathrm{SN}}(t+\tau )`$. Here, $`\tau =0`$ refers to GRB 980425/SN 1998bw (Iwamoto et al. 1998). If $`\tau <0`$ the SN preceded the onset of the GRB.

We describe the afterglow light curve by a broken power-law (Beuermann et al. 1999; Rhoads & Fruchter 2001),

$$F_\nu ^{\mathrm{AG}}(t)=\text{const}[(t/t_b)^{\alpha _1n}+(t/t_b)^{\alpha _2n}]^{1/n},$$
(A2)
with const=$`2^{1/n}\mathrm{\hspace{0.17em}10}^{0.4m(t_b)}.`$ Here $`t`$ is the time after the burst (in the observer frame), $`\alpha _1`$ is the pre-break decay slope of the afterglow light curve, $`\alpha _2`$ is the post-break decay slope, and $`t_b`$ is the break time. The parameter $`n`$ characterizes the sharpness of the break; a larger $`n`$ implies a sharper break. In most cases the parameter $`n`$ (Eq. A2) had to be fixed, otherwise the iteration did not converge. The reason was the usually too small number of data points around the break time. In these cases we set $`n=10`$, producing a relatively sharp break in the light curve. However, this procedure did not strongly affect the deduced supernova parameters $`k`$ and $`s`$ (Eq. A1).

In the observer frame the flux density of the time-dependent supernova light is given by (cf. Dado et al. 2002a)

$$F_\nu ^{\mathrm{SN}}(t)=\frac{1+z_{\mathrm{SN}}}{1+z_{\mathrm{bw}}}\frac{d_{L,\mathrm{bw}}^2}{d_{L,\mathrm{SN}}^2}F_{\mathrm{bw}}(\nu \frac{1+z_{\mathrm{SN}}}{1+z_{\mathrm{bw}}};t\frac{1+z_{\mathrm{bw}}}{1+z_{\mathrm{SN}}}).$$
(A3)
Here ’SN’ stands for the GRB supernova under consideration, and ’bw’ represents SN 1998bw ($`z`$=0.0085; Tinney, Stathakis, & Cannon 1998). We calculated the luminosity distance, $`d_L`$, assuming a flat universe with matter density $`\mathrm{\Omega }_M=`$ 0.3, cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=`$ 0.7, and Hubble-constant $`H_0=65`$ km s-1 Mpc-1.

We always fitted photometric magnitudes. After manipulating Eqs. (A1, A2), the apparent magnitude of the OT in a given photometric band is given as

$$m_{\mathrm{OT}}(t)=2.5\mathrm{log}\{10^{0.4m_c}[(t/t_b)^{\alpha _1n}+(t/t_b)^{\alpha _2n}]^{1/n}+k\mathrm{\hspace{0.17em}10}^{0.4m_{\mathrm{SN}}(t/s)}+10^{0.4m_{\mathrm{host}}}\}.$$
(A4)
Again, $`t/s`$ was replaced by $`t+\tau `$ when we allow for a delay between SN and GRB. Equation (A4) has eight free parameters: $`\alpha _1,\alpha _2,n,t_b,k,s,m_{\mathrm{host}}`$, and the constant $`m_c`$, which absorbs the constant of Eq. (A2) and corresponds to the magnitude of the fitted light curve for the case $`n=\mathrm{}`$ at the break time $`t_b`$. If there is no break in the light curve, then Eq. (A4) reduces to

$$m_{\mathrm{OT}}(t)=2.5\mathrm{log}\{10^{0.4m_1}t^\alpha +k\mathrm{\hspace{0.17em}10}^{0.4m_{\mathrm{SN}}(t/s)}+10^{0.4m_{\mathrm{host}}}\},$$
(A5)
where $`m_1`$ is the brightness of the afterglow at $`t=1`$ day after the burst (if $`t`$ is measured in days).

### A.2 Redshifting the SN 1998bw light curves

Equation (A1) requires as an input the function $`F_\nu ^{\mathrm{SN}}(t)`$ for arbitrarily frequencies in the optical bands. Spectra from SN 1998bw are available in the literature, but the time coverage of published broad-band photometry is much better. Therefore, we constructed $`F_\nu ^{\mathrm{SN}}(t)`$ based on published $`UBVRI`$ light curves (Galama et al. 1998), assuming that we can smoothly interpolate between adjacent photometric bands. Thereby, we have taken into account that various broad-band features are inherent to the spectral energy distribution of SN 1998bw that develop with time (e.g., Patat et al. 2001; Stathakis et al. 2000). Therefore, for different photometric bands SN 1998bw light curves peak at different times $`t_\nu ^{\mathrm{max}}`$, at different flux densities, $`F_\nu ^{\mathrm{max}}=F_\nu (t_\nu ^{\mathrm{max}})`$, and have different shapes.222For reasons of clarity, in this section we omit the index ’SN’ at $`F_\nu `$; all flux densities refer to SN 1998bw. In the following, we demonstrate our numerical approach using the $`U`$ and $`B`$ bands as an example (Fig. 7).

Let $`\nu _U`$ be the central frequency of the $`U`$ band, $`\nu _B`$ be the central frequency of the $`B`$ band, and $`\epsilon `$ be defined as $`0\epsilon 1`$. Let $`\nu ^{}`$ be the frequency for which the light curve $`F_\nu (t)`$ is required, then the relation $`\nu ^{}=\nu _U+\epsilon (\nu _B\nu _U)`$ defines the value of $`\epsilon `$. Assuming a smooth behavior of $`F_\nu (t)`$ between the $`U`$ band and the $`B`$ band, for the frequency-dependent peak flux of the light curve at the frequency $`\nu ^{}`$, we assume

$$\mathrm{log}F_\nu ^{}^{\mathrm{max}}=\mathrm{log}F_U^{\mathrm{max}}+\epsilon (\mathrm{log}F_B^{\mathrm{max}}\mathrm{log}F_U^{\mathrm{max}}).$$

Similarly, for the frequency-dependent peak time of the supernova light curve at a fixed frequency, $`\nu ^{}`$, we write

$$t_\nu ^{}^{\mathrm{max}}=t_U^{\mathrm{max}}+\epsilon (t_B^{\mathrm{max}}t_U^{\mathrm{max}}).$$
(A6)
Finally, in order to model the frequency-dependent shape of the SN 1998bw light curves, we normalize them to their peak flux and peak time. Then, at a given frequency we have $`F_\nu (t)=\eta _\nu F_\nu ^{\mathrm{max}},`$ where $`\eta `$ is a function of the ratio $`t/t_\nu ^{\mathrm{max}}`$ and $`0\eta 1`$. Correspondingly, our ansatz for the shape function $`\eta _\nu ^{}`$ for a redshifted SN 1998bw is

$$\mathrm{log}\eta _\nu ^{}(x)=\mathrm{log}\eta _U(x)+\epsilon [\mathrm{log}\eta _B(x)\mathrm{log}\eta _U(x)],$$
(A7)
where $`x=t_{\mathrm{host}}/t_\nu ^{}^{\mathrm{max}}`$ and

$$t_{\mathrm{host}}=t\frac{1+z_{\mathrm{bw}}}{1+z_{\mathrm{SN}}}$$
(A8)
is measured in the host frame (symbols follow Eq. A3).

Once an ensemble of functions $`F_\nu (t)`$ has been calculated, the apparent magnitude of the redshifted SN 1998bw in a given photometric band (Eqs. A4, A5) is obtained by integrating over the flux density per unit wavelength ($`F_\lambda (t)=\nu ^2/cF_\nu (t)`$), multiplied by the corresponding filter response function $`S_\lambda `$. For $`S_\lambda `$ we used the transmission curves for Bessel filters provided on the internet pages of the European Southern Observatory for VLT-FORS1 with reference to Bessel (1979). For the transformation between photometric magnitudes and flux densities we used the calibration constants provided by Fukugita, Shimasaku, & Ichikawa (1995; their table 9) and Zombeck (1990). The so calculated broad-band light curves of a redshifted SN 1998bw were then used as an input function for Eqs. (A4, A5).

The results of our numerical procedure were compared with corresponding results published by Dado et al. (2002a) and Bloom et al. (2002), and we found close agreement. We used our procedure to correctly predict the color evolution of GRB-SN 030329 (Zeh, Klose, & Greiner 2003), and have performed a very good numerical fit for the light curves of GRB-SN 011121 (Greiner et al. 2003b). The limit of our procedure is given by the chosen photometric band in combination with the redshift of the burster. Once we can no longer interpolate in between the $`UBVRI`$ bands, but have to extrapolate into the UV domain (cf. Bloom et al. 1999), results become less accurate.

## Appendix B Notes on individual bursts with detected supernova bump

GRB 970228:  The light curves were constructed based on the data compiled by Galama et al. (2000). The $`R_c`$ band light curve shows evidence for extra light appearing $``$2 weeks after the burst, which can be attributed to an underlying SN 1998bw component at the redshift of the burster (Galama et al. 2000; Reichart 1999). There is no evidence for a break, but its presence cannot be excluded due to the rather sparse data set. In the $`I_c`$ band the SN bump is most visible, but the light curve suffers from a lack of early-time data. In the $`V`$ band the significance for extra light is even lower, again due to the lack of observational data.

GRB 980703:  The search of an SN component in the afterglow of GRB 980703 is affected by the relatively bright host. Most data were taken from Bloom et al. (1998), Castro-Tirado et al. (1999), Holland et al. (2001), and Vreeswijk et al. (1999). Evidence for a late-time bump is rather weak.

GRB 990712:  This burst had a relatively bright host galaxy hampering the long-term study of its afterglow. We used the data presented by Fruchter et al. (2000a), Hjorth et al. (2000), and Sahu et al. (2000) to analyze the light curves. Although the $`R_c`$ band light curve is well-sampled, the bright host may have hidden the detectability of a break at later times. Our numerical procedure finds evidence for extra-light, confirming the finding by Björnsson et al. (2001).

GRB 991208:  We used the compilation of data by Castro-Tirado et al. (2001), with additional data from Dodonov et al. (1999), Halpern et al. (1999), Garnavich et al. (1999), and Fruchter et al. (2000b), including our late-time observation of the host in early 2003 to analyze the light curves. The $`R_c`$ band data can be fitted with or without the inclusion of a break. In the former case the break is mainly due to a single data point at $`t7`$ days. Most likely, the afterglow was discovered after a break had already occurred in the light curve. The afterglow parameters were determined in the $`R_c`$ band. These parameters fit well in the $`V`$ band. Unfortunately, in the $`I_c`$ band no data were obtained during the time of the SN bump.

GRB 000911:  This was a long-lasting burst with a duration of $``$500 seconds (Hurley et al. 2000; Price et al. 2002b). The optical afterglow was in detail observed by Price et al. (2002b) and Lazzati et al. (2001). We confirm the finding by Lazzati et al. that the published data show evidence for a bump in $`VRI`$ at later times.

GRB 010921:  This burst occurred in a rather crowded stellar field and had a relatively large error box (Hurley et al. 2001), which hampered the early detection of its afterglow (Price et al. 2001). The light curve is therefore not well-sampled. Using the data published by Park et al. (2002) and Price et al. (2002a) combined with our late-time observations of the host, our numerical procedure finds evidence for extra light with its peak time $`2`$ weeks after the burst. This result is obtained when we adopted a single power-law decay, where the decay slope $`\alpha `$ was deduced from the $`r^{}`$ band light curve. Our procedure finds a SN component with $`k=0.68\pm 0.48`$ and $`s=0.68\pm 0.28`$. Within the given uncertainties this is not in conflict with the upper limit reported by Price et al. (2003).

GRB 011121:  This was the nearest known burst at the time of its discovery (excluding GRB 980425/SN 1998bw). It showed clear evidence for an underlying SN component in several photometric bands (Bloom et al. 2002; Dado, Dar, & de Rújula 2002b; Garnavich et al. 2003; Greiner et al. 2003b). In our fit we included late-time Hubble Space Telescope data (Bloom et al. 2002). A break is apparent in the light curve at $`t1`$ day (see Greiner et al. 2003b; note that in Greiner et al. we assumed a Galactic extinction towards SN 1998bw of 0 mag).

GRB 020405:  This is the second burst with known redshift and a well-observed bump in its late-time afterglow. We used host-subtracted data provided by N. Masetti (private communication) to analyze the light curves. Extra light apparent in the late-time light curve can be attributed to an underlying SN component, as already noted by Masetti et al. (2003). Our procedure also detects a break in the light curves at $`t2`$ days.

GRB 021211:  For the fit we included data published by Della Valle et al. (2003), Fox et al. (2003), Li et al. (2003), and Pandey et al. (2003). A weak bump is apparent at late times, which is most likely due to an underlying SN component given its (weak) spectral confirmation (Della Valle et al. 2003).


format: text

text:
# To appear in The Astrophysical Journal Letters, the 2003-November-20th issue Modeling the Infrared Emission from the ϵ Eridani Disk

## 1 Introduction

Together with Vega, Fomalhaut ($`\alpha `$ PsA), and $`\beta `$ Pictoris, the nearby (distance to the Earth $`d3.2\mathrm{pc}`$) low-mass ($`m_{}0.8m_{}`$) K2V main-sequence (MS) star $`ϵ`$ Eridani constitutes the “Big-Four” family of prototypical “Vega-excess” stars, which are MS stars with infrared (IR) radiation in excess of what is expected from their stellar photospheres (Gillett 1986). The IR excess associated with “Vega-excess” stars is generally attributed to the thermal emission of dust grains orbiting the central star in the form of a disk or ring and heated by its stellar radiation. The dust is thought to be “second generational” in nature rather than “primordial”; i.e., the dust is considered as debris generated by collisions of larger objects such as cometary and/or asteroidal bodies rather than remnants left over from the star formation process (Backman & Paresce 1993; Lagrange, Backman, & Artymowicz 2000; Zuckerman 2001).

In recent years, $`ϵ\mathrm{Eri}`$ has received much attention since it is thought to be analogous to the solar system at an age about 10–20% of that of the Sun: (1) in comparison with the other three “Big-Four” prototypes which all are much more luminous A-type stars, $`ϵ\mathrm{Eri}`$ is closer in type, mass, and luminosity ($`L_{}0.33L_{}`$; Soderblom & Däppen 1989) to the Sun (see Greaves & Holland 2000); (2) $`ϵ\mathrm{Eri}`$ has a ring-like structure of $``$55$`\mathrm{AU}`$ in radius (which is very similar to the Kuiper Belt in size) and a central cavity inside $``$30$`\mathrm{AU}`$ (which is equivalent in scale to Neptune’s orbit) as revealed by its 850$`\mu \mathrm{m}`$ image obtained using the SCUBA (Submillimetre Common-User Bolometer Array) camera at the James Clerk Maxwell Telescope (JCMT; Greaves et al. 1998). (3) as the oldest (with an age of $`800\mathrm{Myr}`$; Henry et al. 1996) “Big-Four” Vega-type star, the formation of planets in its disk may have completed since planet formation is thought to be complete by about 10 to 100$`\mathrm{Myr}`$ (Lissauer 1993); (4) a roughly Jupiter-mass planet with a period of $``$6.9$`\mathrm{yr}`$ appears to be present inside the central cleared region of the $`ϵ\mathrm{Eri}`$ disk, as suggested by the radial velocity measurements spanning the years 1980.8–2000.0 (Hatzes et al. 2000).

To better understand the origin and evolution of protoplanetary dust disks and the formation process of planetary systems, we need to know the physical, chemical, and dynamical properties of their constituent grains. In previous papers (Li & Greenberg 1998, Li & Lunine 2003a,b), we have shown that a simple porous dust model is successful in reproducing the spectral energy distributions from the near-IR to millimeter wavelengths of the disks around the Herbig Ae/Be stars $`\mathrm{HR}\mathrm{\hspace{0.17em}4796}\mathrm{A}`$ and $`\mathrm{HD}\mathrm{\hspace{0.17em}141569}\mathrm{A}`$ and the Vega-type star $`\beta `$ Pictoris, including the 9.7$`\mu \mathrm{m}`$ amorphous and the 11.3$`\mu \mathrm{m}`$ crystalline silicate features or the 7.7$`\mu \mathrm{m}`$ and 11.3$`\mu \mathrm{m}`$ polycyclic aromatic hydrocarbon (PAH) features. In this Letter we will study the dust IR to submillimeter emission of the disk around the more evolved star $`ϵ\mathrm{Eri}`$. The objectives of this Letter are four-fold: (1) we wish to infer the physical and chemical properties of the dust in the $`ϵ\mathrm{Eri}`$ disk and to infer the disk structure; (2) we wish to gain insights into the evolution of dust mass, size, morphology, and mineralogy at different stellar ages; (3) we wish to know how widely applicable is the porous dust model; if it can be shown that the porous dust model is “universally” valid for disks at different evolutionary stages and with different geometrical structures, it will be a valuable guide for interpreting future data sets obtained by the upcoming Space Infrared Telescope Facility (SIRTF); (4) we wish to make mid-IR spectral, imaging and broadband photometry predictions for the $`ϵ\mathrm{Eri}`$ disk; these predictions can be compared to future SIRTF observations in order to test further the validity of the porous dust model.

## 2 The Porous Dust Model

We have discussed in detail in Li & Greenberg (1998) and Li & Lunine (2003a,b) that the dust in protoplanetary disks most likely consists of porous aggregates of protostellar materials of interstellar origin, although there is no consensus regarding the degree to which the original interstellar materials have been processed in protostellar nebulae. Unless there are IR signatures (e.g. the 11.3$`\mu \mathrm{m}`$ crystalline silicate feature) suggestive of the thermal history of the dust available, it is difficult to determine the composition of the dust as a priori. Therefore, in modeling the IR emission of dust disks we usually invoke two extreme dust types: “cold-coagulation” dust made of coagulated but otherwise unaltered protostellar interstellar grains, and “hot-nebula” dust made of grains that are highly-processed in protostellar nebulae where silicate dust is annealed and carbon dust is destroyed by oxidization (see Li & Lunine 2003a,b for details). Since the $`ϵ\mathrm{Eri}`$ disk is too cold to emit at the mineralogical “fingerprinting” mid-IR features (e.g. the 9.7 and 11.3$`\mu \mathrm{m}`$ silicate bands and the 3.3, 6.2, 7.7, 8.6, and 11.3$`\mu \mathrm{m}`$ PAH bands), we will just adopt the “cold-coagulation” dust model. It has been shown in Li & Lunine (2003a,b) that the major differences between the IR emission from these two dust models lie in their mid-IR spectral features, while their continuum emission are very similar.

For the dust properties, we take the porosity (or fluffiness; the fractional volume of vacuum in a grain) to be $`P=0.90`$ for 3 reasons: (i) dust models with $`P0.90`$ are successful in modeling the IR emission from the disks around $`\mathrm{HR}\mathrm{\hspace{0.17em}4796}\mathrm{A}`$ (Li & Lunine 2003a), $`\mathrm{HD}\mathrm{\hspace{0.17em}141569}\mathrm{A}`$ (Li & Lunine 2003b), and $`\beta `$ Pictoris (Li & Greenberg 1998); (ii) a porosity in the range of $`0.80P0.90`$ is expected for dust aggregates formed through coagulation as demonstrated both theoretically (Cameron & Schneck 1965) and experimentally (Blum, Schräpler, & Kozasa 2003); (iii) a porosity of $`P0.90`$ for the “cold-coagulation” dust is consistent with the mean mass density of cometary nuclei for which the ice-coated “cold-coagulation” dust aggregates are plausible building blocks (see Greenberg & Li 1999). We assume a power-law dust size distribution $`dn(a)/daa^\alpha `$ which is characterized by a lower-cutoff $`a_{\mathrm{min}}`$, upper-cutoff $`a_{\mathrm{max}}`$ and power-law index $`\alpha `$; we take $`a_{\mathrm{min}}=1\mu \mathrm{m}`$ and $`a_{\mathrm{max}}=1\mathrm{cm}`$ (see Li & Lunine 2003a,b).111We assume all grains are spherical in shape; the grain size $`a`$ is defined as the radius of the sphere encompassing the entire aggregate. The choice of $`a_{\mathrm{min}}`$ and $`a_{\mathrm{max}}`$, as extensively discussed in Li & Lunine (2003a,b) for the disks around $`\mathrm{HR}\mathrm{\hspace{0.17em}4796}\mathrm{A}`$ and $`\mathrm{HD}\mathrm{\hspace{0.17em}141569}\mathrm{A}`$, will not affect our conclusion.  As far as the dust composition is concerned, we take the “cold-coagulation” dust model – the dust is assumed to be composed of amorphous silicate and carbonaceous materials (and H2O-dominated ices in regions colder than $``$110–120$`\mathrm{K}`$); the mixing mass ratios for the silicate, carbon and ice constituent grains are approximately derived in Li & Lunine (2003a), by assuming cosmic abundance, to be $`m_{\mathrm{carb}}/m_{\mathrm{sil}}0.7`$ and $`m_{\mathrm{ice}}/(m_{\mathrm{sil}}+m_{\mathrm{carb}})0.8`$ where $`m_{\mathrm{sil}}`$, $`m_{\mathrm{carb}}`$, and $`m_{\mathrm{ice}}`$ are respectively the total mass of the silicate, carbon, and ice subgrains.

Guided by the 850$`\mu \mathrm{m}`$ SCUBA/JCMT image and the 850$`\mu \mathrm{m}`$ flux radial profile (see Figs. 1,2 of Greaves et al. 1998), we will adopt a Gaussian functional formula for the dust spatial surface density distribution $`\sigma (r)=\sigma _\mathrm{p}\mathrm{exp}[4\mathrm{ln}2\{(rr_\mathrm{p})/\mathrm{\Delta }\}^2]`$ which describes the dust ring peaking at $``$55$`\mathrm{AU}`$ from the star. The extent of the $`ϵ\mathrm{Eri}`$ disk is cut off at an inner boundary of $`r_{\mathrm{in}}=0.05\mathrm{AU}`$ where small grains are heated to temperatures $`1500\mathrm{K}`$, and an outer boundary of $`r_{\mathrm{out}}=200\mathrm{AU}`$ outside of which there is little dust emission.

We use Mie theory together with the Bruggman effective medium theory (Krügel 2003) to calculate the absorption cross sections of the fluffy heterogeneous dust aggregates. Dielectric functions are taken from Draine & Lee (1984) for amorphous silicate dust; Li & Greenberg (1997) for carbonaceous dust; and Li & Greenberg (1998) for H2O-dominated “dirty” ice. The $`ϵ\mathrm{Eri}`$ stellar radiation is approximated by the Kurucz (1979) model atmosphere spectrum for K2V stars with an effective temperature of $`T_{\mathrm{eff}}=5000\mathrm{K}`$, a surface gravity of $`\mathrm{lg}`$$`g`$=4.0, and a solar metallicity $`Z/Z_{}=1`$. The calculated dust IR emission spectrum is compared with the IRAS flux at 12, 25, 60, and 100$`\mu \mathrm{m}`$ and the SCUBA flux at 450 and 850$`\mu \mathrm{m}`$ (see Table 1 in Greaves et al. 1998). There also exist 1300$`\mu \mathrm{m}`$ photometric measurements for the $`ϵ\mathrm{Eri}`$ disk (e.g., Chini, Krügel, & Kreysa 1990, Chini et al. 1991), but we will not include these data since these measurements were all made using single-element bolometers with telescope beams smaller than the disk extent so that the flux obtained by these measurements may be underestimated (Zuckerman & Becklin 1993, Weintraub & Stein 1994).

## 3 Model Results

The model has four free parameters: $`\alpha `$ – the dust size distribution power-law exponent, $`r_\mathrm{p}`$ – the disk radial location where the dust spatial distribution peaks; $`\mathrm{\Delta }`$ – the FWHM of the dust spatial distribution; $`\sigma _\mathrm{p}`$ – the mid-plane surface density at $`r=r_\mathrm{p}`$. Observational constraints include the observed IR spectral energy distribution and the SCUBA 850$`\mu \mathrm{m}`$ image. The latter is particularly important for constraining the spatial distribution of the dust in the $`ϵ\mathrm{Eri}`$ disk. Instead of using the smoothed 850$`\mu \mathrm{m}`$ image presented in Greaves et al. (1998), we constructed an unsmoothed 850$`\mu \mathrm{m}`$ image from all 1997–2000 data in the SCUBA archive using the data reduction procedure outlined in Bendo et al. (2003). However, both the image and the radial profile obtained from our analysis show an overall good agreement with those of Greaves et al. (1998).

As shown in Figure 1, the $`P=0.90`$ porous dust model222In fact the porosity is reduced to $`P^{}0.73`$ after ices fill in some of the vacuum (see Appendix B in Li & Lunine 2003a) which occurs at $`r3\mathrm{AU}`$.  with $`\alpha 3.1`$, $`r_\mathrm{p}55\mathrm{AU}`$, $`\mathrm{\Delta }30\mathrm{AU}`$, and $`\sigma _\mathrm{p}966\mathrm{cm}^2`$ provides an almost perfect fit to the observed IR to submillimeter emission (Fig. 1a) and to the azimuthally averaged flux radial distribution at 850$`\mu \mathrm{m}`$ (Fig. 1b) except that the model is deficient in emitting at 850$`\mu \mathrm{m}`$ inside $``$28$`\mathrm{AU}`$ (see below for further discussion).333The model-predicted 850$`\mu \mathrm{m}`$ radial profile is convolved with the SCUBA $`14^{\prime \prime }`$ beam using the actual SCUBA point spread function (PSF).  This model requires a total dust mass of $`m_\mathrm{d}4.22\times 10^{25}\mathrm{g}7.1\times 10^3m_{}`$. Let the goodness of fit be $`\chi ^2/N_{\mathrm{obs}}\frac{1}{N_{\mathrm{obs}}}_{i=1}^{N_{\mathrm{obs}}}\left\{\left(\left[\lambda F_\lambda \right]_{\mathrm{mod}}\left[\lambda F_\lambda \right]_{\mathrm{obs}}\right)/\left[\lambda \mathrm{\Delta }F_\lambda \right]_{\mathrm{obs}}\right\}^2`$ where $`\left[\lambda F_\lambda \right]_{\mathrm{obs}}`$ and $`\left[\lambda F_\lambda \right]_{\mathrm{mod}}`$ are respectively the observed and model-predicted flux densities, $`\left[\lambda \mathrm{\Delta }F_\lambda \right]_{\mathrm{obs}}`$ is the observational flux uncertainty, $`N_{\mathrm{obs}}6`$ is the number of data points to be fit. This model gives $`\chi ^2/N_{\mathrm{obs}}0.54`$.

The apparent rise in flux density within $`r`$$``$$`4.5^{\prime \prime }`$ ($``$15$`\mathrm{AU}`$) was also seen in the smoothed SCUBA 850$`\mu \mathrm{m}`$ radial profile of Greaves et al. (1998; see their Fig. 2). This rise may be associated with the “small peak $`5^{\prime \prime }`$$`10^{\prime \prime }`$ south of star” seen in the SCUBA 850$`\mu \mathrm{m}`$ image (see Fig. 1 of Greaves et al. 1998). If this is real, we then need to resort to an additional dust component lying close to the star. This inner warm component would resemble the zodiacal dust cloud in our solar system. However, the SCUBA images are quite noisy in the inner region, so the anomalous emission inside $`r4.5^{\prime \prime }`$ may be instrumental noise or image processing artifacts. We therefore exclude this “inner component” emission in our model.

## 4 Discussion

We have seen in §3 that the very same porous dust model shown to be successful in modeling the IR emission from the dust disks around the Herbig Ae/Be stars $`\mathrm{HR}\mathrm{\hspace{0.17em}4796}\mathrm{A}`$ (Li & Lunine 2003a) and $`\mathrm{HD}\mathrm{\hspace{0.17em}141569}\mathrm{A}`$ (Li & Lunine 2003b) and the young MS star $`\beta `$ Pictoris (Li & Greenberg 1998) is also successful in reproducing the IR to submillimeter SED as well as the 850$`\mu \mathrm{m}`$ flux density radial profile of the more evolved MS star $`ϵ\mathrm{Eri}`$ ($``$$`800\mathrm{Myr}`$ old), although the model appears unable to provide sufficient 850$`\mu \mathrm{m}`$ emission inside $``$28$`\mathrm{AU}`$ (see Fig. 1b). It is very likely that this inner component emission is not real (see §3); but if one really wants to push us to account for this emission, we find that it could be explained by invoking an inner “zodiacal” dust cloud of $``$$`5.2\times 10^4m_{}`$ mass located at $``$2$`\mathrm{AU}`$ from the star.

SIRTF will be capable of sensitive imaging using the Infrared Array Camera (IRAC) at 3.6, 4.5, 5.8, and 8.0$`\mu \mathrm{m}`$, and using the Multiband Imaging Photometer (MIPS) at 24, 70, and 160$`\mu \mathrm{m}`$. SIRTF will also be able to perform low-resolution 5–40$`\mu \mathrm{m}`$ and high-resolution 10–37$`\mu \mathrm{m}`$ spectroscopic observations using the Infrared Spectrograph (IRS) instrument. However, the dust in the $`ϵ\mathrm{Eri}`$ disk seems to be too cold to emit in an appreciable quantity in the IRAC bands (see Fig. 1a). Neither can this dust be heated sufficiently to emit at the characteristic mid-IR spectral features. But the MIPS imaging will provide powerful constraints on the $`ϵ\mathrm{Eri}`$ dust spatial distribution. We have therefore calculated the MIPS 24, 70, and 160$`\mu \mathrm{m}`$ flux density radial distributions for the best-fit porous dust model (see Fig. 1c). All profiles are convolved with the MIPS model PSFs. In Table 1 we show the model-predicted band-averaged intensities.

The dust in the $`ϵ\mathrm{Eri}`$ disk is stable against the radiation pressure. As shown in Figure 2a, $`\beta _{\mathrm{RP}}`$, the ratio of radiative pressure (RP) force to gravitational force is found to be smaller than 0.3 for all grains. But the Poynting-Robertson (PR) drag would remove grains smaller than $``$270$`\mu \mathrm{m}`$ at $`r=55\mathrm{AU}`$ in a timescale ($`\tau _{\mathrm{PR}}`$) shorter than the age of the $`ϵ\mathrm{Eri}`$ system (see Fig. 2b).444Closer to the star, the threshold grain size below which grains will be ejected increases since $`\tau _{\mathrm{PR}}r^2`$ (see Eq.22 in Li & Lunine 2003b).  By integrating the PR dust removal rate $`\left(4\pi /3\right)\rho a^3/\tau _{\mathrm{PR}}(a,r)`$ \[where $`\rho `$ is the mass density of the dust aggregate: $`\rho 0.25(0.46)\mathrm{g}\mathrm{cm}^3`$ for the $`P=0.90`$ ($`P^{}=0.73`$ ice-coated) cold-coagulation type dust; see Appendix B in Li & Lunine 2003a\] over the whole size range and over the entire disk, we estimate the PR dust mass loss rate to be $`4.3\times 10^{11}m_{}\mathrm{yr}^1`$ for the porous dust model. Over the life span of $`ϵ\mathrm{Eri}`$, roughly $`0.034m_{}`$ of dust is lost by the PR drag. Therefore, the observed grains need to be continuously replenished, very likely from cascade collisions and evaporation of larger bodies such as comets and/or asteroids.

The inner “zodiacal” cloud, if it indeed exists, could be maintained by grains released from the sublimation of bombarding comets near the star.555In contrast, the cold outer ring is more likely to be collisionally replenished, although the sublimation of cometary CO volatiles (with a sublimation temperature of $``$20$`\mathrm{K}`$) occurring inside $`r`$$`<`$110$`\mathrm{AU}`$ does generate dust outflow. But dust outflow only becomes significant when comets are close to the star, say, at $`r3\mathrm{AU}`$ where H2O ice starts to sublimate.  Taking the total dust mass loss $``$$`3\times 10^{16}\mathrm{g}`$ in the 1997 apparition of comet Hale-Bopp (Jewitt & Matthews 1999) as a representative dust supply rate (per comet per apparition), we estimate a cometary infalling rate of $``$7$`\mathrm{yr}^1`$, roughly consistent with the estimation of $`<`$15$`\mathrm{yr}^1`$ derived by Dent et al. (1995) from the observed CO upper limit assuming sublimating comets as a CO replenishment source in the $`ϵ\mathrm{Eri}`$ disk.666If we adopt the Halley dust loss rate $``$$`1\times 10^{14}\mathrm{g}`$ per orbit (Lisse 2002), we would expect a comet infall rate $``$300 times more frequent. This is much higher than what is observed ($``$20–30 comets per year) for the present solar system. But we note that the number of cometary bodies and the bombardment rate were much higher at the early stage of our solar system. The solar system zodiacal cloud has a mass loss rate of $``$$`5\times 10^{14}m_{}\mathrm{yr}^1`$; Lisse (2002) argued that the short-period comets alone are enough to supply the dust. In the $`\beta `$ Pictoris disk, the number of observed cometary bodies (crossing the line of sight) is $``$300–400 per year (Beust & Morbidelli 1996); however, the actual number would be much higher since many objects do not cross the line of sight and are therefore not visible in spectroscopy.

The grain-grain collision timescale is only $`0.1\mathrm{Myr}`$ for dust at $`r_\mathrm{p}=55\mathrm{AU}`$, significantly shorter than the PR timescale $`\tau _{\mathrm{PR}}`$. But it is unlikely for grain-grain collisions to effectively remove the dust from the disk; instead, their major role is to re-distribute the dust over different size bins through fragmentation. The partially cleared central cavity is most likely created by the accumulation of gas and grains into planetesimals, as supported by the non-detection of circumstellar gas in this disk (Dent et al. 1995, Liseau 1999); and by the fact that the PR drag alone is not sufficient to clear the dust out to large radii up to $``$30$`\mathrm{AU}`$ (Jura 1990, Greaves et al. 1998, Dent et al. 2000).

Dent et al. (2000) modeled the $`ϵ\mathrm{Eri}`$ SED in terms of compact grains, assuming the dust absorption efficiency to be $`Q_{\mathrm{abs}}=1`$ for $`\lambda `$$`<`$$`30\mu \mathrm{m}`$ and $`Q_{\mathrm{abs}}=\left(\lambda /30\mu \mathrm{m}\right)^{0.8}`$ for $`\lambda `$$``$$`30`$$`\mu \mathrm{m}`$, and adopting a uniform surface density distribution \[i.e. $`\sigma (r)\mathrm{constant}`$\] in the range of 50$``$$`r`$$``$80$`\mathrm{AU}`$. Satisfactory fit to the observed SED was achievable. Using emissivities calculated from real dust mixtures, however, Sheret, Dent, & Wyatt (2003) found that compact grains are not able to fit the observed SED.

The model presented in this Letter is somewhat simplified in the sense that we have assumed a symmetrical morphology for the $`ϵ\mathrm{Eri}`$ disk. Asymmetries and discrete flux enhancements (“clumps” or “bright spots”) are actually seen in its 850$`\mu \mathrm{m}`$ SCUBA images (Greaves et al. 1998). As suggested by Ozernoy et al. (2000) and Quillen & Thorndike (2002), these features may be caused by the capture of grains into the mean motion resonances with a planet of Jupiter-mass embedded in the disk. This simplification would not affect the general conclusions of this Letter since the IR emission is mainly determined by the general distribution of dust in the disk, that is, a narrow ring.

## 5 Conclusion

The porous dust model, previously developed for the extended wedge-like disk around the $``$15$`\mathrm{Myr}`$-old MS star $`\beta `$ Pictoris and the narrowly confined ring-like disks around the Herbig Ae/Be stars $`\mathrm{HR}\mathrm{\hspace{0.17em}4796}\mathrm{A}`$ ($``$8$`\mathrm{Myr}`$ old) and $`\mathrm{HD}\mathrm{\hspace{0.17em}141569}\mathrm{A}`$ ($``$5$`\mathrm{Myr}`$ old), is shown also applicable to the dust disk around the $``$800$`\mathrm{Myr}`$-old MS star $`ϵ\mathrm{Eri}`$, the oldest “Big-Four” Vega-excess star. Modeled as a ring peaking at $``$55$`\mathrm{AU}`$ from the star with a FWHM $``$30$`\mathrm{AU}`$, analogous to the young solar system Kuiper Belt, the observed IR to submillimeter emission as well as the 850$`\mu \mathrm{m}`$ flux density radial distribution of the $`ϵ\mathrm{Eri}`$ disk is closely reproduced by porous grains of $``$$`7.1\times 10^3m_{}`$ mass consisting of $``$90% vacuum in volume. Predictions are made for future SIRTF/MIPS imaging observations that will offer a useful test of the porous dust model.

We thank J.S. Greaves, E.K. Holmes, A. Lecavelier des Etangs, R. Malhotra, K.A. Misselt, C. Papovich, W.T. Reach, G. Schneider, Z. Sekanina, M.D. Silverstone, and the anonymous referee for helpful discussions and/or suggestions. A. Li thanks the University of Arizona for the “Arizona Prize Postdoctoral Fellowship in Theoretical Astrophysics”. This research was supported in part by a grant from the NASA origins research and analysis program. G.J. Bendo is a Guest User of Canadian Astronomy Data Centre, which is operated by the Dominion Astrophysical Observatory for the National Research Council of Canada’s Herzberg Institute of Astrophysics.