text:
# Optical polarization observations of NGC 6231: evidence for a past SN fingerprintBased on observations obtained at Complejo Astronómico El Leoncito (CASLEO), operated under agreement between the CONICET and the National Universities of La Plata, Córdoba, and San Juan, Argentina.

## 1 Introduction

The polarimetric techniques are a very useful tool to find important information ($`P_{\lambda max}`$,$`\lambda _{max}`$, magnetic field direction, etc) from the dust located in front of a star (or a luminous object), or from dust located inside a stellar cluster. Young open clusters are very good candidates to carry out polarimetric observations, because previous photometric and spectroscopic studies of these clusters give detailed information of the main sequence, so we can characterize the physical parameters of the member stars (age, distance, memberships, etc), in order to study the extinction of the dust in the direction to the cluster and within the cluster. For example, from the study of the $`P_{\lambda max}`$ versus the $`E_{BV}`$ of each star of a cluster it is possible to find the polarimetric components of the extinction produced by Galactic dust located between the Sun and the cluster (e.g. Tr 27, Feinstein et al. (2000)). By subtracting the effect of this dust on the line of sight from the data it is possible to study the component associated with the internal extinction of the cluster. In addition to this valuable information, stars with intrinsic polarimetric properties can be isolated. The angle of the polarization vector is an estimator of the direction of the Galactic magnetic field. By observing members and nonmembers stars (foreground and background) the direction of the magnetic field towards the cluster at several distances can be inferred.

In a typical open cluster the polarimetric data show, for “normal” stars (that is, not considering stars with intrinsic polarization), a distribution of angles that can be described by a mean fixed value which characterizes the cluster. The polarization angles of the cluster members show some scatter ($`10^{}`$ or less) over this mean angle, mainly because of density inhomogeneities in the dust along the line of sight to the cluster, inhomogeneities of the dust inside the cluster, to the intrinsic polarimetric properties of the stars and also due to more complex magnetic field structures. For example, in Tr 27 the polarimetric data show a mean angle of $`\theta 30^{}`$ for the cluster with a FWHM of $`10^{}`$ (Feinstein et al., 2000), in Ara OB1 the polarimetric observations give a FWHM of $`9^{}.9`$ for NGC 6167, $`13^{}.1`$ for NGC 6193 and $`9^{}.5`$ for NGC 6204 (Waldhausen et al., 1999). In the case of Stock 16 (Feinstein et al., 2003) the FWHM is 6.7, which seems to be slightly smaller than the other clusters.

NGC 6231 (l= 345.5 , b=1.2) is known in the literature to be one of the youngest open clusters in the Galaxy and one of the brightest in the Southern Hemisphere. The cluster is embedded in an extended system of O and early-B stars, namely the very young association Sco OB1, that contains the ring-shaped HII region IC4678, centered on NGC 6231. Usually, NGC 6231 is considered as the nucleus of Sco OB1 (Perry, Hill & Christodoulou, 1991). Its age determinations range between $`3\times 10^6`$ years (Mermilliot & Maeder, 1986) and $`6.9\times 10^6`$ years (Perry, Hill & Christodoulou, 1991). It is known as one of the richest open clusters, containing WR stars, $`\beta `$ Cepheid variables and more than 100 O and B stars (Shobbrook, 1983). According to photometric CCD observations of Sung et al. (1998) and of Baume, Vázquez & Feinstein (1999) there is evidence of a population of pre-main sequence stars (PMS). The distance obtained by deep CCD photometry by Baume, Vázquez & Feinstein (1999) locates the cluster NGC 6231 at $`1990\pm 200`$ pc from the Sun.

The average visual extinction toward the cluster is about $`A_V1.4`$ (Balona & Laney, 1995; Baume, Vázquez & Feinstein, 1999). According to Raboud, Cramer & Bernasconi (1997) in the northern part of the cluster the absorption is constant, but it increases from the central part towards the South. A similar trend, but less pronounced, is observed from West to East in the southern part. Part of the absorption suffered by the light from the cluster stars is due to matter within the distance interval 100-1300 pc from the Sun (Perry, Hill & Christodoulou, 1991), but the field shows a relatively good transparency between about 400 and 1000 pc at a constant value averaging of $`A_V=0.75`$ mag (Raboud et al. 1997).

## 2 Observations

The bulk of the observations were obtained on two observing runs. The first one was performed during 1995 (May 30 to June 3). The second one was on 3 nights (August 17-20) in 2001. Two more runs on May 2002 and on August 2002 were used for re-observing and checking some stars. All runs were carried out using the five-channel photopolarimeter of the Torino Astronomical Observatory (Scaltriti et al., 1989) attached to the 2.15 m telescope at the Complejo Astronómico El Leoncito (San Juan, Argentina). Polarimetric measurements have been made for 36 stars to the cluster NGC 6231. Most of them were observed through the Johnson broad band UBVRI filters ($`\lambda _{Ueff}=0.360\mu m,\lambda _{Beff}=0.440\mu m,\lambda _{Veff}=0.530\mu m,\lambda _{Reff}=0.690\mu m,\lambda _{Ieff}=0.830\mu m`$).

Our results are listed in Table 1 which shows, in self explanatory format, the stellar identification as given by Seggewiss (1968), the average of the percentage of polarization ($`P_\lambda `$) and the position angle of the electric vector ($`\theta _\lambda `$) observed through each filter, with their respective mean errors.

The Torino polarimeter collects data for all the filters simultaneously (UVBRI), so all the wavebands have the same exposure times, and the S/N varies between filters. Therefore, the values from different filters may be of different quality. Observations whose values are not over the 3$`\sigma `$ error level were discarded, and they are not reported in Table 1 neither used in this work.

Several standard stars for null polarization and for the zero point of the polarization position angle were measured for calibration purposes. At least three standard of angle stars were measured every night at different times. Errors were handled as described by Maronna, Feinstein and Clocchiatti (1992).

In order to verify the lack of systematic differences between both main runs we have re-observed a common group of stars in the second run.

