text:
# HST Snapshot Study of Variable Stars In Globular Clusters: Inner Region of NGC 64411footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., (AURA), under NASA Contract NAS 5-26555.

## 1 Introduction

Globular clusters (GCs) are tracers of the formation of the galaxy to which they belong. One key way to analyze the stellar populations of a GC is to determine its RR Lyrae (RRL) content, which is tied to the horizontal branch (HB) morphology. Oosterhoff (1939) discovered a fundamental division of Galactic GCs into two groups based on the mean periods of their RRL stars. It was later shown that this division was associated with GC metallicity: Metal-rich GCs tend to have shorter period RRL stars compared to metal-poor GCs (Arp 1955b; Kinman 1959). Understanding the Oosterhoff effect and its effect on RRL stars is key for determining accurate distances and ages for the Galactic GCs, and ultimately has consequences on how we view the formation of our Galaxy.

Sandage & Wildey (1967) and van den Bergh (1967) found that a number of metal-poor GCs have very different HB morphologies even though they have very similar values of \[Fe/H\]. It was concluded that at least one parameter, in addition to metallicity, affects the HB. They also pointed out that the most distant GCs exhibited strong second parameter effects, having red HBs despite their low metallicities. Since then it has been shown that a GC’s galactocentric distance plays a role in the strength of the second parameter in the sense that at small galactocentric distances the HB morphology and \[Fe/H\] are closely tied to one another, while going out to larger distances one finds that there is more scatter in the HB morphologies, with the HB type becoming redder for a given \[Fe/H\] (Searle & Zinn 1978; Zinn 1980; Lee, Demarque, & Zinn 1994). A difference in cluster ages has been proposed as one explanation for the second parameter effect, but there are cases where this may not be the only explanation for unusual HB morphologies (e.g., Fusi Pecci, Bellazzini, & Ferraro 1996; Fusi Pecci & Bellazzini 1997). Any second-parameter affecting the HB morphology will also affect any RRL stars found in that cluster. Thus, as pointed out by Catelan, Sweigart, & Borissova (1998), RRL stars may also provide a probe into the second-parameter phenomenon, which in turn may be directly related to the origin of the Oosterhoff dichotomy itself (Castellani 1983; Renzini 1983; Lee & Zinn 1990).

While the second-parameter effect has primarily been found in more metal-poor GCs, two metal-rich systems, NGC 6388 and NGC 6441, have been found to have very unusual HB morphologies for their metallicities (Rich et al. 1997). With metallicities of $`[\mathrm{Fe}/\mathrm{H}]0.6`$ these clusters not only have strong red HBs similar to the canonical metal-rich GC 47 Tucanae, they also have components to their HBs that extend blueward across the instability strip and terminate in blue tails. In addition, the red HBs in both clusters are sloped toward increasing brightness as one moves blueward along the HB, which cannot be explained by differential reddening alone. To make these GCs even more interesting, Sweigart & Catelan (1998) have claimed that non-canonical effects are needed to explain the sloped nature of the HBs. Rich et al. speculated that the high central concentrations of NGC 6388 and NGC 6441 may have led to higher interactions in the inner regions of the clusters, resulting in the bluer HBs—without however providing an explanation for the sloping nature of the HBs. However, both Rich et al. and Layden et al. (1999) found that this hypothesis did not hold because the blue HB stars were not centrally concentrated.

In attempting to model the HBs in NGC 6388 and NGC 6441, Sweigart & Catelan (1998) proposed three alternatives, (i) a high cluster helium abundance, (ii) an increase in helium core mass due to rotation on the red giant branch, and (iii) enhancement of the envelope helium abundance due to deep mixing. The high helium abundance scenario has subsequently been ruled out on observational grounds (Layden et al. 1999; Raimondo et al. 2002). Other explanations for these unusual HB morphologies have ranged from increased mass loss along the red giant branch, to a spread in metallicity, to the clusters being dwarf galaxy remnants (Piotto et al. 1997; Sweigart 2002; Ree et al. 2002).

All three scenarios introduced by Sweigart & Catelan (1998) require that the HBs of NGC 6388 and NGC 6441 be unusually bright for the metal abundances of the clusters. As a consequence, any RRL stars found in these clusters would also be unusually bright and would therefore have unusually long periods for their metallicity. Surveys of the variable star content of NGC 6388 (Silbermann et al. 1994; Pritzl et al. 2002) and NGC 6441 (Layden et al. 1999; Pritzl et al. 2001) found that their RRL stars did indeed have periods that are unusually long for the cluster metallicity. On the other hand, the gravities (for the hot HB stars) measured by Moehler, Sweigart, & Catelan (1999) imply that the blue HB stars in these clusters are not anomalously bright, contrary to the Sweigart & Catelan (1998) predictions.

Due to the compact nature of these clusters \[$`\mathrm{log}\rho (L_{}\mathrm{pc}^3)=5.34`$ (NGC 6388), 5.25 (NGC 6441)\], it is very difficult to survey the inner regions of the clusters, where most of the cluster variables are expected to be found, from the ground. Making use of the Hubble Space Telescope (HST) snapshot program, we have obtained images of the inner region of NGC 6441 to survey its variable star content with the Wide Field Planetary Camera 2 (WFPC2). We present in this paper the results of our survey and compare the properties of the variable stars found in the inner region of the cluster to those variables found in the outer regions from the ground-based surveys. We also examine the nature of the Population II Cepheids (P2Cs) found in NGC 6441 and NGC 6388 and compare them with those in other GCs.

## 2 Observations and Reductions

Eighteen sets of observations in $`B`$, $`V`$, and $`I`$ were taken by the HST with the Wide-Field Planetary Camera 2 starting on 1999 March 16 and ending on 2000 June 30 as part of the SNAP 8251 program. The complete list of observations is given in Table 1, including NGC 6441 HST archive observations which were added to our reductions (GO 5667, GO 6095, and GO 6780). Exposure times for each filter are also listed in Table 1. NGC 6441 was centered on the WFC3 chip for our program. No restrictions were made on the orientation of the observations. Therefore, rotations between the images, in addition to translations, produced partial and varying overlap among the images obtained on the other three chips.

The observations were reduced using the daophot/ allstar/ allframe routines (e.g., Stetson 1987, 1994). We used standard point-spread functions that had been prepared for each chip and filter in the course of the Hubble Space Telescope Key Project on the Extragalactic Distance Scale (e.g., Freedman et al. 2001; Stetson et al. 1998) and our previous research on outer-halo globulars (e.g., Harris et al. 1997, Stetson et al. 1999). After the geometric relationships relating the various images to a common coordinate system and a master star list for the field had been generated, allframe was used to obtain final stellar positions on the sky, and a magnitude measurement for each star at each epoch. Aperture corrections were determined from comparatively bright, isolated stars in each image to relate them to a well-defined, common instrumental zero-point. Charge-transfer efficiency corrections and fundamental photometric zero-points similar to those of Stetson (1998), but based on a much larger database, and the color corrections of Holtzman et al. (1995) were applied in the standard fashion.

Figure 1 shows the photometric uncertainty in $`V`$ at the level of the HB as a function of clustercentric radius. We assumed $`m_{V,\mathrm{RR}}=17.54`$ mag and chose a range of magnitudes around this value of $`\pm 0.25`$ mag. All candidate variable stars were removed. The $`\sigma _V`$ relationship is nearly independent of radius. This suggests that at these limits our photometric errors are dominated by photon statistics and fitting errors and not by crowding for most stars. The increased photometric uncertainty at large radii occurs because, with the differing roll angles of the spacecraft, not all stars could be measured at all epochs.

Variable stars were found by a variety of methods including a variant of the Welch & Stetson (1993) algorithm (see Stetson 1996). Finding charts for the variables can be seen in Figure 2. We used the period finding routines described in Layden et al. (1999) and Layden & Sarajedini (2000) to determine the periods of the variable stars. This technique involves folding each star’s magnitude-time data by a sequence of periods, and at each period, fitting the resulting light curve with a set of template light curves. If the period being tested is near the star’s true period, the template will fit the data well and produce a small value of $`\chi ^2`$. Probable periods are represented by $`\chi ^2`$ minima in the array of periods and templates considered. We searched for periods in the range appropriate for each class of variable; e.g., for RRL stars the range was 0.2 to 1.0 day. Once a period was found, the search was narrowed to about 0.0001 day around the candidate period to determine the best period. In almost all cases, we found a single, dominant $`\chi ^2`$ minimum, suggesting that aliasing was not a problem despite the irregular sampling and long time interval over which the data were gathered. The few exceptions involved a small number of first-overtone (RRc) stars which produced light curves with significant scatter, and long period variables (LPVs) whose entire light cycle was not covered by our observations. We derived pulsation amplitudes and intensity-weighted magnitudes from the fitted templates. We estimate the uncertainty in these magnitudes to be 0.02–0.03 mag for the RRL stars and Population II Cepheid (P2C) stars. Typical uncertainties in the RRL and P2C periods are estimated to be 0.00002 and 0.002 days. The magnitude and period uncertainties are considerably larger for the LPV stars. In part, this is because the templates we used are not optimal for LPV stars, but a more significant discrepancy arises because many of the LPVs may not demonstrate pulsations with a regular period. The analysis here is intended only to classify the type of pulsation and to establish approximate mean-light parameters so the stars may be placed in the CMD. A comprehensive analysis of the LPV stars in NGC 6441 is reserved for a future paper.

Figure 3 shows the positions of the variables in the color-magnitude diagrams for NGC 6441 in $`V`$ versus (a) $`\mathrm{B V}`$ and (b) $`VI`$, where variables are plotted according to their intensity-weighted $`V`$ magnitudes and $`BV`$ and $`VI`$ colors. Figures 3c and 3d give an expanded view of the HB. The variable stars along the HB are clearly RRL stars; those brighter than the HB and located toward the red but blueward of the red giant branch are P2Cs; those found after the turnover of the red giant branch are LPVs. However, as noted above, a number of the LPVs have poor periods due to aliasing and have larger uncertainties in their mean magnitudes. We also found one $`\delta `$ Scuti star (V95), but given its magnitude it is likely a member of the field. Due to the small number of data points and the gaps in the phase coverage, we did not attempt to do a Fourier analysis of the light curves.

In the next sections we present the detailed results of our survey for variable stars and discuss their properties. Comparisons are also made to the variables found in ground-based surveys to determine any similarities and differences in their properties as compared to those found in this study.

## 3 RR Lyrae Stars

We were able to find 38 candidate RRL stars in the inner region of NGC 6441, with 24 of them being new. Previous ground-based surveys (Layden et al. 1999; Pritzl et al. 2001) of the cluster found another 25 RRL stars which are not included in the field coverage of the present study. We also confirmed the variability of three stars whose variablility had been suspected, but not confirmed. Thus, the total number of RRL stars known in NGC 6441 with definite classification and membership is now 63. Figure 4 shows the light curves of the 57 variables we identified where Layden’s template-fitting routines were used to fit the data. Tables 2, 3, and 4 list the $`B,V,I`$ photometric data for the variable stars. The mean properties of the variable stars are given in Table 5. The mean periods for the 26 fundamental (RRab) and 12 RRc stars found in this survey are 0.753 d and 0.365 d, respectively. These values match up very well with those found from the ground-based surveys, which were $`P_{ab}=0.75`$ d and $`P_c=0.38`$ d for 24 and 9 stars, respectively (Pritzl et al. 2001). The variables listed in Table 5 from V55 to V102, were included in both the ground-based and HST surveys, so that the two samples are not completely independent. Combining the periods of the RRL stars with good light curves, accurate classifications, and certain membership found in this survey and the ground-based survey, we find the average periods for 42 RRab and 21 RRc stars in NGC 6441 to be 0.759 d and 0.375 d, respectively. The ratio of RRc to RRL stars is 0.33. As noted in Pritzl et al. (2000, 2001), the mean period of the RRab stars is higher than that found in the typical Oosterhoff II cluster. The mean period of the RRc stars is similar to that for Oosterhoff II clusters, while the $`N_c/N_{\mathrm{RR}}`$ ratio is somewhere between the values found for Oosterhoff I and II clusters.

