text:
# A 250 GHz Survey of High Redshift QSOs from the Sloan Digital Sky Survey

## 1 Introduction

The existence of massive black holes at the centers of galaxies has long been predicated on consideration of the energetics of active galactic nuclei (AGN): accretion of matter onto a massive black hole is an order of magnitude more efficient at converting mass into energy than is stellar nucleosynthesis (Begelman, Blandford, & Rees 1984). The last few years have seen an explosion of dynamical evidence for massive black holes at the center of galaxies, and in the direct measurement of black hole masses from the dynamics of circumnuclear regions (Richstone et al. 1998; Miyoshi et al. 1996; Ghez et al. 1998; Genzel et al. 2000; Tanaka et al. 1995). These observations indicate that the overwhelming majority of spheroidal galaxies in the nearby universe contain massive black holes (Kormendy 2000). Moreover, a clear correlation has been found between the black hole mass and the velocity dispersion of the stars in the spheroid (Gebhardt et al. 2000, Ferrarese & Merritt 2000). This remarkably tight correlation suggests a ‘causal connection between the formation and evolution of the black hole and the bulge’ (Gebhardt et al. 2000; Kormendy 2000; Richstone et al. 1998).

Support for this general idea comes from: (i) the observed space density of high redshift Quasi-Stellar Objects (QSOs), which would require a massive black hole in most low $`z`$ spheroidal galaxies (Kormendy & Ho 2000), (ii) the observation of similar rapid increases in the space densities of AGN and starburst galaxies from $`z=0`$ to $`z=2`$ (Boyle & Terlevich 1998, Blain et al. 1999), and (iii) the observation of co-eval AGN and star formation activity in some nearby nuclear starburst galaxies (Rigopoulou et al. 1999; Carilli, Wrobel, & Ulvestad 1999). On the other hand, the QSO population shows an abrupt turn-over in the co-moving number density at $`z3`$ (Schmidt, Schneider, & Gunn 1995a; Kennefick, Djogovski, & de Carvahlo 1995; Fan et al. 2001b). It is unclear whether such a turnover exists in the starburst population (Blain et al. 1999, Dunlop 2000). Also, questions remain about the relative timescales for the commencement and duration of the starburst and AGN processes (Sanders & Mirabel 1996), and the origin at very high redshift of the ‘seed’ black hole (Richstone et al. 1998).

A particularly intriguing observation that has fueled the debate over co-eval starbursts and AGN at high redshift is the recent detection at 250 GHz of copious thermal dust emission from high redshift QSOs by Omont et al. (1996a). Omont et al. observed a sample of $`z>4`$ QSOs from the Automatic Plate Measuring (APM) survey and found that 6 of 16 sources show dust emission with 250 GHz flux densities of 3 mJy or greater, with implied Far Infrared (FIR) luminosities, $`\mathrm{L}_{\mathrm{FIR}}>10^{12}`$ L, and dust masses $`10^8`$ M. Follow-up observations of three of these dust-emitting QSOs revealed CO emission as well, with implied molecular gas masses $``$ few $`10^{10}`$ M (Guilloteau et al. 1997, 1999; Ohta et al. 1996; Omont et al. 1996b; Carilli, Menten, & Yun 1999).

Two different mechanisms have been proposed for heating the dust in high redshift QSOs. Given the large dust and gas masses, Omont et al. (1996a) make the circumstantial argument that star formation is inevitable, and hence that the dominant dust heating mechanism may be star formation. Supporting evidence comes from deep radio observations at 1.4 GHz, which show that the ratio of the radio continuum to submillimeter continuum luminosity from some of these sources is consistent with the well established radio-to-far IR correlation for low redshift star forming galaxies (Yun et al. 2000). And Omont et al. (2001) make the simple point that copious star formation is required in order to produce the dust. The implied star formation rates are so extreme, $``$ 103 M year-1, that a significant fraction of the stars in the QSO host galaxy can be formed in $`<10^9`$ years.

Alternatively, the dust may be heated by radiation from the AGN. Sanders et al. (1989) argue that radiation from the AGN is the dominant dust heating mechanism, since in most cases it requires the absorption of only a small fraction ($`20\%`$) of the AGN UV luminosity by dust, even for the high redshift sources (Carilli et al. 2000). An important constraint on AGN dust heating models is provide by the minimum size of about 1 kpc for the dust emission region one derives from the observed luminosity and dust temperature (see section 3). Hence, while the optical ‘big blue bump’ emission is thought to arise in a hot ($`10^{4.5}`$ K) accretion disk on pc-scales, the rest-frame far IR emission must be due to warm dust on kpc-scales. In the Sanders et al. (1989) model for AGN-powered IR emission from QSOs, dust heating on kpc scales is accomplished by assuming that the dust is distributed in a kpc-scale warped disk, thereby allowing UV radiation from the AGN to illuminate the outer regions of the disk. An alternate mechanism for large-scale dust heating has been proposed by Maloney, Hollenbach, & Tielens (1996) in which hard X-rays penetrate the very high column density gas in the vicinity of the AGN.

In this paper we present extensive observations at cm and mm wavelengths of a sample of high redshift QSOs from the Sloan Digital Sky Survey (SDSS; York et al. 2000). Our observations include sensitive radio continuum imaging at 1.4 and 43 GHz with the Very Large Array (VLA), and photometry at 250 GHz using the Max-Planck Millimeter Bolometer array (MAMBO) at the IRAM 30m telescope. These observations are a factor three more sensitive than previous studies of high redshift QSOs at these frequencies (Schmidt et al. 1995b; Omont et al. 1996a; Stern et al. 2000). Details of the multifrequency radio continuum imaging will be given elsewhere (Rupen et al. 2001, in preparation). In this paper we present the radio results that are relevant to the analysis and interpretation of the dust continuum emission. We compare the dust emission properties of these sources with their optical and radio continuum and spectral properties, and with those of lower redshift sources. We assume $`\mathrm{H}_\mathrm{o}=50`$ km s-1 Mpc-1, $`\mathrm{q}_\mathrm{o}=0.5`$, $`\mathrm{\Lambda }=0`$, and we define the spectral index, $`\alpha `$, as a function of frequency, $`\nu `$, as: $`\mathrm{S}_\nu \nu ^\alpha `$.

## 2 The Sample

The QSO sample is the result of optical spectroscopy of objects of unusual color from the northern Galactic Cap and the Southern Equatorial Stripe of the SDSS. The survey has yielded over 100 QSOs with $`z3.6`$, including eight of the ten highest redshift QSOs known (Fan et al. 1999, 2000b, 2001a,b; Schneider et al. 2000a,b, 2001; Zheng et al. 2000).

We observed 41 of these high redshift QSOs, taken from Fan et al. (1999), Fan et al. (2000a), and Schneider et al. (2000a). The QSO properties are listed in Table 1. The observed sources span a range of $`\mathrm{M}_\mathrm{B}=26.1\mathrm{to}28.8`$, and a redshift range of $`z=`$ 3.6 to 5.0. Comparative numbers for the APM sample of 16 QSOs observed by Omont et al. (1996a) are $`\mathrm{M}_\mathrm{B}=26.8\mathrm{to}28.5`$, and $`z=`$ 4.0 to 4.7.

The sources in Table 1 have a mean optical spectral power law index, $`\alpha =0.67`$, with a dispersion of $`0.29`$, a mean readhift, $`z=4.20`$, with a dispersion of $`0.44`$, and a mean absolute blue magnitude, $`\mathrm{M}_\mathrm{B}=27.07`$, with a dispersion of $`0.57`$. We can compare these values to those of the 39 QSOs at $`z>3.6`$ from the SDSS Fall Equatorial stripe presented in Fan et al. (2001b). The Fall Equatorial stripe sample is a statistically complete sample of high-$`z`$ QSOs selected according to the same criteria as the sources listed in Table 1. The Fall Equatorial stripe sample has a mean optical spectral power law index, $`\alpha =0.79`$, with a dispersion of $`0.34`$, a mean readhift, $`z=4.06`$, with a dispersion of $`0.39`$, and a mean absolute blue magnitude, $`\mathrm{M}_\mathrm{B}=26.82`$, with a dispersion of $`0.58`$. Hence, within the scatter of the distributions, the sources in Table 1 can be considered representative of the complete high-$`z`$ QSO sample from the Fall Equatorial stripe.

## 3 Observations

Observations were made using MAMBO (Kreysa et al. 1999) at the IRAM 30m telescope on Pico Veleta in Spain, in December 1999 and February 2000. MAMBO is a 37 element bolometer array sensitive between 190 and 315 GHz. The half-power sensitivity range is 210 to 290 GHz, but the overall profile is asymmetric, with a sharp rise in sensitivity at low frequency, then a gradual decrease in sensitivity to higher frequency. Convolving the frequency response of the bolometer array with the atmospheric transmission curve (which also decreases with increasing frequency), and with the rising spectrum of a typical dust emitting source at high redshift, leads to an effective central frequency of 250 GHz. The beam for the feed horn of each bolometer is matched to the telescope beam of 10.6′′, and the bolometers are arranged in a hexagonal pattern with a beam separation of 22′′. Observations were made in standard on-off mode, with 2 Hz chopping of the secondary by 50′′ in azimuth. The data were reduced using the MOPSI software package of the Max-Planck-Institut für Radioastronomie (Zylka 1998). Pointing was monitored every hour, and was found to be repeatable to within 2′′. The sky opacity was monitored every hour; zenith optical depths ranged between 0.23 and 0.36. Gain calibration was performed using observations of Mars, Uranus, and Ceres. We estimate a 20$`\%`$ uncertainty in absolute flux density calibration based on these observations. The target sources were centered on the central bolometer in the array, and the temporally correlated variations of the sky signal (sky noise) detected in the remaining bolometers were subtracted from all the bolometer signals. The total on-target plus off-target observing times for the sources are listed in Table 1. The typical rms sensitivity in 1 hour was 0.5 mJy for the on-source bolometer beam.

The 1.4 GHz VLA observations were made in the A (30 km) and BnA (mixed 30 km and 10 km) configurations on July 8, August 14, September 30, and October 8, 1999, using a total bandwidth of 100 MHz with two orthogonal polarizations. Each source was observed for a total of 2 hours, with short scans made over a large range in hour angle to improve Fourier spacing coverage. Standard phase and amplitude calibration were applied, as well as self-calibration using background sources in the telescope beam. The absolute flux density scale was set using observations of 3C 48. The final images were generated using the wide field imaging and deconvolution capabilities of the AIPS task IMAGR. The theoretical rms noise value is between 20 and 30 $`\mu `$Jy beam-1, depending on the range of telescope elevations over which the source was observed, and for about half the sources the measured noise values are in this range. For the remaining sources, the noise levels are significantly higher due to side-lobe confusion by bright sources. The Gaussian restoring CLEAN beams were between 1.5′′ and 3′′ FWHM.

Some of the sources had previous VLA detections as part of the FIRST survey (SDSS J1053-0016, SDSS J1235-0003, SDSS J1412-0101; Becker, White, & Helfand 1995). We re-observed these sources, and obtained the same flux densities to within 15$`\%`$ of the FIRST values in all cases.