Three other bright stars have been observed previously by Serkowski et al. (1975): these are HD 152234 (star # 290), HD 152235 and HD 152248. Although in that paper no observation errors were given, their data (polarization and angle) are very similar to those reported in the present paper (see Table 2).

To look for nearby objects, the observed stars were checked with the Hipparcos/Tycho Catalog Data. Stars #80,110,112,220,290 and HD 152248 were included in these catalogs but are so far away from the Sun that no useful parallax measures could be obtained.

## 3 Results

Figure 1 shows the sky projection of the V-band polarization for the stars in NGC 6231. The dot-dashed line indicates the Galactic parallel $`b=1^{}.1`$. The first striking feature of this figure is the abnormal distribution of the polarization angles. Commonly in a typical open cluster the polarization vectors have a similar mean orientation of the polarization angle, and the individual stars show little scatter over this fixed direction. Previous studies for 5 clusters indicate that this scatter is equal or less than 10 degrees of FWHM (Waldhausen et al., 1999; Feinstein et al., 2000, 2003). As seen in Fig. 1, this is not the case in NGC 6231.

It is possible to classify the stars of NGC 6231 in three groups according to their location and the angle of the polarization vectors as shown in Fig. 1. Stars in the northern part of the image in Fig. 1 form the first group, which has vectors orientated in the direction of the Galactic disk, see for example stars #1,6,16,184 and 189. To the South three more stars #34,102 and 266 also have vectors orientated with the same direction. Stars #16,102, 184 and 189 are considered as nonmembers of the cluster (Baume, Vázquez & Feinstein, 1999) and this fact can explain the alignment of the polarization vectors with the Galactic plane. We believe that all this group probably consists of nonmembers or less reddened members which are closer to the Sun (“front-side stars”) and we are seing them without the effect of the intra-cluster dust.

It is very noticeable that most of the polarization vectors of the stars in the inner core of the cluster seem to form part of a semi-circular structure and none of them appears aligned with the direction of the Galactic disk. We consider these stars as the second group.

The third group consist of stars with their polarimetric vectors at right angles to the direction of the Galactic disk. These are stars #70,73,80,105,110,112,248,253,254 to the South.

The second group is a very unusual finding. In the next section we will discuss some probable explanations based on our polarimetric observations and on data from other authors obtained with different astronomical techniques (photometry, spectroscopy, etc).

## 4 Analysis and discussion

### 4.1 Serkowski law

To analyze the data, the polarimetric observations were fitted using the Serkowski’s law of interstellar polarization (Serkowski et al., 1975). This is:

$`P_\lambda /P_{\lambda max}=e^{Kln^2(\lambda _{max}/\lambda )}[1]`$

If the polarization is produced by aligned interstellar dust particles, the observed data (in terms of wavelength, UBVRI) will follow , where each star is characterized by $`P_{\lambda max}`$ and $`\lambda _{max}`$.

Adopting $`K=1.66\lambda _{max}+0.01`$ (Whittet et al., 1992), we have fitted our observations and computed the $`\sigma _1`$ parameter (the unit weight error of the fit) in order to quantify the departure of our data from the “theoretical curve” of Serkowski’s law. A $`\sigma _1`$ value larger than 1.5 is considered as an indication of the presence of a component of intrinsic stellar polarization. Another criterion of intrinsic stellar polarization is to compute the dispersion of the position angle for each star normalized by the average of the position angle errors ($`\overline{ϵ}`$, Marraco et al. 1993).

The $`\lambda _{max}`$ values can also be used to test the origin of the polarization. In fact, since the average value of $`\lambda _{max}`$ for the interstellar medium is 0.545 $`\mu m`$ (Serkowski et al., 1975), objects showing $`\lambda _{max}`$ rather lower than this value are also candidates for having an intrinsic component of polarization (Orsatti et al., 1998). The values we have obtained for $`P_{\lambda _{max}}`$, the $`\sigma _1`$ parameter, $`\lambda _{max}`$, and $`\overline{ϵ}`$ together with the identification of stars are listed in Table 3.

Figure 2 shows the observed $`P_\lambda `$ and $`\theta _\lambda `$ vs $`\lambda `$ of those stars which are candidates for having an intrinsic component of polarization (stars #70,73 and 254). For comparison purposes, the best Serkowski’s law fit has been plotted as a continuous line. Star #220 also shows an intrinsic component of polarization on Table 2, but it is already known as polarimetric variable. First Luna (1982) and then St-Louis et al. (1987) have made an extensive polarimetric study of this double-line binary (WC7 + O5-8) system.

### 4.2 Polarization efficiency

Fig. 3 shows the plot of $`P_{\lambda max}`$ vs. $`E_{BV}`$ (Sung et al., 1998; Raboud, Cramer & Bernasconi, 1997; Baume, Vázquez & Feinstein, 1999) for each of the three groups. This plot is very useful for study the polarization efficiency, defined as the ratio $`P_{\lambda max}/E_{BV}`$ which indicates how much polarization is obtained for a certain amount of extinction.

To understand the behavior of the polarimetry and the extinction a detailed discussion of the distribution of the dust on the line of sight to the cluster is necessary. Neckel & Klare (1980) have studied the extinction values of the ISM and distances in our Galaxy using more that 11.000 stars. Their Figure 6j(211) (which is the plot in the direction of NGC 6231) shows two steeps in the absorption increase: the first one is very near to the Sun and associated with the Lupus cloud (Crawford, 2000) at 150 pcs ($`A_V=0.8`$) and the second step is about 700 to 1000 pc and is also detected in the diffuse interstellar lines (Crawford, 2001). This second step of absorption is responsible for $`A_V0.6`$. Baume, Vázquez & Feinstein (1999) have reconstructed Neckel & Klare’s plot with new CCD data, and found a more complex structure obtaining a similar result for the distribution of dust between the Sun and the cluster. Nevertheless, Figure 7 of Baume, Vázquez & Feinstein (1999) shows that NGC 6231 has also intracluster dust.

Most of the studies show that the cluster has differential reddening. For example, Sung et al. (1998) have made a bidimensional map of the distribution of $`E_{BV}`$ over the NGC 6231 cluster (see Fig. 5 of their paper). They show that the core of the cluster is less reddened and seems to lie in a hole in the reddening material ($`E_{BV}0.45`$). To the South they find that the extinction increases to $`E_{BV}0.65`$. The differential reddening of the cluster is extended discussed in Raboud, Cramer & Bernasconi (1997) and Baume, Vázquez & Feinstein (1999).

Fig. 3 is the plot of $`P_{\lambda max}`$ vs $`E_{BV}`$ for each group of stars. The upper panel are stars of the first group, the middle panel is for second group of stars and the bottom panel are data of the third group of stars. The solid line in the first panel of Fig. 3 is the empirical upper limit relation for the polarization efficiency by the interstellar polarization $`P_{\lambda max}=RA_v=9E_{BV}`$ (for normal dust, $`R=3.2`$) given by Serkowski et al. (1975). The dashed line ($`P_{\lambda max}/E_{BV}=5`$) represents the observed normal efficiency of the polarizing properties of the dust (Serkowski et al., 1975).

Nevertheless, the $`E_{BV}`$ data were collected from different sources and probably some of these values have large errors, some qualitative properties emerge from the figure. It is noticiable that all the groups have a polarizing efficiency lower than the $`P_{\lambda max}/E_{BV}5`$ expected. Because for an extinction of $`E_{BV}=0.43`$ (Baume, Vázquez & Feinstein, 1999) and the average polarization efficiency, the typical value should be around $`P_{\lambda max}2\%`$. This value is much higher than the one observed in the cluster (see Table 3). Therefore, in all the groups we have an important depolarization and this efect is notorious for the stars of the semi-circular structure (group 2) that have a much lower polarization in comparison with the first and third groups (see Fig. 3 middle panel). There are not stars with polarization greater than 0.95% in this group.

There is a question that remains to be answered. Why we can see this semicircular structure behind two layers of dust which orientate the polarization vectors in different directions. It is pretty obvious that just having only one dust layer on the line of sight would compose and change the direction of the polarization vectors, washing out this semicircular structure. But on the other hand, the observed semicircular structure seems quite perfect. As we have discussed in the previous paragraph, in the particular case of NGC 6231 there are two layer of dust, a first one at Lupus (hereafter the first component of extinction) and a second one in front of the cluster locate at 1 kpc, (Raboud, Cramer & Bernasconi, 1997; Baume, Vázquez & Feinstein, 1999; Neckel & Klare, 1980), hereafter the second component. It is very important to notice that the two components are orthogonal and one of them can compensate or diminish the other. The first component is dominated by dust in the Galactic plane, which has an angle of polarization of $`45`$ degrees, the second component seems to be at right angles to the Galactic plane at $``$ 135 degrees. Fig. 3 shows a huge depolarization for the second group of stars. We think that in the northern part the first component dominates the observed polarization, while in the South the second component is the strong one. But in the core, where the semicircular structure lays, both components (first and second) cancel each other, therefore making a polarization “window” where we can observe inside the cluster. This interpretation also explains the depolarization of the stars locate in the semi-circular structure (second group).

We tested this interpretation of the extinction toward NGC 6231, by analyzing the polarization of the nonmembers stars. For example, in the first group, #102 is considered a F0III (Perry, Hill & Christodoulou, 1991) star so it must be located at about 379 pc from the Sun and star # 189 is a B9-8V (Raboud, Cramer & Bernasconi, 1997) being located at 1 kpc. Both stars are behind the Lupus cloud and have polarization angles similar to the direction of the Galactic disk. Therefore, we can deduce that the Lupus cloud is responsible for the polarization of the light in the direction of the Galactic disk. This cloud is reddening the stars behind by $`E_{BV}=0.25`$ in average (Neckel & Klare, 1980). For larger distances Raboud, Cramer & Bernasconi (1997) have found that the interstellar medium shows a relatively good transparency between 400 and 1000 pc.

Star # 105 is also a nonmember star but its polarization (and location in the cluster) is related with the third group of stars. Assuming that this star is a main sequence object reddened by normal dust it is possible to found an unique solution from the color-color diagram that gives an $`E_{BV}=0.52`$. This value is larger than the estimate for the Lupus cloud and shows that this object is located behind another layer of dust which is polarizing the stellar light at right angles to the direction of the Galactic disk. Therefore, this star is probably located behind the second cloud at a distance of more 1 kpc from the Sun.

### 4.3  The semicircular polarization pattern

As despicted in Fig. 1, the polarization vectors of the stars belonging to our second group (at the core of the cluster) displayed a semicircular pattern. This is very unlikely to be the result of a distribution of polarization angles at random, because the pattern is well organized. We think we are in the presence of an energetic phenomenon that has changed the magnetic field orientation, resulting in a realignment of the dust grains, which manifests itself in a rotation of the polarization vectors. We propose that a supernova event has occurred some time ago and produced a shock over the ISM. It has pushed the gas and dust (and the associated magnetic field), causing the dust grains to re-oriented to the new direction of the magnetic field. We will discuss in the next paragraphs the relevant questions associated with the proposed supernova event.

Where did this phenomenon occur ?. Tracing lines perpendicular to the most conspicuous polarimetric vectors it is possible to obtain an approximate center at the crossing point where these lines overlap each other. This center seems to be at $`\alpha =16^h54^m22^s`$, $`\delta =41^{}46^{^{}}30^\mathrm{"}`$, J2000.0 (Fig. 4). This semi-circle is very regular and the vectors converge clearly to a center. We think this is a solid argument about the existence of center of this pattern, which means that the energetic phenomenon was produced at a certain fixed location in the cluster. This supports the idea that a supernova exploded in the near past of the cluster, located at the center of the crossing lines as shown in Fig. 4.

Could a supernova event have happened in NGC 6231 ? The earliest star on the main sequence, not considering binaries, is at present of spectral type O9 (#292, HD 326329), and there is evidence of more massive stars that had evolve out of the main sequence (Baume, Vázquez & Feinstein, 1999). Therefore, it is very probable that a star more massive than HD 326329 had evolved to a supernova stage in the past. Also, Baume, Vázquez & Feinstein (1999) have computed the initial mass function of NGC 6231 (hereafter IMF) with binary correction and they have obtained a slope of $`x=1.14\pm 0.10`$, for $`x=log(dN/\delta M)/log(M)`$. This slope of the distribution allows the existence of a former star of about 70 $`M_o`$ (an O3-O4) without any change in the parameters of the observed IMF. This means that one or few (because we are dealing with small number statistics) very massive stars could have lived in the past with no contradiction with the actual IMF. So, the evolutionary status of the cluster population in NGC 6231 indicates that some massive star have evolved and probably became supernovae.

Is there other evidence on the ISM of this supernova event ? Crawford and collaborators (e.g. Crawford et al. (1989), Crawford (1989, 1992, 2001)) have made a very detailed study of the interstellar lines (NaI, CaII, KI) towards Sco OB1. Although these studies are related to the whole Sco OB1 association, these series of papers include in the analysis some of the brightest stars of NGC 6231. They have found several components and classified them into three major groups (Crawford, 2001). The first one consists of several components with radial velocities in the range $`20v_{helio}6kms^1`$ which are consistent with the galactic rotation along the 1900 pc to Sco OB1. All the clouds in this component are presumably foreground, and unrelated with Sco OB1. For example components in the velocity range $`6v_{helio}6kms^1`$ are probably associated with the Lupus molecular cloud complex at a distance of 150 pc.

The second group of interstellar lines found by Crawford (2001, and references therein) is the “shell component” which has velocities more negative than -20 $`kms^1`$ (-22 $`kms^1`$ is the velocity expected for the galactic rotation model at a distance of 1900 pc) and these clouds of interstellar matter seem to be related with the shells surrounding the Sco OB1 association. In a previous work Crawford et al. (1989) found a very low NaI/CaII relation in this component. This finding was explained by the removal of the Ca atoms from the grain surfaces by shock velocities of $`20kms^1`$. New data (Crawford, 2001) indicate that the observed shell component has an upper limit for the temperature of the order of $`T=450^{}K`$ (as observed in the spectrum of HD 152235, a member of NGC 6231), which is low for a post-shock temperature ($`T10^4`$). But as Crawford (2001) noted, the cooling time of a shock for the regular interestellar medium ($`v_s=20kms^1,n_H=10cm^3`$) gives $`t_{cool}=6300yr`$, meanwhile the depletion time scale for the Ca to become part of the interestellar grains is $`t_{dep}=1.6\times 10^7yr`$. Therefore, shocks in the near past can account for the low NaI/CaII observed in the interstellar lines.

The third component of interestellar line has velocities more positive than those predicted by the galactic rotation models. In NGC 6231, HD 152249 exhibits an interstellar absorption component at a velocity of +19 $`kms^1`$. Note that HD 152249 is located in the core of the NGC 6231, and 1’ to the South of the semicircular polarization pattern. This component has a low NaI/CaII ratios ($`0.4`$), which (as described above) is an indication of shocks. This component could be located near the sun in the “Upper-Centarus Lupus” or it could be an expanding interestellar structure related to NGC 6231.

Both second and third interstellar components exhibit the low NaI/CaII ratio and in our opinion it is a strong support for shocks in the region and gives clues over a turbulent scenario in the near past of NGC 6231. These interstellar line components are probably located in the cluster, as discussed by Crawford in his series of papers, and were produced by shocks in the region. This point of view of Crawford is consistent with the interpretation of our polarimetric observations.

Sung et al. (1998), studying the Pre Main Sequence (PMS) stars in NGC 6231, have found a discontinuity in the number of stars near the MS band with ages between 30 Myr and 12 Myr. They have suggested that probably some PMS are not $`H\alpha `$ emitters and they are not detected by the typical $`RH\alpha `$ vs. V-I diagram. They explained the lack of $`H\alpha `$ emission suggesting that all the material surrounding the PMS stars, material that would normally produce strong H$`\alpha `$ emission, could be swept away by stellar winds from massive stars, Wolf-Rayet stars or past supernova explosions. Their result is also compatible and a strong support to our polarimetric findings.

Other interesting result showing some clues to support the idea of energetic events on the cluster NGC 6231 came from Manchanda et al. (1996). To explain the COS B diffuse emission of gamma-rays observations, Manchanda et al. (1996) propose that strong winds in the cluster interact with the interstellar medium creating the suitable condition for cosmic ray enhancement.

## 5 Other possible explanations of the semicircular polarization pattern

We have proposed a supernova event to explain our observations, but we have analyzed other possible interpretations. For example, a loop in the magnetic field would lead to similar observations, however several papers (Neckel & Klare, 1980; Baume, Vázquez & Feinstein, 1999) show very clearly that the dust along the line of sight is distributed in layers and there is no dust between these layers. Thus, there is not a continuum of dust to be aligned by the magnetic field. Therefore, we have discarded this interpretation of the semicircular polarization pattern.

Other explanation for the observed polarization pattern might not need a supernova. Same observed bubbles in the Milky Way are supposed to be made by the strong winds of Wolf-Rayet stars over a long period of time. In some of these interstellar bubbles, like WR 101 (Cappa et al., 2002), the Wolf-Rayet star appears to be locate in an eccentric position close to the densest part of the nebula. This phenomenon could be explained by an inhomogeneous structure of the ISM near the Wolf-Rayet star or from the motion of the star inside the shell (Cappa et al., 2002). This would be the case if the Wolf-Rayet star #220 (HD 152270) is responsible for the bubble in NGC 6231. The location of this star is eccentric to the semicircular pattern.

## 6 Summary

We have observed multicolor (UBVRI) linear polarization ($`P`$ and $`\theta `$) for 35 of stars on the direction of NGC 6231. We have found a very unusual distribution of the angles of the polarization vectors. To the North of the cluster the polarization angles are dominated by the direction of the Galactic disk; to the South are perpendicular to this orientation. Close to the core of the cluster, the polarization vectors are well organized in a regular pattern, showing a semicircular morphology which we interpreted as a re-orientation of the dust particles by the magnetic field.

We suppose that an energetic event (very probably a supernova) has occurred some time ago and produced a shock pattern on the local ISM. A supernova is compatible with the evolutionary status of the cluster. We discuss in this paper independent confirmations of this event from studies of the interstellar lines (Crawford, 2001) to the photometric results of the PMS stars given by Sung et al. (1998).

We wish to acknowledge the technical support of CASLEO during our observing runs. Also we want to acknowledge the useful discussions with Ana M. Orsatti, Virpi Niemela, Cristina Cappa, Paula Benaglia which are greatly appreciated. We also aknowledge the use of the Torino Photopolarimeter built at Osservatorio Astronomico di Torino (Italy) and operated under agreement between Complejo Astronomico El Leoncito and Osservatorio Astronomico di Torino. We also want to thank to Dr. Michael Corcoran for provide us the Rosat images of the region.


format: text

text:
# Limits on Very High Energy Emission from Gamma-Ray Bursts with the Milagro Observatory

## 1. Introduction

Very high energy gamma-ray burst (GRB) observations have the potential to constrain theoretical models of both the prompt and extended phases of GRB emission. Models based on both internal and external shocks have predicted VHE fluence comparable to, or in certain situations stronger than, the keV – MeV radiation, with durations ranging from shorter than the keV – MeV burst to extended TeV afterglows (Dermer & Chiang, 2000; Pilla & Loeb, 1998; Zhang & Mézáros, 2001). GRB emission above 100 GeV is particularly sensitive to the Lorentz factor and the photon density of the emitting material — and thus the distance of the radiating shock from the source — due to $`\gamma \gamma e^+e^{}`$ absorption in the emission region.

EGRET observed several GRBs at GeV energies and observed no evidence of a high energy rollover in the GRB spectrum (Dingus, 1995), and recent results by Gonzalez (2003) indicate that the spectrum of some GRBs contains a very hard, luminous, long duration component. At energies above a few 100 GeV several observation attempts have been made, including: coincidence observations between the Tibet air-shower array and BATSE (Amenomori, 1996), rapid follow-up observations by the Whipple Air-Cherenkov telescope (Connaughton et al., 1997), and coincident and monitoring studies by AIROBICC (Padilla et al., 1998), Whipple (Connaughton et al., 1998), and the Milagro prototype (Atkins et al., 2000a; Leonor, 2000; McCullough, 2001). While several of these studies have reported evidence for VHE GRB emission, no clear picture of the VHE emission from gamma-ray bursts has emerged.

In recent years the rate of possible GRB coincident detections has been relatively low due to the absence of sensitive, large field of view GRB detectors like BATSE. However, Milagro is a very sensitive detector with a low $``$100 GeV energy threshold and wide $``$ 2 sr field of view, and is fully capable of autonomous identification of a VHE burst of emission. The low energy threshold is particularly important because of the reduced attenuation by extragalactic background light near 100 GeV (Jelley, 1966; Primack et al., 1999; Stecker & de Jager, 1998), which dramatically increases the volume of space observed. This paper details the monitoring of the northern sky for VHE transients of 40 s – 3 hours duration with the Milagro observatory.

## 2. The 40 s to 3 hour Transient Search in Milagro

Very high energy gamma rays incident on the earth interact in the upper atmosphere and produce extensive air showers (EAS) which propagate to lower altitudes. Milagro uses the water Cherenkov technique to detect EASs by converting the wave front of relativistic particles in the EAS into a front of Cherenkov light, detected by photomultiplier tubes (PMTs) (Atkins et al., 2000b). The amplitude and arrival time of the PMT signals are analyzed to determine the direction and characteristics of the $``$1,800 EAS per second detected by Milagro, most of which are initiated by cosmic ray nuclei. The water Cherenkov technique enables Milagro to operate continuously.

The precision of the Milagro event reconstruction depends strongly on the shower characteristics, leading to a point spread function (PSF) for the initiating particle direction which varies considerably from event to event. The arrival times measured by the PMTs are fit to a plane, after making corrections based on the pulse amplitude in each PMT and the distance of the PMT from the reconstructed shower core. The shower fit proceeds through several iterations, at each stage removing times with large residuals to the fit shower plane. Showers reconstructed within $`45^{}`$ of zenith are kept for further analysis. The number of PMT times surviving in the final iteration is called $`n_{\mathrm{fit}}`$, which is used in conjunction with the reduced $`\chi ^2`$ to characterize the quality of the reconstruction. The showers are then separated into 13 classes, based on the values of $`n_{\mathrm{fit}}`$ and $`\chi ^2`$, such that events within each class have similar angular PSFs (see Morales (2002) for details). In general, the width of the PSF decreases as $`n_{\mathrm{fit}}`$ increases or $`\chi ^2`$ decreases. The class with the worst PSF is removed from further analysis, at a negligible cost in sensitivity, because the PSF was so broad that it could not be readily handled by the code used to accumulate the sky maps.

The characteristics of the light recorded by PMTs at the bottom of the reservoir can be used to reject background cosmic ray showers in favor of gamma rays using a quantity we call “compactness” (Atkins et al., 2003). Events are assigned a weight, between 0 and 1, based on the relative probability of the event being a background cosmic ray or a gamma ray, given the compactness value and to which of the 12 remaining PSF classes the event belongs. To speed the computation, events with weight less than 0.5 are rejected.

This search implements an analysis method developed by Morales et al. (2004) which efficiently includes the event-by-event PSF and photon probability (i.e. background rejection) information, enhancing the sensitivity of the observations. The technique uses the characteristics described above to determine a spatial photon probability density — the expected PSF of the event multiplied by the probability the shower was initiated by a gamma ray — which is then added to a sky map to form the estimated total photon density. Therefore, rather than considering each photon as a point at each reconstructed event location, the new technique instead places a smooth distribution at that location, with the shape of the distribution mirroring the reconstruction accuracy of the detector. The probability of the background, which primarily consists of cosmic ray events, producing an observed peak in the sky map is then given by the expected statistical distribution of the total photon density at that location on the sky map. The expected statistical distributions, which are typically not Gaussian in shape, can be calculated directly from the PSF characteristics of the detector and placed into a lookup table to enable very fast searches for unknown gamma-ray sources. This technique is similar to maximum likelihood analyses in sensitivity, but is computationally fast enough to perform a real time transient search in a high data-rate detector like Milagro .

The 40 s – 3 h transient search is performed by first dividing the data into 20 second intervals and creating sky maps, both for the actual data (“on source”) and for the expected background in each time interval. From these initial maps, 40 s sky maps are formed by summing all pairs of successive 20 s maps. The resulting set of sky maps is temporally oversampled by a factor of two, with every 20 s map used twice so that a 40 s burst is guaranteed to be largely contained in one of the overlapping 40 s maps. A search is performed on each of these maps for emission using the test statistic described below, then pairs of the 40 s maps are combined to make 80 s maps ($`1/2`$ of the 40 s maps are discarded to maintain the factor of two oversampling), and so on. The nine specific time scales that were searched are 40, 80, 160, 320, 640, 1280, 2560, 5120, and 10240 seconds. The total observation days for each time scale is shown in the second column of Table 1. The longer time scales have slightly less observation time due to gaps from occasional data dropouts. Full details of the transient search can be found in Morales (2002).

The significance of excesses in the on-source sky maps, compared to the background sky maps, is evaluated using a test statistic. The test statistic is very closely related to the probability of an observation being produced by a background fluctuation, with approximations used to speed up the analysis breaking the exact one-to-one relationship and leading to a slight flattening of the test statistic distributions (Morales et al., 2004). Since the value of the test statistic decreases with increasing transient flux (lower probability of the observation being produced by a background fluctuation), a discovery search simply looks for values of the test statistic below that expected from the number of independent observations. If a transient is identified, the statistical significance can be calculated from the test statistic.

The data used for this analysis were taken between 2001 May $`2^{\mathrm{nd}}`$ and 2002 May $`22^{\mathrm{nd}}`$. Example distributions of the test statistic for 4 time scales are shown in Figure 1, with no observed location for any of the nine search intervals having a test statistic below the $``$$`10^{12}`$ value expected from the number of independent observations. There is no evidence for transient VHE emission of 40 s to 3 hours duration in the Milagro data between 2001 May $`2^{\mathrm{nd}}`$ and 2002 May $`22^{\mathrm{nd}}`$.

## 3. Limits on VHE GRB Emission

Two sets of constraints are presented here — upper limits on the VHE flux at the earth, and sensitivity distributions which constrain the absolute VHE luminosity of gamma-ray bursts. The flux limits have the advantage of being independent of cosmological models and directly comparable to previous VHE GRB limits. However, due to the energy-dependent absorption of VHE photons by extragalactic background light (EBL), it is difficult to compare flux limits at the earth to theoretical calculations. The sensitivity distributions place the observations into a theoretical context by correcting for absorption and cosmological effects. In conjunction with model predictions these can be used to determine the maximum VHE luminosity at the source which is consistent with the observations presented in this paper.

Since no observation in any of the nine time scales had a test statistic below $`10^{12}`$, we can use this as a threshold to set a conservative limit on the maximum flux consistent with these observations. The relationship between the gamma ray flux and the test statistic was determined using a combination of Monte Carlo generated transient VHE signals and real data background. The simulated detector signal was modeled from 100 GeV to 21 TeV as a function of zenith angle and observed spectrum, propagated through the Milagro reconstruction code, added to representative data, and fed to the 40 s to 3 hour transient detection program. If the test statistic for the simulated signal fell below the $`10^{12}`$ threshold of this study, the signal would have been stronger than any observed fluctuation on any time scale.

The flux limits represent the VHE transient signals which would have produced a signal with a test statistic below $`10^{12}`$ $`90\%`$ of the time (signal from a single burst which is excluded at the $`90\%`$ confidence level). For this paper $`E^{2.0}`$ power-law spectra extending from 100 GeV to hard spectral cutoffs at 300 GeV, 1 TeV, or 21 TeV (as observed local to the earth) were chosen. These spectra serve as approximate models for an inverse Compton spectral bump at TeV energies or a hard VHE extension of the observed GRB spectrum past the multi-GeV observations of EGRET (Dingus, 1995), with the spectral cutoff approximating either an intrinsic cutoff for nearby sources or EBL absorption for more distant sources. The limits are presented in Figure 2 as flux density as a function of zenith angle and burst duration.111The full set of limits for $`E^{2.0}`$ and $`E^{2.4}`$ spectra without cutoffs are shown in Morales (2002). Given the $`0.8`$ years of observation time and the $`2`$ sr field of view, these flux limits provide the strongest constraints to date on the the number of nearby bursts which emit TeV photons. For example, at a redshift of 0.024 (the 100 Mpc distance of the nearby bursts identified by Norris (2002), with the cosmology used in that paper) our 90% confidence limit for a spectrum which extends to 21 TeV and burst emission which lasts 40 s at a 10.5 degree zenith angle is $`1.1\times 10^{46}`$ ergs/s.

However, at larger distances the usefulness of the limits from Figure 2 is complicated by the absorption of VHE photons by extragalactic background light (Jelley, 1966; Primack et al., 1999; Stecker & de Jager, 1998). The observed spectrum is a convolution of the emitted spectrum and the extragalactic background light absorption — neither of which are well determined. Limits on the absolute VHE luminosity of GRBs are thus necessarily model dependent. In setting the limits which follow, Milagro is more sensitive than might initially be expected because bursts above $``$1051 erg/s are sufficiently bright to be detected even after extinction of a few optical depths.

In an effort to place the current observations in context, a set of assumptions about the emitted spectrum, EBL absorption, and cosmology have been chosen and sensitivity distributions calculated within this theoretical framework. For a set of source luminosities and durations, simulated GRBs were created with an $`E^{2.0}`$ emission spectrum from 100 GeV to 21 TeV, an isotropic sky position, and following the star formation rate in z. The EBL absorption was determined by Bullock et al. (2004) (similar to Primack et al. (1999)) using recent results from semi-analytic modeling, and includes collisional starburst effects and a Kennicut initial mass function (Kennicut, 1983).222The magnitude of the EBL absorption is quite uncertain, and the models by Stecker & de Jager (1998) have considerably less attenuation. Consequently the limits presented here may be conservative. A $`\mathrm{\Lambda }`$CDM cosmology ($`\mathrm{\Omega }_\mathrm{M}=.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=.7`$, $`\mathrm{h}=0.65`$) was used for both the EBL calculation and the cosmological corrections. The results are rather insensitive to the emitted spectrum and cutoff energy due to the strong EBL absorption which eliminates nearly all emission above 300 GeV at redshift $`0.3`$, and limits the distance to which Milagro can observe VHE emission to redshift $`0.7`$. The star formation rate (SFR) was modeled as $`10^{1.45z1.52}`$ M/Mpc3/year for $`z<0.9`$ and 0.61 M/Mpc3/year for $`z0.9`$ within each redshift bin, following the discussion in Hopkins et al. (2001).

The results of the sensitivity calculation are presented in Table 1 as the probability that Milagro would have observed a VHE GRB (i.e. test statistic $`<10^{12}`$) as a function of the absolute luminosity,333The luminosity quoted is the isotropic luminosity assuming no beaming at the source, and is used to facilitate comparison with MeV to GeV measurements. distance, and duration of the source. Due to increased atmospheric depth away from detector zenith, the low energy threshold of Milagro increases with zenith angle, as can be seen in Figure 2 by the decreased sensitivity at large zenith angles when the emission is limited to lower energies. As GRBs become increasingly distant the EBL absorption effectively truncates the observed spectrum above a few hundred GeV. These effects lead to the effective field of view of Milagro becoming narrower with increased redshift, producing the gradual decline in detection probability seen in Table 1.

Given theoretical expectations for the GRB distance, luminosity, and duration distributions, the probabilities in Table 1 can be used to calculate the upper limits provided by these observations. For a 90% confidence upper limit, the maximum GRB rate consistent with these observations is given by multiplying the expected number of GRBs within each redshift bin by the detection probability in Table 1, and adjusting the total number of GRBs until the mean number of events seen by Milagro is 2.3. If the model GRB distance distribution does not follow the SFR, the probabilities in Table 1 can still be used as approximate values, or corrections can be applied to adjust for the distribution of events within each bin. The results are typically quoted as the number of GRBs/$`\mathrm{M}_{}`$ of star formation or the average number of GRBs/$`\mathrm{Gpc}^3/\mathrm{year}`$.

As an example, consider a model which predicts that GRBs follow the star formation rate, with all GRBs emitting a characteristic 80 s pulse of VHE emission. The resulting 90% confidence upper limits for this model are $`6.2\times 10^8`$ GRBs/$`\mathrm{M}_{}`$ of star formation (an average of $`4.8`$ GRBs/$`\mathrm{Gpc}^3/\mathrm{year}`$ over 0$`<`$z$`<`$0.5) for an isotropic luminosity of $`10^{51}`$ ergs/s, or $`1.1\times 10^8`$ GRBs/$`\mathrm{M}_{}`$ of star formation (an average of 0.8 GRBs/$`\mathrm{Gpc}^3/\mathrm{year}`$ over 0$`<`$z$`<`$0.5) for a luminosity of $`10^{52}`$ ergs/s, with significantly tighter constraints if the GRB distribution trails the star formation rate (i.e. there are more low redshift GRBs). Of the thirty-six GRBs with known distances, five have a redshift below 0.5.444One can use this fraction to roughly estimate that if all bursts were 1051 erg/s, they would yield a flux brighter than the sensitivity in Figure 2 ($`10^7`$ cm-2 s-1) out to redshift of about 0.5 for zenith angles less than $`30^{}`$, and approximately $`(5/36)\times 0.84/(4\pi )=1\%`$ of all bursts would be detectible by Milagro. If the GRBs detected by BATSE follow the same distance distribution, a rough estimate of the observed GRB rate yields $``$$`2.6`$ GRBs/$`\mathrm{Gpc}^3/\mathrm{year}`$, or $``$$`3.4\times 10^8`$ GRBs/$`\mathrm{M}_{}`$ of star formation if they follow the SFR. While detailed model calculations are needed to convert the probabilities in Table 1 into meaningful upper limits, comparison with the limits from this simple model suggests that if GRBs follow the star formation rate, the typical luminosity of 40 s – 3 h VHE GRB counterparts is constrained to be similar to or less than the prompt keV – MeV emission.

## Acknowledgments

We acknowledge Scott Delay and Michael Schneider for their dedicated efforts in the construction and maintenance of the Milagro experiment. This work has been supported by the National Science Foundation (under grants PHY-0070927, -0070933, -0075326, -0096256, -0097315, -0206656, -0302000, and ATM-0002744) the US Department of Energy (Office of High-Energy Physics and Office of Nuclear Physics), Los Alamos National Laboratory, the University of California, and the Institute of Geophysics and Planetary Physics. MFM was a NASA Graduate Student Researcher.


text:
# Spontaneous squeezing of a vortex in an optical lattice

## I Introduction

Vortices play a crucial role in explaining rotational and dissipative properties of superfluids. In the superfluid systems studied until now, such as liquid Helium Donnelly (1991), superconductors Blatter et al. (1994), and Bose-Einstein condensates Matthews et al. (1999); Madison et al. (2000); Abo-Shaeer et al. (2001); Hodby et al. (2002); Butts and Rokshar (1999); Williams and Holland (1999); Dodd et al. (1997); Zambelli and Stringari (1998); Fetter and Svidzinsky (2001), a vortex line behaves as a classical object. However, placing a Bose-Einstein condensate in an optical lattice Cataliotti et al. (2001); Burger et al. (2001); Orzel et al. (2001); Greiner et al. (2002a, b) leads to an unprecedented control over the system parameters, which has already enabled experimental studies on quantum phase-transitions Greiner et al. (2002a); Jaksch et al. (1998); van Oosten et al. (2001), superfluidity Cataliotti et al. (2001); Burger et al. (2001), and number squeezing Orzel et al. (2001). Furthermore, a Bose-Einstein condensate in a one-dimensional optical lattice has a similar layered structure as the cuprate superconductors and is a promising way to achieve the quantum Hall regime in a Bose-Einstein condensed gas Cor . This suggests that a careful examination of the quantum properties of a vortex line in a one-dimensional optical lattice is warranted.

A characteristic feature of all superfluids is their ability to support persistent currents. Quantized vortices play a crucial role in understanding the decay of these currents, and also the superfluid’s response to rotation or to an external magnetic field in the case of superconductor. In pioneering work Fetter developed the quantum theory of vortices in liquid Helium Fetter (1967a) and so-called type-$`\mathrm{II}`$ superconductors Fetter (1967b). In these systems quantum fluctuations turned out to play such a small role that they are experimentally inaccessible. Here we show, however, that truly quantum mechanical behavior of the macroscopic vortex occurs in a Bose-Einstein condensate in a one-dimensional optical lattice. This is due to the reduced dimensionality as well as the reduced number of atoms in every well of the optical lattice. Remarkably, it turns out that the vortex can spontaneously become strongly squeezed, which is reflected in the quantum-mechanical probability distribution of the vortex position. Unlike coherent states, which are described by Gross-Pitaevskii theory, squeezed states are highly nonclassical.

## II Bose-Hubbard Hamiltonian

We consider a Bose-Einstein condensate in a one-dimensional optical lattice with lattice spacing $`\lambda /2`$, where $`\lambda `$ is the wavelength of the laser beams creating the standing wave of the lattice. The optical lattice splits the Bose-Einstein condensate into a stack of $`N_s`$ weakly-coupled pancake condensates, each containing $`N`$ atoms. For concreteness we use in the following always $`{}_{}{}^{87}\mathrm{Rb}`$ atoms. The depth of the lattice $`V_L`$ can be easily changed and controls the tunneling of an atom from one pancake condensate to the next, i.e., the strength of the interlayer Josephson coupling. The condensate also experiences a harmonic trapping potential in the radial direction with an oscillator frequency $`\omega _r`$. The longitudinal trapping potential along the direction of the laser beams is assumed to be so weak that it can be neglected. While the lattice is taken to be deep enough to allow us to use a tight-binding approximation and to include only the weak nearest-neighbor Josephson coupling, it is also taken to be shallow enough to support a superfluid state as opposed to the Mott-insulator state van Oosten et al. (2001).

Furthermore, we consider a vortex line that pierces at the position $`(x_n,y_n)`$ through each layer of the stack labelled by the index $`n`$. In each layer the density of the pancake condensate varies as a function of the radial distance. As a result the vortex experiences an effective potential that depends on its distance from the origin Fetter and Svidzinsky (2001). Without rotation of the gas this potential is approximately an inverted parabola. Hence, the vortex is energetically unstable and will tend to spiral out of the system. However, the vortex is stabilized if we rotate the Bose-Einstein condensate with a rotation frequency $`\mathrm{\Omega }`$ that is larger than a critical frequency $`\mathrm{\Omega }_c`$. In the rest of this paper we consider the regime where the vortex is energetically stable and $`\mathrm{\Omega }>\mathrm{\Omega }_c`$.

In the tight-binding approximation the attraction between the nearest-neighbor parts of the vortex turns out to be harmonic with respect to their separation. We denote the typical strength of this attraction by $`J_V`$. Physically the attractive interaction is due to the energy cost for phase differences between the two layers. The theory of the resulting coupled harmonic oscillators can be quantized by introducing the bosonic annihilation operators for the eigenmodes of the vortex line Bretin et al. (2003); Mizushima et al. (2003); Martikainen and Stoof (2003). These modes are the Kelvin modes (or kelvons) and correspond physically to a wiggling of the vortex line with a wavelength $`2\pi /k`$. The operators $`\widehat{v}_k`$ in momentum space are related to the coordinate space kelvon operators $`\widehat{v}_n`$ and to the vortex positions through

$$\widehat{x}_n=\frac{R}{2\sqrt{N}}(\widehat{v}_n^{}+\widehat{v}_n)=\frac{R}{2\sqrt{NN_s}}\underset{k}{}e^{ikn\lambda /2}\left(\widehat{v}_k^{}+\widehat{v}_k\right)$$
(1)
and

$$\widehat{y}_n=\frac{iR}{2\sqrt{N}}(\widehat{v}_n^{}\widehat{v}_n),$$
(2)
where $`R`$ is the typical radial size of each pancake condensate. In this way we obtain the dispersion relation $`\mathrm{}\omega _K(k)=\mathrm{}\omega _K(0)+J_V\left[1\mathrm{cos}\left(k\lambda /2\right)\right]`$ of the kelvons Martikainen and Stoof (2003). The energetic stability of the vortex is reflected in a positive value of the dispersion at zero momentum. Physically the frequency at zero momentum is the precession frequency of a slightly displaced straight vortex line around the center of the Bose-Einstein condensate.

It is crucial for our purposes to realize that, due to the Euler dynamics of the vortex, the coordinates $`\widehat{x}_n`$ and $`\widehat{y}_n`$ are canonically conjugate variables and therefore obey the Heisenberg uncertainty relation $`[\widehat{x}_n,\widehat{y}_n]=iR^2/2N`$. This means that the position of the vortex is always smeared out by quantum fluctuations. For a vortex in a lattice these fluctuations can become large due to the reduced particle numbers in every layer, as opposed to the total number of particles, and the radial spreading of the pancake condensate wave function as the lattice depth is increased.

The kelvon dispersion relation was determined by taking into account only the lowest order expansion of the effective potential experienced by the vortex. However, expanding up to fourth order in the vortex displacements results in an interaction $`\left(V_0/2\right)_n\widehat{v}_n^{}\widehat{v}_n^{}\widehat{v}_n\widehat{v}_n`$ for the kelvons. In the regime where the vortex is stable $`V_0<0`$ and this interaction is attractive. The physics of the vortex line is thus described by a one-dimensional Bose-Hubbard model with a negative interaction strength. The corresponding Hamiltonian is given by

$$\widehat{H}=\underset{n}{}\left[\mathrm{}\omega _K(0)+J_V\right]\widehat{v}_n^{}\widehat{v}_n\frac{J_V}{2}\underset{n,m}{}\widehat{v}_m^{}\widehat{v}_n+\frac{V_0}{2}\underset{n}{}\widehat{v}_n^{}\widehat{v}_n^{}\widehat{v}_n\widehat{v}_n,$$
(3)
where $`n,m`$ denotes the nearest-neighbor layers.

We calculate the parameters of this Bose-Hubbard model using a variational ansatz we used in our earlier work Martikainen and Stoof (2003), where the density of the condensate wave function in each layer was proportional to $`\mathrm{exp}\left[\left(x^2+y^2\right)/R^2\right]`$. In this way we find that the strength of the nearest-neighbor coupling is given by

$$J_V=\frac{\mathrm{}\omega _r}{4\pi ^2}\mathrm{\Gamma }[0,\frac{l_r^4}{R^4}]\left(\frac{\omega _L\lambda }{\omega _rl_r}\right)^2\left(\frac{\pi ^2}{4}1\right)\mathrm{exp}\left(\frac{\lambda ^2m\omega _L}{4\mathrm{}}\right)$$

and the interaction strength is given by

$$V_0=\frac{2\mathrm{}\omega _r\left(l_r/R\right)^2\mathrm{\Gamma }[0,l_r^4/R^4]3\mathrm{}\omega _r\left(l_r/R\right)^24\mathrm{}\mathrm{\Omega }}{4N},$$

where $`\omega _L=\sqrt{8\pi ^2V_L/m\lambda ^2}`$ is the oscillator frequency of the optical lattice, $`l_r=\sqrt{\mathrm{}/m\omega _r}`$, and $`\mathrm{\Gamma }[a,z]`$ is the incomplete gamma function com . Furthermore, a straight and slightly displaced vortex precesses around the condensate center with the frequency

$$\omega _K(k=0)=\left(\omega _rl_r^2/2R^2\right)\left(1\mathrm{\Gamma }[0,l_r^4/R^4]\right)+\mathrm{\Omega }.$$
(4)
The precession frequency of the vortex changes from negative to positive at the critical rotation frequency

$$\mathrm{\Omega }_c=\left(\omega _rl_r^2/2R^2\right)\left(\mathrm{\Gamma }[0,l_r^4/R^4]1\right).$$
(5)
When $`\mathrm{\Omega }>\mathrm{\Omega }_c`$ the vortex is locally energetically stable. The condensate size $`R`$ is given by the solution of the transcendental equation $`\left(l_r/R\right)^4\left[1\gamma +2Na\sqrt{m\omega _L/2\pi \mathrm{}}4\mathrm{ln}\left(l_r/R\right)\right]=1`$, where $`\gamma `$ is the Euler-Mascheroni constant and $`a`$ is the three-dimensional scattering length. Without the kelvon interaction the ground state of the vortex corresponds to a vortex line along the symmetry axis of the Bose-Einstein condensate and with a quantum mechanical uncertainty such that $`\widehat{x}_n^2=\widehat{y}_n^2`$. An attractive interaction will lead to squeezing with $`\widehat{x}_n^2\widehat{y}_n^2`$ as we show next.

## III Vortex squeezing

First we have to identify the correct order parameter for the vortex line. Since $`\mathrm{\Omega }>\mathrm{\Omega }_c`$ there occurs no Bose-Einstein condensation of the kelvons and the appropriate order parameter for the vortex line is not $`\widehat{v}_n`$, which would signal the above mentioned energetic instability and the tendency of the vortex to move away from the symmetry axis. However, the correct order parameter can be identified with $`\mathrm{\Delta }=V_0\widehat{v}_n\widehat{v}_n`$. To arrive at the associated mean-field theory, we quadratically expand the Hamiltonian around this order parameter and obtain apart from a constant Stoof (1994)

$$\widehat{H}=\underset{n}{}\left[\mathrm{}\omega _K(0)+J_V\right]\widehat{v}_n^{}\widehat{v}_n\frac{J_V}{2}\underset{n,m}{}\widehat{v}_m^{}\widehat{v}_n+\underset{n}{}\left[\frac{\mathrm{\Delta }}{2}\widehat{v}_n^{}\widehat{v}_n^{}+\frac{\mathrm{\Delta }^{}}{2}\widehat{v}_n\widehat{v}_n\right].$$
(6)
We diagonalize this result by means of a Bogoliubov transformation. It is most convenient to work in the momentum space where we define new creation and annihilation operators with the help of of the Bogoliubov amplitudes $`u_k`$ and $`v_k`$ as $`\widehat{b}_k=u_k\widehat{v}_k^{}v_k\widehat{v}_k`$ and $`\widehat{b}_k^{}=u_k^{}\widehat{v}_kv_k^{}\widehat{v}_k^{}`$. Furthermore, to ensure bosonic commutation relations we must have $`|u_k|^2|v_k|^2=1`$. The Hamiltonian in Eq. (6) is then diagonalized with the choise

$$|v_k|^2=\frac{|\mathrm{\Delta }|^2}{\left(E(k)+\omega _K(k)\right)^2|\mathrm{\Delta }|^2}$$
(7)
and

$$|u_k|^2=\frac{\left(E(k)+\omega _K(k)\right)^2}{\left(E(k)+\omega _K(k)\right)^2|\mathrm{\Delta }|^2},$$
(8)
where $`E(k)=\sqrt{\omega _K(k)^2|\mathrm{\Delta }|^2}`$ is the dispersion of the Bogoliubov quasiparticles.

Requiring self-consistency of the approach leads to a gap equation

$$\frac{1}{V_0}=\frac{1}{N_s}\underset{k}{}\frac{1+2N_k}{2E(k)},$$
(9)
where $`N_k=1/(e^{\beta E(k)}1)`$ is the Bose distribution. We solve the gap equation for the order parameter, analogous to the BCS theory of superconductors Schrieffer (1964). In Fig. 1 we present the behavior of the order parameter as a function of temperature and rotation frequency for typical parameter values. We see from Fig. 1 that there is a critical temperature below which the order parameter becomes nonzero and then increases monotonically with decreasing temperature. Furthermore, the absolute value of the order parameter at zero temperature is usually close to the precession frequency. Only very close to the critical rotation frequency $`\mathrm{\Omega }_c`$ the deviation becomes large.

The complex order parameter $`\mathrm{\Delta }=|\mathrm{\Delta }|e^{i\varphi }`$ has an interesting physical interpretation in terms of the quantum mechanical uncertainty of the vortex position. It turns out that in a coordinate system rotated by an angle $`\theta `$ we have

$$\widehat{y}_n^2\widehat{x}_n^2=\frac{|\mathrm{\Delta }|R^2\mathrm{cos}\left(\varphi 2\theta \right)}{|V_0|NN_s},$$
(10)
where $`\varphi `$ is the phase of the order parameter. This implies that a nonzero value of the order parameter is reflected in the squeezing of the vortex position distribution. The main axis of the uncertainty ellipse is at an angle $`\theta =\varphi /2`$ with the x-axis. In equilibrium the uncertainty ellipse of the vortex position distribution is independent of the layer index $`n`$. Therefore, the measurement of the vortex positions in different layers samples the same distribution and provides a signature for the expected transition into a squeezed state. In Fig. 2 we plot the relative squeezing as a function of rotation frequency at fixed temperature. As can be seen from this figure the size of the uncertainty ellipse can easily be several times the radial trap length and also becomes very strongly deformed as the rotation frequency is increased. Furthermore, the long axis of the uncertainty ellipse can be much larger than the vortex core size. The results in Fig. 2 were obtained with only one set of typical parameters. By a careful choice of, for example, lattice depth, number of particles, number of layers, or temperature a considerable degree of tunability is possible. For example, a shallower lattice implies a larger relative squeezing.

## IV Phase fluctuations

For our theory to be valid, constraints have to be set for the temperature. Since we assume a pure condensate, we are ignoring the influence of the noncondensate atoms on the vortex line. Therefore, the temperature should be much lower than the critical temperature for Bose-Einstein condensation in the layers. This condition is relatively easy to satisfy. However, a more stringent condition is set by the phase fluctuations of the order parameter. In an infinite one-dimensional system phase fluctuations destroy the long-range order. In the finite system we are considering here, phase fluctuations can only be excited if the temperature is high enough Petrov et al. (2000); Andersen et al. (2002); Khawaja et al. (2002). In order to calculate the energy cost for a phase gradient of the order parameter we must determine the associated stiffness or superfluid density $`\rho _S(T)`$. Since the lattice breaks the Galilean invariance this calculation is not entirely standard.

The superfluid density is defined by the systems response to an imposed phase gradient Rey et al. (2003). Therefore, we want to study the long-wavelength effects of the transformation $`\widehat{v}_n\widehat{v}_ne^{i\varphi _n/2}`$. The phase gradients enter only in the nearest-neighbor coupling and result in an interaction energy

$$\widehat{H}_I=\frac{J_V}{2}\underset{n,m}{}\widehat{v}_m^{}\widehat{v}_n\left[i\left(\varphi _n\varphi _m\right)/2\left(\varphi _n\varphi _m\right)^2/8\right],$$

which can be considered as a disturbance in the long-wavelength limit. With this interaction energy we calculate in second-order perturbation theory the contribution to the free energy of the phase gradients. Up to second order the effective imaginary time action for the phase fluctuations is

$$S_\varphi =_0^\mathrm{}\beta 𝑑\tau H_I\left(\tau \right)\frac{1}{2\mathrm{}}_0^\mathrm{}\beta _0^\mathrm{}\beta 𝑑\tau 𝑑\tau ^{}H_I\left(\tau \right)H_I\left(\tau ^{}\right).$$

Furthermore, it is sufficient to keep only the terms up to second order in the phase differences between layers. In this way we find the action for the phase gradients as

$$S_\varphi =_0^\mathrm{}\beta 𝑑\tau \frac{J_V}{8N_s}\underset{k^{}}{}\left\{\mathrm{cos}\left(k^{}\lambda /2\right)\left[|u_k^{}|^2N_k^{}+|v_k^{}|^2\left(N_k^{}+1\right)\right]J_V\beta \mathrm{sin}^2\left(k^{}\lambda /2\right)N_k^{}\left(N_k^{}+1\right)\right\}\times \underset{k}{}k^2\varphi _k\varphi _k.$$
(11)
This result should be compared with the standard (Galilean invariant) Landau Schrieffer (1964) result

$$S_L=_0^\mathrm{}\beta 𝑑\tau 𝑑z\frac{\mathrm{}^2}{2m_\mathrm{\Delta }}\rho _s(T)|\varphi |^2,$$
(12)
where the superfluid density at temperature $`T`$ is

$$\rho _s(T)=\frac{1}{N_s\lambda }\underset{k}{}\left[|u_k|^2N_k+|v_k|^2\left(N_k+1\right)\frac{\mathrm{}^2\beta }{m_\mathrm{\Delta }}k^2N_k\left(N_k+1\right)\right].$$
(13)
When we take the lattice spacing to zero while keeping $`J_V\lambda ^2`$ constant, we can see that Eq. (11) recovers the standard result when the mass of the pairing field is defined as $`m_\mathrm{\Delta }=2m_K`$, where

$$m_K=\frac{4\mathrm{}^2}{J_V\lambda ^2}.$$
(14)
is the effective mass of the kelvons. Furthermore, from this we see that the superfluid density in our case is defined by

$$\rho _s(T)=\frac{1}{N_s\lambda }\underset{k}{}\left\{\mathrm{cos}\left(k\lambda /2\right)\left[|u_k|^2N_k+|v_k|^2\left(N_k+1\right)\right]J_V\beta \mathrm{sin}^2\left(k\lambda /2\right)N_k\left(N_k+1\right)\right\}.$$
(15)
Equating the energy cost due to a $`4\pi /N_s\lambda `$ phase gradient with the thermal energy gives us an estimate for the temperature scale $`T_\varphi =J_V\rho _S(T)\pi ^2\lambda /2N_sk_B`$ of the phase fluctuations. In order to avoid phase fluctuations, we should be well below this temperature scale.

At temperatures higher than $`T_\varphi `$ the vortex can still be squeezed, but the phase fluctuations result in different main axis of the vortex position uncertainty ellipse in different layers. In Fig. 3 we plot the critical temperature for Bose-Einstein condensation, the critical temperature for the transition into a squeezed vortex, and the temperature scale of the phase fluctuations for typical parameter values. This figure indicates that the vortex spontaneously squeezes at a temperature that is easily accessible experimentally. Below this temperature there is a region in which the vortex is squeezed, but with fluctuating main axes. At even lower temperatures the phase fluctuations become negligible and long-range order over the size of the system is established.

## V Summary and conclusions

We studied the equilibrium squeezing of the vortex line in an optical lattice and predicted that strong squeezing is indeed possible in the experimentally realistic parameter regime. Although the kelvon interaction induced by the vortex displacement is the most important one, other mechanisms for kelvon interactions also exist. In principle, the kelvons are also coupled to the collective modes of the Bose-Einstein condensate, for example, to the quadrupole mode Martikainen and Stoof (2003, 2004). Such processes typically induce an effective attractive kelvon interaction and thus renormalize the value of the interaction strength, but do not make it positive. Therefore, we expect that the squeezing transition is robust with respect to coupling to the collective modes of the condensate.

###### Acknowledgements.

This work is supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM) and by the Nederlandse Organisatie voor Wetenschaplijk Onderzoek (NWO).


text:
# Polarization effects in metallic films perforated with a bidimensional array of subwavelength rectangular holes

## Introduction

For several years, the motivation for investigating surface plasmon properties takes a crucial importance in the context of plasmonics. In addition, the consideration of metallic gratings have been renewed in the context of the results report by Ebbesen et al . They reported on particularly attractive optical transmission experiments on periodical arrays of subwavelength cylindrical holes in a thin metallic layer deposited on glass, a dielectric substrate. Two specific characteristics of their results are the transmission which is higher than the simple sum of individual holes contributions and the interesting pattern of the zeroth order transmission versus wavelength .

During the last few years, these experimental results have received considerable attention and they have implied several theoretical and experimental works \[1-17\]. Many assumptions were suggested to explain the transmission properties of these devices. For instance, it has been suggested that these phenomena can be described in terms of short range diffraction of evanescent waves , or in terms of dynamical diffraction . Another explanation suggest the role of cavity resonance into the holes to explain the transmission enhancement .

However, the main results tend to prove that surface plasmons (SPs) play a key role in the features observed in the transmission curves \[1-7,9,10\]. Nevertheless, it remains many questions which must be answered to clarify the processes involved in these experiments. For instance, the exact role of SPs was not clearly assessed. Indeed, the observed transmission presents a set of peaks and minima. Many authors suggested that transmission peaks correspond to SPs resonance \[1-7,9,10\] whereas others suggested that SPs resonance were related to minima .

In a recent work , we have presented many results which tend to qualify these two hypothesies. We have shown that the concept of resonant Wood anomalies  can be invoked to interpret the role of SPs in the Ebbesen experiments. We have shown that the transmission pattern is better described by Fano’s profiles  correlated with interferences between non resonant processes and the resonant response of SPs coupled with nonhomogeneous diffraction orders . We have shown that each maximum of transmission (preceded by a minimum) corresponds to one maximum of a Fano’s profile (preceded by the related minimum) . Moreover, whereas such Fano’s profile in transmission is related to a resonant process, the location of its maximum (or minimum) does not necessarily correspond to the location of the resonance .

In such devices made of subwavelength hole arrays, polarization properties start only to be underlined in some experimental works . In the original experiments, for square gratings made of circular holes in chromium films for instance , it has been underlined that the zeroth order transmission pattern does not depend on the polarization of the incident light at normal incidence . Obviously, for a square grating made of rectangular holes for instance, one can expect a polarization dependent transmission. Indeed, for a square grating with circular holes the main axies of the grating can be inverted whereas for a square grating with rectangular holes these axies cannot be inverted without break the grating symmetry.

In the present paper, we will show the polarization dependency of the transmission for rectangular holes arrays contrary to circular holes arrays. We will present a simple analytical interpretation of our numerical results. However, we will underline the limits of this simple interpretation as one must not forgot the real complexity of the problem which yet requires a numerical study. Despite a simple appearance, the present problem is not trivial. Our results will lead us to suggest a new kind of polarizer which presents the advantages brought by the original Ebbesen devices .

## Numerical approach and studied device

In the following, our simulations rest on a coupled modes method associated with the use of the scattering matrix (S matrix) formalism. Taking into account the periodicity of the device, the permittivity is first described by a Fourier series. Then, the electromagnetic field is described by Bloch’s waves which can also be described by a Fourier series. In this context, Maxwell’s equations take the form of a matricial first order differential equation along the $`z`$ axis perpendicular to the $`x`$ and $`y`$ axies where the permittivity is periodic . The heart of the method is to solve this equation. One approach deals with the propagation of the solution step by step by using the S matrix formalism. More explicitly, we numerically divide the grating along the $`z`$ axis into many thick layers for which we calculate the S matrix. The whole S matrix of the system is obtained by using a special combination law applied step by step to each S matrices along the $`z`$ axis. Indeed, it is well know that S matrices and their combinations are much better conditionned than transfert matrices . In this way, we can calculate the amplitudes of the reflected and transmitted fields, for each diffracted order (which corresponds to a vector $`𝐠`$ of the reciprocal lattice) according to their polarization ($`s`$ or $`p`$) . Note that our algorithm has been accurately compared with other methods such as FDTD or KKR . In the present work the convergence is obtained from two harmonics only, i.e. for 25 vectors of the reciprocal lattice. Furthermore, there is no convergence problem associated with discontinuities here, and we do not need to use Li’s method .

As a continuation of our previous works, we turn our attention to the case of chromium films . The values of chromium and glass permittivities are those obtained from experiments . For instance, glass permittivity is approximately equal to $`2.10`$ for wavelengths ranging from $`1000`$ to $`1500`$ $`nm`$. Here, the thickness of the metallic film is taking equal to $`h=100`$ $`nm`$. The film is perforated with rectangular holes of sides $`a=350`$ $`nm`$ and $`b=700`$ $`nm`$ respectively. The shorter side is along the $`x`$ axis and the longer side is along the $`y`$ axis. The holes shape a square grating of parameter $`c=1000`$ $`nm`$ (fig.1a). Incident light is normal to the interface i.e. normal to the $`Oxy`$ plane. The incident light is linearly polarized and the orientation of the incident electric field $`𝐄_{inc}`$ is given by the angle $`\theta `$ between $`𝐄_{inc}`$ and the $`y`$ axis (fig.1a) in the $`Oxy`$ plane. Because the transmitted zeroth order is normal to the $`Oxy`$ plane, its polarization is described in the $`Oxy`$ plane too. The transmitted zeroth order electric field will be described by the complex electric field components $`E_x`$ and $`E_y`$ along the $`x`$ and $`y`$ axies respectively. The polarization of the transmitted zeroth order is then described by the modulus of the amplitudes $`\left|E_x\right|`$ and $`\left|E_y\right|`$ and the dephasing $`\delta =\mathrm{arg}\left\{E_y/E_x\right\}.`$ The amplitude of the incident light is equal to $`1`$ $`V.m^1`$.

## Results and analysis

First, we represent in fig.1b a schematic view of the reciprocal lattice related to the square grating of paramater $`c`$ in reduced coordinates. Each point corresponds to a reciprocal lattice node, i.e. to one diffraction order $`𝐠=\frac{2\pi }{c}(i,j)`$. Solid lines represent a contours representation of the modulus of the Fourier transform of the rectangular hole. For the most common diffraction orders we give the values of the modulus of the Fourier transform corresponding to the rectangular hole. Such values are representative of the typical order of magnitude for the relatives amplitudes of the diffraction orders. For instance the orders $`(\pm 1,0)`$ are associated with a relative amplitude of $`0.810`$ whereas it is of $`0.318`$ for the orders $`(0,\pm 1)`$, i.e. $`2.5`$ times less. Obviously, for a circular hole, the Fourier transform would be symmetric if we substitute $`i`$ by $`j`$ and vice versa. Then, the relative amplitude of the orders $`(\pm 1,0)`$ and $`(0,\pm 1)`$ would be equal. This difference between circular and rectangular holes is crucial. Indeed, in our recent work  we have specified the role of SPs in the Ebbesen experiment through the concept of resonant Wood anomalies . As explained above, the transmission pattern is described by Fano’s profiles correlated with interferences between non resonant processes and the resonant response of the surface plasmons coupled with nonhomogeneous diffraction orders. More precisely, the key process in those device appears as follows. We suppose that the incident light is polarized along the $`[1,0]`$ axis of the grating. The incident light then diffracts against the grating and generates nonhomogeneous resonant diffraction orders e.g. $`(1,0)`$ which is $`p`$ polarized. Such a resonant order is coupled with a surface plasmon (obviously $`(0,1)`$ is also generated but it is $`s`$ polarized and then it is not a resonant order, i.e. it is not coupled with surface plasmons). It becomes possible to excite this surface plasmon which leads to a feedback reaction on the order $`(1,0)`$. Then this diffraction order can diffract too and generate a contribution to the homogenous zeroth diffraction order $`(0,0)`$. In fact $`s`$ polarized diffraction orders give no contribution to the enhancement of the transmission, contrary to $`p`$ polarized diffraction orders . Obviously, the contribution of the order $`(1,0)`$ depends on the relative amplitude related to the Fourier transform of the hole. Now, if the incident light is polarized along the $`[0,1]`$ axis of the grating, the order $`(0,1)`$ becomes the resonant $`p`$ polarized order whereas $`(1,0)`$ becomes the non-resonant $`s`$ polarized order. Then two cases occur. If the Fourier transform is symmetric, the substitution of the $`p`$ polarized order $`(1,0)`$ by the $`p`$ polarized order $`(0,1)`$ let invariant the transmission. If the Fourier transform is not symmetric, then the contribution to zeroth order of both $`p`$ polarized orders will be not the same according to the incident polarization. In our case, if the incident light is polarized along the $`[1,0]`$ axis of the grating, the $`p`$ polarized order $`(1,0)`$ (which is associated with a relative amplitude of $`0.810`$ of the Fourier transform) gives an important contribution to the zeroth order transmission. If the incident light is polarized along to the $`[0,1]`$ axis of the grating, the $`p`$ polarized order $`(0,1)`$ (which is associated with a relative amplitude of $`0.318`$ of the Fourier transform) gives a weak contribution to the zeroth order transmission. Then the transmission $`t_x`$ for light polarized along the $`x`$ axis will be greater than the transmission $`t_y`$ for light polarized along the $`y`$ axis. This let to some issues. As explained above, the electric field of the incident light can be written as

$$𝐄_{inc}=E_{inc}\left[\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right]$$
(1)
Assumming the fact that both components of the incident electric field are not transmitted in the same way, the zeroth order transmitted field is written as

$$𝐄=\left[\begin{array}{c}E_x=t_xE_{inc}\mathrm{sin}\theta \\ E_y=t_yE_{inc}\mathrm{cos}\theta \end{array}\right]$$
(2)
where $`t_x`$ and $`t_y`$ respectively are the transmission coefficients along the $`x`$ and $`y`$ axis respectively. The modulus of the amplitudes are then

$$\{\begin{array}{c}\left|E_x\right|=\left|t_xE_{inc}\right|\mathrm{sin}\theta \\ \left|E_y\right|=\left|t_yE_{inc}\right|\mathrm{cos}\theta \end{array}$$
(3)
As the zeroth order transmission can be written as

$$T=\sqrt{\frac{\epsilon _s}{\epsilon _v}}\frac{1}{\left|E_{inc}\right|^2}\left(\left|E_x\right|^2+\left|E_y\right|^2\right)$$
(4)
it is easier to show that one can write $`T`$ as

$$T=T_{90^\text{o}}+\left(T_{0^\text{o}}T_{90^\text{o}}\right)\mathrm{cos}^2\theta $$
(5)
where $`T_{0^\text{o}}=\sqrt{\frac{\epsilon _s}{\epsilon _v}}\left|t_y\right|^2`$ and $`T_{90^\text{o}}=\sqrt{\frac{\epsilon _s}{\epsilon _v}}\left|t_x\right|^2`$ ($`\epsilon _s`$ and $`\epsilon _v`$ are the permittivities of the substrate and the vaccum respectively). At last, note that the dephasing $`\delta `$ is such that

$$\delta =\mathrm{arg}\left\{\frac{t_y}{t_x}\right\}$$
(6)
One shows that eq.5 is similar to the well known Malus Law for polarizer. To illustrate our assumptions, fig.2 shows the zeroth order transmission against wavelength for different $`\theta `$ values. We clearly show the dependence of the transmission pattern on the incident polarization. For $`\theta =90^\text{o},`$ i.e. the incident electric field is along the $`x`$ axis, the transmission is maximal and identical to what one observes in the Ebbesen experiments . It is shown that the transmission increases with the wavelength, and that it is characterized by sudden changes in the transmission marked 1 to 2 on the figure. These values are shifted toward larger wavelengths when the grating size increases . Moreover, for $`\theta =0^\text{o}`$, i.e. when the incident electric field is along the $`y`$ axis, the transmission is minimal. This result is, at least qualitatively, in agreement with our hypothesis. Note that wavelengths 1 and 2 are close to Rayleigh’s wavelengths. In diffraction gratings such wavelengths are those for which a diffracted order becomes tangent to the plane of the grating. At normal incidence, and for a diffraction order $`𝐠=\frac{2\pi }{c}(i,j)`$, Rayleigh’s wavelength is defined as

$$\lambda _R^{u,i,j}=c\sqrt{\epsilon _u}(i^2+j^2)^{\frac{1}{2}}$$
(7)
where $`\epsilon _u`$ represents either the permittivity of the vacuum ($`\epsilon _v`$), or of the dielectric substrate ($`\epsilon _d`$). For the diffraction orders such that $`(\pm 1,0)`$ or $`(0,\pm 1),`$ the values of Rayleigh’s wavelengths in the present case are $`1000`$ $`nm`$ for the vacuum/metal interface, and $`1445.29`$ $`nm`$ for the substrate/metal interface.

Fig.3 shows two zoomed plots around wavelengths 1 (fig.3a) and 2 (fig.3b). In fig.3a, one shows that the transmission minima are localized at the same wavelength whatever $`\theta .`$ Surprisingly, in this narrow wavelength domain, transmission for $`\theta =0^\text{o}`$ is greater than for $`\theta =90^\text{o}`$ contrary to previously. The minimum which corresponds to the wavelength 1, is located at $`1001.44`$ $`nm`$ just after the related Rayleigh’s wavelength. Fig.3b shows that transmission minima are weakly shifted toward large wavelengths as $`\theta `$ increases (minima locations for each $`\theta `$ value, are marked by vertical black dashs). For $`\theta =0^\text{o}`$, the minimum is located at $`1445.29`$ $`nm`$, i.e. the wavelength 2 which corresponds to the related Rayleigh’s wavelength. For $`\theta `$ greater than zero, minimum is shifted. Moreover, one notes a sudden change in the transmission at the Rayleigh’s wavelength whatever $`\theta `$. The minima observed here correspond to minima of Fano’s profiles . In fig.3a, the minima are not of the same kind of the minima in Fig.3b. This is not true minima of the Fano’s profile. All occur as if the minimum of the Fano’s profile disappears behind the Rayleigh’s wavelength towards low wavelength. In other words, the minima in Fig.3a come from the cutoff and the discontinuity introduced in Fano’s profiles at the Rayleigh’s wavelength. This had been widely explained in ref. .

For more details, in fig.4, we have reported the wavelengths $`\lambda _m`$ related to the transmission minima (related to fig.3b) against $`\theta `$ (black dots). For comparison, and by analogy with eq.5, the solid line is related to the following equation

$$\lambda _{m,\theta }=\lambda _{m,90^\text{o}}+(\lambda _{m,0^\text{o}}\lambda _{m,90^\text{o}})\mathrm{cos}^2\theta $$
(8)
Quantitatively, fig.5 shows the variation of the transmission against the angle of polarization $`\theta `$ for both wavelengths $`1300`$ $`nm`$ and $`1900`$ $`nm`$ (fig.5a) and both wavelengths $`1147.64`$ $`nm`$ and $`1001.44`$ $`nm`$ (fig.5b) in the neighborhood of the above mentionned minima. The numerical computations give results indicated by the dots on the figure. Using the numerical values of the transmission for $`\theta =90^\text{o}`$ and $`\theta =0^\text{o}`$ respectively, i.e. $`T_{90^\text{o}}`$ and $`T_{0^\text{o}}`$, one plots the transmission against $`\theta `$ for both wavelength by using eq.5 (solid lines). One observes then a quantitative agreement between the pattern from eq.5 and the results obtained numerically. Note the weakly decreasing transmission at $`1001.44`$ $`nm`$ as already mentionned, which is also in full agreement with eq.5.

In order to more detail our approach, we present on fig.6a the amplitudes $`\left|E_x\right|`$ and $`\left|E_y\right|`$ of the zeroth transmitted order against wavelength for two different values of $`\theta .`$ One shows that the $`\left|E_y\right|`$ amplitude does not exhibit high convex regions in transmission contrary to the $`\left|E_x\right|`$ amplitude. As explained, the convex regions in the zeroth order transmission (i.e. the high transmission domains between minima $`1`$ and $`2,`$ and just after minimum $`2`$) are induced by the scattering of non-homogenous resonant diffraction orders such as $`(\pm 1,0)`$ or $`(0,\pm 1)`$. For the $`\left|E_x\right|`$ amplitude, the orders $`(\pm 1,0)`$ are associated with a relative amplitude of $`0.810`$ of the Fourier transform of the hole profile. This large contribution of these orders enables an enhanced transmission via multi-scattering. This effect appears as the reason of the existence of the convex region of high transmission as shown in previous work . On the contrary, for the $`\left|E_y\right|`$ amplitude the resonant orders $`(0,\pm 1)`$ are associated with a relative amplitude of $`0.318`$ of the Fourier transform of the hole profile. This contribution is not large enough to obtain a significant enhanced transmission. For $`\theta =30^\text{o}`$, the $`\left|E_y\right|`$ amplitude contribution to the whole transmission is greater than the $`\left|E_x\right|`$ contribution. Although the transmission along the $`x`$ axis is greater than along the $`y`$ axis, as $`𝐄_{inc}`$ components depend on $`\theta `$ (see eq.1) the $`x`$ axis component of the incident electric field is small enough to lead $`\left|E_x\right|`$ to be lower than $`\left|E_y\right|`$ (see eq.2). On the contrary for $`\theta =60^\text{o}`$, the $`\left|E_x\right|`$ amplitude contribution to the whole transmission is greater than the $`\left|E_y\right|`$ contribution. Indeed, the transmission along the $`x`$ axis is obviously greater than along the $`y`$ axis, but now the $`x`$ axis component of the incident electric field is also greater than the one along the $`y`$ axis. Then, the amplitude $`\left|𝐄\right|`$ pattern for $`\theta =60^\text{o}`$ is essentially defined by the pattern of the amplitude $`\left|E_x\right|`$ whereas the amplitude $`\left|𝐄\right|`$ pattern for $`\theta =30^\text{o}`$ is essentially defined by the pattern of the amplitude $`\left|E_y\right|`$. In addition, the whole zeroth order transmission would be greater for $`\theta =60^\text{o}`$ than for $`\theta =30^\text{o}.`$ Again, these results are in qualitative agreement with our assumptions. Quantitatively, fig.6b shows the variation of the amplitudes $`\left|E_x\right|`$ and $`\left|E_y\right|`$ against the angle of polarization $`\theta `$ for both wavelength $`1300`$ $`nm`$ and $`1900`$ $`nm`$. The numerical computations give the results indicated by the dots on the figure. Using the numerical values of the amplitudes for $`\theta =90^\text{o}`$ and $`\theta =0^\text{o}`$, one plot the transmission against $`\theta `$ for both wavelength by using eq.3 (solid lines). As previously, one observes a quantitative agreement between the pattern from eq.3 and the results obtained numerically.

In the possibility of later experimental results, we wish to add some details for comparison. Fig.7a gives the dephasing $`\delta `$ against wavelength for both angle $`\theta =60^\text{o}`$ and $`\theta =30^\text{o}.`$ One shows the peaks close to the minima of transmission, i.e. close to the surface plasmon resonance according to our previous work . Despite the absence of an explicit dependence against $`\theta `$ for $`\delta `$ according to eq.6, one notes that such a dependence exist after all. So, fig.7b shows the dephasing against $`\theta `$ for both wavelength $`1300`$ $`nm`$ and $`1900`$ $`nm`$. One shows for $`\theta =0`$ that the dephasing is equal to zero, and the polarization of the zeroth transmitted order is thus the same as the incident field. Though the dephasing tends to be different from zero when $`\theta =90^\text{o}`$, as $`\left|E_y\right|`$ tends to zero (see fig.6) the polarization is the same as the incident field. By contrast, the transmission polarization will be elliptical for others values of $`\theta `$. Such a dependence of $`\delta `$ against $`\theta `$ and wavelength probably express a dependence of $`t_x`$ and $`t_y`$ against $`\theta `$. This maybe corresponds to a propagating constant against $`z`$-axis that is different according to the $`x`$-axis or $`y`$-axis. In addition, as shown in ref., $`t_x`$ and $`t_y`$ result from many contributions such as SPs via multi-scattering. Because SPs exhibit a polarization dependency, this could explain the $`t_x`$ and $`t_y`$ behaviour. Nevertheless, as one shows, the dependence of $`\delta `$ against $`\theta `$ and wavelength is not well understood and we hope to clarify this point in a later work. Moreover, it will be necessary to search for the effects induced by a parameters change.

At last, it is striking that our numerical results coincide with our simple and intuitive explanation which leads to predict that a square grating with rectangular holes can behave like a polarizer. Indeed, the whole multiscattering processes are complex and there are numerous contributions to the zeroth order transmission as explained above. It was then not necessarily obvious, that our intuitive model match up with our numerical results, all the more that the results shown in fig.7 underline the limits of our simple interpretation.

## Conclusion

In the present paper, we have shown that contrary to the original experiments (which used a grating made of circular holes) , the zeroth order transmission pattern depends on the polarization of the incident light when using rectangular holes. The results presently obtained match up with our previous work . In this context, we obtain a device with properties globally similar to those of the original Ebbesen experiments, but which exibits in addition, some properties of a polarizer. Assuming the recent technological interest for this kind of metallic gratings, we think that the kind of device introduced in the present paper could lead to new technological applications. At last, note that, during the submission process, Gordon et al  have published experimental results concerning transmission through a bidimensional array of elliptical holes in metallic film. Their results clearly establish the polarization dependence of the transmission with asymmetric holes. These experimental results corroborate the present theoretical approach and in this way, the experimental and theoretical works supply many complementary results.

## Acknowledgments

We acknowledge the use of Namur Scientific Computing Facility (Namur-SCF), a common project between the FNRS, IBM Belgium, and the Facultés Universitaires Notre-Dame de la Paix (FUNDP).

This work was carried out with support from EU5 Centre of Excellence ICAI-CT-2000-70029 and from the Inter-University Attraction Pole (IUAP P5/1) on ”Quantum-size effects in nanostructured materials” of the Belgian Office for Scientific, Technical, and Cultural Affairs.

## References

 T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff, Nature (London) 391, 667 (1998)

 H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec, Phys. Rev. B, 58, 6779 (1998)

 U. Schröter, D. Heitmann, Phys. Rev. B, 58, 15419 (1998)

 T. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen, JOSA B, 16, 1743 (1999)

 J.A. Porto, F.J. Garcia-Vidal, J.B. Pendry, Phys. Rev. Lett., 83, 2845 (1999)

 D.E. Grupp, H.J. Lezec, T.W. Ebbesen, K.M. Pellerin, T. Thio, Appl. Phys. Lett. 77 (11) 1569 (2000)

 T. Thio, H.J. Lezec, T.W. Ebbesen, Physica B 279, 90 (2000)

 E. Popov, M. Nevière, S. Enoch, and R. Reinisch, Phys. Rev. B 62, 16100-16108 (2000)

 A. Krishnan, T. Thio, T. J. Kim, H. J. Lezec, T. W. Ebbesen, P.A. Wolff, J. Pendry, L. Martin-Moreno, F. J. Garcia-Vidal, Opt. Commun., 200, 1 (2001)

 L. Martin-Moreno, F.J. Garcia-Vidal, H.J. Lezec, K.M. Pellerin, T. Thio, J.B. Pendry, T.W. Ebbesen, Phys. Rev. Lett., 86, 1114 (2001)

 J.-M. Vigoureux, Optics Comm., 198, 4-6, 257 (2001)

 Q. Cao, Ph. Lalanne, Phys. Rev. Lett., 88 (5), 057403 (2002)

 M.M.J. Treacy, Phys. Rev. B 66, 195105, (2002)

 M. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux, Phys. Rev. B, 67, 085415 (2003)

 Z. Sun, Y.S. Jung, H.K. Kim, Appl. Phys. Lett. 83 (15) 3021 (2003)

 E. Altewischer, M.P. van Exter, J.P. Woerdman, JOSA B, 20, 1927 (2003)

 R. Gordon, A. G. Brolo, A. McKinnon, A. Rajora, B. Leathem, and K. L. Kavanagh, Phys. Rev. Lett. 92, 037401 (2004)

 R.W. Wood, Phys. Rev. 48, 928 (1935)

 V.U. Fano, Ann. Phys. 32, 393 (1938)

 J.P. Vigneron, F. Forati, D. André, A. Castiaux, I. Derycke, A. Dereux, Ultramicroscopy, 61, 21 (1995)

 J.B. Pendry, P.M. Bell, NATO ASI Series E Vol. 315 (1995)

 V. Lousse, K. Ohtaka, Private Communication (2001)

 L. Li, JOSA A 13 (9) 1870 (1996)

 L. Li, JOSA A 14 2758 (1997)

 D.W. Lynch, W.R. Hunter, in Handbook of Optical Constants of Solids II, edited by E.D. Palik (Academic Press, New York, 1991)

## Captions

Figure 1 : (a) Diagrammatic view of the system unders study, (b) Simplified view of the Fourier transform of the holes grating.

Figure 2 : zeroth order transmission against wavelength for various polarization, i.e. for various values of $`\theta .`$

Figure 3 : close view of zeroth order transmission against wavelength for various polarization, i.e. for various values of $`\theta .`$ (a) in the vicinity of wavelength (1). (b) in the vicinity of wavelength (2).

Figure 4 : Wavelengths related to the transmission minima against $`\theta `$ in the vicinity of wavelength (2) (see fig.3b).

Figure 5 : (a), zeroth order transmission against $`\theta `$ for both wavelength $`1300`$ $`nm`$ and $`1900`$ $`nm`$. (b), zeroth order transmission against $`\theta `$ for both wavelength $`1001.44`$ $`nm`$ and $`1447.64`$ $`nm`$.

Figure 6 : (a) $`\left|E_x\right|`$ and $`\left|E_y\right|`$ amplitudes against wavelength for $`\theta =30^{}`$ and $`\theta =60^{}.`$ (b) $`\left|E_x\right|`$ and $`\left|E_y\right|`$ amplitudes against $`\theta `$ for both wavelength $`1300`$ $`nm`$ and $`1900`$ $`nm`$.

Figure 7 : (a) dephasing $`\delta `$ against wavelength for $`\theta =30^{}`$ and $`\theta =60^{}.`$ (b) dephasing $`\delta `$ against $`\theta `$ for both wavelength $`1300`$ $`nm`$ and $`1900`$ $`nm`$.