Considering only the RRab stars newly discovered in this survey, we find a mean period of 0.765 d. This is slightly larger than one might expect given the result from the ground-based surveys (Layden et al. 1999; Pritzl et al. 2001). Since the HST observations allowed us to better survey the inner region of NGC 6441, it may be that there is a difference between the mean properties of the RRL in the inner and outer regions. To make a more quantitative argument, we find the RRab stars within 20 arcsec of the cluster center to have a mean period of 0.766 d (17 RRab stars). Those outside this assumed radius have a mean period of 0.754 d (25 stars). The radius of 20 arcsec was chosen as the best to illustrate the difference in mean period. For radii less than this we begin to enter small-number statistics and for more than this the mean periods progressively get closer to matching.

Another way to view the shift in period is shown in Figure 5. The histogram in Fig. 5a illustrates that the distribution of RRL stars in NGC 6441 is closer to that of an Oosterhoff II cluster than that of an Oosterhoff I cluster, as discussed in Pritzl et al. (2000, 2001). A comparison of Fig. 5b and 5c shows an excess of longer period RRab stars within 20 arcsec of the cluster center. Some caution must be exercised, however, before one can conclude that the difference is real. A two-sample K-S test run on the inner and outer RRab results in a $`D`$ statistic of 0.2259, which means that the two distributions are the same at a 62.1% confidence level. A Student $`t`$-test, which specifically tests if two distributions are different based on their means, was also run and resulted in a $`t`$ value of 0.411, which means that the two distributions are indistinguishable at a 68.3% confidence level. This confirms the K-S test, thus suggesting that the two samples are not significantly different, in a statistical sense. While our HST observations have been able to survey the inner regions of NGC 6441, crowding in the regions just outside this survey area may have prevented low amplitude, long period RRab stars from being detected in the ground-based surveys. As an example, adding five stars with periods of 0.85 d outside 20 arcsec would change the mean period to match that found for the RRab inside 20 arcsec. We also note that if V118 is truly a P2C (see §4) and not an evolved RRL star, the mean period for the RRab stars within 20 arcsec would be 0.753 d. This matches up very well with the mean period for the RRab stars outside 20 arcsec and argues that there may be no significant difference between the RRab populations inside and outside 20 arcsec. Further observations of NGC 6441 with better resolution at larger radii may help clear up any distinctions between the inner and outer RRab stars. Alternately, application of the ISIS image subtraction package (Alard 2000) may improve the detection of low amplitude variables in the crowded regions of NGC 6441 observed from the ground. We did not discuss the RRc stars given their relatively low numbers.

We derive the mean magnitude for the RRL stars to be $`m_{V,\mathrm{RR}}=17.54\pm 0.02`$ mag, where the error is the standard deviation of the mean. This matches up well with what was found in the ground-based survey of Pritzl et al. (2001; $`m_{V,\mathrm{RR}}=17.51\pm 0.02`$ mag). As seen in Figure 3, a number of the candidate RRL stars are found at magnitudes brighter than the majority of the RRL stars. These variable stars have periods of 0.90582 d (V56), 0.76867 d (V110), 0.61419 d (V112), 0.97923 d (V118), 0.80573 d (V136), and 0.55581 (V145). It may be that differential reddening is playing a role given that there is a small number of RRL stars fainter than the majority of the RRL stars. It has already been shown from the RRL star light curves that differential reddening affects the ground-based observations of NGC 6441 (Layden et al. 1999; Pritzl et al. 2001). We did not attempt to determine the reddenings from the RRab stars in this survey via their $`\mathrm{B V}`$ color at the minimum (Blanco 1992) due to the small number of data points. It may also be the case that some of these stars have evolved further from the zero-age HB and are thus more luminous. The periods for two of the stars, V118 (and, to a lesser extent V56) and V145, are long for RRab and RRc stars, respectively. In fact, given the magnitude and the period of V118, it may be the faintest P2C in the cluster (see §4).

### 3.1 Period-Amplitude Diagram

Figure 6 shows the period-amplitude diagrams for the newly discovered RRL stars and those previously found (Pritzl et al. 2001). In general, the RRL stars from the two data sets lie in the same regions of the diagrams. There is more of a spread among the RRab in Fig. 6c because the $`I`$-band light curves have more scatter than the $`B`$ and $`V`$ light curves. From Figure 6 we see that the general properties of the RRL stars found in the inner region of the cluster are the same as those found in the outer regions.

It was noted in Layden et al. (1999) and Pritzl et al. (2001) that some of the RRab stars were likely blended with other stars, producing lower than expected amplitudes for their periods. Such blends are likely given the compact nature of this cluster. The candidate blended stars from the ground-based surveys can be seen in Figure 6 as those RRab stars which are scattered toward lower amplitudes. In Table 6, we show the comparison between the variable stars found in both this survey and the ground-based survey of Pritzl et al. There are clearly a number of stars which appeared significantly brighter in the ground-based survey and had lower amplitudes than found in the present survey. We also examined the images for stars neighboring those listed in Table 6. All of the stars listed which showed significant differences in their magnitudes and periods were found to have nearby stars that likely affected the photometry of the ground-based studies. The higher resolution of the WFPC2 indicates that the previous conclusion of the “low amplitude” RRab stars in the ground-based survey being blends was correct. There is also the possibility that some of these stars may be experiencing the Blazhko effect which also affects the amplitudes of RRL stars.

Layden et al. (1999) and Pritzl et al. (2000, 2001) found that the RRab stars in NGC 6441 lie at longer periods in the period-amplitude diagram compared to field RRab stars of similar metallicity. This confirms that the RRL stars, and therefore the HB, are unusually bright for their metallicity. Pritzl et al. (2001) also showed that the RRab stars in NGC 6441 fall near a line in the period-amplitude diagram established by Clement (2000) as defining the location of RRab stars in Oosterhoff type II clusters, a result first noted by Clement and reconfirmed here in Fig. 6b. As noted above, this is consistent with the mean period of the RRab actually being more like that of an Oosterhoff type II cluster than an Oosterhoff type I cluster (see also Walker 2000). As discussed in detail in Pritzl et al. (2002), it is difficult to model such GCs as NGC 6388 and NGC 6441 as Oosterhoff II systems, under the hypothesis that their variables are evolved from a position on the blue zero-age HB, due to the small number of progenitors on the blue HB (see also Fig. 3).

We also find one RRc star of unusually long period in the inner region of NGC 6441 (V145). Two similarly long period RRc stars were found in the ground-based studies (Layden et al. 1999; Pritzl et al. 2001). In Figure 6 we see that most of the RRc stars group together, while the longer period ones are seen not to partake in the more usual, inverted parabola-like (e.g., Bono et al. 1997a) locus of RRc stars in the period-amplitude plane, for reasons which are still unclear. As noted for NGC 6441 in Pritzl et al. (2001), and for the long period RRc stars in NGC 6388 in Pritzl et al. (2002), there are very few GCs that produce this type of long period RRc star.

In Figure 6 we can also see that there are two candidate RRab stars with longer periods than expected for their amplitudes (V56 and V118). As noted above, these are two of the stars which are brighter than the majority of the other RRL stars. This prompts the question as to whether these stars are evolved RRL stars or possibly P2Cs. We will discuss the issue further in the following section.

### 3.2 $`I`$-band Period-Luminosity Relations

RRL stars provide excellent distance indicators given their nearly uniform absolute magnitudes. However, while the HB is flat at the RRL level in $`B`$ and $`V`$, thus destroying any significant period-luminosity (P-L) relation that might otherwise be present due to the important range in $`T_{\mathrm{eff}}`$ compared to $`L`$, the same may not apply in redder (and bluer) bandpasses, where the HB becomes significantly non-horizontal at the RRL level. Many studies have been done to determine the $`K`$-band P-L relations for RRL stars (Longmore et al. 1990; Nemec, Nemec, & Lutz 1994; Bono et al. 2001). The RRL P-L relations in the $`I`$-band have, on the other hand, been neglected. We searched through the literature and found only Layden & Sarajedini (2003) derived a P-L relation in the I-band for the RRL in NGC 3201. Here we derive the $`I`$-band P-L relations for the RRL in NGC 6441 using the large number of RRL stars we found in this cluster.

As discussed here, in Layden et al. (1999), and in Pritzl et al. (2001, 2002), the RRL stars in NGC 6388 and NGC 6441 are unusually bright for their metallicity. Given that the general properties of the RRL stars in both of these clusters are more like Oosterhoff II clusters than Oosterhoff I, we assume that their RRL stars have absolute magnitudes similar to those of Oosterhoff II clusters such as M15. Using Eqn. 7 in Lee, Demarque, & Zinn (1990) and adopting a metallicity of $`[\mathrm{Fe}/\mathrm{H}]=2.0`$, we find an absolute magnitude for the RRL stars in NGC 6441 to be $`M_V=+0.48`$ mag. Our adopted magnitude for RRL stars in Oosterhoff Type II clusters is broadly consistent with the recent parallax determination for RRL stars by Benedict et al. (2002). They obtained $`M_V=+0.61`$ for RR Lyrae itself, which, at $`[\mathrm{Fe}/\mathrm{H}]=1.4`$, is more metal-rich and might be slightly fainter than RRL stars in a more metal-poor system. In calculating the distance modulus for NGC 6441, we adopt $`m_{V,\mathrm{RR}}=17.54`$ mag, $`E(\mathrm{B V})=0.51`$, and $`R_V=3.1`$ to find $`\mu _0=15.48`$ mag. To derive the individual absolute magnitudes we used $`A_I=1.85E(\mathrm{B V})`$. Figure 7 shows the P-L relations for the RRL stars. Removing outlying stars, we determine the P-L relations to be,

$$M_{I,\mathrm{ab}}=0.25(\pm 0.04)1.78(\pm 0.12)\mathrm{log}P$$
(1)
$$M_{I,\mathrm{c}}=0.61(\pm 0.06)1.78(\pm 0.14)\mathrm{log}P,$$
(2)
with correlation coefficients $`r=0.792`$ and 0.947, respectively. The root-mean-square deviations for the RRab and RRc relations are 0.06 and 0.04, respectively. As noted in §3, it is uncertain whether the outlying stars, which arebrighter than the majority of RRL stars of similar period, are evolved or not. Still, it is clear that the $`I`$-band provides a good filter to derive P-L relations for the RRL stars. We note that the zero-point for the above relations is uncertain due to the fact that we have assumed that the RRL stars in NGC 6441 have the same luminosity as Oosterhoff II RRL stars with $`[\mathrm{Fe}/\mathrm{H}]=2.0`$. The precise luminosity of the NGC 6441 RRL stars is not conclusively known given the unusual nature of its HB.