We chose four of the dust emitting sources, (SDSS J0150+0041, SDSS J0255+0048, SDSS J0338+0021, SDSS J1112+0050) as a sub-sample for 43 GHz observations, in order to search for high frequency, flat spectrum emission that might be synchrotron self-absorbed at 1.4 GHz. Observations were made on Feb 20, 2000, using the VLA in the B (10 km) configuration in standard continuum mode. Fast switching phase calibration was employed (Carilli & Holdaway 1999), and observing times were 0.5 hr at 43 GHz, resulting in an rms noise of 0.35 mJy beam-1.

## 4 Results and Analysis

The results of our 250 GHz survey of high redshift QSOs are given in Table 1. These include: (i) flux densities at 250 GHz (S250; column 3), (ii) flux densities at 1450 $`\mathrm{\AA }`$ (S1450; column 4), (iii) flux densities at 1.4 GHz at the position of the optical QSO (S1.4; column 5), (iv) absolute blue magnitudes (column 6; Fan et al. 1999, 2000a), (v) observing time at 250 GHz, and (vi) notes on source properties (column 8). In the following analysis we consider a source with a measured flux density three times greater than the rms noise (3$`\sigma `$) to be a valid detection. Such sources are marked with an asterisk in Table 1. We should also point out that the mm, cm, and optical observations were not simultaneous, and hence that variability could affect the observed spectra at some level.

We detect 16 of the 41 sources at S$`{}_{250}{}^{}`$ 1.4 mJy. Only four of the sources have S$`{}_{250}{}^{}3`$ mJy. A higher detection rate at the limit of S$`{}_{250}{}^{}3`$ mJy was found in the survey of the APM QSO sample by Omont et al. (1996), in which they detected 6 of 16 sources. This difference in detection rates appears to be statistically significant: if there were a true universal fraction of 10$`\%`$ of high-$`z`$ QSOs with S$`{}_{250}{}^{}3`$ mJy, then there is about a 1$`\%`$ chance that in 16 objects 6 or more would meet this criterion. The optical/near IR color selection criteria for the APM sample are similar to those of the SDSS sample. The one significant difference is that the rest frame blue luminosities for the APM are somewhat higher on average (by about 0.5 mag), than those of the SDSS sources in Table 1 (see section 2). However, given the lack of a strong correlation between MB and S250 (see Figure 2 below), the cause for the possibly higher detection rate at high S250 for the APM sample relative to the SDSS sample remains unknown. Further observations of the mm properties of high $`z`$ QSOs that are in progress may clarify this issue (Omont et al. 2001).

Considering the radio properties of the sources in Table 1, nine sources in Table 1 are detected at 1.4 GHz at $`3\sigma `$. For the nine radio-detected sources in Table 1, eight are unresolved, with size limits (FWHM) of about 1′′, while SDSS J0232-0000 is possibly extended, with a size $`1.3^{\prime \prime }`$. Miller et al. (1990) suggested a division between radio loud and radio quiet QSOs at a rest frame 5 GHz spectral luminosity of $`10^{33}`$ erg s-1 Hz-1. This corresponds to a flux density of S1.4 = 1 mJy for a source at $`z=4.2`$ assuming a radio spectral index of –0.8. According to this criterion, 2 of 16 SDSS QSOs detected at 250 GHz can be considered radio-loud (SDSS J1235–0003 and SDSS J1412–0101), while 2 of 25 sources not detected at 250 GHz can be considered radio-loud (SDSS J0153–0011 and SDSS J1053–0016). Both fractions are consistent with the value of 10$`\%`$ of optically selected high redshift QSOs being radio loud according to this criterion (Schmidt et al. 1995b; Stern et al. 2000).

### 4.1 Evidence for thermal dust emission at 250 GHz

We first address the question of whether the 250 GHz emission is thermal dust emission, or non-thermal synchrotron radiation. For most of the sources detected at 250 GHz, the 1.4 GHz flux density (or upper limit) is an order of magnitude or more below the 250 GHz flux density. Two of the sources have 1.4 GHz flux densities comparable to, or larger than, the 250 GHz flux density. The source SDSS J1412–0101 has S$`{}_{1.4}{}^{}=3.9`$ mJy and S$`{}_{250}{}^{}=4.5`$ mJy. Observations of this source at 8.4 GHz revealed a 1.1 mJy source (Rupen et al. 2001), such that the cm source has a falling spectrum of index –0.7, making it unlikely that the mm continuum is a continuation of the synchrotron emission spectrum from the AGN. The source SDSS J1235–0003 has S$`{}_{250}{}^{}=1.6`$ mJy, and S$`{}_{1.4}{}^{}=4.5`$ mJy. Observations on the same day as those at 1.4 GHz found a 5 GHz flux density of S$`{}_{5}{}^{}=17`$ mJy, and a 15 GHz flux density of S$`{}_{15}{}^{}=7`$ mJy. The source was unresolved at all frequencies. In this case, it remains possible that the 250 GHz emission is non-thermal.

It is possible that emission from a compact AGN could be synchrotron self-absorbed at low frequency. To check this idea, we performed 43 GHz observations of 4 of the 250 GHz detected sources, as listed in section 3; none of the sources were detected. Assuming a 2$`\sigma `$ upper limit of 0.7 mJy at 43 GHz implies a rising spectrum between 43 and 250 GHz of index $`\alpha >0.5`$. While this limit is well below the index expected for thermal dust emission ($`\alpha _{43}^{250}3`$ for a source at $`z=4.2`$), it does argue against compact, synchrotron self-absorbed sources as the origin of the 250 GHz emission, as such sources typically show mm spectral indices $`\alpha _{\mathrm{mm}}0\pm 0.3`$ (Sanders et al. 1989). In the following discussion we will assume that the 250 GHz emission from the sources in Table 1 is thermal dust emission, with the possible exception of SDSS J1235–0003. These data imply that in mm surveys of high redshift QSOs, the fraction of flat spectrum, radio loud AGN with synchrotron emission extending into the mm-regime is at most a few percent. Studies of lower redshift QSOs show that the fraction of flat spectrum radio loud sources in optically selected samples is also only a few percent (Hopkins et al. 1999).

A source with S$`{}_{250}{}^{}=1`$ mJy at $`z=4.2`$ and a dust spectrum of the type seen in the low redshift starburst galaxy Arp 220 (corresponding roughly to a modified black body spectrum of dust emissivity index = 1 and temperature = 50 K) has an FIR luminosity, LFIR, of $`1.1\times 10^{12}`$ L, where $`\mathrm{L}_{\mathrm{FIR}}`$ corresponds to the integrated luminosity between 42 and 122 $`\mu `$m (Condon 1992). If we assume that star-formation is the mechanism giving rise to the dust emission, we can use the relations in Omont et al. (2001) and Carilli et al. (2001) to derive the dust mass, $`\mathrm{M}_\mathrm{D}`$, and total star formation rate, SFR, from $`\mathrm{S}_{250}`$. The relations are: $`\mathrm{M}_\mathrm{D}=1.1\times 10^8\times \mathrm{S}_{250}`$ M, and $`\mathrm{SFR}=360\times \mathrm{S}_{250}\mathrm{M}_{}\mathrm{year}^1`$, with S250 in mJy. This again assumes $`\mathrm{T}_\mathrm{D}=50`$ K.

A conservative lower limit to the solid angle, $`\mathrm{\Omega }_s`$, of a dust emitting source can be derived from the observed flux density making the extreme assumption of optically thick emission (Downes et al. 1999). Assuming $`\mathrm{T}_\mathrm{D}50\mathrm{K}`$ (Benford et al. 1999) and $`z4`$, leads to: $`\mathrm{\Omega }_\mathrm{s}>0.0043\mathrm{S}_{250}\mathrm{arcsec}^2,`$ where S250 is the 250 GHz flux density in mJy. This implies absolute lower limits to the angular diameters of the sources in Table 1 between $`0.08^{\prime \prime }`$ and $`0.15^{\prime \prime }`$, corresponding to physical sizes of 0.5 to 1 kpc. Given the likelihood that the observed 250 GHz emission is from an optically thin grey body, it is almost certain that the sources are larger than these lower limits. Also, this size corresponds to the total emitting area. Multiple smaller regions within the 10.6′′ telescope beam are certainly possible (see section 5). Observations at mm wavelengths with sub-arcsecond resolution are required to constrain the source sizes.

### 4.2 Trends with Redshift and MB

Figure 1 shows values of S250 versus redshift for the 41 QSOs in Table 1. The solid line is the expected flux density of Arp 220. This curve shows clearly the effect of the large ‘inverse $`K`$’ correction for thermal dust emission at mm and submm wavelengths, with the observed flux density of a source such as Arp 220 remaining roughly constant over the observed redshift range (Blain and Longair 1993). There is an interesting trend for a lower detection rate for the SDSS QSOs with increasing redshift: for $`z<4.4`$ we detect 14 of 28 sources while for $`z>4.4`$ we detect 2 of 13 sources. A larger sample of sources, over a larger redshift range, is required to verify this trend.

Figure 2 shows values of LFIR versus MB for the 41 QSOs in Table 1. At the limit of L$`{}_{\mathrm{FIR}}{}^{}>1.8\times 10^{12}`$ L, we detect 9 of 23 sources with M$`{}_{\mathrm{B}}{}^{}>27`$, and 7 of 17 sources with M$`{}_{\mathrm{B}}{}^{}<27`$. This is consistent with the lack of a strong correlation between MB and LFIR in the lower redshift, lower optical luminosity sample of Palomar-Green (PG) QSOs presented by Sanders et al. (1989; see also Chini, Kreysa, & Biermann 1989), and is indicative of the large scatter in the dust properties of QSOs in general (see Figure 4 below), and the relatively narrow range in MB in Figure 2.

### 4.3 Comparison to the Radio-to-FIR Correlation for Star Forming Galaxies

Figure 3 shows the relationship between redshift and the 250-to-1.4 GHz spectral index for a star forming galaxy taken from the model of Carilli & Yun (2000) based on 17 low redshift galaxies (roughly equivalent to a modified black body spectrum with $`\mathrm{T}_\mathrm{D}=50\mathrm{K}`$ and dust emissivity index = 1.5). The relationship relies on the tight radio-to-FIR correlation seen for star forming galaxies in the nearby universe (Condon 1992). The dotted curve gives the rms scatter for the 17 galaxies. The solid symbols are the results for the SDSS QSOs detected at 250 GHz, while the open symbols are for the APM QSO sample from Omont et al. (1996) and Yun et al. (2000). The arrows indicate non-detections in the radio, and hence lower limits to the spectral index. On this diagram, a point located below the curve would indicate a source that is radio-loud relative to the standard radio-to-FIR correlation for star forming galaxies, while a point located above the curve would indicate a source which is radio-quiet relative to this relationship (Condon 1992). The lower limits are below the curve, and hence are consistent with a star forming galaxy spectrum, although a factor two better sensitivity at 1.4 GHz is required to provide more stringent constraints in this regard.