## 4 Population II Cepheids

In the ground-based photometry we were able to find four candidate P2Cs in NGC 6388 (Pritzl et al. 2002), but none in NGC 6441 (Layden et al. 1999; Pritzl et al. 2001). At the time it seemed strange that in two clusters with such similar HB morphology and RRL star populations, one would contain P2Cs and the other would not. In this HST snapshot survey we have found six definite P2C candidates in NGC 6441. This makes NGC 6441, along with NGC 6388, the most metal-rich GC known to contain P2Cs. Interestingly, the P2Cs in both clusters have been found near the inner regions of the clusters, though this may reflect just the comparatively small population of P2C stars compared to RRL variables. The discovery of P2Cs in these two clusters is consistent with the idea that these types of variables are found only in clusters with blue HB components (Wallerstein 1970; Smith & Wehlau 1985). <sup>2</sup><sup>2</sup>2Palomar 3 is a notable exception to this otherwise seemingly general rule (Borissova, Ivanov, & Catelan 2000).

Of the six P2Cs, five are of W Virginis type (W Vir; $`10<P<20`$ d) and one is of BL Herculis type (BL Her; $`P<10`$ d). The difference between these types is not only defined by their periods, but also by their origins. As reviewed by Wallerstein (2002, but see also Bono et al. 1997b), BL Her stars are believed to derive from stars which have evolved from the HB bluer than the RRL gap, and which are now evolving redward through the instability strip at luminosities brighter than those of RRL stars. By contrast, W Vir stars, while also considered to be the progeny of blue HB stars (e.g., Gingold 1985), are believed to be undergoing helium shell flashes leading them along blueward loops into the instability strip as they move up the asymptotic giant branch (AGB). This is consistent with studies of period changes in P2Cs in GCs. Wehlau & Bohlender (1982) found that of 12 BL Her stars in Galactic GCs, nine had increasing periods, three had unchanging periods, and zero had decreasing periods. Increasing periods are expected for variables evolving redward from the blue HB toward the AGB. For the W Vir stars the picture is not so clear, but Clement, Hogg, & Yee (1988) found that a number of W Vir stars showed both increases and decreases in period, and that at least two variables showed long term period decreases consistent with blueward evolution. The period change observations for W Vir stars are thus at least consistent with evolution of W Vir stars into and out of the instability strip during blueward loops.

### 4.1 Period-Luminosity Relation for Population II Cepheids

Given the relatively large number of P2Cs found in both NGC 6441 and NGC 6388, we would like to determine the P-L relations for these stars. We use the distance modulus derived in §3.2 for NGC 6441, $`\mu _0=15.48`$ mag, to derive the absolute magnitudes for the P2Cs found in this survey. For NGC 6388, we follow the prescription used for NGC 6441 to derive its distance modulus. Since we assume that the RRL stars in NGC 6388 are of a similar nature as those in NGC 6441, we adopt the same absolute magnitude for its RRL stars, $`M_V=+0.48`$ mag.

In Table 7 we list the data for the P2Cs in NGC 6388 and NGC 6441. For our discussions here, we will include V118 with the P2Cs, although its precise classification is uncertain. We have chosen not to include V56 since its period and magnitude are not unreasonably different from the other RRL stars. We have assumed $`R_V=3.1`$, $`R_B=4.1`$, and $`A_I=1.85E(\mathrm{B V})`$, and have adopted mean reddenings of $`E(\mathrm{B V})=0.40\pm 0.03`$ and $`0.51\pm 0.02`$ for NGC 6388 and NGC 6441, respectively (Pritzl et al. 2001, 2002). Figure 8 shows the absolute magnitudes of the P2Cs in NGC 6388 and NGC 6441 for the various filters. In each panel of Figure 8, we see that the P2Cs form a surprisingly good correlation from shorter to longer periods, with the lone exception of one star in NGC 6388. V29 of NGC 6388 was shown to be unusually bright for its period in Pritzl et al. (2002) leading to the idea that it may be blended with another star. Leaving that star out, we determined the following P-L relations,

$$M_V=0.05(\pm 0.05)1.64(\pm 0.05)\mathrm{log}P$$
(3)
$$M_B=0.31(\pm 0.09)1.23(\pm 0.09)\mathrm{log}P$$
(4)
$$M_I=0.36(\pm 0.01)2.03(\pm 0.03)\mathrm{log}P,$$
(5)
with correlation coefficients $`r=0.996`$, 0.985, and 0.998, respectively. The root-mean-square deviation for each relation is 0.07, 0.10, and 0.06 for $`V`$, $`B`$, and $`I`$, respectively. Both the slopes for $`M_V`$ and $`M_B`$ match up quite well with those found by McNamara (1995, Eqns. 11 and 12) for P2Cs with periods less than ten days, but are very different compared to the results of Alcock et al. (1998, Eqn. 5) for longer period P2Cs. The tight fits to the relations also indicate that all of the P2Cs in NGC 6388 and NGC 6441 are likely to be pulsating in one mode, presumably the fundamental mode. Following the suggestion by Alcock et al. that adding a color term results in better relations, we calculated the P-L-C relations for the NGC 6441 and NGC 6388 P2Cs and found no significant improvement in the quality of the fit. This may be explained by the already excellent correlation of the P2Cs in the P-L relations and the fact that Alcock et al. found the RV Tauri (RV Tau; $`P>20`$ days) benefited the greatest in adding in the color term.

### 4.2 Comparison of Period-Luminosity Relations

Defining a clear P-L relation is important in the use of Cepheid variable stars in distance determinations. While other Cepheid-type variables such as the classical and anomalous types have their P-L relations well-defined, the relations for the P2Cs in general are still somewhat uncertain. Arp (1955a) argued that P2Cs follow at least two parallel P-L relations, one each for fundamental and overtone mode pulsators. Others later favored a single P-L relation to describe the P2Cs (e.g., Dickens & Carey 1967; Demers & Harris 1974). Harris (1985) also assumed a single P-L relation, but with a break toward a steeper slope for P2Cs with periods greater than 12.5 days. The theoretical $`B`$-band P-L relations of Bono et al. (1997b) show a break in slope near $`P=15`$ days at the faint boundary of the domain in which P2Cs might exist, but no break in slope at the brighter boundary.

In their paper on Population II variable stars, Nemec, Nemec & Lutz (1994) revived Arp’s model stating that P2Cs are pulsating in either fundamental or first overtone modes, which followed their theme that all pulsating variable stars can be found to have both pulsation modes. However, McNamara (1995), using the data from Nemec et al., stated that the observations were better fit by one line describing the BL Her stars ($`P<10`$ d) and another fitted to the W Vir stars ($`P>10`$ d), similar to what was argued by Harris (1985). Recently, Alcock et al. (1998) used a sample of 33 P2Cs in the Large Magellanic Cloud (LMC) in the period range of $`0.9<\mathrm{log}P<1.75`$ to compare the various P-L relations. While noting that the Galactic GC P-L relations provided adequate fits to the shorter period LMC P2Cs, they found that those P-L relations did not provide a good fit for the longer period RV Tau stars, which appear more luminous than the P-L relations predicted. Combining the Galactic GC and LMC data sets, Alcock et al. found that the shorter period Galactic GC P2Cs provided a good linear extension to the LMC data with the RV Tau stars turning toward brighter magnitudes. To account for this, Alcock et al. included a color term to their P-L relations for $`\mathrm{log}P>0.9`$, which greatly reduced the scatter in the data and the up-turn of the RV Tau stars. They concluded that the RV Tau stars are direct extensions of the other P2C types. Unfortunately, the Alcock et al. data were taken with $`V`$ and $`R`$ filters, while most of the Galactic GC data is in $`B`$ and $`V`$. This makes the application of their P-L-C relation difficult to apply to other data sets.

To help clarify the issue of the differing slopes found by various studies, we investigated the P2Cs in other GCs and incorporated them in our analysis. Table 8 lists the properties of the Galactic GC P2Cs summarizing the data found in the literature. We estimated the absolute magnitudes for each star using Eqn. 7 of Lee et al. (1990) and the metallicities, reddenings, and $`V_{\mathrm{HB}}`$ given in Harris (1996). The format is the same as in Table 7. We plot the best available data for the Galactic GC P2Cs in the top row of Figure 9, including those found by Alcock et al. (1998) in the LMC where they adopted a distance modulus of 18.5 mag. The only P2Cs not included in the plot were those from NGC 4372, which are clearly too faint. The periods listed for these stars may be correct, but the faint magnitudes may be due to the fact that Kaluzny & Krzeminski (1993) were searching for contact binaries and SX Phoenicis stars resulting in the candidate P2Cs being saturated. There is some scatter among the P2Cs in Fig. 9a, which is to be expected since it is known that there is some width to the instability strip. Fig. 9b shows a lot of scatter which may possibly be attributed to the presence of older photographic data, as compared to the $`I`$ data in Fig. 9c which is only from CCD data. We plot in the bottom row of Figure 9 only those P2Cs with CCD photometry, except again for NGC 4372. This helps to reduce the scatter, as seen in Fig. 9d, but some still remains. Without any color information we are not able to tell if this scatter may be intrinsic, possibly due to the width of the P2C instability strip, or if it is photometric. It should be noted that the scatter seen at shorter periods in Fig. 9d can be attributed to the $`\omega `$ Cen P2Cs, a number of which have poor light curves (Kaluzny et al. 1997). For $`B`$ and $`I`$, the CCD data are mostly taken from this study, but it is encouraging that the other data are consistent with what we have found.

Taking a closer look at the $`B`$ and $`I`$ data, it does not look like there is a change in the slope between the BL Her and W Vir stars. However, this is a tentative conclusion given the fact that most of the data comes from this survey alone. If we examine the $`V`$ data in Fig. 9d, where more points are available, we can make some stronger arguments. There is clearly a change in the slope toward brighter magnitudes at longer periods. The question is whether this change is found at or around 10 or 12 days as suggested by Harris (1985) and McNamara (1995) or at a longer period. If we consider only the BL Her and W Vir stars, there is statistically no change in the slope between the two groups of P2Cs. On the other hand, if we consider the RV Tau stars, they clearly have a different slope from that of the shorter period Cepheids. It appears, therefore, that a single P-L relationship can be used to describe the P2Cs until periods greater than about 20 days are reached, when variables with RV Tau characteristics begin to enter the mix. Unfortunately, given the scarcity of data, we are not able to test the conclusion of Alcock et al. (1998) that adding in the color term to the RV Tau stars brings them more into line with the other P2C types.

Another possibility is raised by Fig. 9d. Alcock et al. (1998) identified a number of W Vir stars found in the LMC that are brighter than the bulk of the W Vir stars. These variables appear to constitute a straight line in combination with the brighter RV Tau stars that appears to be parallel to a line that could be drawn through the BL Her and W Vir stars. Alcock et al (1998) classified their RV Tau stars as those having periods longer than 20 days. It may be that some of their W Vir stars are actually RV Tau stars, and vice versa. There appears to be in reality no sharp dividing line between the W Vir and RV Tau stars, with the longest period W Vir variables gradually taking on RV Tau characteristics (Alcock et al. 1998).

At this point, it is uncertain how well the relations created from the NGC 6441 and NGC 6388 P2Cs apply to other GCs given the small amount of CCD photometry. Two of the more important questions are how good are the zeropoints we are using to derive the absolute magnitudes and is there any metallicity dependence in the relations which is currently not being taken into account in our relations. Clearly, more CCD photometry needs to be taken of P2Cs in order to clarify these issues, especially in the $`B`$ and $`I`$ filters, before any firm conclusions can be drawn.

## 5 Summary and Conclusions

We have found 57 variable stars in the inner regions of NGC 6441 using the HST snapshot observing program, 41 of which were previously undiscovered. Of the RRL stars found in this survey, 26 are RRab stars and 12 are RRc stars. There are no clear differences in the properties of the RRL stars in the inner and outer regions of the cluster. This is illustrated in the period-amplitude diagram where the RRL stars in this survey are located in the same regions as those found in the ground-based surveys (Layden et al. 1999; Pritzl et al. 2001). Combining the RRL stars in this survey and those found in ground-based surveys, we find the mean periods of the RRab and RRc stars to be 0.759 d and 0.375 d. This reaffirms that NGC 6441 is unusual when compared to other Galactic GCs in that, whereas it is a metal-rich GC, its RRL stars have mean periods that are similar to, if not greater than, those of very metal-poor GCs. This reconfirms the conclusions of Pritzl et al. (2001) that the RRL stars in NGC 6441 are unusually bright for their metallicity. This breaks the general metallicity-luminosity relation for the RRL stars and has consequences on their use as distance indicators. A possible trend of the inner RRab stars ($`r<20`$ arcsec) having a longer mean period than the remaining RRab stars was also detected, but further observations will be needed to make definitive conclusions. We also used the $`I`$-band photometry to derive P-L relations for the RRL stars.

We also found six P2Cs in NGC 6441 making it, along with its “sister” cluster NGC 6388, the most metal-rich GC to contain this type of variable star. Five of the P2Cs are of W Vir type and one is a BL Her star. Also, V118 may be an additional BL Her variable given its bright magnitude and long period compared to the RRL stars. Under the assumption that the RRL stars in NGC 6388 and NGC 6441 have absolute magnitudes similar to those of metal-poor GCs, we have derived the $`B`$, $`V`$, and $`I`$ period-luminosity relations for their P2Cs. These P-L relations do not show the change in slope around $`P=10`$ day which some have suggested and show a surprisingly tight correlation. We also find that only one P-L relation is needed to describe the data, i.e., it appears that the P2Cs are pulsating in one mode. An examination of the P2Cs in GCs for which $`V`$ CCD photometry is available indicates that the break in the slope of the P2C P-L relation does not come at periods of 10 or 12 days, but at longer periods, where the P2Cs begin to take on RV Tau characteristics.

Support for Proposal number SNAP 8251 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Support for M.C. was provided by Proyecto de Inicio DIPUC 2002-04E and Proyecto FONDECYT Regular 1030954.


format: text

text:
# Survival Probability of a Ballistic Tracer Particle in the Presence of Diffusing Traps

## Abstract

We calculate the survival probability $`P_S(t)`$ up to time $`t`$ of a tracer particle moving along a deterministic trajectory in a continuous $`d`$-dimensional space in the presence of diffusing but mutually noninteracting traps. In particular, for a tracer particle moving ballistically with a constant velocity $`c`$, we obtain an exact expression for $`P_S(t)`$, valid for all $`t`$, for $`d<2`$. For $`d2`$, we obtain the leading asymptotic behavior of $`P_S(t)`$ for large $`t`$. In all cases, $`P_S(t)`$ decays exponentially for large $`t`$, $`P_S(t)\mathrm{exp}(\theta t)`$. We provide an explicit exact expression of the exponent $`\theta `$ in dimensions $`d2`$ and for the physically relevant case, $`d=3`$, as a function of the system parameters.

PACS numbers: 05.70.Ln, 05.40.+j, 02.50.-r

The calculation of the survival probability of a tracer particle moving in the presence of diffusing traps is a problem of long standing interest as it appears, in various guises, in a wide variety of contexts such as reaction-diffusion systems, chemical kinetics, predator-prey models and ‘walker persistence’ problems. The tracer particle dies instantly upon meeting any of the diffusing traps. Perhaps the simplest of all these problems is the case when the diffusing traps are noninteracting and the motion of the tracer particle is governed by its own intrinsic dynamics that depends on the specific problem. For example, if the tracer particle is static, this problem is known as the target annihilation problem. Of particular interest is the case when the tracer particle itself has a diffusive motion, a problem that was first studied by Bramson and Lebowitz and has recently seen a flurry of activity. It is, however, somewhat frustrating that, despite various new developments, this diffusive target annihilation problem has defied a direct exact solution. In contrast, we show in this paper that the ballistic target annihilation problem, where the tracer particle moves ballistically with a constant velocity, is exactly solvable.

A precise definition of the general problem is as follows. Consider a set of particles initially (at time $`t=0`$) distributed randomly in a continuous $`d`$-dimensional space with average density $`\rho `$. Each of these particles subsequently undergoes independent diffusive motion with the same diffusion constant $`D`$. A tracer particle is introduced into the system at $`t=0`$ at the origin and subsequently moves according to its own prescribed equation of motion. This motion can be either deterministic or stochastic, depending on the problem. For a given trajectory $`\stackrel{}{R}_0(t)`$ of the tracer particle, we ask: what is the probability $`P_S(t)`$ that none of the random walkers hits the tracer particle up to time $`t`$? Evidently $`P_S(t)`$ depends implicitly on the trajectory $`\stackrel{}{R}_0(t)`$. For a deterministic motion of the tracer particle where the trajectory $`\stackrel{}{R}_0(t)`$ is prescribed, its survival probability is precisely $`P_S(t)`$. On the other hand, for stochastic motion of the tracer particle the survival probability is obtained by subsequently averaging $`P_S(t)`$ over all possible trajectories of the tracer particle. Thus the basic step, in either case, is to compute $`P_S(t)`$ for a given deterministic trajectory $`\stackrel{}{R}_0(t)`$. In this paper we limit ourselves to the study of this general deterministic trajectory problem and, in particular, present explicit exact results when the trajectory is ballistic (see Fig. 1).

For a given fixed trajectory $`\stackrel{}{R}_0(t)`$, this problem can be formally reduced to a single-particle problem as follows. Let $`Q(\stackrel{}{r}_i,t)`$ be the probability that the $`i`$-th random walker, starting from the initial position $`\stackrel{}{r}_i`$, does not meet the trajectory $`\stackrel{}{R}_0(t)`$ up to time $`t`$. Note that $`Q(\stackrel{}{r}_i,t)`$ depends implicitly on the trajectory $`\stackrel{}{R}_0(t)`$. The initial position $`\stackrel{}{r}_i`$ of the $`i`$-th walker is a random variable distributed over a volume $`V`$ of a $`d`$-dimensional space with uniform probability density $`1/V`$. We assume that there are $`N`$ such walkers. Eventually we will take the limit $`N\mathrm{}`$, $`V\mathrm{}`$ keeping the density $`\rho =N/V`$ fixed. Since the random walkers are independent, one immediately gets $`P_S(t)=_{i=1}^NQ(\stackrel{}{r}_i,t)`$ where the angle brackets denote the average over the initial positions $`\stackrel{}{r}_i`$ of the random walkers. We next write $`Q(\stackrel{}{r}_i,t)=1P(\stackrel{}{r}_i,t)`$ and use the product measure of the initial condition to write, $`P_S(t)=\left[1\frac{1}{V}P(\stackrel{}{r},t)𝑑\stackrel{}{r}\right]^N`$. Taking the $`V\mathrm{}`$ limit at fixed $`\rho `$ gives

$$P_S(t)=\mathrm{exp}\left[\rho P(\stackrel{}{r},t)𝑑\stackrel{}{r}\right]=\mathrm{exp}[\mu (t)],$$
(1)
where $`\mu (t)=\rho P(\stackrel{}{r},t)𝑑\stackrel{}{r}`$ and $`P(\stackrel{}{r},t)`$ is the probability that a single random walker starting at $`\stackrel{}{r}`$ hits the trajectory $`\stackrel{}{R}_0(t)`$ before time $`t`$. Note that the calculation of $`P(\stackrel{}{r},t)`$ and hence that of $`\mu (t)`$ is, in general, very hard even for this single particle problem due to the moving boundary $`\stackrel{}{R}_0(t)`$. Such problems, generally termed Stefan problems, are well known to be formidable to solve. In some simple cases one can derive the asymptotic behaviors of $`P(\stackrel{}{r},t)`$ for large $`t`$. Unfortunately, this asymptotic form of $`P(\stackrel{}{r},t)`$ can not be used to perform the integral over $`\stackrel{}{r}`$ in Eq. (2) since this integral usually diverges. One therefore needs an alternative approach.

Fortunately, a technique for computing $`\mu (t)`$ for a general trajectory $`\stackrel{}{R}_0(t)`$ has recently been developed . It was shown that, for $`d<2`$, $`\mu (t)`$ satisfies the exact integral equation

$$\rho =_0^t𝑑t^{}\dot{\mu }(t^{})G(\stackrel{}{R}_0(t),t|\stackrel{}{R}_0(t^{}),t^{}),$$
(2)
where $`\dot{\mu }=d\mu /dt`$ and $`G(\stackrel{}{R}_0(t),t|\stackrel{}{R}_0(t^{}),t^{})=[4\pi D(tt^{})]^{d/2}\mathrm{exp}\{[\stackrel{}{R}_0(t)\stackrel{}{R}_0(t^{})]^2/4D(tt^{})\}`$ is the standard $`d`$-dimensional diffusion propagator. For the marginal dimension $`d=2`$, Eq. (2) is still valid, though one has to introduce an ultraviolet cut-off in the diffusion propagator. The derivation of this integral equation has been detailed in Ref. . Note that for continuous space with $`d>2`$, if the tracer is a point particle then the random walkers will never meet the point particle trajectory. This is, however, not true on a lattice. So, if one sticks to a continuous space, the problem is sensible for $`d>2`$ provided the tracer particle has a finite size. The formula in Eq. (2) assumes a point particle trajectory and hence is valid only for $`d2`$. For $`d>2`$, there is no equivalent formula and one has to use other methods (see later in the paper).

For a general trajectory $`\stackrel{}{R}_0(t)`$ it is not easy to invert the integral equation (2) to obtain an explicit expression for $`\mu (t)`$. However, the advantage of Eq. (2) is that in many cases it can be exploited to derive exact asymptotic results. For example, for the ‘diffusive target annihilation’ problem, Eq. (2) has been used to derive asymptotically exact bounds for the survival probability. Note, however, that in the diffusive case, one had to average over all trajectories of the tracer particle weighted with the Wiener measure. Another solvable case is in $`d=1`$ when the tracer particle moves deterministically as $`R_0(t)=c\sqrt{t}`$. For this case an exact expression for $`\mu (t)`$ for large $`t`$ was obtained.