The radio-detected QSOs fall below the star forming galaxy curve in Figure 3. Two of the sources (SDSS J1235–0003 and SDSS J1412–0101) are clearly radio loud AGN (see section 4.0). The fainter radio detections may indicate radio emission from the AGN, but at a level below that required for the source to be defined as radio loud according to Miller et al. (1990). An alternative explanation for these sources is a warmer dust component: increasing the dust temperature by a factor two or so could give rise to the observed offset with respect to the model in Figure 3 (Blain 1999).

### 4.4 Spectral Energy Distributions

Figure 4 shows the mean radio-through-optical Spectral Energy Distributions (SEDs) for radio-loud and radio-quiet QSOs derived from observations of a large sample of PG QSOs by Sanders et al. (1989). The solid curve shows the SED for radio loud QSOs, while the dashed curve shows the SED for radio quiet QSOs. The hatched regions indicate the scatter in the measured values for the PG sample. Note that the QSOs in the PG sample are typically an order of magnitude less luminous in the rest-frame UV than the high redshift sources, with most of the PG sources in the range: $`\mathrm{M}_\mathrm{B}=23\mathrm{to}27`$.

The data points in Figure 4 show the results for the mm detections in the SDSS and APM QSO samples, normalized by their blue spectral luminosities. The normalized mm detections all fall within the range defined by the PG sample. The upper limits for the non-detected sources at 250 GHz would fall at the low end of the range defined by the Sanders et al. (1989) SEDs.

The two radio loud SDSS QSOs according to the Miller et al. (1990) definition detected at 250 GHz are clearly evident in the normalized SEDs in Figure 4, while the normalized upper limits for the non-detected sources at 1.4 GHz are consistent with the Sanders et al. (1989) radio quiet QSO SED. The remaining three radio-detected SDSS QSOs in Figure 4 fall in-between the radio quiet and radio loud SEDs, as do the QSOs from the APM sample. The large gap between the radio loud and radio quiet SEDs from Sanders et al. (1989), as reproduced in Figure 4, represents the possible bi-modality of the radio properties of QSOs, as suggested by Stocke et al. (1992). The reality of this bi-modality has been called into question recently by White et al. (2000), and Stocke et al. (1992) show that the bi-modality is less clearly delineated for high blue luminosity QSOs relative to lower luminosity sources. Further analysis of the radio detections within this sample, including studies of radio spectra and spatial structure, will be given in Rupen et al. (2001).

### 4.5 Broad Absorption Lines

Omont et al. (1996a) suggest a possible correlation between the Broad Absorption Line (BAL) phenomenon and thermal dust emission in QSOs. In the APM sample of Omont et al. (1996a), two of six sources detected at 250 GHz have BAL features. The BAL fraction for high $`z`$ QSOs in general is about 10$`\%`$ (Fan et al. 2001a).

Our data on the SDSS sample provide only marginal support for this idea, with 4 of the 16 sources detected at S250 GHz showing evidence for BAL features, as compared to 3 sources with BALs out of the 25 non-detected sources at 250 GHz.

## 5 Discussion

We detect 16 of 41 high redshift QSOs from the SDSS survey with S$`{}_{250}{}^{}1.4`$ mJy. Assuming $`\mathrm{T}_\mathrm{D}=50\mathrm{K}`$, we show that the 250 GHz emission is most likely thermal emission from warm dust on a scale $`>1`$ kpc in the host galaxy of the QSO. The implied dust masses for the sources are in the range of $`1.5\mathrm{to}5.9\times 10^8`$ M. If the dust is heated by star formation, the star formation rates are in the range 500 to 2000 M year-1. Sensitive radio observations show a radio-loud fraction (according to the definition of Miller et al. 1990) for the sample as a whole of 10$`\%`$, consistent with previous observations of high redshift QSOs, with no clear relationship between dust emission and radio-loud QSOs.

The blue-normalized optical-through-radio SEDs for these sources are within the broad range delineated by the PG sample of lower redshift, lower luminosity QSOs (Figure 4). This result suggests that for the sources detected at 250 GHz, we are not seeing a major new broad-band spectral emission ‘feature’ in high redshift QSOs relative to the SEDs of lower redshift QSOs. Conversely, the blue-normalized limits for the non-detections at 250 GHz in the SDSS high $`z`$ QSO sample are at the low end of the distribution defined by low $`z`$ QSOs, leaving open the possibility of a population of high $`z`$ sources which are under-luminous in the rest frame submm relative to the SED of low $`z`$ QSOs. While this result could be used to argue against models involving a major burst of star formation in the host galaxies of high redshift QSOs, the large scatter in the mm SEDs could effectively mask such a phenomenon. Also, Sanders et al. (1989) point out that from cm to submm wavelengths, the radio quiet SED model based on the PG QSO sample is consistent with the observed SEDs for star forming galaxies.

A number of studies have addressed the question of the dust heating mechanism in high $`z`$ AGN by SED modeling, with mixed results. Radiative transfer models by Rowan-Robinson (1999) suggest that the observed SEDs at rest frame wavelengths $`>50\mu `$m are consistent with an extreme starburst model for most high redshift dust emitting QSOs. On the other hand, models by Andreani et al. (1999) and Willott et al. (1999) suggest that the dust emission spectra from 3$`\mu `$m to 30$`\mu `$m can be explained by dust heated by the AGN. The cm-to-mm SEDs may also indicate star formation in a few sources (Yun et al. 2000), however the radio limits in most sources are inadequate to constrain such a model (Figure 3). Mid-infrared spectroscopy using the Space Infrared Telescope Facility of the PAH features and other spectral lines may provide a star formation vs. AGN diagnostic (Genzel et al. 1998).

Spatially resolved observations may be critical in addressing this issue, since a dust emitting region spatially distinct from the optical QSO can be used to argue for at least some dust heating by star formation. Such an argument has been used convincingly in the case of CO and dust emission in the vicinity of a number of high redshift radio galaxies and AGN (Papadopoulos et al. 2000, 2001; Ivison et al. 2000). To date, the only high redshift QSO in which such a spatial segregation has been observed is the $`z=4.7`$ QSO BR1202-0725, in which the cm and mm continuum, and the CO line emission, show two sources of roughly equal strength separated by 4′′, while the optical line and continuum emission is restricted to a single source (Hu, McMahon, & Egami 1996; Omont et al. 1996b; Ohta et al. 1996; Yun et al. 2000; Kohno et al. 2000, in preparation).

The question of the dominant dust heating mechanism, ie. AGN versus star formation, is difficult to answer for low redshift sources (Rigopoulou et al. 1999). Addressing the question at high redshift is just that much more difficult. The sources detected at 250 GHz in Table 1 present an ideal opportunity for addressing this issue in detail. A number of follow-up observations of these sources are required to properly address this question, including: (i) more sensitive radio continuum observations to better characterize the cm-to-mm SEDs, (ii) submm observations to determine the dust temperature, (iii) high resolution (sub)mm observations to determine the spatial distribution of the dust, and (iv) CO observations to better constrain the gas mass and kinematics. Such observations are at the limit of what is possible with current instrumentation, but will become routine with the expanded VLA and the Atacama Large Millimeter Array.

The VLA is a facility of the National Radio Astronomy Observatory (NRAO), which is operated by Associated Universities, Inc. under a cooperative agreement with the National Science Foundation. This work was based on observations carried out with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research made use of the NASA/IPAC Extragalactic Data Base (NED) which is operated by the Jet Propulsion Lab, Caltech, under contract with NASA. CC acknowledges support from the Alexander von Humboldt Society. DPS acknowledges support from National Science Foundation Grant AST99-00703. XF and MAS acknowledge support from the Research Corporation, NSF grants AST96-16901 and AST00-71091 and an Advisory Council Scholarship.


format: text

text:
# ONE- AND TWO-DIMENSIONAL SUBWAVELENGTH SOLITONS IN SATURABLE MEDIA

## I Introduction

The standard approach to description of spatial solitons in nonlinear optical waveguides is based on the use of the nonlinear Schrödinger (NLS) equation, which replaces the Maxwell’s equations (ME) in the paraxial approximation  (see also Refs. ). As is known, this approximation can be insufficient for the description of very narrow, subwavelength (subwavelength) solitons, with a size $`{}_{}{}^{<}{}_{}{}^{}`$ $`\lambda `$, where $`\lambda `$ is the carrier wavelength in vacuum . Terms which are neglected when deriving the NLS equation from ME couple to the propagation constant, affecting dynamics of very narrow solitons with arbitrary polarization . Particular cases of the TM (transverse-magnetic)- and TE (transverse–electric) polarized subwavelength solitons in media with purely cubic (Kerr) nonlinearity were considered in detail in Refs. . It was found that, in the TM case, the size of both bright and dark solitons cannot be essentially smaller than $`\lambda /2`$. In the TE case, the size of the dark soliton is also limited from below, while the bright soliton may formally be arbitrarily narrow; however, narrow TE solitons are subject to a strong instability.

The analysis performed in those works was formal in the sense that subwavelength soliton solutions in the model with the Kerr nonlinearity imply unrealistically large values of the nonlinear correction $`\mathrm{\Delta }n`$ to the refractive index, $`\mathrm{\Delta }n_{}^>\mathrm{\hspace{0.17em}\hspace{0.17em}1}`$. As is known , the quadratic (in the amplitude of the electromagnetic field) correction $`\mathrm{\Delta }n`$, which defines the Kerr nonlinearity, is, as a matter of fact, only the first term of the expansion in powers of the squared field. A detailed analysis shows that results of the formal consideration, based on dropping the higher-order corrections and using the ensuing truncated model in the range where the quadratic correction $`\mathrm{\Delta }n`$is allowed to be large, are ambiguous: they strongly depend on a particular stage of the analysis at which the higher-order terms are omitted. This ambiguity is especially conspicuous in the case of the TM narrow solitons: depending on the choice of the truncation stage, one arrives at a conclusion that the width of the bright TM soliton remains limited from below, or may become arbitrarily small.

In fact, the physical problem of the existence of subwavelength solitons is unresolvable within the framework of the Kerr model, and the only possibility to obtain a definite result is to adopt a more realistic nonlinearity, with saturation of $`\mathrm{\Delta }n`$ in some form for strong fields. While an exact form of the saturation in real materials is not usually known, there are two commonly used models, viz., the saturable nonlinearity proper , and the cubic-quintic (CQ)  one. The former nonlinearity is realistic for alkali metal vapors , and the latter nonlinearity, which combines self-focusing cubic and self-defocusing quintic terms, is well documented to describe nonlinear optical properties of the PTS crystal . Moreover, formation of (2+1)D solitons in this material, consistent with the predictions of the CQ model, was recently observed in the experiment reported in Ref. . We also note that the (2+1)D and (3+1)D versions of the CQ model have recently attracted considerable attention as they give rise to stable or almost stable multidimensional (spatial or spatiotemporal) solitons .