The purpose of this paper is to present another solvable case where the tracer particle moves ballistically through the system. For $`d2`$, we will consider the tracer to be a point particle with a trajectory $`\stackrel{}{R}_0(t)=ct\widehat{𝐳}`$, where $`c`$ denotes the velocity of the particle and $`\widehat{𝐳}`$ denotes the unit vector along the direction of motion, which we choose to be the $`z`$ axis. For $`d>2`$, we will consider the tracer particle to be a ball of finite radius $`a`$. For $`d2`$, we show how Eq. (2) can be exploited to derive an explicit exact result for $`P_S(t)`$ for all $`t`$. For $`d>2`$, we use a different method. We will present an explicit result for $`\mu (t)`$ in the physically relevant dimension $`d=3`$, though our method can, in principle, be used for any $`d>2`$. Note that this problem has an alternative physical description. In terms of the relative coordinates $`\stackrel{}{R}_i=\stackrel{}{r}_i\stackrel{}{R}_0`$, it represents a system of noninteracting particles diffusing in the presence of a constant drift velocity in the negative $`z`$ direction. The physical origin of this drift could be, for example, an external field such as gravity or an electric field. The survival probability $`P_S(t)`$ of the tracer particle, in this alternative formulation, is simply the probability that none of the particles hits the origin up to time $`t`$.

We consider first the case $`d<2`$. The substitution $`\stackrel{}{R}_0(t)=ct\widehat{𝐳}`$ in Eq. (2) reduces it to a convolution form,

$$\rho =_0^t𝑑t^{}\dot{\mu }(t^{})\frac{\mathrm{exp}\left[c^2(tt^{})/4D\right]}{\left[4\pi D(tt^{})\right]^{d/2}},$$
(3)
which can, subsequently, be solved by the Laplace transform method. We denote $`\alpha =c^2/4D`$ and define the Laplace transform $`\stackrel{~}{\mu }(s)=_0^{\mathrm{}}\mu (t)e^{st}𝑑t`$. Taking Laplace transforms on both sides of Eq. (3) and using $`\mu (0)=0`$, we obtain $`\stackrel{~}{\mu }(s)=A\frac{(\alpha +s)^{1d/2}}{s^2}`$, where $`A=\rho (4\pi D)^{d/2}/\mathrm{\Gamma }(1d/2)`$ is a constant. To invert the Laplace transform, we write $`\stackrel{~}{\mu }(s)=A(1/s+\alpha /s^2)(\alpha +s)^{d/2}`$. The inverse Laplace transform of the first factor, $`(1/s+\alpha /s^2)`$, is simply $`(1+\alpha t)`$ and that of the second factor, $`(\alpha +s)^{d/2}`$, can be easily found to be $`e^{\alpha t}t^{d/21}/\mathrm{\Gamma }(d/2)`$. We then use the convolution theorem again to write $`\mu (t)`$ as

$$\mu (t)=\frac{A}{\mathrm{\Gamma }(d/2)}_0^t\left[1+\alpha (tt_1)\right]e^{\alpha t_1}t_{1}^{}{}_{}{}^{d/21}𝑑t_1.$$
(4)
The integral on the right-hand side can be expressed in closed form as

$$\mu (t)=B\left[(1+\alpha t)\gamma (d/2,\alpha t)\gamma (d/2+1,\alpha t)\right],$$
(5)
where $`B=\rho (4\pi D)^{d/2}\mathrm{sin}(\pi d/2)/[\pi \alpha ^{d/2}]`$ and $`\gamma (\nu ,x)=_0^x\mathrm{exp}(y)y^{\nu 1}𝑑y`$ is the incomplete Gamma function. Note that the exact expression for $`\mu (t)`$ in Eq. (5) is valid for all $`t`$ and $`d<2`$.

Let us consider some limiting cases of the general result in Eq. (5). In the limit when the velocity $`c0`$ at fixed $`t`$, our problem reduces to the static ‘target annihilation’ problem for which exact results are already available. Expanding Eq. (5) for small $`\alpha =c^2/4D`$ we find, to leading order, $`\mu (t)=A_1t^{d/2}`$ where $`A_1=2\rho (4\pi D)^{d/2}\mathrm{sin}(\pi d/2)/\pi d`$. This indicates a stretched exponential decay for the survival probability $`P_S(t)=\mathrm{exp}(A_1t^{d/2})`$ and the constant $`A_1`$ matches exactly that of the static case derived earlier. Note that the limit $`c0`$ for fixed $`t`$ is equivalent to the limit $`t0`$ with $`c`$ fixed, since in Eq. (5) the time $`t`$ appears only in the scaling combination $`\alpha t`$. Thus the result $`P_S(t)=\mathrm{exp}(A_1t^{d/2})`$ also holds for small $`t`$ with $`c`$ fixed. We now consider the opposite limit of $`t\mathrm{}`$ at fixed $`c`$. In this case we find, from Eq. (5), $`\mu (t)A_2t`$ as $`t\mathrm{}`$ where $`A_2=\rho (4\pi D)^{d/2}\mathrm{sin}(\pi d/2)\mathrm{\Gamma }(d/2)/[\pi \alpha ^{(d2)/2}]`$. This indicates an exponential decay for the survival probability at late times, $`P_S(t)\mathrm{exp}(\theta t)`$ where the decay exponent $`\theta =A_2`$ is given by the following exact expression

$$\theta =\rho \pi ^{d/21}(4D)^{d1}\mathrm{sin}(\pi d/2)\mathrm{\Gamma }(d/2)c^{2d},$$
(6)
as a function of the $`4`$ physical parameters $`\rho `$, $`c`$, $`D`$ and $`d`$.

The marginal dimension $`d=2`$ is a special case. The integral equation (2) is still valid provided one introduces an ultraviolet cut-off reflecting the necessity of a lattice structure. Alternatively, a short time cut-off $`t_0`$ can be introduced in the diffusion propagator: $`G(\stackrel{}{R}_0(t),t|\stackrel{}{R}_0(t^{}),t^{})=\mathrm{exp}\{[\stackrel{}{R}_0(t)\stackrel{}{R}(t^{})]^2/4D(tt^{}+t_0)\}/[4\pi D(tt^{}+t_0)]`$, where we consider $`t_0`$ to be small. To extract the leading asymptotic behavior, one can put $`t_0=0`$ inside the exponential and retain it only in the denominator of the propagator. Using this propagator in Eq. (2), substituting $`\stackrel{}{R}_0(t)=ct\widehat{𝐳}`$ and then taking the Laplace transform as for $`d<2`$ gives $`\stackrel{~}{\mu }(s)=4\pi \rho D/[s^2\stackrel{~}{g}(s)]`$ where $`\stackrel{~}{g}(s)=_0^{\mathrm{}}𝑑t\mathrm{exp}[(\alpha +s)t]/(t+t_0)`$. Unlike the $`d<2`$ case, it is now difficult to invert the Laplace transform. However, the large-$`t`$ behavior of $`\mu (t)`$ can be easily extracted from the $`s0`$ behavior of the Laplace transform. As $`s0`$, $`\stackrel{~}{g}(s)F(\alpha t_0)`$ where $`F(x)=_0^{\mathrm{}}𝑑u\mathrm{exp}(u)/(u+x)`$. Inverting $`\stackrel{~}{\mu }(s)`$, one then gets $`\mu (t)[4\pi \rho D/F(\alpha t_0)]t`$ as $`t\mathrm{}`$. In the limit $`t_01/\alpha `$, one gets $`F(\alpha t_0)\mathrm{log}(\alpha t_0)`$. Thus, the survival probability again decays exponentially for large $`t`$, $`P_S(t)=\mathrm{exp}(\theta t)`$, where $`\theta `$ is now nonuniversal, $`\theta =4\pi \rho D/[\mathrm{log}(\alpha t_0)]`$.

We turn now to the case $`d>2`$. We consider a spherical tracer particle of radius $`a`$ moving ballistically with constant velocity $`c`$ in the $`z`$ direction. Unfortunately, for $`d>2`$ we do not yet have an analog of Eq. (2). Thus one has to resort to the original single particle formulation in Eq. (1). Fortunately, for the ballistic case in $`d>2`$, the exact asymptotic behavior of $`P_S(t)`$ for large $`t`$ can be derived even within this formulation. It turns out to be advantageous in this case to consider the alternative description of the problem in terms of the relative coordinates $`R_i(t)=r_i(t)\stackrel{}{R}_0(t)`$ where the traps diffuse independently in the presence of an external drift along the negative $`z`$ direction. Upon shifting to this relative coordinate, one finds from Eq. (1),

$$P_S(t)=\mathrm{exp}\left[\rho P(\stackrel{}{R},t)𝑑\stackrel{}{R}\right]=\mathrm{exp}[\mu (t)],$$
(7)
where $`P(\stackrel{}{R},t)=1Q(\stackrel{}{R},t)`$ and $`Q(\stackrel{}{R},t)`$ is the probability that a trap, diffusing in the presence of an external drift and starting at the initial position $`\stackrel{}{R}`$ outside the sphere of radius $`a`$, does not hit the surface of the sphere before time $`t`$. The integral in Eq. (7) is now restricted to the region $`|\stackrel{}{R}|a`$.

It is easy to see that the probability $`Q(\stackrel{}{R},t)`$ satisfies a backward Fokker-Planck equation,

$$\frac{Q}{t}=D^2Qc\widehat{𝐳}\stackrel{}{}Q,$$
(8)
where the first term on the right-hand side represents the diffusion and the second term represents the drift. This is a backward equation since we are varying the initial position $`\stackrel{}{R}`$ of the trap. Eq. (8) holds in the region $`|\stackrel{}{R}|a`$ with the boundary conditions $`Q(\stackrel{}{R},t)=0`$ for $`|\stackrel{}{R}|=a`$ and $`Q(\stackrel{}{R},t)1`$ as $`|\stackrel{}{R}|\mathrm{}`$. This is because if the particle starts at the surface of the sphere at $`t=0`$, the probability that it does not hit the surface before time $`t`$ vanishes for all $`t>0`$. Similarly, if the particle starts at infinity, with probability $`1`$ it will not hit the sphere in any finite time $`t`$. Evidently the probability $`P(\stackrel{}{R},t)=1Q(\stackrel{}{R},t)`$ also satisfies the same backward Fokker-Planck equation in (8) (with $`Q`$ replaced by $`P`$) but with a reversal of boundary conditions, $`P(\stackrel{}{R},t)=1`$ for $`|\stackrel{}{R}|=a`$ and $`P(\stackrel{}{R},t)0`$ as $`|\stackrel{}{R}|\mathrm{}`$.

Note that this Fokker-Planck equation can be written as a continuity equation, $`_tP+\stackrel{}{}\stackrel{}{J}=0`$ with a current $`\stackrel{}{J}=D\stackrel{}{}P+c\widehat{𝐳}P`$. Before solving this equation, we first make a simple observation. We have, from Eq. (7), $`\mu (t)=\rho _{|\stackrel{}{R}|a}P(\stackrel{}{R},t)𝑑\stackrel{}{R}`$. Therefore, $`\dot{\mu }(t)=_{|\stackrel{}{R}|a}_tPd\stackrel{}{R}`$. Using the continuity equation and Gauss’s divergence theorem, we obtain

$$\dot{\mu }(t)=\rho D_{|\stackrel{}{R}|=a}\stackrel{}{}P𝑑\stackrel{}{S},$$
(9)
where the surface integral is over the sphere of radius $`a`$. In deriving Eq. (9) we have used the boundary conditions on $`P`$ and the identity $`_{|\stackrel{}{R}|=a}\widehat{𝐳}𝑑\stackrel{}{S}=0`$.

It turns out that, for $`d>2`$, Eq. (8) has a stationary solution. This is because the particle always has a finite probability to escape the sphere for $`d>2`$. Thus, in order to extract the leading asymptotic behavior of $`\mu (t)`$ for $`t\mathrm{}`$, one can replace $`P(\stackrel{}{R},t)`$ by its stationary solution $`P_{st}(\stackrel{}{R})`$ on the right hand side of Eq. (9). A subsequent integration over $`t`$, using $`\mu (0)=0`$, shows that, to leading order for large $`t`$, $`\mu (t)=\theta t`$ where

$$\theta =\rho D_{|\stackrel{}{R}|=a}\stackrel{}{}P_{st}𝑑\stackrel{}{S}.$$
(10)
Note that while Eq. (9) is valid for all $`d`$, the result in Eq. (10) is valid only for $`d>2`$. This is because the trick of replacing $`P(\stackrel{}{R},t)`$ by its stationary solution $`P_{st}(\stackrel{}{R})`$ does not work for $`d2`$ since there is no stationary solution in that case.

Thus for $`d>2`$ the survival probability also decays exponentially for large $`t`$, $`P_S(t)\mathrm{exp}(\theta t)`$ with the exponent $`\theta `$ given by the general formula in Eq. (10). Obtaining an explicit expression for $`\theta `$ requires knowledge of the stationary solution of Eq. (8) which we now provide for the physically relevant dimension $`d=3`$. The stationary solution satisfies the equation

$$D^2P_{st}=c\widehat{𝐳}\stackrel{}{}P_{st},$$
(11)
with the boundary conditions $`P_{st}(\stackrel{}{R})=1`$ for $`|\stackrel{}{R}|=a`$ and $`P_{st}(\stackrel{}{R})0`$ as $`|\stackrel{}{R}|\mathrm{}`$. This problem does not have a radial symmetry. Fortunately, the substitution $`P_{st}(\stackrel{}{R})=\mathrm{exp}(\beta z)\psi (\stackrel{}{R})`$, with $`\beta =c/2D`$, restores the radial symmetry since $`\psi (\stackrel{}{R})`$ satisfies the Poisson equation $`^2\psi =\beta ^2\psi `$. The general solution satisfying the boundary condition at infinity can be obtained by standard techniques, to give

$$P_{st}(\stackrel{}{R})=R^{1/2}e^{\beta R\mathrm{cos}\varphi }\underset{l=0}{\overset{\mathrm{}}{}}b_lP_l(\mathrm{cos}\varphi )K_{l+\frac{1}{2}}(\beta R),$$
(12)
where $`R=|\stackrel{}{R}|`$, $`\varphi `$ is the angle between $`\stackrel{}{R}`$ and the $`z`$ axis, $`P_l(x)`$ is the Legendre polynomial of degree $`l`$ and $`K_\nu (x)`$ is the modified Bessel function of index $`\nu `$. The unknown coefficients $`b_l`$ are determined from the boundary condition $`P_{st}(R=a)=1`$. Substituting $`R=a`$ in Eq. (12) and using the orthogonality properties of the functions $`P_l(x)`$ gives, after some algebra,

$$P_{st}(\stackrel{}{R})=\sqrt{\frac{a}{R}}e^{\beta R\mathrm{cos}\varphi }\underset{l=0}{\overset{\mathrm{}}{}}(l+\frac{1}{2})a_lP_l(\mathrm{cos}\varphi )\frac{K_{l+\frac{1}{2}}(\beta R)}{K_{l+\frac{1}{2}}(\beta a)},$$
(13)
where $`a_l=_1^1P_l(x)e^{\beta ax}𝑑x=(1)^l\sqrt{2\pi /\beta a}I_{l+1/2}(\beta a)`$.

Substituting the stationary solution, Eq. (13), into Eq. (10) and performing the surface integral, we finally obtain, after a few steps of algebra, the following rather nontrivial expression for $`\theta `$ in terms of the physical parameters $`\rho `$, $`D`$, $`a`$ and $`c`$:

$$\theta =2\pi a\rho D\left[12\pi \underset{l=0}{\overset{\mathrm{}}{}}(1)^l(l+\frac{1}{2})\frac{K_{l+\frac{1}{2}}^{}(\beta a)}{K_{l+\frac{1}{2}}(\beta a)}I_{l+\frac{1}{2}}^2(\beta a)\right],$$
(14)
where $`K_\nu ^{}(x)=dK_\nu (x)/dx`$ and $`\beta =c/2D`$. The series in Eq. (14) can be summed numerically. If $`H(\beta a)`$ is the function in the square brackets, then $`H(x)`$ is monotically increasing, with $`H(0)=2`$ and $`H(x)/x1`$ for $`x\mathrm{}`$. The former result, corresponding to $`c=0`$, recovers the known result, $`\theta =4\pi a\rho D`$, for a static target . The latter, corresponding to $`c\mathrm{}`$, can be understood by noting that in this limit the probability that the sphere has not been hit by a trap is given by $`\mathrm{exp}(\rho V)`$, where $`V=\pi a^2ct`$ is the volume swept out by the sphere in time $`t`$ and we require that this volume initially contains no traps. Hence $`\theta \pi a^2c\rho `$ in this limit, corresponding to $`H(x)x`$.

In summary, we have studied the general problem of calculating the survival probability of a tracer particle moving along a deterministic trajectory in the presence of diffusing traps. In particular, when the tracer particle moves ballistically we have shown that its survival probability $`P_S(t)\mathrm{exp}(\theta t)`$ for large $`t`$ in all dimensions. We have derived exact expressions for the exponent $`\theta `$ in terms of the system parameters for $`d2`$ and for $`d=3`$.


text:
# Untitled Document

The paper is being withdrawn since the results are incorporated in paper arxiv.org/abs/math.AG/0306195.


text:
# Polarization of the prompt 𝛾-ray emission from the 𝛾-ray burst of 6 December 2002

(Nature, 2003 May 22, 423, 415-417)

Space Sciences Laboratory, and
Department of Phyisics, University of California at Berkeley,
Berkeley, CA, 94720, USA

Observations of the afterglows of $`\gamma `$-ray bursts (GRBs) have revealed that they lie at cosmological distances, and so correspond to the release of an enormous amount of energy1; 2. The nature of the central engine that powers these events and the prompt $`\gamma `$-ray emission mechanism itself remain enigmatic because, once a relativistic fireball is created, the physics of the afterglow is insensitive to the nature of the progenitor. Here we report the discovery of linear polarization in the prompt $`\gamma `$-ray emission from GRB021206, which indicates that it is synchrotron emission from relativistic electrons in a strong magnetic field. The polarization is at the theoretical maximum, which requires a uniform, large-scale magnetic field over the $`\gamma `$-ray emission region. A large-scale magnetic field constrains possible progenitors to those either having or producing organized fields. We suggest that the large magnetic energy densities in the progenitor environment (comparable to the kinetic energy densities of the fireball), combined with the large-scale structure of the field, indicate that magnetic fields drive the GRB explosion.

We used the *Reuven Ramaty High Energy Solar Spectroscopic Imager* (*RHESSI*)3 to make these$`\gamma `$-ray observations of GRB021206. *RHESSI* has an array of nine large-volume (300 cm3) coaxial germanium detectors with high spectral resolution, designed to study solar X-ray and$`\gamma `$-ray emission (3 keV–17 MeV). *RHESSI* has high angular resolution (2′′) in the $``$1 field of view of its optics; however, the focal plane detectors are unshielded, open to the whole sky. Thus, while the chances are small that *RHESSI* will see a GRB in its imaging field of view, it measures them frequently in the focal plane detectors themselves. These observations provide highresolution spectra, individual photon times and energies, as well as the potential for polarization measurements. *RHESSI* is not optimized to act as a $`\gamma `$-ray polarimeter, but several aspects of its design make it the most sensitive instrument to date for measuring astrophysical $`\gamma `$-ray polarization.

In the soft $`\gamma `$-ray range of $``$0.15–2.0 MeV, the dominant photon interaction in the *RHESSI* detectors is Compton scattering. A small fraction of incident photons will undergo a single scatter in one detector before being scattered and/or photoabsorbed in a second separate detector, which are the events sensitive to the incident $`\gamma `$-ray polarization. Linearly polarized $`\gamma `$-rays preferentially scatter in directions perpendicular to their polarization vector. In *RHESSI* this scattering property can be used to measure the intrinsic polarization of astrophysical sources. The sensitivity of an instrument to polarization is determined by its effective area to scatter events, and the average value of the polarimetric modulation factor4; 5, $`\mu (\theta ,E)`$, which is the maximum variation in azimuthal scattering probability for polarized photons. This factor is given by $`\mu =(d\sigma _{}d\sigma _{})/(d\sigma _{}+d\sigma _{})`$, where $`d\sigma _{}`$,$`d\sigma _{}`$ are the Klein-Nishina differential cross-sections for Compton scattering perpendicular and parallel to the polarization direction, respectively, which is a function of the incident photon energy $`E_\gamma `$, and the Compton scatter angle $`\theta `$ between the incident photon direction and the scattered photon direction. For a source of count rate $`S`$ and fractional polarization $`\mathrm{\Pi }_s`$, the expected azimuthal scatter angle distribution is $`dS/d\varphi =(S/2\pi )[1\mu _m\mathrm{\Pi }_s\mathrm{cos}(2(\varphi \eta ))]`$, where $`\varphi `$ is the azimuthal scatter angle, $`\eta `$ is the direction of the polarization vector, and $`\mu _m`$ is the average value of the polarimetric modulation factor for the instrument. Although *RHESSI* has a small effective area ($``$20 cm2) for events that scatter between detectors, it has a relatively large modulation factor in the 0.15–2.0 MeV range, $`\mu _m0.2`$, as determined by Monte Carlo simulations described below.

In comparison with other $`\gamma `$-ray instruments (COMPTEL, BATSE) that have attempted to measure polarization in the past5; 6, *RHESSI* has the major advantage of quickly rotating around its focal axis (centred on the Sun) with a 4 s period. Rotation averages out the effects of asymmetries in the detectors and passive materials that could be mistaken for a modulation. Because polarimetric modulations repeat every 180, any source lasting more than half a rotation (2 s) will be relatively insensitive to the systematic uncertainties that typically plague polarization measurements. Finally, although the *RHESSI* detectors have no positioning information themselves, they are relatively loosely grouped on the spacecraft, allowing the azimuthal angle for a given scatter to be determined to within $`\mathrm{\Delta }\varphi =13^{}`$ r.m.s. This angular uncertainty will decrease potential modulations by a factor of 0.95, which is included in our calculated modulation factor.