Our objective is to consider subwavelength spatial solitons of both TM and TE types parallel in the CQ and saturable models, in order to arrive at general conclusions concerning their width and other fundamental properties. Moreover, both (1+1)- and (2+1)-dimensional \[(1+1)D and (2+1)D\] cases will be considered, in the latter case the spatial soliton being a cylindrical beam. We will conclude that, in both models and for both dimensions, there always exists a finite minimum value of the soliton’s width, and fundamental properties of the subwavelength solitons are fairly similar in these two models. On the basis of these results, one may infer that, in fact, qualitative results obtained in this work (first of all, the existence of the minimum width) apply to narrow spatial solitons in any non-Kerr model.

The rest of the paper is organized as follows. In section 2, we derive, directly from ME, general equations that describe narrow spatial TM solitons in both (1+1)D and (2+1)D geometries. In section 3, we consider a much simpler TE case. In fact, for this case soliton solutions are well known (in particular, they can be found in an exact form in the (1+1)D geometry). This allows us to rigorously prove that the relative soliton’s width, divided by the carrier wavelength, cannot be smaller than a certain minimum value. In section 4, combining analytical considerations and direct numerical solutions, we find (1+1)D and (2+1)D fundamental TM solitons and arrive at the conclusion that their width is also limited from below.

## II Equations for the transverse-magnetic case

We start the analysis with the most nontrivial case of the TM polarization. In this case, irrespective of the particular nonlinearity, the real (physical) vectorial electric $``$ and magnetic $``$ fields can be taken as

$$_y=H(x)\mathrm{cos}\mathrm{\Phi },_x=_z=0,$$
(1)
$$_x=E(x)\mathrm{cos}\mathrm{\Phi },_z=E_z(x)\mathrm{sin}\mathrm{\Phi },_y=0$$
(2)
in the (1+1)D case, and in the (2+1)D case,

$$_x=(y/r)H(r)\mathrm{cos}\mathrm{\Phi },_y=(x/r)H(r)\mathrm{cos}\mathrm{\Phi },_z=0,$$
(3)
$$_x=(x/r)E(r)\mathrm{cos}\mathrm{\Phi },_y=(y/r)E(r)\mathrm{cos}\mathrm{\Phi },_z=E_z(r)\mathrm{sin}\mathrm{\Phi };$$
(4)
note a phase shift $`\pi /2`$ between the transverse and longitudinal components in these expressions. Here, $`x,y`$ and $`z`$ are the transverse and propagation coordinates, $`r^2x^2+y^2`$, $`t`$ is time, and the common phase of all the fields is $`\mathrm{\Phi }=\beta z\omega t`$, $`\beta `$ and $`\omega `$ being, respectively, the propagation constant and frequency of the carrier wave. In the (2+1)D case, the expressions (3) and (2) correspond to the standard TM$`_{n1}`$mode in a cylindrical waveguide .

The simplest saturable model assumes an isotropic material characterized by the relation

$$𝒟_{\mathrm{sat}}=\epsilon _0\left[1+\frac{\left(\epsilon _2/\epsilon _0\right)^2}{1+\left(\epsilon _4/\epsilon _2\right)^2}\right]$$
(5)
between the electric induction and strength, with positive nonlinear permeabilities $`\epsilon _2`$ and $`\epsilon _4`$. Expansion and truncation of Eq. (5) yields the CQ model in the form

$$𝒟_{\mathrm{CQ}}=\left(\epsilon _0+\epsilon _2^2\epsilon _4^4\right).$$
(6)
A difference between the two models appears in the case when the truncation leading from $`𝒟_{\mathrm{sat}}`$ to $`𝒟_{\mathrm{CQ}}`$ is no longer valid. Below, we present analysis in a detailed form for the CQ model. For the saturable one, it is quite similar, but formulas are more cumbersome. Final results will be displayed for both models together.

Following the usual rotating-wave approximation, the next step is to substitute Eqs. (1) and (2) or (3) and (4) into the relation (6) (or (5), in the case of the saturable nonlinearity), and collect all the contributions to the fundamental harmonics $`\mathrm{sin}\mathrm{\Phi }`$ and $`\mathrm{cos}\mathrm{\Phi }`$, neglecting higher-order harmonics. This yields direct relations between the electric-field induction and strength,

$$𝒟_{x,y}=\epsilon _t_{x,y},𝒟_z=\epsilon _l_z,$$
(7)
where, in the case of the CQ nonlinearity, effective nonlinear susceptibilities are

$$\epsilon _t\epsilon _0+\frac{1}{4}\epsilon _2(3E^2+E_z^2)\frac{1}{8}\epsilon _4(5E^4+2E^2E_z^2+E_z^4),$$
(8)
$$\epsilon _l\epsilon _0+\frac{1}{4}\epsilon _2(E^2+3E_z^2)\frac{1}{8}\epsilon _4(E^4+2E^2E_z^2+5E_z^4).$$
(9)
Insertion of the above expressions (1) and (2) or (3) for the magnetic and electric field and (4) into the Maxwell’s vectorial equation for the electric field,

$$\times =\frac{1}{c}\frac{}{t},$$

where $`c`$ is the light velocity in vacuum, we arrive at a scalar ODE

$$E_z^{}+\beta E=(\omega /c)H,$$
(10)
where the prime stands for $`d/dx`$ or $`d/dr`$ in the (1+1)D and (2+1)D cases, respectively. Further, we substitute the (1+1)D or (2+1)D expressions (1) and (2) or (3) and (4), in combination with the relations (7) and (8), (9) which define the electric induction, into the Maxwell’s vectorial equation for the magnetic field,

$$\times =\frac{1}{c}\frac{𝒟}{t}.$$

This yields two more equations, one of which is again an ODE, while the other one is just an algebraic relation,

$`\epsilon _tr^{1D}(r^{D1}H)^{}`$ $`=`$ $`(\omega /c)E_z\epsilon _l,`$ (11)
$`\beta H`$ $`=`$ $`(\omega /c)E,`$ (12)
where $`D=1,2`$ is the transverse dimension.

Eliminating the magnetic field $`H`$ from Eqs. (10) and (11) by means of Eq. (12), we then obtain a system of two equations for the fields $`E`$ and $`E_z`$,

$`E_z^{}`$ $`=`$ $`E(1P),`$ (13)
$`r^{1D}\left(r^{D1}Ev_1\right)^{}`$ $`=`$ $`E_zv_2.`$ (14)
Here, the variables have been rescaled as $`(x,y,r)\beta \sqrt{\gamma }(x,y,r),[E,E_z][E,(E_z\beta /\sqrt{\gamma })](\omega /\beta \sqrt{\gamma })\sqrt{3\epsilon _2}/2`$, with $`\gamma 1\left(\beta _0/\beta \right)^2`$, where $`\beta _0\sqrt{\epsilon _0}\omega /c`$ is the propagation constant in the linear regime, and

$`P`$ $``$ $`E^2+{\displaystyle \frac{\gamma }{3}}E_z^2+\sigma \left(E^4+{\displaystyle \frac{2\gamma }{5}}E^2E_z^2+{\displaystyle \frac{\gamma ^2}{5}}E_z^4\right),`$ (15)
$`v_1`$ $``$ $`1+{\displaystyle \frac{\gamma }{1\gamma }}P,`$ (16)
$`v_2`$ $``$ $`1+{\displaystyle \frac{\gamma }{1\gamma }}\left[{\displaystyle \frac{1}{3}}E^2+\gamma E_z^2+\sigma \left({\displaystyle \frac{1}{5}}E^4+{\displaystyle \frac{2\gamma }{5}}E^2E_z^2+\gamma ^2E_z^4\right)\right],`$ (17)
$`\sigma `$ $``$ $`\gamma (1\gamma )^1s,`$ (18)
$`s`$ $``$ $`\left(10/9\right)\epsilon _0\epsilon _4\epsilon _2^2.`$ (19)
These equations will be used below to study TE and TM solitons in the CQ model, which will be paralleled by the same analysis for the model with the saturable nonlinearity. Before proceeding to that, it may be relevant to revisit the case of the TM soliton in the Kerr (rather than CQ) (1+1)D medium. In that case, it is necessary to bear in mind the pure cubic (Kerr) nonlinearity applies as long as the nonlinear correction to the refractive index is much smaller than the linear index. In the present notation, this condition can be shown to amount to inequalities

$$E^2(x)1\text{and}\left(E^{}\right)^21.$$
(20)
Taking these into regard and performing the corresponding expansions in the above equations (13) through (17), one can eliminate the longitudinal electric field $`E_z`$ and derive an eventual equation for the transverse field,

$$E^{\prime \prime }\gamma E+\left[E^2+\frac{1}{3}\left(E^{}\right)^2\right]E=0,$$
(21)
where $`\gamma `$ was defined by Eq. (LABEL:gamma). In comparison with the traditional paraxial (NLS) approximation, the only new term in Eq. (21) is $`\left(1/3\right)\left(E^{}\right)^2E`$, which, as a matter of fact, is a contribution from the longitudinal component $`E_z`$ in the present TM case.

Obviously, Eq. (21) can give rise to bright solitons only in the case $`\gamma >0`$, which will be assumed to hold hereafter. Note that the classical broad solitons correspond to the case $`\gamma 1`$. In that case, the broad soliton with the first correction produced by the extra term can be easily found:

$$E_{\mathrm{sol}}(x)=\sqrt{2\gamma }\left[\left(1+\frac{\gamma }{18}\right)\mathrm{sech}\left(\sqrt{\gamma }\xi \right)\frac{\gamma }{9}\mathrm{sech}^3\left(\sqrt{\gamma }\xi \right)\right].$$
(22)
Numerical solution of Eq. (21) shows that it does give rise to solitary-wave solutions with an arbitrary small width; however, when the soliton becomes too narrow, the solution violates the applicability conditions (20), which makes it necessary to modify the nonlinearity, i.e., to consider the model with the CQ or saturable nonlinearity, which will be done below.

## III The transverse-electric case

The above equations, derived for the CQ nonlinearity, also contain an essentially simpler case of the TE polarization, which can be obtained setting formally $`\gamma =0`$ but keeping $`\sigma 0`$. In this case, $`v_1=v_2=1`$,$`P=E^2+\sigma E^4`$, and in the (1+1)D geometry an eventual equation for $`E(x)`$ takes a well-known form

$$E^{\prime \prime }=EE^3\sigma E^5,$$
(23)
a commonly known exact soliton solution to which can be written as

$$4/E^2(x)=1+\sqrt{1+16\sigma /3}\mathrm{cosh}(2x).$$
(24)
It exists provided that $`\sigma >3/16`$, or, with regard to Eq. (19),

$$(\beta /\beta _0)^2<(\beta /\beta _0)_{\mathrm{cr}}^21+3/(16s).$$
(25)
As $`(\beta /\beta _0)^2`$ is approaching the value $`(\beta /\beta _0)_{\mathrm{cr}}^2`$, the average width of the soliton (24), defined as per Eq. (30), diverges as $`\left|\mathrm{ln}\left[1+3/(16s)(\beta /\beta _0)^2\right]\right|`$, while the soliton’s amplitude remains finite, $`E_{\mathrm{max}}=2`$.

The main characteristic of the spatial soliton is its width. For the (1+1)D soliton (24), the FWHM  width $`\mathrm{\Delta }`$ can be easily found,