Prompt $`\gamma `$-ray emission from GRB021206 was detected with *RHESSI* on 6 December 2002 at 22:49 UT (Fig. 1). This GRB was also observed7 with the Interplanetary Network (IPN), which reported a 25–100 keV fluence of $`1.6\times 10^4`$ erg cm-2, and a peak flux of $`2.9\times 10^5`$ erg cm-2 s-1, making this an extremely bright GRB. The IPN localized8 GRB021206 to a 57 square-arcminute error box located 18 from the Sun. For *RHESSI* we analysed photons in the energy range 0.15–2.0 MeV that scattered between two, and only two, detectors for the 5 s integration period shown in Fig. 1. Scattered photons constitute roughly 10% of the total 0.15–2.0 MeV light-curve events. Counts were binned by the centre-to-centre azimuthal angle between the two detectors around the *RHESSI* roll axis, corrected for the rotation of the spacecraft at the time of the photon event. This azimuthal distribution is plotted in Fig. 2. The top panel shows the raw data, as well as the expected variation for an unpolarized GRB due only to the light-curve variability. The bottom panel shows the residual of the data once this unpolarized distribution is subtracted. The residual shows a large modulation, which we interpret as a linear polarization of $`\mathrm{\Pi }_m=(80\pm 20)`$%. This observation is the first astrophysical polarization measurement at $`\gamma `$-ray energies. Note that the uncertainty on this polarization amplitude reflects in large part our uncertainty in the modulation factor. The fact that we have measured a polarization somewhere in this range has a significance of $`<10^8`$ (or a confidence $`>5.7\sigma `$).

Linear polarization is generally considered a clear indication of synchrotron emission. For electrons with an energy spectrum characterized by a power-law distribution with spectral index $`p`$, synchrotron photons are emitted with a linear polarization9 of $`\mathrm{\Pi }=(p+1)/(p+7/3)`$. For shock acceleration10, typical values of $`p=`$2–3 correspond to linear polarizations of 70–75%. For unresolved sources, polarizations from many directions generally add together to produce net polarizations that are a fraction of this maximum value. While the source is unresolved in our observation, if the source electrons are moving with a bulk Doppler factor of $`\mathrm{\Gamma }`$, we are only viewing the source over a solid angle $`\mathrm{\Omega }_\gamma 1/\mathrm{\Gamma }^2`$. Thus, for typical values of $`\mathrm{\Gamma }300`$ that have been implied from GRB afterglow observations, we are effectively resolving a source region of solid angle $`\mathrm{\Omega }_\gamma 10^5`$ sr. Our measurement of this high polarization is consistent with synchrotron origin for the initial GRB from a region of nearly uniform magnetic field. For this emission process to be radiatively efficient as is implied by many afterglow observations11; 12, the magnetic energy densities must be close to equipartition (comparable to the kinetic energy densities of the fireball)13; 14.

This polarization from the prompt $`\gamma `$-ray emission is significantly higher than optical polarizations of 1–3% typically measured from afterglows15; 16; 17, as well as the optical polarization of 10% recently reported18 from the afterglow of GRB020405. The afterglow emission implies a strong magnetic field behind the shocks19, although the field energy density can be well below equipartition20. However, the implied magnetic field strengths are too large to have been due to a progenitor field being dragged along by the expanding fireball21, or compression of the ISM magnetic fields in external shocks22. Therefore, the magnetic fields responsible for the afterglow synchrotron emission are probably turbulent fields that have built up behind the shocks21; 23, which is consistent with the relatively small optical afterglow polarizations. This locally generated field invoked for the afterglow has influenced many researchers to consider the magnetic field responsible for the prompt burst of $`\gamma `$-rays as a turbulent, fireball-induced field as well. In the standard internal shock model24, the prompt $`\gamma `$-rays are produced by synchrotron emission of electrons accelerated to relativistic energies by shock acceleration, requiring magnetic energy densities near equipartition in the progenitor environment. However, late-time turbulent fields have no direct implications on whether the fields responsible for the prompt $`\gamma `$-ray emission are predominantly turbulent or organized21; 25. Our observation conclusively shows that the engine driving the GRB has a strong, large-scale magnetic field.

Another potential source of polarization would be as follows: if unpolarized $`\gamma `$-rays are initially beamed into a small-angle jet, and then scatter at an angle $`\theta `$ into our line of sight. The polarization induced by this Compton scattering9 is $`\mathrm{\Pi }=(1\mathrm{cos}^2\theta )/(1+\mathrm{cos}^2\theta )`$ for photon energies below $``$0.1 MeV, decreasing as the photon energy approaches and exceeds the electron rest mass, 0.511 MeV. For our energy band, we estimate a maximum Compton- induced polarization of 70% for scatter angles near 90. However, this process would be inefficient, requiring a much larger $`\gamma `$-ray energy budget than the synchrotron case. First, the distribution of scattered $`\gamma `$-rays would be nearly isotropic, requiring an energy budget $`4\pi /\mathrm{\Omega }_j`$ larger than if the observed $`\gamma `$-rays originated from a collimated jet of opening solid angle $`\mathrm{\Omega }_j`$. Second, to maintain such a high polarization as we observed, the $`\gamma `$-rays must undergo only a single scatter into our direction because secondary scatters will erase the induced polarization from the initial scatter. This condition requires that the scatter medium be optically thin, $`\tau <<1`$, and that the luminosity of the unseen initial $`\gamma `$-ray beam be larger than the observed scattered $`\gamma `$-ray emission by a factor, $`1/\tau `$. Because Compton-scattering-induced polarization requires a total $`\gamma `$-ray luminosity several orders of magnitude larger than that implied by synchrotron emission, and elaborate source geometries, the synchrotron origin for the polarization is preferred.

We suggest that our observation is evidence that the magnetic fields are actually powering the GRB explosion itself. It has been argued that a “passive” magnetic field–that is, a field dragged from the surface of central object with a magnetic dipole moment, but not driving the GRB–could not be strong enough to produce the prompt $`\gamma `$-ray emission without an additional locally generated turbulent magnetic field25. If this conclusion holds, our observation is consistent with models of a magnetically driven GRB fireball, where the driving magnetic field was generated by extracting the rotational energy of an accretion disk around a central compact object through differential rotation26; 27, by directly extracting the spin energy of a black hole threaded by magnetic field lines28, or by extracting the spin energy of a highly magnetized neutron star29. Alternatively, our observation of a large-scale magnetic field could support models of dynamos in the post-shock flows30 if the shocks can be shown to be unstable on large size scales and on timescales comparable to the prompt $`\gamma `$-ray emission.

Methods

The top panel of Fig. 2 shows the azimuthal scatter distribution of 0.15–2.0 MeV events for the 5 s integration period shown in Fig. 1, corrected for the rotation of the spacecraft at the time of each photon event. This distribution is the sum of three components: the GRB scatter event rate which averages $`820\pm 8`$ counts per bin, the chance coincidence rate which averages $`374\pm 6`$ counts per bin, and the background scatter event rate of $`49\pm 2`$ counts per bin. These rates were determined and confirmed using event rates before and after the burst, combined with studies of the readout times of multiple detectors during scattered events, and were verified independently using our Monte Carlo simulations. Although the rotation will average out systematic variations in the scatter angle distribution, we still have to correct for the complex time profile of the burst itself, which will cause variations for an unpolarized source owing to the finite number of potential scatter angles *RHESSI* can measure at any given instant. We modelled this effect by using the 0.15–2.0 MeV total count rate in the *RHESSI* instrument (Fig. 1) as the time dependent flux template for a photon transport Monte Carlo simulation, and using the time-averaged GRB photon spectrum as measured by *RHESSI* for our input spectrum. This simulation used the detailed *RHESSI* mass model that has been developed under CERN’s GEANT package, allowing us to model the instrument response to a GRB at the IPN8 sky coordinates for each rotation angle and instantaneous flux, assuming an unpolarized source. This distribution is also presented in the top panel of Fig. 2. We are looking for a modulation signal relative to this variation induced by the GRB time profile.

In the bottom panel of Fig. 2 we show the residual of the measured distribution once we have subtracted away the simulated response for an unpolarized GRB, showing our absolute modulation signal. For an unpolarized source we would expect this distribution to be flat, which we can rule out at an extremely high confidence level ($`\chi ^2=83.5`$, 11 degrees of freedom, d.f.). When we fit this with a modulation curve, the fit improves significantly ($`\chi ^2=16.9`$, 9 d.f.), with an amplitude of $`128\pm 16`$ counts per bin. Statistically, this is a reasonable fit to the data (the probability of $`\chi ^2>16.9`$ is 5%), but could be improved with even more detailed Monte Carlo simulations including timedependent spectral variability. Using the count rates given above and the simulated distribution for an unpolarized GRB, we performed further numerical simulations to determine the probability that an unpolarized GRB could produce a modulation as large as the one we measure due to random Poisson counting statistics. We found this probability to be very low, $`<10^8`$, which translates to a confidence that we have measured a polarization at a level $`>5.7\sigma `$. Finally, we estimated the modulation factor to be $`\mu _m=0.19\pm 0.04`$, using both a separate photon transport code which fully treats polarization in scattering and uses a simplified mass model, as well as analytical estimates based on the GEANT simulation with the full *RHESSI* mass model. Combining the modulation amplitude, the total source scatter event rate, and the *RHESSI* modulation factor, we derive a measured polarization $`\mathrm{\Pi }_m=(80\pm 20)`$%.

A number of tests were performed to check that themeasured modulation is real. First, we verified that the simulated variation induced by the GRB light curve is accurate by comparing it to an angular distribution of events that were chance coincidences in two detectors. These interactions are nearly simultaneous, but separated by enough time to distinguish them as chance coincidences, not real scattered photons. This distribution should exhibit the same variations owing to the GRB light curve, but no polarization. When we subtracted the simulated distribution from the chance-coincident distribution, we found no evidence for a residual modulation. We have performed a number of independent checks to make sure we do not see modulations from other sources as well. We have verified that extended *RHESSI* background observations show no sign of modulations. In addition, we have done a preliminary analysis of a strong solar $`\gamma `$-ray flare observed on 23 July 2002, where we see some evidence for a modulation, but corresponding to a polarization $`<10`$%. Therefore we feel confident that we have characterized the systematic effects in *RHESSI* to below the 10% polarization level.

Received 7 February; accepted 12 March 2003

Acknowledgements We thank D. Smith for help in learning *RHESSI* data analysis and providing simulation support, K. Hurley for IPN data and references, R. Lin, E. Quataert, J. Arons, C. Matzner and I. Fisk for discussions, and especially the *RHESSI* team for making all of their data immediately available to the public at http://rhessidatacenter.ssl.berkeley.edu.

Correspondence and requests for materials should be addressed to S.E.B. (boggs@ssl.berkeley.edu) or W.C. (wcoburn@ssl.berkeley.edu)


text:
# General unit-disk representation for periodic multilayers

## Abstract

We suggest a geometrical framework to discuss periodic layered structures in the unit disk. The band gaps appear when the point representing the system approaches the unit circle. We show that the trace of the matrix describing the basic period allows for a classification in three families of orbits with quite different properties. The laws of convergence of the iterates to the unit circle can be then considered as universal features of the reflection.

Photonic band gap structures dowling  can be dealt with by computing the band structure (using e.g. Bloch theory) or from the point of view of scattering Cen00 . The essential difference is that a scattering experiment always involves a finite structure, while Bloch waves imply an endlessly repetition of the basic period.

In the context of electromagnetic optics, photonic crystals (that is, one-dimensional periodic layered structures) have attracted recently a lot of attention because the striking property of acting as omnidirectional reflectors: they can reflect light at any polarization, any incidence angle, and over a wide range of wavelengths Fin98 ; Dow98 ; Yab98 ; Chi99 ; Lek00 .

The appearance of strong reflection (stop bands) depends on the properties of the basic period. Different theoretical approaches, involving equivalent medium theories, group velocity, and other quantities Gra00 ; Mot00 , have been developed to account for the detailed structure of these stop bands and their edges. Each one of these models emphasizes some aspects of the problem but, at the same time, has some drawbacks.

In the present Letter we introduce a geometrical setting that allows for a deeper understanding of periodic systems. Our treatment is quite general and only assumes linearity: it applies to any physical system whose transfer matrix belongs to the group SU(1, 1). The key point for our purposes is the fact that the multilayer transfer function induces a bilinear transformation in the unit disk Yon02 ; Mon02 . Since perfect mirrors are represented precisely by the unit circle, the route to a stop band can be understood as the convergence of the point representing the action of the system to the unit circle.

We start by examining the basic period of the system, which consists of a stack of plane-parallel layers sandwiched between two semi-infinite ambient ($`a`$) and substrate ($`s`$) media that we shall assume to be identical, since this is the common experimental case. Hereafter all the media are supposed to be lossless, homogeneous, and isotropic.

A monochromatic linearly polarized plane wave falls from the ambient making an angle $`\theta _0`$ with the normal to the first interface and with an amplitude $`E_a^{(+)}`$. We consider as well another plane wave of the same frequency and polarization, and with amplitude $`E_s^{()}`$, incident from the substrate at the same angle $`\theta _0`$. The output fields in the ambient and the substrate will be denoted $`E_a^{()}`$ and $`E_s^{(+)}`$, respectively.

The field amplitudes at each side of the multilayer are related by the linear relation

$$\left(\begin{array}{c}E_a^{(+)}\\ E_a^{()}\end{array}\right)=𝖬_{as}\left(\begin{array}{c}E_s^{(+)}\\ E_s^{()}\end{array}\right),$$
(1)
where the multilayer transfer matrix $`𝖬_{as}`$ can be shown to be Mon99a ; Mon99b

$$𝖬_{as}=\left[\begin{array}{cc}1/T_{as}& R_{as}^{}/T_{as}^{}\\ R_{as}/T_{as}& 1/T_{as}^{}\end{array}\right]\left[\begin{array}{cc}\alpha & \beta \\ \beta ^{}& \alpha ^{}\end{array}\right].$$
(2)
Here the complex numbers $`R_{as}`$ and $`T_{as}`$ are, respectively, the overall reflection and transmission coefficients for a wave incident from the ambient. Because $`|R_{as}|^2+|T_{as}|^2=1`$, we have $`det𝖬_{as}=+1`$ and then the set of lossless multilayer matrices reduces to the group SU(1, 1).

We are often interested in the transformation properties of field quotients rather than the fields themselves. Therefore, we introduce the complex numbers

$$z=\frac{E^{()}}{E^{(+)}},$$
(3)
for both ambient and substrate. Equation (1) defines then a transformation on the complex plane $``$, mapping the point $`z_s`$ into the point $`z_a`$, according to

$$z_a=\mathrm{\Phi }[𝖬_{as},z_s]=\frac{\beta ^{}+\alpha ^{}z_s}{\alpha +\beta z_s},$$
(4)
which is a bilinear (or Möbius) transformation. This action can be seen as a function $`z_a=f(z_s)`$ that will be called the multilayer transfer function. One can check that the unit disk, the external region and the unit circle remain invariant under the multilayer action Yon02 . Note that $`|z_a|=|z_s|=1`$ for a perfect mirror with light incident from both the ambient and the substrate, so the transformation (4) relates then points on the unit circle. When no light strikes from the substrate $`z_s=0`$ and $`|z_a|=1`$, so a mirror maps the origin into a point on the unit circle.

In what follows, the idea of fixed points of the transformation will prove to be essential. These points can be defined as the field configurations such that $`z_a=z_sz_f`$ in Eq. (4); i.e., $`z_f=\mathrm{\Phi }[𝖬_{as},z_f]`$, whose solutions are

$$z_{f\pm }=\frac{1}{2\beta }\left\{2i\mathrm{Im}(\alpha )\pm \sqrt{[\mathrm{Tr}(𝖬_{as})]^24}\right\}.$$
(5)
The trace of $`𝖬_{as}`$ provides then a suitable tool for the classification of multilayers San01 .

When $`[\mathrm{Tr}(𝖬_{as})]^2<4`$ the multilayer action is elliptic and it has only one fixed point inside the unit disk, while the other lies outside. When $`[\mathrm{Tr}(𝖬_{as})]^2>4`$ the action is hyperbolic and it has two fixed points on the unit circle. Finally, when $`[\mathrm{Tr}(𝖬_{as})]^2=4`$ the multilayer action is parabolic and it has only one (double) fixed point on the unit circle.

To proceed further let us note that by taking the conjugate of $`𝖬_{as}`$ with any matrix $`𝖢`$ SU(1, 1), that is

$$\widehat{𝖬}_{as}=𝖢𝖬_{as}𝖢^1,$$
(6)
we obtain another matrix of the same type, since $`\mathrm{Tr}(\widehat{𝖬}_{as})=\mathrm{Tr}(𝖬_{as})`$. Conversely, if two multilayer matrices have the same trace, a matrix $`𝖢`$ satisfying Eq. (6) can be always found.

The fixed points of $`\widehat{𝖬}_{as}`$ are then the image by $`𝖢`$ of the fixed points of $`𝖬_{as}`$. In consequence, given any matrix $`𝖬_{as}`$, it can always be reduced to a unique $`\widehat{𝖬}_{as}`$ with one of the following canonical forms:

$`\widehat{𝖪}(\phi )`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{exp}(i\phi /2)& 0\\ 0& \mathrm{exp}(i\phi /2)\end{array}\right],`$ (9)
$`\widehat{𝖠}(\chi )`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{cosh}(\chi /2)& i\mathrm{sinh}(\chi /2)\\ i\mathrm{sinh}(\chi /2)& \mathrm{cosh}(\chi /2)\end{array}\right],`$ (12)
$`\widehat{𝖭}(\eta )`$ $`=`$ $`\left[\begin{array}{cc}1i\eta /2& \eta /2\\ \eta /2& 1+i\eta /2\end{array}\right],`$ (15)
that have as fixed points the origin (elliptic), $`+i`$ and $`i`$ (hyperbolic) and $`+i`$ (parabolic), and whose physical significance has been studied before Mon01 .

The concept of orbit is especially appropriate for obtaining a picture of these actions. Given a point $`z`$, its orbit is the set of points $`z^{}`$ obtained from $`z`$ by the action of all the elements of the group. In Fig. 1.a we have plotted typical orbits for each one of the canonical forms. For matrices $`\widehat{𝖪}(\phi )`$ the orbits are circumferences centered at the origin. For $`\widehat{𝖠}(\chi )`$, they are arcs of circumference going from the point $`+i`$ to the point $`i`$ through $`z`$. Finally, for $`\widehat{𝖭}(\eta )`$ the orbits are circumferences passing through the point $`+i`$ and joining the points $`z`$ and $`z^{}`$.

For an arbitrary matrix $`𝖬_{as}`$ the orbits can be obtained by transforming the orbits described before with the appropriate matrix $`𝖢`$ . In Fig. 1.b we have plotted typical examples of such orbits for elliptic, hyperbolic, and parabolic actions. We stress that once the fixed points of the matrix $`𝖬_{as}`$ are known, one can ensure that $`z_a`$ will lie in the orbit associated to $`z_s`$.

Assume now that we have a finite periodic structure obtained by repeating $`N`$ times the basic period represented by $`𝖬_{as}`$. The overall transfer matrix for this system is $`𝖬_{as}^N`$. In the unit-disk picture, the transformed field by the $`N`$-period structure is represented by the point

$$z_N=\mathrm{\Phi }[𝖬_{as},z_{N1}]=\mathrm{\Phi }[𝖬_{as}^N,z_0],$$
(16)
where $`z_0`$ is the initial point $`z_s`$. The idea that maps iterates could be applied to this problem was recognized before, though with a somewhat different approach Fel00 ; Fel01 .

Henceforth, we shall take $`z_0=0`$, which is not a serious restriction and corresponds to the case in which no light incides from the substrate $`[E_s^{()}=0]`$, as it happens usually. Note also that all the points $`z_N`$ lie in the orbit associated to the initial point $`z_0`$ by the single period, which is determined by its fixed points: the character of these fixed points determine thus the behavior of the periodic structure.

To illustrate how this geometrical approach works in practice, we take the single period as a Fabry–Perot-like system formed by two identical plates separated by a spacer of phase thickness $`\delta _2`$. This is a symmetric system for which $`R_{as}`$ and $`T_{as}`$ can be easily computed Mon02 . By varying $`\delta _2`$ we can choose to work in the elliptic, the hyperbolic, or the parabolic case. In Fig. 2 we have plotted the sequence of successive iterates obtained numerically for these three regimes.

In the elliptic case, it is clear that the points $`z_N`$ revolve in the orbit centered at the fixed point and the system never reaches the unit circle.

On the contrary, for the hyperbolic and parabolic cases the iterates converge to one of the fixed points on the unit circle, although with different laws, which correspond to the band stop and band edges of the system, respectively Lek00 .

To gain further insights into these behaviors, we compute explicitly the $`N`$th iterate. This can be easily done for the canonical forms in Eq. (9) and then, conjugating as in (6) we obtain, after some lengthy calculations, that for a hyperbolic action one has

$$z_N=\frac{1\xi ^N}{1\xi ^N(z_{f+}/z_f)}z_{f+},$$
(17)
where $`\xi =(\alpha +\beta z_f)/(\alpha +\beta z_{f+})`$ is a complex number satisfying $`|\xi |<1`$. Analogously, for the parabolic case we have

$$z_N=\frac{N\beta z_f^2}{N\beta z_f1},$$
(18)
where $`z_f`$ is the (double) fixed point. It is quite obvious that in both cases $`z_N`$ converges to one of the fixed points on the unit circle, so $`|z_N|1`$ when $`N`$ increases, a typical behavior of perfect mirror. In the mathematical literature this limit point is referred to as the Denjoy-Wolff point of the map Kap99 .

To characterize the convergence of $`z_N`$ we note that, because $`z_0=z_s=0`$, this initial point is transformed by the single period into $`z_a=R_{as}`$. Therefore, $`z_N`$ represents the reflection coefficient of the overall periodic structure $`R_{as}^{(N)}`$, which is obviously different from $`(R_{as})^N`$. One can then compute that for the hyperbolic case Mon03

$$|z_N|^2=\frac{|\beta |^2}{|\beta |^2+[\mathrm{sinh}(\chi )/\mathrm{sinh}(N\chi )]^2},$$
(19)
that approaches the unit circle exponentially with $`N`$, as one could expect from a band stop, while for the parabolic action

$$|z_N|^2=\frac{|\beta |^2}{|\beta |^2+(1/N)^2},$$
(20)
that goes to unity with a typical behavior $`O(N^2)`$. This is universal in the physics of reflection, as put forward in a different framework by Yeh Yeh88  and Lekner Lek87 .

To conclude, we expect that the geometrical scenario presented here could provide an appropriate tool for analyzing and classifying the performance of periodic multilayers in an elegant and concise way that, additionally, is wider enough to accommodate other periodic systems appearing in physics.

We acknowledge José F. Cariñena and José M. Montesinos for illuminating discussions.