$$\mathrm{cosh}\mathrm{\Delta }=\left(1+16\sigma /3\right)^{1/2}+2.$$

In unnormalized units, the FWHM width is

$$W=(4\lambda /3\pi n_0)\sqrt{s}F(16\sigma /3),$$

where

$$F(y)\left|y\right|^{1/2}\mathrm{ln}\left[(1+y)^{1/2}+2+\sqrt{3+4(1+y)^{1/2}+(1+y)^1}\right].$$
(26)
The function $`F`$, and hence the width of the TE soliton, attain a minimum value at $`\sigma =0.1493`$. The ratio of the corresponding minimum width of the TE soliton, normalized to the carrier wavelength, is

$$W_{\mathrm{min}}^{(\mathrm{TE})}/\lambda =0.825\sqrt{\epsilon _4}/\epsilon _2.$$
(27)
In the next section, we will also use a different (average) definition of the width, based on Eq. (30). It is easy to find that the minimum value of the average width (30) differs from the minimum FWHM value (27) by less than $`8\%`$.

Consideration of TE solitons in the (2+1)D geometry leads to an equation

$$\frac{d^2E}{d^2r}+\frac{1}{r}\frac{dE}{dr}\frac{1}{r^2}E=EE^3\sigma E^5$$
(28)
($`r`$ again being the radial variable), cf. Eq. (23). In fact, Eq. (28) is exactly the same equation which describes known solitons with an internal spin (vorticity) $`s=1`$ in the CQ model in the (2+1)D geometry . The origin of an effective “vorticity” term, $`E/r^2`$, in Eq. (28) is the structure of the (2+1)D ansätze (3) and (4) for the vectorial fields.

Equation (28) has been studied in detail numerically elsewhere . Although the minimum width of a family of solitons generated by this equation was not specially considered, a straightforward consequence of results presented in Refs.  is that the minimum width is always finite.

Thus, it is possible to rigorously prove the existence of a finite lower bound for the width of spatial TE solitons in the (1+1)D and (2+1)D models with the CQ nonlinearity. Quite similarly, the same results can be obtained for TE solitons in the model with the saturable nonlinearity (we do not display technical details here, as they do not contain anything essentially novel).

## IV Transverse-magnetic narrow spatial solitons

The most interesting case is that with the TM polarization, as it does not reduce to a single-component equation. As it was mentioned above, the necessary condition for the existence of localized solutions to the corresponding equations (13) and (14), supplemented by Eqs. (15) through (19), is $`\beta _0^2<\beta ^2`$. Below, we will use a relative propagation constant, $`\beta /\beta _0(1\gamma )^{1/2}`$, as a measure of the departure from the paraxial approximation: in the limit of a very broad soliton, one has $`\beta /\beta _01`$, while in the opposite limit of an infinitely narrow soliton (if any), $`\beta /\beta _0\mathrm{}`$.

A straightforward analysis of Eqs. (13) and (14) makes it possible to prove that, in the (1+1)D case for both CQ and saturable nonlinearities, localized solutions with no zeros of $`E(x)`$ and $`E_z(x)`$ at $`x0`$ do not exist if $`E(0)=0`$. Further analysis, details of which are not displayed here, shows that, in the (1+1)D geometry, a localized fundamental-soliton solution has an even monotonically decreasing transverse component $`E(|x|)`$, while the longitudinal one $`E_z(x)`$ is odd, with $`E_z(0)=0`$ and a single extremum at finite $`\left|x\right|`$, see Fig. 1a (in this and next figures, the results are simultaneously shown for the CQ and saturable models). In the (2+1)D case, fundamental solitons feature a different structure: they have $`E(0)=0`$, and both $`E(r)`$ and $`E_z(r)`$ may have zeros at finite $`r`$, see Fig. 1b.

Proceeding to the width of the TM soliton, we use the electromagnetic energy density,

$$𝒬=_x𝒟_x+_y𝒟_y+_z𝒟_z+_x_x+_y_y,$$
(29)
to define the soliton’s average half-width as

$$W=_0^{\mathrm{}}𝒬r^D𝑑r/_0^{\mathrm{}}𝒬r^{D1}𝑑r.$$
(30)
The FWHM definition of $`W`$, which was used in the previous section for the TE soliton, is ambiguous for TM solitons, as they have three different components. The results for the average width, which were obtained from numerical solutions of Eqs. (13) and (14), are summarized in Fig. 2, showing $`W`$ vs. the relative propagation constant at different values of the material constant $`s`$ defined by Eq. (19). Basic features of the dependence are the same as those following from the above analytical expressions for the TE soliton in the (1+1)D case: there is a finite minimum of the soliton’s width, and the solitons do not exist beyond a critical value of $`\beta /\beta _0`$, at which the soliton’s width diverges while its amplitude remains finite, cf. Eq. (25). Naturally, the minimum value of the width depends on the dimension of space, type of the nonlinearity, and, for a fixed nonlinearity, on values of material constants. As well as in the case of the TE solitons, the minimum width tends to zero and the maximum value of $`\beta /\beta _0`$ diverges as $`\epsilon _40`$.

Note that the soliton shapes shown in Fig. 1 pertain to the fixed value of the relative propagation constant $`\beta /\beta _0=1.05`$, which is rather close to the paraxial limit. Comparison with other numerical results demonstrates that, for both CQ and saturable models, the shapes of the (1+1)D and (2+1)D TM solitons remain fairly similar to those shown in Fig. 1 with the increase of $`\beta /\beta _0`$ up to the point where the minimum width is attained. With the subsequent increase of $`\beta /\beta _0`$ up to the value at which the solitons cease to exist, the shape changes much more.

The above considerations did not tackle the stability problem. While a consistent stability analysis requires very tedious simulations of the Maxwell’s equations with the CQ or saturable nonlinearity, which is beyond the scope of this work, we tried to test the stability by means of the simple Vakhitov-Kolokolov (VK) criterion  (although its applicability to the present model is not obvious): the propagation constant must have a positive slope as a function of the beam’s power (norm of the solution),

$$N=2^{2D}\left(2\pi \right)^{D1}_0^{\mathrm{}}\left(_x^2+_y^2+_z^2\right)r^{D1}𝑑r$$

(cf. Eq. (30)), which is a dynamical invariant of the model. Our numerical calculations have shown that all the TM solitons, in both the (1+1)D and (2+1)D geometries, satisfy the VK criterion.

## V Conclusion

In this work, we have studied very narrow (1+1)D and (2+1)D spatial TE and TM solitons supported by the quintic-cubic and saturable nonlinearities, starting the full Maxwell’s equations. For TE solitons, we have obtained a single equation, which is equivalent to the stationary version of the nonlinear Schrödinger equation with the corresponding nonlinearity, while for TM solitons we end up with a system of equations for the transverse and longitudinal components of the electric field, rather than with a single one. It is also noteworthy that TE solitons in the (2+1)D case are described by the same equation as vortex solitons with spin $`s=1`$ in the (2+1)D nonlinear Schrödinger equation.

Detailed analysis was displayed for the cubic-quintic model, while final results were given for the saturable model as well. Full qualitative similarity of the results obtained for both models (although exact scales may be quite different, see caption to Fig. 2) makes it very plausible that the situation is basically the same for any physically relevant non-Kerr medium. The most important conclusions are that, in all the cases, there is a finite minimum of the soliton’s width, and the solitons cease to exist at a critical value of the propagation constant, at which their width diverges.

## Figure captions

Fig. 1. Typical examples of the fundamental spatial TM soliton in the (1+1)D (a) and (2+1)D (b) geometries. The relative propagation constant is $`\beta /\beta _0=1.05`$, and the material constant (19) is $`s=1`$. In this and next figures, the solid and dotted curves pertain, respectively, to the cubic-quintic and saturable models.

Fig. 2. The width of the (1+1)D (a) and (2+1)D (b) fundamental spatial TM solitons vs. the relative propagation constant $`\beta /\beta _0`$ at various fixed values of the material constant $`s`$. For the saturable model, the deviation of the relative propagation constant from $`1`$ is five times that shown on the horizontal axis for the cubic-quintic model.


text:
# A Catalog of OB Associations in the spiral galaxy NGC 300 Based in part on observations obtained with the 1.3 m Warsaw telescope at the Las Campanas Observatory of the Carnegie Institution of Washington, and on observations obtained at the European Southern Observatory, La Silla, Chile

## 1 Introduction

OB associations constitute physical groups of young, massive, gravitationally unbound stars, which have formed in the same molecular cloud. Observations of these short lived objects provide much valuable information about regions of recent and/or current star formation, parent galaxies and star formation processes in general. In particular, a comparison of the properties of systems of associations in different galaxies may help us to better understand the influence of the environment on mechanisms of star formation.

A lot of work has been dedicated in the past to the detection and analysis of properties of associations in nearby galaxies. However, as thoroughly discussed by Hodge (1986), the identification of OB associations is a very difficult task. Indeed, due to very low surface density of stars it is not possible to detect them just based on the density contrast with the surrounding field. For a long time, associations were identified from visual searches of photographic plates. Subjective identification together with nonuniform observational material and poor resolution often led to disagreement between different authors. As a result of these limitations it was not possible to perform unambiguous comparisons of the properties of stellar associations found in different galaxies. To overcome these difficulties new, automatic and objective techniques of detecting OB associations were developed (Battinelli 1991, Wilson 1991). All one needs in order to use these methods are the positions of the sample of OB stars. A successful application of these techniques to search for OB associations in several nearby galaxies (Battinelli 1991, Magnier et al. 1993, Wilson 1992) reflects that they are very useful tools for building uniform catalogues of OB associations in galaxies well suited for comparison purposes.

In this contribution we present the results of a search for OB associations in the nearly face-on spiral galaxy NGC 300. This galaxy, which is a member of the Sculptor Group, presents massive recent star formation, yet no systematic work to discover OB associations in this stellar system has been undertaken. To our knowledge only five stellar groups assumed to be associations were observed by Breysacher et al. (1997) during their search for WR stars in stellar associations. A systematic search with the aim to uncover the complete system of OB associations in this important nearby galaxy is therefore long overdue.

Our paper is organized as follows. In Sect. 2 we briefly describe the observations and data reduction techniques applied in this work. The algorithm used for searching for potential associations is described in Sect. 3. Results are presented in Sect. 4. Sect. 5 contains some final remarks on our findings.

## 2 Observations

The observational data have been collected in the course of a long term project whose principal goal is to discover a large number of Cepheid variables in NGC 300 in order to calibrate the effect of metallicity on Cepheid luminosities. Up to now several hundred images in $`BVRI`$ bands were collected with the Wide Field Camera attached to the ESO/MPI 2.2 m telescope on La Silla. This camera has a pixel scale of 0.238 arcsec/pixel, and its large field of view of 33 x 34 arcmin allows to cover the whole galaxy with one pointing. For the search of OB associations we used 10 $`B`$ and 10 $`V`$ images obtained on two photometric nights, at low airmass and under 0.8 arcsec seeing conditions. Exposure time was set to 360 sec for all images in each filter. As can be seen from Fig. 1 the limiting magnitude was about 22.5 mag in $`B`$ and about 23 mag in $`V`$. After debiasing and flatfielding the images with the IRAF 111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF. package, each of the eight 2k x 4k chips was divided into 16 slightly overlapping sub-frames. Profile photometry with $`DAOPHOT`$ and $`ALLSTAR`$ programs was performed for all stars detected in all subframes. Then the photometry on all subframes was tied together using stars located in the overlapping regions.

In order to calibrate our data we performed follow–up observations with the Warsaw 1.3 m, photometric telescope located at Las Campanas Observatory, Chile. The telescope was equipped with a 2048 x 2048 CCD camera. The scale is 0.41 arcsec/pixel, which corresponds to a field of view of about 14.5 x 14.5 arcmin. More details about the instrumental system can be found in Udalski et al. (1997). Observations were conducted through $`BVRI`$ filters, during 2 photometric nights. We monitored four fields, covering about 25 x 25 arcmin of the central part of NGC 300. During each night a large number of Landolt standard stars (Landolt 1992) covering a wide range of colors and airmasses, were also observed in order to tie our new secondary standards in NGC 300 to the Landolt system. Preliminary reductions (i.e. debiasing and flatfielding) were performed with the IRAF package while the PSF photometry was done with the $`DAOPHOT`$/$`ALLSTAR`$ package.

The accuracy of the transformation of the instrumental magnitudes and colors to the standard system is about 0.04 mag for the $`B`$ and $`V`$ bands, and slightly better for the $`(BV)`$ color. As a next step we carefully selected objects with reliable photometry and used them to calibrate the photometry from the WFI camera. In this process, we eliminated background galaxies, severely blended and saturated stars by visual examination.

Having established a sequence of about 300 photometric secondary standard stars in the field of NGC 300 (results to be published in a forthcoming paper), we transformed the WFI photometry from each chip using a least squares method. Usually more than 20 stars common between a given chip and our catalog of secondary standards were identified, and the accuracy of internal transformations was better than 0.01 mag (rms).

We are continuing the observations of secondary standards stars in NGC 300, and the accuracy of their brightness calibration will soon be significantly improved. A more detailed description of the project, observations, and reduction procedures will be presented in a followup paper.

## 3 Search for Associations

### 3.1 Identification Method

The Path Linkage Criterion (PLC) technique (Battinelli 1991) was adopted in our search for OB associations in NGC 300. For a detailed description of this technique the reader is referred to Battinelli (1991). In short, this method assumes that any given two stars in an ensemble of OB stars belong to the same association if and only if it is possible to connect these two stars, by successively linking OB stars located between them, separated from each other by no more than a certain fixed distance parameter, or search radius, called $`d_\mathrm{s}`$. The advantage of the PLC method is that $`d_\mathrm{s}`$ can be unambiguously derived based on the given catalog of early type stars using the function $`f_\mathrm{p}(d)`$. This function describes the number of groups containing at least p stars, for any given value of the distance parameter d. The behaviour of this function is as follows: a steep rise until a maximum value is reached, followed by a gradual decrease towards a limiting value of 1, which corresponds to the situation that the value of d is so large that all stars can be connected using it (i.e. all stars are assigned to just one group). The optimum value of the distance parameter $`d_s`$ is defined as the value of d corresponding to the maximum of the function $`f_\mathrm{p}(d)`$.

### 3.2 Application of the PLC technique to NGC 300

Before we can start to search for OB associations we need to select the catalog of blue stars. In the case of the photometric data it means that we need to adopt brightness and color cutoffs. The sample of blue stars was selected from the photometric catalog of stars in NGC 300 using the following criteria: $`V<`$ 22 mag and -0.6 $`<BV<`$ 0.4 mag. We note that NGC 300 has a very low foreground reddening, in the order of $`E(BV)=0.02`$ mag, so the color cutoff for the unreddened $`BV`$ color is also approximately 0.4 mag. Altogether 4016 stars satisfied our selection criteria and entered our list of blue stars. Unfortunately the distance to NGC 300 is rather uncertain at the present time; distance moduli obtained by different authors differ quite significantly. Based on an analysis of a variety of distance indicators, van den Bergh (1992) found a distance modulus of 26.0 mag, while Freedman et al. (1992) obtained a value of 26.66 mag, using $`BVRI`$ CCD observations of Cepheids in this galaxy. If we adopt the latter value, our brightness cutoff corresponds to an absolute magnitude of $`M_v=4.7`$ mag. In Fig. 1 we present the $`BV`$ color - magnitude diagram (CMD) from about 14000 stars observed in the field of NGC 300, with an indication of the location of the selected blue stars.

The next step was to derive the value of $`d_s`$ based on our catalog of blue stars. To do this we apply the PLC technique with several values for the distance parameter and several values for the minimum number of stars to define an association, p. The results are presented in Fig. 2. It can be seen that the position of the maximum of the function $`f_p(d)`$ depends slightly on p and is located in the interval from about 35 to 40 microradians. In order to check on the influence of our choice for the value of $`d_s`$ on the number of associations we find, we used three search distances, namely $`d_\mathrm{s}`$= 35, 37.5 and 40 microradians which correspond to linear distances of 72, 78 and 83 pc, respectively ( assuming a distance modulus of 26.66 mag), and then compared the results. We found that all but three objects were common in all trials. This is because in most cases we can just add a few (10 - 15 %) new stars to already detected associations if we use a slightly larger value of $`d_\mathrm{s}`$. As a conclusion, the number of detected OB associations is very insensitive to the adopted value of the search radius, in the range suggested by Fig. 2. We adopted $`d_\mathrm{s}=37.5`$ microradians as the most appropriate search radius.

The last parameter that needs to be specified is the minimum number of stars p any potential association should possess. It is clear that too small a number of stars could produce many spurious detections. On the other hand, if we demand an unreasonably large number of stars many small association may be lost. We are then at the tradeoff between completeness and reliability of our sample of associations. In order to check how many detections are caused by chance coincidence, as a function of p, a statistical test was performed. One hundred random distributions of a number of stars equal to the number of our sample of blue stars, and distributed over the same area, were created and the PLC technique with our adopted distance parameter value was applied to search for potential groups. This experiment was repeated for different minimum numbers of stars. The plot of the resulting fraction of the spurious detections to the total number of detected associations, as a function of the minimum number of stars, is presented in Fig. 3. We can see that the spurious detections constitute less than 10 % if we adopt 6 stars as a minimum population of OB stars of potential associations. This number seems to be a good compromise and we therefore adopted it.

## 4 Results

### 4.1 The Catalog

The application of the PLC technique, with a threshold of number of stars equal to 6 and the distance parameter value of 37.5 microradians as discussed above, resulted in the detection of 117 groups of blue stars. Table 1 contains their description. Columns 1, 2 and 3 give the designation, and the equatorial coordinates of the detected associations. The number of blue stars belonging to the group, its size in both microradians and parsecs (distance modulus of 26.66 assumed), cross-identification with the catalog of HII regions in NGC 300 (Deharveng et al. 1988) and objects studied by Breysacher et al. (1997), (arabic and roman numbers respectively) are given in the following columns.

The coordinates of the OB associations were determined as the mean coordinates of their separate member stars. Following Battinelli (1991), the real diameter of associations is defined as $`0.5(\mathrm{\Phi }_{R.A.}cos(\mathrm{\Phi }_{Dec})+\mathrm{\Phi }_{Dec})d_{\mathrm{NGC}300}`$, where $`\mathrm{\Phi }_{R.A.}`$ and $`\mathrm{\Phi }_{Dec}`$ correspond to the diameter in right ascention and declination respectively and $`d_{NGC300}`$ stands for distance to the galaxy. We assumed $`d_{\mathrm{NGC}300}=2.1Mpc`$ (distance modulus of 26.66 mag). Fig. 4 displays the spatial distribution of the detected associations over the disc of NGC 300. Maps of six typical associations are presented in the appendix.

We would like to stress that our spatial resolution allows detection of associations as small as about 15 pc. We therefore expect that our catalog is near complete even down to potentially very compact associations.

### 4.2 Size distribution

In Fig. 5 we show the size distribution of the detected associations. There is a sharp peak located at about 60 microradians. This value depends somewhat on the assumed value of $`d_\mathrm{s}`$. Using $`d_\mathrm{s}=35`$ and $`d_\mathrm{s}=40`$ microradians, respectively, sizes are about 15 % different (smaller and bigger, respectively), but the shape of the distribution remains the same. Table 2 gives numerical values for the location of the peak of the size distribution, for a range of distance parameters, and for the two different distance moduli to NGC 300 discussed before.

Until a more precise distance to NGC 300 is known (which is one of the main purposes of our NGC 300 project), it is not possible to carry out an unambiguous comparison of our size results with those obtained for OB associations in other galaxies.

### 4.3 Correlation with HII regions

In most cases the ionization of hydrogen is caused by UV photons emitted from nearby hot, massive stars, so one can suspect a correlation of the spatial distributions of OB associations and H II regions. To check this we compared the spatial distribution of OB associations obtained in this paper with the spatial distribution of H II regions in NGC 300 cataloged by Deharveng et al. (1988). The result is that about 60 % of the detected OB associations do clearly overlap with H II regions (see Tables 1 and 3 for detailed cross-identification) , which indicates that there is indeed the expected strong correlation. Thus both OB associations and H II regions appear to be about equally well suited tracers of recent or ongoing star formation.

### 4.4 Stellar Complexes

We can see from Fig. 5 that besides the stellar groups with sizes typical for OB associations (of about 100 pc) there are also some groups having much bigger dimensions. Regarding their size they resemble Galactic star complexes, as described by Efremov et al. (1987).

By visual inspection of the objects identified as having sizes larger than 300 pc (potential star complexes) we found that most of them appear to consist of a number of smaller, well defined associations with sizes close to 100 pc (see the map of one of the stellar complexes presented in the appendix). These were not detected as individual stellar groupings by the PLC technique due to the high density of blue stars with respect to the surrounding regions in these parts of the galaxy.

To check on this in an objective way we selected blue stars from each of these regions and applied the PLC method for a second time. In 14 out of 18 cases we detected that the apparently large single associations do indeed consist of several smaller associations. Their description is given in Table 3. The names of new objects were designated using the name of the corresponding associations detected during the first application of the PLC technique, followed by consecutive roman letters. (e.g. AS\_102a, AS\_102b, etc.) Their parameters were derived in the same manner as described in Sect. 4.1. As a conclusion, very large single OB associations (like AS\_008, AS\_031, AS\_084 and AS\_096) seem to be very rare.

It is interesting to note that 14 out of 18 potential star complexes are located in the western part of NGC 300, while the distribution of detected OB associations seems to be fairly uniform (see Fig. 4). This fact suggests that for some unknown reasons the physical conditions for the formation of large objects have been more favorable in this part of the galaxy.

## 5 Summary

We have presented the results of a search for OB associations in the spiral galaxy NGC 300, using photometric data obtained within a broader project aiming at a precise distance determination to this galaxy from its Cepheid variables. Magnitudes and colors were obtained for about 14000 stars in a 0.5 x 0.5 square degree field centered on the nucleus of NGC 300. From this sample, we selected blue stars using the constraints $`V<22`$ mag and $`0.6<BV<0.4`$ mag. Application of the PLC method resulted in the detection of 117 groups of blue stars. Statistical considerations show that the expected contamination of this sample by spurious detections is less than 10 %. This suggests that our sample constitutes a homogeneous group which is well suited for comparison with systems of associations found in other galaxies. We detected the expected strong correlation between the spatial distributions of associations and H II regions in NGC 300, a fact which lends additional strong support to the relability of our OB association searching procedure.

Most of the detected OB associations have sizes comparable to typical associations from the Galaxy and other Local Group galaxies. However, a few objects have sizes up to five times larger, corresponding rather to stellar complexes than associations. The PLC technique was applied once more to stars from the regions of each of the potential stellar complexes. This revealed that most of these large complexes are composed of well-defined smaller subgroups having sizes corresponding to typical OB associations. This suggests a hierarchical arrangement of OB stars groupings, already found previously in other galaxies (Efremov et al. 1987, Battinelli et al. 1996).

Once we will have an accurate distance determination to NGC 300 from our project, we will use the results of this paper to compare the system of OB associations of NGC 300 to those of other galaxies for which such studies have been performed.

###### Acknowledgements.

We would like to thank the OGLE team, especially Drs: Andrzej Udalski, Michal Szymański and Marcin Kubiak for their kind help in observations and data reduction. We are very grateful to the European Southern Observatory for large amounts of observing time with the Wide Field Camera on the 2.2 m telescope. It is a pleasure to thank the 2p2 team for expert support and for doing several of the observations in service mode. It is also a pleasure to express our gratitude to Dr. Paolo Battinelli for providing us with his computer programs. WG acknowledges financial support through Fondecyt grant No. 1000330. We would also like to thank Dr. Azzopardi for valuable comments.
Figure Captions

Fig. 1. CMD diagram for about 14000 stars observed in the field of NGC 300. The rectangular region contains blue stars used for searching for OB associations.

Fig. 2. Behaviour of the function $`f_\mathrm{p}(d)`$ for a set of p (minimum number of stars in an association) values.

Fig. 3. The percantage of expected spurious detections of associations, versus p, the minimum number of stars in an association.

Fig. 4. Spatial distribution of detected OB associations (filled circles) over the disc of NGC 300. Dots are individual blue stars. The center of the figure corresponds to $`\alpha _{2000}`$ = 0h54m53.s4, $`\delta _{2000}`$ = -374103′′.

Fig. 5. Size distribution of associations in NGC 300.


text:
# The Peculiar Type Ia Supernova 1999by: Spectroscopy at Early Epochs1footnote 11footnote 1Based on observations obtained at David Dunlap Observatory, University of Toronto, Canada

## 1 Introduction

In the last decade it became evident that some of the observed Type Ia supernovae are intrinsically sub-luminous compared to the majority of so-called “normal” SNe Ia. SN 1991bg (Filippenko et al., 1992) is now considered as their prototype, and the sample of well-observed events includes SN 1992K (Hamuy et al., 1994), SN 1997cn (Turatto et al., 1998) and SN 1998de (Modjaz et al., 2000). Beside being about 1 magnitude fainter than normal SNe Ia at maximum, sub-luminous SNe Ia also exhibit peculiar colors and spectra. For example, the intrinsic (unreddened) $`(BV)_0`$ color around maximum is 0.7 $``$ 0.8 magnitude redder than normal SNe Ia, reaching $`(BV)_0=1.5`$ mag at 10 days after maximum light. The spectroscopic peculiarity is represented by $`i)`$ the appearance of strong Ti II band in the blue region and $`ii)`$ the increase of the ratio of the depths of Si II absorption troughs at 5800 Å and 6150 Å ($`R`$(SiII)) (e.g. Nugent et al., 1995). Moreover, sub-luminous SNe Ia have steeply declining light curves. Their decline rate parameter is $`\mathrm{\Delta }m_{15}(B)>1.9`$ instead of 1.0 $``$ 1.5 as for normal SNe Ia (e.g. Phillips et al., 1999). Because both $`R`$(SiII) and $`\mathrm{\Delta }m_{15}(B)`$ correlate with the peak luminosity, the observational data indicate that sub-luminous SNe Ia have less energetic explosions and less ejected Ni masses than normal SNe Ia, although the correlation between the different quantities cannot be described with a single parameter (Hamuy et al., 1996, Hatano et al., 2000).

In this paper we present medium resolution spectra of the sub-luminous Ia-type supernova SN 1999by occured in NGC 2841. Chronicles of the discovery and the earliest observations can be found in IAU Circulars No. 7156$``$7159. Gerardy & Fesen (1999) pointed out very soon that the spectrum of SN 1999by is of Type Ia, while Garnavich et al. (1999) noted that SN 1999by belongs to the group of sub-luminous SNe Ia based on the $`R`$(SiII) line depth ratio. Early-phase spectra of such SNe are available only for SN 1991bg and SN 1997cn. Thus, the spectroscopic observations made around maximum light may supply interesting details about peculiar SNe Ia. Very recently Howell et al. (2001) presented spectropolarimetric evidence for the strong ($`20\%`$) asphericity of the ejecta of SN 1999by, suggesting a connection between the observed asymmetry and the mechanism producing sub-luminous SNe Ia.

Our new observations are described in Sect.2, while their interpretation is given in Sect.3. Sect.4 summarizes the results of this paper.

## 2 Spectroscopic observations

SN 1999by was observed with the 1.88 m telescope and Cassegrain spectrograph of the David Dunlap Observatory on three nights between May 10 and 21, 1999. The 150A grating in 2nd order with an order separation filter in the stellar beam was applied, resulting in a reciprocal dispersion of 1.22 Å/pixel on the CCD detector. The recorded spectra cover about 1250 Å centered on 6000 Å.

Reduction was performed by standard routines in $`IRAF`$333$`IRAF`$ is distributed by NOAO which is operated by the Association of Universities for Research in Astronomy (AURA) Inc. under cooperative agreement with the National Science Foundation. The wavelength scale was calibrated with FeAr spectral lamp exposures. The background light due to the host galaxy was approximated with a quadratic polynomial perpendicular to the dispersion axis and subtracted from the SN aperture at each wavelength.

Because of the rarity of photometric nights at the DDO site (in the proximity of Toronto), precise flux calibration of the spectra was not possible. However, an approximate correction for the CCD spectral response has been determined by obtaining the spectra of three nearby field stars (#2, #3, #4 of Skiff, 1999; see Fig.1 of Tóth & Szabó, 2000). Their $`(BV)`$ indices (0.59, 0.75, 0.60 mag, respectively) corresponds to $`T_{eff}55006000`$ K according to the tables of Schmidt-Kaler (e.g. Carroll & Ostlie, 1996). The measured continua of these stars could be approximated by a smooth polynomial. The real continuum fluxes (assuming $`\mathrm{log}g=4`$) were estimated from model spectra of Buser & Kurucz (1992). Then, the wavelength-dependent sensitivity of the used CCD spectrograph was calculated by dividing the measured continua with the chosen model continuum flux values. The resulting function turned out to be smooth and slightly decreasing to the red. We have looked for any sudden drop in the sensitivity function that might distort the depths of spectral lines and found nothing.

The SN spectra were then corrected by dividing the measured spectra with the approximate sensitivity function described above. Because the used comparison stars were not spectroscopic standards, the corrected “continuum” slopes of the SN spectra may still be slightly distorted. Nevertheless, the SN spectra calibrated in such a way are acceptable for a comparative analysis described in Sect.3.1.

The reduced spectra are collected and plotted in Fig.1. The epochs of the observations are indicated on the right side of each spectrum. The phases relative to B-maximum were determined by adopting the result of Bonanos, et al. (1999) that the maximum in $`B`$ occured on JD 2451308.5 (May 9/10, 1999), which was also confirmed by photometry of Tóth & Szabó (2000) (note that the JD of maximum light was misprinted in both papers).

Although the spectral coverage is rather small, the appearance of the strong Si II line around 6150 Å unambigously identifies a Type Ia SN. The other strong feature, the absorption trough at 5800 Å is also attributed to Si II. The relative strength of this line to that of the 6150 Å line $`R`$(SiII) (Nugent et al., 1995) suggests that SN 1999by is a sub-luminous event, as mentioned in Sect.1. This is illustrated in Fig.2, where the maximum-light spectrum of SN 1999by is plotted together with a pre-maximum spectrum of SN 1998aq (Vinkó et al., 1999) obtained with the same instrument and setup (unfortunately, no spectrum of SN 1998aq closer to its maximum light was available for us). The spectra in Fig.2 were normalized to the local “continuum” for illustrative purpose. It is visible that the 5800 Å line is much weaker in SN 1998aq than in SN 1999by. This line gets stronger in pre-maximum spectra of SNe Ia when approaching the maximum, however, in “normal” SNe the increase of $`R`$(SiII) is only 0.04$``$0.06 (Riess et al., 1998b). Thus, the 5800 Å feature was not much stronger in the spectrum of SN 1998aq at maximum.

## 3 Data analysis and interpretation

### 3.1 Line depths and radial velocities

We have calculated the value of $`R`$(SiII) of SN 1999by following the definition of Nugent et al. (1995). $`R`$(SiII) = 0.63 was obtained for the May 11 spectrum, while $`R`$ = 0.66 was determined for the May 14 one, in accord with the result of Garnavich et al. (1999) (see also Hatano et al., 2000). This is more than a factor of 2 larger than its usual value of $`R=0.25\pm 0.1`$ for normal SNe Ia (Riess et al, 1998b, Hatano et al., 2000).

Very recently Howell et al. (2001) published three spectra of SN 1999by made just before our observations. From their Fig.1 we estimated $`R`$(SiII) = 0.68, 0.70, 0.66 at epochs $`\tau =2,1,0`$, respectively. These agree very well (within a few hundredths) of the results derived from our spectra, suggesting that our observations probably do not contain wavelength-dependent systematic errors distorting the line-depth ratio.

$`R`$(SiII) correlates with the luminosity parameter $`\mathrm{\Delta }`$ introduced by Riess et al. (1996) in the Multi-Color Light Curve Shape (MLCS) method. A calibration of the linear dependence of $`R`$(SiII) on $`\mathrm{\Delta }`$ was given by Riess et al. (1998b), resulting in $`\mathrm{\Delta }=2.33R0.55(\pm 0.22)`$ for $`+1`$ day phase, while for $`+4`$ days the slope and zero point changes to 2.06 and $``$0.46, respectively. Substituting the $`R`$(SiII) values given above one can get $`\mathrm{\Delta }=0.90\pm 0.23`$ for the luminosity parameter of SN 1999by (the given error takes into account the quoted uncertainty of the $`R`$ \- $`\mathrm{\Delta }`$ calibration). For “normal” SNe Ia $`\mathrm{\Delta }0`$, therefore, this high value suggests that SN 1999by is one of the most sub-luminous SNe Ia (Bonanos et al., 1999). Note that this spectroscopic $`\mathrm{\Delta }`$ is not suitable for describing the multi-color light curve of SN 1999by, which requires $`\mathrm{\Delta }>1.0`$. The light curve decline parameter, $`\mathrm{\Delta }m_{15}(B)=1.87`$, is also indicative of the low luminosity (Bonanos et al., 1999).

Radial velocities were inferred from the absorption minima of the 6150 Å Si II trough (fitted by a low-order polynomial around minimum). Table 1 lists the resulted radial velocities, corrected for the host galaxy center recession velocity $`v_{gx}=638`$ kms-1 (as given in the NASA/IPAC Extragalactic Database). As the anonymous referee pointed out, SN 1999by was closer to the edge of NGC 2841 and could have significant orbital velocity in addition to the center-of-galaxy radial velocity. The rotation velocity of NGC 2841 at distances greater than 10 kpc from its center is $`300`$ kms-1 (Sofue et al., 1999). Thus, the uncertainty of the radial velocities is at least $`\pm `$ 300 kms-1 due to the possible orbital velocity as well as the large width of the Si II line. Nevertheless, this uncertainty is only a few percent of the Si II expansion velocity.

Fig.3 compares the radial velocities of SN 1999by with those of SN 1998aq (Vinkó et al.,1999) and SN 1994D (Patat et al.,1996). The decline of the radial velocities of SN 1999by is much stronger than in the other two cases. SN 1991bg showed similar behavior (Filippenko et al.,1992), while the other sub-luminous SNe Ia were not well covered spectroscopically around maximum. The stronger velocity decline is also in good agreement with the recent result of Hatano et al. (2000) that sub-luminous events have lower velocities at 10 days after maximum than do normal SNe Ia. This may support the hypothesis that sub-luminous SNe Ia arise from sub-Chandrasekhar mass progenitors, thus have less energetic explosion than normal Type Ia events (see Modjaz et al., 2000 for recent review of explosion models).

### 3.2 Absolute magnitude and distance

Application of the Multi-Color Light Curve Shape (MLCS) method for SN 1999by has been published very recently by Tóth & Szabó (2000). They used only the $`V`$ light curve for template fitting and $`(BV)`$ data for the reddening determination, and revealed $`\mathrm{\Delta }>1.6`$ mag. This may indicate that SN 1999by is extremely sub-luminous. For example, SN 1991bg, the prototype, had $`\mathrm{\Delta }=1.44\pm 0.1`$ (Riess, Press & Kirshner, 1996) which is known as one of the least luminous SNe Ia to date.

We re-analysed the $`BVRI`$ photometric dataset published by Tóth & Szabó (2000) supplemented by the data of Hanzl (1999) using the original MLCS template curves of Riess, Press & Kirshner (1996). Because the photometric data are inhomogeneous, systematic errors may distort the parameters derived from them. Thus, the results below can be considered only as crude estimates.

Although Riess et al. (1998a) developed a second version of MLCS (MLCS-2, based on an extended list of calibrating SNe), the updated template vectors have not been published yet. However, recently Vinkó et al. (2001) demonstrated that for nearby SNe, even with moderate reddening, the MLCS-1 and MLCS-2 methods produce consistent results, if the maximum magnitude of the fiducial ($`\mathrm{\Delta }=0`$) light curve is set properly. We adopted $`M_V^{max}=19.46`$ mag for the fiducial $`V`$ light curve (Riess, personal communication), while the original MLCS-1 method was based on $`M_V^{max}=19.36`$ mag. The difference has been taken into account by adding $`+0.1`$ mag to the distance modulus computed with MLCS-1.

For SN 1999by, $`E(BV)=0.01`$ was adopted, consistent with the all-sky reddening map of Schlegel et al. (1998) (Tóth & Szabó derived $`E(BV)=0.05\pm 0.02`$ from their SN photometry). The template vectors were fitted simultaneously to all $`BVRI`$ data, i.e. the residual of all data were combined in a single $`\chi ^2`$-function. The fitting was restricted for $`\tau >2`$ days, because the available pre-maximum data of SN 1999by could not be adequately fitted together with the post-maximum data (see also Tóth & Szabó, 2000). Also, the minimum photometric error was increased to $`\pm 0.15`$ mag taking into account some possible systematic uncertainties in the used dataset (due to e.g. errors in the standard transformation).

The simultaneous fitting of the template curves resulted in $`\mathrm{\Delta }=1.4\pm 0.1`$ and $`\mu _0=31.16\pm 0.15`$ mag for the luminosity parameter and the reddening-free distance modulus, respectively. The given uncertainty of the distance modulus is only the rms error of the fitting, disregarding the possible systematic errors in the absolute magnitudes of the template light curves. The quality of the fit can be judged from Fig.4. It is seen that the fitting is not perfect, however, we believe that this is close the optimal result that can be achieved from the present datasets. The minimum of the reduced chi-squared that could be reached was $`\chi ^21.5`$. Still, the solution may be systematically incorrect, because for such under-luminous SN the template vectors in MLCS-1 are based only on the light curves of SN 1991bg, and the light variation of SN 1999by may differ substantially from that. The large deviation of the pre-maximum points may be due to such inconsistence. Also, photometric errors higher than we adopted ($`\pm 0.15`$ mag) cannot be ruled out.

The photometric $`\mathrm{\Delta }`$ is significantly higher than the value found from spectroscopy ($`\mathrm{\Delta }=0.90`$, Sect.3). It was not possible to reach an adequate fit with $`\mathrm{\Delta }=0.9`$ given by spectroscopy. Thus, the current spectroscopic determination of $`\mathrm{\Delta }`$ does not give a suitable result for SN 1999by, which is probably true for other extreme sub-luminous SNe Ia as well. An improved $`R`$(SiII) $``$ $`\mathrm{\Delta }`$ calibration is presented in the next section.

Although sub-luminous SNe Ia are usually eliminated from the sample of “good” distance indicators, the MLCS method gives information on their distances as well. The new distance modulus of SN 1999by given above corresponds to $`17.1\pm 1.2`$ Mpc geometric distance, which is significantly higher than $`d=1012`$ Mpc given by Tóth & Szabó (2000) and Tully (1988). Note that the template vectors of the MLCS-1 method were calibrated from a training set of SNe adopting Tully-Fisher- and surface-brightness-fluctuation distances, but the final absolute magnitudes were shifted to match the Cepheid-based distance scale (Riess, Press & Kirshner, 1996). Thus, the higher distance given above simply reflects the difference between various distance scales, and may still contain a systematic error. Both Cepheid- and Tully-Fisher distances undergo significant revisions nowadays (e.g. Gibson et al., 2000, Sakai et al., 2000), so this question can be answered only when the problems of the different distance scales are solved.

### 3.3 Comparison with other SNe Ia

The absolute magnitudes of SNe Ia correlate with $`R`$(SiII) as was presented very clearly by Nugent et al. (1995). Although Hatano et al. (2000) demonstrated that this correlation is not a one-parameter sequence, the dependence of $`M_V^{max}`$ on $`R`$(SiII) may still be a somewhat useful relation for predicting SNe luminosities from spectra observed around maximum (at least as a first approximation). Now the correlation can be slightly improved using the new data of sub-luminous SNe Ia.

We have collected values of the light curve parameter $`\mathrm{\Delta }`$ that were derived via the MLCS method and the values of $`R`$(SiII), both at $`V`$maximum. These are listed in Table 2. Because $`\mathrm{\Delta }`$ measures the deviation of the $`V`$ light curve of a particular SN from the fiducial SN Ia light curve at maximum, it is also directly proportional to the maximum absolute magnitude $`M_V^{max}`$. However, at present, $`\mathrm{\Delta }`$ seems to be a more useful parameter than $`M_V^{max}`$, because it is determined directly from the shape of the (multi-color) light curves, therefore it is not affected by the adopted zero-point of the distance scale.

The result is plotted in Fig.5, where the continuous line represents the approximate relation

$`\mathrm{\Delta }=`$ $`3.5R\text{(SiII)}`$ $`0.8`$ (1)
$`\pm 0.44`$ $`\pm 0.17`$
while the dashed line shows the relation given by Riess et al. (1998b) (see Sect.3.1). It is seen that the Riess et al. calibration gives systematically lower $`\mathrm{\Delta }`$ for the most sub-luminous SNe Ia. The new relation fits these SNe better, thus it may be used to predict $`\mathrm{\Delta }`$ from the measurement of $`R`$(SiII) around maximum. Such prediction may be useful, for example, in the “snapshot” distance estimate method (Riess et al., 1998b) when only a few nights’ data are available for a particular SN.

From Fig.5, it may also be interesting that the most sub-luminous SN 1991bg, SN 1997cn and SN 1999by represent a small group with very consistent line-depth ratio and $`\mathrm{\Delta }`$ (or absolute magnitude). It may mean that sub-luminous SNe Ia might be better distance indicators than previously thought, if their sub-luminous status is properly identified from spectroscopy and their difference from normal SNe Ia is taken into account. It is known (e.g. Hamuy et al., 1994) that such sub-luminous SNe Ia may be more frequent events than their discovery rate suggests, thus, it is possible that many more such objects will be detected in the extensive supernova-search programs. However, the present sample is obviously too small to draw a certain conclusion.

## 4 Summary

The results of this paper can be summarized as follows:

Three spectra of the sub-luminous Type Ia SN 1999by around maximum were obtained at DDO, centered on the Si II trough at 6150 Å. The $`R`$(SiII) line depth ratio (Nugent et al., 1995) is determined to be 0.63 at maximum, indicative of an intrinsically sub-luminous SNe Ia.

Radial velocities of the 6150 Å Si II line were determined from measuring the Doppler-shift of the line core. It is shown that SN 1999by exhibited stronger decline of its radial velocity curve than other, spectroscopically normal SNe Ia. This may support the low-energy, sub-Chandrasekhar mass explosion proposed for such SNe.

Absolute magnitude and distance of SN 1999by have been derived via the MLCS method. The assumption that the reddening is $`E(BV)=0.01`$ led to $`M_V^{max}`$=$``$18.06$`\pm 0.1`$ mag and $`d=17.1\pm 1.2`$ Mpc. The maximum absolute magnitude as well as the line depth ratio $`R`$(SiII) of SN 1999by are in very good agreement with those of other extreme sub-luminous SNe Ia, SN 1991bg and SN 1997cn. An approximative linear relation has been determined to predict the light curve parameter $`\mathrm{\Delta }`$ of SNe Ia over a wide range of observed $`R`$(SiII).

This research has been supported by Hungarian OTKA Grant No.T032258, the Magyary Postdoctoral Fellowship to JV from Foundation for Hungarian Education and Science (AMFK), the “Bolyai János” Research Scholarship to LLK. Thanks are due to the David Dunlap Observatory for allocating the necessary telescope time. The NASA Astrophysics Data System, the SIMBAD database and the Canadian Astronomy Data Centre were used to access data and references. The availability of these services are gratefully acknowledged.