text:
# Giant Molecular Clouds in M33 I – BIMA All-Disk Survey

## 1. Introduction

This paper presents the results of an interferometric survey in 12CO ($`J=10`$) of the entire H$`\alpha `$ disk of M33. Observations were made with the Berkeley-Illinois-Maryland Association array. The proximity and low inclination of M33, an SA(s)cd spiral in the Local Group, permits direct and detailed comparisons of molecular, atomic, and ionized hydrogen gas on the scale of individual GMCs. The angular resolution of the BIMA survey is $`13^{\prime \prime }`$, which corresponds to a linear resolution of 50 pc at the distance of M33 (0.85 Mpc, Lee, Freedman, & Madore, 1993). This resolution is adequate to separate, but not to resolve, individual giant molecular clouds (GMCs) if they are similar in size to those in the Milky Way (Blitz, 1993). A follow-up paper by (Rosolowsky, et al., 2003, Paper II) describes a study of rotation properties, virial masses, and morphology of 36 GMCs from this survey which were reobserved at BIMA with 20 pc linear resolution. The BIMA array is well-suited for this extensive mapping project for three reasons. (1) The 6-meter diameter antennas have a large ($`100^{\prime \prime }`$ FWHM) field of view at 115 GHz; this keeps the number of pointing centers to a manageable level. (2) With 10 telescopes, the array affords good snapshot imaging capability; visiting each field only a few times per night yields a clean synthesized beam. (3) The SIS receivers (Engargiola & Plambeck, 1998) have excellent sensitivity; double sideband receiver temperatures at 115 GHz are typically 40 K.

Of the numerous mechanisms proposed for molecular cloud formation (Elmegreen, 1990), it is uncertain which, if any, is responsible for most GMCs in spiral galaxies. Validating any cloud formation model requires an unbiased survey of molecular clouds in a galaxy. Prior to the present study of M33, there existed no complete catalog of giant molecular clouds in an external spiral galaxy. Even cloud catalogs in the Milky Way are incomplete because velocity blending along the line-of-sight can prevent unambiguous separation into individual clouds (e.g. Issa, MacLaren, & Wolfendale, 1990). Only in the Large Magellenic Cloud, a Local Group irregular, has CO been surveyed with sufficient resolution to measure accurately the GMC mass spectrum (Fukui et al., 2001).

A decade ago, Wilson & Scoville (1990) undertook a pioneering survey of extragalactic GMCs in the nuclear region of M33. With the Owens Valley interferometer they mapped 12CO ($`J=10`$) in 19 fields with a $`1^{}`$ field-of-view and a $`7^{\prime \prime }`$ synthesized beam. Their data analysis yielded a catalog of 38 GMCs with $`(0.2\text{ — }4)\times 10^5M_{\mathrm{}}`$ and a mass distribution index similar to that of the Milky Way. Although their observations covered a significant fraction of the central 2 kpc, they positioned most fields to include areas of high optical extinction detected in CO emission with the NRAO 12 m telescope, possibly biasing the cloud sample. This paper attempts to address this bias by presenting a complete survey of GMCs in M33 over most of the optical disk.

## 2. Observations

Observations were made with the BIMA array in its most compact (“D”) configuration at various times between 1999 August and 2000 October. As shown in Figure 1, M33 was mosaiced using a grid of 759 pointing centers, spaced by $`78^{\prime \prime }`$ in a hexagonal pattern. Circles indicate the FWHM primary beamwidth of the BIMA antennas, which is $`100^{\prime \prime }`$ at 115.27 GHz. The grid spacing was chosen to cover the galaxy with as few pointing centers as possible while maintaining a point source sensitivity that varies by no more than 10% across the mosaic.

The CO($`10`$) line was observed in the upper sideband of the SIS mixers. Single sideband system temperatures, scaled to the top of the atmosphere, were typically 300 to 500 K at 115.27 GHz. The correlator covered the LSR velocity range $`414`$ to $`+59`$ km s-1 with 1.015 km s-1 channel spacing, on each of the 45 baselines.

Each night we cycled through observations of about 50 grid points, using an integration time of 25 seconds per point. The slew and settle time of the telescope between integrations was 5 seconds, so a complete loop through 50 grid points required about 25 minutes. Observations of the calibrator, 0136$`+`$478, were interleaved with complete loops. For a typical 6-hour track, each pointing center was observed 10–12 times. With 45 baselines, this yields about 500 measurements of the source visibility across the $`(u,v)`$ plane, resulting in a nearly round synthesized beam with $`10`$% sidelobes.

A total of 45 observing tracks was obtained, with durations of a few hours to 10 hours. The total useful integration time accumulated on M33 was 132 hours, which corresponds to an average of 10.5 minutes per pointing center.

### 2.1. Flux Scale

For each observing track the amplitude gain of each antenna was scaled by a constant factor to obtain the correct visibility amplitude on the calibrator, 0136$`+`$478. The flux density of this quasar was variable over the course of the observations (2.3 Jy in 1999 September, 3.0 Jy in 2000 May, 3.7 Jy in 2000 July, and 3.4 Jy in 2000 Oct). We monitored the calibrator flux density by making a brief observation of the compact H II region W3OH once per night. The flux density of W3OH was measured to be $`3.8\pm 0.2`$ Jy by comparison with primary flux calibrators Mars and Uranus. All comparisons with W3OH were made using the lower sideband continuum channel of the receivers at 112 GHz, since measurements at 115 GHz are contaminated by CO emission from the extended molecular cloud associated with W3OH.

The flux densities derived for 0136$`+`$478 were generally reproducible to $`\pm 0.2`$ Jy from night to night. Allowing for possible systematic pointing errors and decorrelation by atmospheric phase fluctuations, we estimate that the final amplitude scale is uncertain by approximately $`\pm 15`$%.

### 2.2. Mapping

After calibrating the visibility data, we generated mosaiced spectral line maps using the Miriad software package111http://www.atnf.csiro.au/computing/software/miriad. The maps were deconvolved with a CLEAN algorithm in order to remove negative bowls around the strongest sources. The synthesized beam varies slightly across the mosaic because individual subgrids have slightly different $`(u,v)`$ coverage. To handle this, Miriad creates a cube of synthesized beams, one for each pointing center, and interpolates among them when deconvolving the map. The maps were restored with a $`13^{\prime \prime }`$ Gaussian.

Most of our analysis was performed on maps with $`2.03`$ km s-1 velocity resolution. These were produced by Hanning smoothing the (1.015 km s-1 spacing) visibility spectra, then generating maps of alternate channels. The resulting channel maps are expected to be nearly statistically independent.

To evaluate the noise level across the mosaic, we computed the RMS fluctuations in the spectrum corresponding to each map pixel, omitting any channels with flux densities $`S>4\sigma `$. From a spectrum with 222 statistically independent channels, the noise is determined to approximately 5% accuracy. We then smoothed the map of the noise to $`30^{\prime \prime }`$ resolution (Figure 2). The white 99% gain contour denotes the effective outer boundary of the mosaic. Outside this contour intensities are not fully corrected for attenuation by the primary beam. For these 2 km s-1 resolution channel maps the RMS noise level over the central half of the mosaic ranges from 0.34 Jy beam-1 to 0.51 Jy beam-1; the average value is 0.44 Jy beam-1, or 0.24 K. At the worst spots near the edge of the mapped region (but inside the 99% gain contour) the noise level is sometimes as high as 0.65 Jy beam-1.

## 3. The GMC catalog

### 3.1. Identifying clouds

Because the sensitivity varies across the mosaic, we searched for molecular clouds by clipping the maps at a constant signal-to-noise ratio instead of a constant flux density threshold. The intensity data cube is converted to a signal-to-noise data cube by dividing through by the noise at each position:

$$s(x,y,\nu )=\frac{I(x,y,\nu )}{\sigma (x,y)f(\nu )},$$
(1)
where $`\sigma (x,y)`$ is the smoothed RMS noise map (Figure 2) and $`f(\nu )`$ is a frequency-dependent noise correction (Figure 3) which arises because of gain variations of $`510\%`$ across the instrumental passband. Real CO emission occupies only a few of the $`20500`$ independent resolution elements in each channel map, thus having little effect on the renormalization. In the normalized, or signal-to-noise, data cube the intensity is expressed as a multiple of the RMS noise at each pixel. For example, $`s(x,y,v)=3`$ represents a $`3\sigma `$ detection of a negative peak.

The maximum and minimum normalized fluxes in each velocity channel are plotted in Figure 4. CO emission is obvious in the range $`260<V_{\mathrm{LSR}}<105`$ km s-1. A histogram of the normalized fluxes in this velocity range (Figure 5) is well-fitted by a normal distribution, apart from the wing of high signal-to-noise pixels attributable to CO detections.

Clearly, one could comfortably accept as a CO detection any pixel with normalized flux $`>5\sigma `$. Unfortunately, such a simple strategy would miss much real molecular emission – Figure 5 shows a clear excess in pixels with flux $`>3.5\sigma `$. Since the velocity maps are nearly statistically independent, we can obtain a higher detection efficiency by searching for neighboring channels with unusually high normalized fluxes. The lowest mass molecular cloud which our survey is likely to detect is of order $`2\times 10^4M_{\mathrm{}}`$. A typical Milky Way GMC with this mass has a velocity dispersion of 1.8 km s-1 (Solomon et al., 1987), which corresponds to a FWHM linewidth of 4.2 km s-1. Therefore, we expect that CO emission from most real clouds should appear in at least two adjacent 2 km s-1 velocity channels. To identify molecular cloud candidates, we wrote an IDL routine to locate pairs of adjacent channels with normalized flux $`>2.7\sigma `$. Each pair, and all contiguous pixels with normalized flux $`>2\sigma `$, are assigned to a candidate molecular cloud. These pixels are then masked to avoid double counting, and a search for another candidate cloud proceeds until all pairs of adjacent $`>2.7\sigma `$ channels have been masked. The assignment of pixels to candidate clouds in this manner does not depend on search order but occasionally links nearby sources which could be regarded as separate clouds. Although we expect from our map statistics to find real CO emission only in the velocity range $`280<V_{\mathrm{LSR}}<85`$ km s-1, the search was conducted over the full 450 km s-1 range of the data cube. A total of 1600 candidate sources were identified in this way.

To evaluate the likelihood that each candidate source is real, we examine the most significant spectrum through the cloud – either the one with the highest summed significance or the one with the greatest number of adjacent channels with flux $`>2\sigma `$. The calculation described in Appendix A evaluates the probability $`P`$ that this spectrum is a noise fluctuation, i.e. that the detection is spurious. As discussed in Appendix A, the likelihood that the cloud is real is then given by $`\mathrm{exp}(NP)1NP`$, where $`N`$ is the number of statistically independent ways in which $`n`$ adjacent channels can be selected from the data cube. For the full data cube $`Nn_{beams}n_{chans}4.5\times 10^6`$. Note that the likelihood computation is based on a single spectrum through each cloud, not on an average over the source. This keeps the calculation simple, since adjacent spatial pixels oversample the synthesized beam and hence are not statistically independent.

Originally, we intended to catalog all candidate sources with $`10\%`$ likelihoods of being noise fluctuations. To our disappointment, however, the algorithm identified several sources which were almost certainly spurious – i.e., well outside the plausible range for CO emission – with false detection likelihoods of $`<0.001`$. There are several possible explanations for this inconsistency. First, there may be outliers caused by a few bad visibility measurements, perhaps from observations obtained in especially poor weather or at low elevation. Second, we may be underestimating the number $`N`$ of statistically independent resolution elements in the data cube. Nyquist sampling the highest spatial frequency visibility data (corresponding to antenna separations of 11.5 k$`\lambda `$) requires a $`9^{\prime \prime }`$ interval, which would imply there are $`50000`$ independent resolution elements in the map. We believe the actual number is close to our original estimate of 20500 because the long baseline data receive little weight in the final image. Third, there may be errors in normalizing the maps by the rms noise. Such errors are particularly likely for the frequency dependent correction (Figure 3) because there were slow changes in the passband gain over the 14 month period in which observations were made.

Because the likelihood calculation is imperfect, we use it only to rank the relative likelihoods of the candidate sources. For convenience we define a ‘likelihood measure’ as $`\mathrm{ln}(NP)`$. Figure 6 plots this quantity vs. LSR velocity for the 1600 candidate sources. Sources with the 60 highest likelihood measures ($`\mathrm{ln}(NP)>13.4`$) fall in the velocity range $`260<V_{\mathrm{LSR}}<105`$ km s-1and are apparent in multiple channel maps, i.e. they are unquestionably real. As the likelihood threshold is lowered, we begin to find a few sources that probably are spurious, since they lie well outside this velocity range. At $`\mathrm{ln}(NP)0`$ the number of obviously spurious sources rises sharply, indicating the noise floor of the data cube.

For our primary GMC catalog, we accept all cloud candidates in the velocity interval $`280<V_{\mathrm{LSR}}<85`$ km s-1 for which $`\mathrm{ln}(NP)>5.83`$. A total of 93 sources fit these criteria. It is straightforward to estimate the number of false detections in this catalog simply by counting the number of sources above the likelihood cutoff which lie outside the 96-channel search range. There are 3 spurious detections in the 48 channels just below the search range, and 3 spurious detections in the 48 channels just above the search range. Hence, we expect that there are of order 6 spurious detections in this catalog.

One can identify fainter sources by narrowing the search range to cover only velocities that lie within the H I emission line profile in each direction. Figure 7 shows the cloud likelihoods from Figure 6 plotted as a function of the H I–CO velocity difference. H I velocities were obtained from the Westerbork map of Deul & van der Hulst (1987), which has $`12\times 24^{\prime \prime }`$ resolution. For the 93 clouds in the primary catalog, the velocity difference exceeds $`\sigma _{\mathrm{HI}}`$ for only 5 sources, where $`\sigma _{\mathrm{HI}}`$ is the velocity dispersion of the H I line. The four sources with velocity differences $`>2\sigma _{\mathrm{HI}}`$ may be false detections.

Narrowing the search range to $`\pm 2\sigma _{\mathrm{HI}}\pm 15`$ km s-1 from the H I velocity, we identified 55 additional clouds with likelihood measure $`\mathrm{ln}(NP)>1.4`$. Figure 8 shows these sources on a likelihood–velocity difference plot. Once again, by counting the number of (presumably) spurious sources outside the velocity search range, we estimate that there are of order 9 false detections in this catalog extension.

It is not useful to catalog individual clouds with likelihoods $`\mathrm{ln}(NP)<1.4`$ because the number of spurious sources is quite large. Nevertheless, in Figure 8 there is an obvious overdensity of GMC candidates with $`V_{\mathrm{LSR}}`$ within $`2\sigma _{HI}`$ of $`V_{\mathrm{HI}}`$, and this may be used to gauge the number of cloud candidates which are falsely rejected. There are 145 candidate sources with $`2<V_{\mathrm{CO}}V_{\mathrm{HI}}/\sigma _{\mathrm{HI}}<4`$ which are presumably spurious. One expects to find an equal number of spurious sources in the range $`V_{\mathrm{CO}}V_{\mathrm{HI}}/\sigma _{\mathrm{HI}}<2`$. Since there are actually 205 candidate sources in this range, we estimate that roughly 60 clouds have been falsely rejected from our GMC catalog in this significance range. This result is not particularly sensitive to how we select boundaries for the overdensity band in Figure 8.

### 3.2. Cloud Positions, Velocities, and Masses

Figure 9 shows the most significant CO spectrum through each cloud, overplotted on the H I spectrum in the same direction. Table 1 lists the positions and properties of the 148 clouds identified in our survey, ranked by likelihood. Sources 1–93 were chosen without regard to the CO–H I velocity difference, and hence are suitable for statistical tests against the H I distribution.

The position and velocity of each cloud listed in Table 1 are intensity-weighted averages over all pixels which constitute the cloud i.e., $`\overline{x}=(Tx)/(T)`$. The velocity width is derived from the intensity-weighted second moment of the velocity distribution,

$$\mathrm{\Delta }v_{\mathrm{FWHM}}=\left[8\mathrm{ln2}\frac{T(v\overline{v})^2}{T}\right]^{0.5}.$$
(2)
To estimate the cloud’s mass, we convert the CO integrated intensity for each pixel to an H2 column density $`N(\text{H}_2)=T(x,y,v)\mathrm{\Delta }vX`$, where $`\mathrm{\Delta }v`$ is the channel velocity width, 2.032 km s-1, and where we assume a conversion factor $`X=2\times 10^{20}\text{cm}^2/(\text{K km s}^1)`$ (Solomon et al., 1987; Strong et al., 1988; Strong & Mattox, 1996; Digel et al., 1999). The higher resolution observations presented in Paper II imply that this conversion factor is consistent with the virial masses of clouds out to radii of 4.5 kpc, despite variations in metallicity. Each pixel corresponds to a physical area $`D^2\mathrm{\Delta }\alpha \mathrm{\Delta }\delta `$, where $`D`$ is the distance to M33, 850 kpc (Lee, Freedman, & Madore, 1993, and references therein), and $`\mathrm{\Delta }\alpha `$ and $`\mathrm{\Delta }\delta `$ are the extent of the pixel in angle, both of which are $`4^{\prime \prime }`$. Thus, the total cloud mass is

$$M=\underset{\mathrm{pixels}}{}1.362m_HXT_{\mathrm{pixel}}\mathrm{\Delta }vD^2\mathrm{\Delta }\alpha \mathrm{\Delta }\delta .$$
(3)
The factor 1.36 is the mass correction for helium with a number fraction of 9% relative to hydrogen nuclei. If expressed in terms of flux density rather than brightness temperature, this is equivalent to

$$M=7.5\times 10^3M_{\mathrm{}}/(\text{Jy km s}^1).$$
(4)
This is 25% lower than the conversion factor adopted by Wilson & Scoville (1990), mostly because of our choice of a lower $`X`$-factor.

### 3.3. Completeness

To gauge the completeness of the survey, we added a regular lattice of artificial test sources to the cube of channel maps and computed the fraction of them recovered by our search algorithm. Each test source was a $`13^{\prime \prime }\times 13^{\prime \prime }`$ Gaussian, corresponding to a GMC unresolved by our synthesized beam, with a Gaussian velocity profile. The velocity width was kept constant at 6 km s-1 FWHM, independent of the flux density (mass) of the test source. Artificial sources were spaced by $`100^{\prime \prime }`$ in right ascension and declination, over the entire mosaic; and by 56.3 km s-1 in velocity, from $`345`$ to $`11`$ km s-1. Sources close to the edge of the mosaic were omitted if more than 5% of their flux lay outside the mosaic’s 99% gain contour. In total, 2625 artificial sources were added to the data cube. Source recovery tests were done with masses ranging from $`3\times 10^4M_{\mathrm{}}`$ to $`10^6M_{\mathrm{}}`$.

Figure 10 shows the fraction of the test sources detected as a function of the local RMS noise level. Over the full mosaic, the survey is complete to $`1.5\times 10^5M_{\mathrm{}}`$. Within the central half of the mosaic, where the noise is typically $`<0.5`$ Jy beam-1, we recover approximately 90% of $`1\times 10^5M_{\mathrm{}}`$ clouds. Outside this region, we recover 70 — 80% of such clouds.

Figure 11 shows the recovered mass versus the actual mass of the test sources. Our algorithm underestimates cloud masses because it integrates the CO luminosity only over pixels inside the 2$`\sigma `$ brightness contour. The underestimates are particularly severe for low mass clouds, for which much of the emission is beneath below the 2$`\sigma `$ noise level. For clouds at the threshold of detectability, $`7\times 10^4M_{\mathrm{}}`$, we recover about half the mass, whereas for clouds more massive than $`2\times 10^5M_{\mathrm{}}`$, we recover more than 80%. The average correction factors from Figure 11 were applied to our catalog sources to obtain more accurate mass estimates; these are given as $`M_{\mathrm{CO}}^{}`$ in column (8) of Table 1.

### 3.4. Extended sources

The completeness estimates discussed above are applicable only to sources that are unresolved by the $`13^{\prime \prime }`$ synthesized beam, i.e. to GMCs $`50`$ pc in diameter. Because an aperture synthesis telescope acts as a spatial filter, attenuating low spatial frequencies, the sensitivity to larger sources is reduced. Essentially, the synthesized beam consists of a narrow positive lobe surrounded by a broad negative bowl. If one convolves an extended source brightness distribution with this beam, the negative contribution from the bowl partially cancels the positive contribution from the central lobe, reducing the observed flux density. For the BIMA D-array, the observed flux density for a $`40^{\prime \prime }`$ FWHM Gaussian source is about half the true flux density. Thus, our survey may not detect a group of molecular clouds with a total mass of $`1\times 10^5M_{\mathrm{}}`$ spread out over a $`40^{\prime \prime }`$ (150 pc) region, though it is quite likely to detect a single GMC with this mass. For sources with radii comparable to the beam size, their surface brightness will be diminished relative to that of a point source by a factor of $`2`$. Fortunately, sources of this size are also significantly more massive ($`10^6M_{}`$, Paper II) than our completeness limit and would be readily detectable.

Single dish observations of M33 with the University of Arizona Steward Observatory 12-m telescope were used to estimate the total fraction of the CO flux recovered by the BIMA survey. The single dish flux calibrations were repeatable within 6%, and agreed with BIMA flux measurements of planetary and spectral line calibrators within 10%; further observational details are given in Paper II. Two sets of single dish observations were made. First, 18 positions with BIMA catalog sources were observed. For these fields the BIMA D-array maps recovered approximately 60% of the single dish CO flux. CO was detected in just 2 of the 18 off-source spectra obtained in these position-switched observations; the 16 non-detections may be used to set an upper limit of 0.3 M pc-2 for the surface density of an extended ($`R4`$ kpc) uniform molecular disk. Second, we observed 15 fields to map a $`7.5^{}`$ cut along the major axis of M33. In Figure 12 we show a composite single dish spectrum derived by shifting the spectrum for each of the 15 fields to a common velocity and averaging. The integrated flux of the of the 12 m spectrum is 1.2 K km s-1 corresponding to $``$5.2 M pc-2. This corresponds to $`1.6\times 10^7M_{\mathrm{}}`$ of molecular gas in the central $`7^{}`$ diameter of the galaxy. Wilson & Scoville (1989) measure $`1.9\times 10^7M_{\mathrm{}}`$ in the same region when converted to our $`X`$-factor, distance, and removing their correction for source-beam coupling.

For comparison between the 12 m and the interferometer, we show composite spectra made from BIMA data within the 12 m fields, convolved to the $`54^{\prime \prime }`$ FWHM 12 m beam, sampled at the centers of the 12 m beams, and averaged together in a fashion identical to the 12 m spectra. Two BIMA spectra were synthesized, one from the catalog source pixels only (thick line), another from the data cube clipped at 2$`\sigma `$ (dotted line). The BIMA catalog spectrum contains 43% of the single dish flux, the clipped spectrum, 54%. In Paper II we argue that most of the missing flux is attributable to complexes of undetected low mass clouds associated with catalog sources.

We can estimate the total molecular mass of M33 from the major axis cut. The integrated single dish molecular mass in this $`7.5^{}`$ long strip is $`2.8\times 10^6M_{\mathrm{}}`$. The corresponding BIMA mass, considering only cataloged clouds with masses greater than $`1.5\times 10^5M_{\mathrm{}}`$, the completeness limit for our survey, is $`9.2\times 10^5`$ M, 33% of the single dish mass. In this computation we have weighted the cloud masses by the appropriate 12-m beam attenuation factors. The total (corrected) mass in the catalog is $`2.3\times 10^7M_{\mathrm{}}`$, of which 1.5$`\times 10^7M_{\mathrm{}}`$ is from clouds more massive than the completeness limit. If the GMC mass spectrum in the central cut is similar to that for the galaxy as a whole, then we infer that the total molecular mass is 3 times this value, or 4.5$`\times 10^7M_{\mathrm{}}`$.

### 3.5. Comparison with Wilson and Scoville 1990

Wilson & Scoville (1990, hereafter WS90) used the Owens Valley array to map 19 fields in the nucleus of M33 with $`7^{\prime \prime }`$ resolution, identifying 38 molecular clouds. Figure 13 is an overlay of the WS90 fields on the central $`10^{}\times 10^{}`$ section of the BIMA survey. Red contours show the outer boundary of clouds identified in the BIMA survey; solid circles mark locations of clouds cataloged by WS90. Frequently, a BIMA cloud corresponds to two or more WS clouds – for example, the complex WS1-WS5 at the center of the galaxy is classified as one cloud in our survey because these sources are connected by low-level emission and lower BIMA resolution blends these clouds together. The BIMA survey detects all of the WS90 clouds with masses $`>2\times 10^5M_{}`$.

To examine a few of our clouds more closely, and to investigate discrepancies between our catalog and the WS90 catalog, we obtained higher resolution BIMA observations of the three fields shown by dashed circles in Figure 13. These fields were observed on 3 nights in 2000 July with the array in its “C” configuration. The integration time was approximately 6 hours per field. The C-array observations were calibrated and combined with the D-array data to generate maps with a $`6.5^{\prime \prime }`$ synthesized beam. The RMS noise level in each 2.03 km s-1 channel map was 0.15 Jy beam-1, approximately 0.32 K, comparable to that of WS90. The C$`+`$D-array maps are compared with the D-array survey images in Figure 14.

Generally there is good agreement between sources in our C+D array images and the WS90 positions, although we fail to confirm some of the less massive WS90 clouds (e.g., WS21, WS22, WS36). One can find several examples showing how the D-array survey blends WS sources (e.g., WS27 and WS28 = catalog source no. 17; WS6 + WS8 = catalog source no. 74).

Figure 15 compares the masses of GMCs which can be clearly identified in both data sets. As indicated by the dashed line, the BIMA masses are expected to be systematically lower by 25% because we used a lower conversion factor from flux density to mass. Generally the BIMA and WS90 masses are in reasonable agreement, although there is a good deal of scatter in the ratio, probably because of the difficulty in establishing the outer boundaries of clouds.

Small clouds detected by the C-array observations are not included in our survey catalog (Table 1) but a further analysis of these clouds is included in Paper II.

## 4. Cloud Properties

### 4.1. GMC Mass Spectrum

Histograms of the mass distribution of GMCs in our catalog are presented on a log-log scale in Figure 16. The dashed line histogram is the observed spectrum of GMC masses uncorrected for signal beneath the $`2\sigma `$ clipping level, as described in §3.3. The solid line histogram is the spectrum of corrected masses. A power law of the form

$$\frac{dN}{dM}=K\left(\frac{M}{M_{\mathrm{}}}\right)^\alpha $$

was fit to the mass bins for which the catalog is complete ($`M>1.5\times 10^5M_{\mathrm{}}`$). The observed and corrected mass distributions are well fit with similar indices ($`\alpha =2.6\pm 0.3`$) but different normalization constants ($`K=1.2`$ vs. $`1.8\times 10^{10}M_{\mathrm{}}^1`$).

For a mass spectrum with $`\alpha <2`$, most of the integrated mass resides in the most massive clouds. For a mass spectrum with $`\alpha >2`$, most of the mass is in the least massive clouds, and the total mass diverges when extended to arbitrarily small mass. Hence, for M33 there must be either a lower mass limit for molecular clouds or a break in the distribution where $`\alpha `$ changes with decreasing cloud mass to less than 2. If the integrated mass

$$M_{total}=_{lowlimit}^{\mathrm{}}M𝑑N=_{lowlimit}^{\mathrm{}}M\left(\frac{dN}{dM}\right)𝑑M$$
(5)
is distributed according to

$$\frac{dN}{dM}=1.8\times 10^{10}\left(\frac{M}{M_{\mathrm{}}}\right)^{2.6}[M_{\mathrm{}}^1]$$
(6)
we compute a strict lower mass limit $`M_L`$ of $`4\times 10^4M_{\mathrm{}}`$ for $`\alpha =2.6`$ and $`M_{total}=4.5\times 10^7M_{}`$. More likely, there is a ”turnover” mass, where the mass spectrum becomes flatter than 2. Since we can put a firm lower limit on the cutoff mass of $`4\times 10^4M_{\mathrm{}}`$ and an upper limit of $`10^5M_{\mathrm{}}`$, where the measured mass distribution is nearly complete, the turnover mass of GMCs is $`7\times 10^4M_{\mathrm{}}`$ with an uncertainty of $`3\times 10^4M_{\mathrm{}}`$. This turnover represents a “characteristic” mass for the molecular clouds in M33 that any theory of GMC formation must be able to reproduce, given the physical conditions in M33.

The highest mass we measure for a GMC in M33 is $`7.3\times 10^5M_{\mathrm{}}`$; there are two GMCs with this mass. Neither is close to the galactic center nor to an especially bright H II region. In total, we find five GMCs more massive than the highest mass cloud ($`4\times 10^5M_{\mathrm{}}`$) reported by Wilson & Scoville (1990). Moreover, they find a lower value of $`\alpha `$ (1.7), which may have resulted from preferentially targeting massive clouds by observing CO-bright single dish fields.

#### 4.1.1 M33, Milky Way, and the LMC

Over a similar mass range, the power law distribution of M33 is steeper than that of the Milky Way (1.6, Solomon et al. 1987; 1.9, Heyer, Carpenter, & Snell 2001). If M33 were as massive as the Milky Way, the normalization constant of its power law distribution would be $`3.6\times 10^{11}M_{\mathrm{}}^1`$ and $`(dN/dM)\mathrm{\Delta }M1.0`$ at $`1.0\times 10^6M_{\mathrm{}}`$. In contrast, the Milky Way is estimated to have more than 100 GMCs with $`M>1.8\times 10^6M_{\mathrm{}}`$ Dame et al. (1986). We conclude that M33, even though it has many bright H II regions (Israel & van der Kruit, 1974), has significantly fewer massive GMCs. As we infer in $`\mathrm{\S }\text{5.6}`$, the majority of H II regions must be (1) fueled by molecular clouds below our completeness threshold or (2) capable of efficiently dispersing their parent GMC complexes (Hartmann, Ballesteros-Paredes, & Bergin, 2001). Extrapolating the observed power law index to a cut-off $`3.0\times 10^4M_{\mathrm{}}`$, we find $`(dN/dM)𝑑M2000`$. This is a sufficient number of molecular clouds to match approximately the observed number of H II regions.

The NANTEN survey of CO in the LMC revealed 168 GMCs (Fukui et al., 2001). The cloud sample has a mass range and completeness limit similar to the BIMA survey. The mass spectrum of GMCs in the LMC has an index of 1.9, similar to that of the outer Milky Way (Heyer, Carpenter, & Snell, 2001), indicating that most of the mass resides in high mass clouds. The physical meaning of the difference between the mass indices of the LMC and M33 is unclear. However, as in M33, molecular clouds in the LMC appear to be tightly correlated with H I overdensities (Mizuno et al., 2001), and H I filaments underlie bright complexes of GMCs (§5.3). Hence, we expect the details of GMC formation to be at least somewhat similar for M33 and the LMC.

#### 4.1.2 Cloud Masses versus Radius

Figure 17 shows the distribution of catalog GMC masses as a function of galactocentric distance. The distribution of clouds with masses greater than our completeness limit (1.5$`\times 10^5M_{\mathrm{}}`$) appears to cut off abruptly beyond 4 kpc. For $`R_{gal}>4.5`$ kpc, we detect only a few low mass clouds despite considerable reserves of H I here (§5.1). However, higher sensitivity measurements are required to determine if we are indeed observing a sharp edge to the molecular disk or a shift to the production of predominantly low mass clouds. For $`R_{gal}<4.5`$ kpc, the mass distribution appears more or less constant to the eye; i.e. the apparent density of points in Figure 17 appears constant with $`R`$ for a given mass.

In a more quantitative comparison, we looked for spatial variations in the GMC mass spectrum by partitioning our catalog into an inner and outer disk sample, each containing an equal mass of molecular gas in catalog GMCs. The radius at which this occurs is 2.0 kpc. This coarse radial binning assures adequate statistics for $`M>1.5\times 10^5M_{\mathrm{}}`$ in two distinct cloud populations. Power law fits give $`\alpha =2.8\pm 0.3`$ for the inner disk and $`\alpha =2.4\pm 0.3`$ for the outer disk. A 2-sided K-S test shows that the two cloud samples are drawn from the same population with a 9% likelihood. Though there is some evidence for a difference in the mass distribution, it is only marginally significant.

## 5. Spatial Comparison of GMCs, H I Structure, and H II Regions

### 5.1. Radial Profiles and Star Formation Rates

For galactic radii $`R<7`$ kpc, we derived the azimuthally averaged mass surface density from our catalog GMCs. The results are presented in Figure 18, along with the surface density of H I mass and H$`\alpha `$ flux. Points represent averages within deprojected radial bins of $`\mathrm{\Delta }R`$ = 0.25 kpc. The H I data are from Corbelli & Salucci (2000) and the H$`\alpha `$ data are derived from the optical H II region catalogs compiled by Hodge, et al. (1999). All three components of the ISM are well fit by an exponential disk $`\mathrm{\Sigma }(R)=\mathrm{\Sigma }_0\mathrm{exp}(R/R_0)`$. Only the fit to the molecular data is shown (dot-dash line) in the figure. The molecular surface density $`\mathrm{\Sigma }_{\mathrm{H2}}`$ (squares) declines exponentially with $`R_0=1.4\pm 0.1`$ kpc, from a peak mass surface density $`2.15M_{\mathrm{}}\text{ pc}^2`$. Corbelli (2003) obtained a total molecular gas mass of $`2\times 10^8M_{}`$ from an FCRAO survey. Accounting for the difference in $`X`$ factors (Corbelli (2003) used $`2.8\times 10^{20}\text{cm}^2/(\text{K km s}^1)`$), her reported mass is 3 times the value we derive. The surface mass densities in the central region are comparable, but the scale length in Corbelli (2003) is 2.5 kpc, resulting in a significantly larger integrated mass. The reason for the difference in the derived scale lengths is unclear and merits further investigation. The H I surface brightness (dashed line) is nearly flat for $`R<6`$ kpc with a scale length $`10_2^{+8}`$ kpc and a peak surface density of $`9.5M_{\mathrm{}}\text{ pc}^2`$. H$`\alpha `$ surface brightness (solid line) has a scale length of $`1.7\pm 0.2`$ kpc and a peak surface brightness of $`1.2\times 10^{32}\text{ergs s}^1\text{pc}^2`$, corresponding to a star formation rate of 9.5 $`\times 10^{10}M_{\mathrm{}}\text{yr}^1\text{ pc}^2`$.

Kennicutt (1989,1998) has argued that there is a simple empirical relation between the SFR and the surface density of gas, $`\mathrm{\Sigma }_{gas}`$, which can be parameterized by a power law SFR $`(\mathrm{\Sigma }_{gas})^n`$. The SFR is most conveniently calculated by scaling the H$`\alpha `$ line emission, powered by OB stars. Similar scale lengths for H$`\alpha `$ and $`\mathrm{\Sigma }_{\mathrm{H2}}`$ in M33 imply $`n=0.9\pm 0.1`$, consistent with a constant SFR per unit molecular mass. Similarly, Wong & Blitz (2002) find for a sample of 7 normal spiral galaxies that $`\mathrm{\Sigma }_{\mathrm{SFR}}`$ strongly correlates with $`\mathrm{\Sigma }_{\mathrm{H2}}`$ only, where the power law index $`n`$ 0.8 — 1.4. Integrating $`L_{\mathrm{H}\alpha }`$ for the catalog of Hodge, et al. (1999) and using Equation 2 of Kennicutt (1998), we compute a star formation rate of 0.24 $`M_{\mathrm{}}\text{yr}^1`$ and a molecular gas depletion time of $`1.9\times 10^8`$ yr for M33. This value is a factor of 1.5 – 4 smaller than the depletion times quoted in Wilson, Scoville, & Rice (1991) who consider individual molecule-rich complexes as opposed to the galaxy as a whole. Our results, along with those of Wong & Blitz (2002), are inconsistent with those of Kennicutt (1989), who finds for disk-averaged observables that the SFR correlates more strongly with the surface density of either H I or H I$`+\mathrm{H}_2`$ than with that of molecular gas alone.

### 5.2. Circular and Radial Motions of Molecular and Atomic Hydrogen

We derived the rotation curve for GMCs cataloged independently of $`V_{\mathrm{HI}}`$ (catalog sources 1 – 93) and for the H I emission along the same line-of-sight from the Westerbork data cube (Deul & van der Hulst, 1987). The observed source positions $`(x,y)`$ were deprojected onto the plane of the galaxy assuming a P.A. of $`21^{}`$ and an inclination angle of $`51^{}`$. These values were varied slightly with galactic radius in accordance with the tilted ring model of Corbelli & Salucci (2000). The rotation center was assumed to coincide with the optical center. At the deprojected positions the rotational velocity

$$v_{rot}(R)=\frac{V_{los}(x,y)V_{sys}}{\mathrm{sin}i\mathrm{cos}\theta },$$
(7)
where $`\theta `$ is the azimuthal angle of the source measured in the plane of the galaxy from the the major axis. Sources within $`45^{}`$ of the minor axis were not included because of the large correction to ($`V_{los}(x,y)V_{sys}`$).

Figure 19 shows the circular velocities for catalog GMCs, averaged in 250 pc bins. The median CO rotation curve (open squares) closely matches the H I rotation curve determined by Corbelli & Salucci (2000) from Westerbork data. The dispersion of GMCs about the H I rotation curve varies between 7.5 and 9.0 km s-1 for $`R<4`$ kpc. There is no apparent trend in the dispersion with galactic radius. These values are comparable to Milky Way values after streaming in the Galaxy is taken into account.

Figure 19 also shows the CO radial velocities

$$v_{rad}(R)=\frac{V_{los}(x,y)V_{sys}}{\mathrm{sin}i\mathrm{sin}\theta },$$
(8)
including only sources within $`\mathrm{\Delta }\theta <30^{}`$ of the minor axis. We find little evidence for large scale radial (non-circular) motions in M33 in agreement with Deul & van der Hulst (1987).

### 5.3. Spatial Comparison of GMC Distribution with H I Filaments

An extraordinary spatial correspondence between GMCs and the distribution of atomic hydrogen is shown in Figure 20. The halftone image is the H I brightness distribution, derived by integrating the Westerbork data cube (Deul & van der Hulst, 1987); the resolution of the Westerbork map is $`24^{\prime \prime }\times 24^{\prime \prime }`$, comparable to the BIMA data. The H I distribution is dominated by bright, narrow filaments in an extended disk. The dark circles mark the positions of the catalog GMCs; cloud masses are proportional to circle areas. Figure 21 shows that nearly all GMCs are found in regions of H I overdensity, but there is no correlation between the magnitude of the overdensity and the cloud mass. Apparently, the amount of ambient H I is independent of the mass of the GMC that forms from it. The correspondence between GMCs and filaments persists down to smaller scales in the central region where the H I brightness map has less dynamic range. The association with the filaments is even more pronounced in the high resolution H I maps of M33 made at the VLA (Thilker, private communication). The average H I mass located within 150 pc of our catalog sources is $`5\times 10^5M_{\mathrm{}}`$, indicating an ample local supply of atomic hydrogen from which observed GMCs could have condensed.

In Figure 22, we see that GMCs form from H I with an column density threshold of $`5\times 10^{20}\text{cm}^2`$, which is comparable to the mean H I column density in the H$`\alpha `$ disk (dashed histogram). The shaded histogram indicates pixels with associated CO emission. The solid line shows all WSRT pixels Deul & van der Hulst (1987). The threshold value is similar to that found by Savage, et al. (1977) and Blitz, Bazell & Desert (1990) in the solar vicinity, and is indicated in Figure 22. Moreover, the H I shows an upper limit in column density of $`3\times 10^{21}\text{ cm}^2,A_V3`$, often at positions associated with molecular gas. This suggests that the atomic gas is providing shielding for the molecular gas, but higher column densities result in more gas in the molecular state.

The close association between GMCs and H I filaments can be used to set an upper limit on the GMC lifetimes. The RMS difference between the CO and H I velocity centroids is 8 km s-1, calculated from clouds 1–93, which were selected independent of the H I velocity. The typical width of an H I filament in the Deul & van der Hulst (1987) map is $`200`$ pc. If cloud lifetimes were significantly longer than 10–20 Myr, clouds would drift off the filaments, destroying the observed association.

### 5.4. Implications for GMC Formation

Elmegreen (1990) reviews many theories of GMC formation, including the agglomeration of small molecular clouds (Scoville & Hersh, 1979). However, there can be little doubt that the GMCs in M33 form out of the atomic gas. If instead the atomic gas resulted from the dissociation of molecular gas, the H I would appear as discrete envelopes around the GMCs rather than as filaments, and there would little atomic gas beyond the apparent edge of the molecular disk at $`R=4`$ kpc. Neither of these features are seen. In addition, there is significantly more atomic gas ($``$ a factor of 5) in the filaments than can be associated with discrete envelopes around GMCs and $`<15\%`$ of this mass can be attributed to photodissociation (Wilson & Scoville, 1991). Furthermore, if GMCs formed by the collisional agglomeration of smaller molecular clouds, the growth time would be comparable to the collision time. This is roughly equal to the collision time of the catalog clouds, which is $`2\times 10^9`$ years based on 140 GMCs inside a radius of 4 kpc, an H2 half thickness of 100 pc and a cloud-cloud velocity dispersion of 10 km s-1. Since the number density of clouds scales roughly as $`1/M_{cloud}`$ and the geometric cross section scales as $`M_{cloud}`$ (owing to GMCs having a constant surface density, see Paper II), this collision time is independent of mass and characterizes the growth time for all self-gravitating clouds. Gravitation effects do not significantly enhance the collision cross section (Scoville & Hersh, 1979; Blitz & Shu, 1980). Additionally, Paper II sets an upper limit for the diffuse H2 surface density of 0.3 M pc-2. Forming a $`10^5`$ $`M_{\mathrm{}}`$ GMC from a diffuse molecular component at this surface density would require almost 108 yr to accumulate material from a 300 pc radius region, provided that the gas velocity dispersion could be converted into ordered motion. Both times are considerably longer than the lifetimes of the GMCs ($`2\times 10^7`$ yr) based on their association with the H I filaments.

It appears that the process of GMC formation in M33 is first and foremost a process of H I filament formation. While a high H I column density is a necessary condition for the formation of a GMC, it clearly is not a sufficient condition, since there are many positions with high H I column densities but no observed CO emission (Figure 22). The filamentary structures extend out to the $`R=8`$ kpc edge of the H I map, well beyond the radius where GMCs are observed, with no significant change in character at the edge of the molecular disk. It appears that some physical parameter, which evidently is a function of galactic radius, determines what fraction of the atomic gas in a filament is converted into molecules. A likely candidate is hydrostatic pressure in the disk, a hypothesis that will be investigated in a subsequent paper.

### 5.5. Are GMCs and H I Holes Physically Linked?

Expanding shells in the interstellar medium may play an important role in the dynamics of the neutral gas, transmitting both mechanical and radiant energy at their boundaries (McKee & Ostriker, 1977). Deul & den Hartog (1990, DdH) cataloged 148 holes (underdensities) in the H I distribution of M33. Ninety-three of the holes (Types 2 and 3 in DdH) show the kinematic evidence of expansion in the Westerbork channel maps. These holes cover $`<20\%`$ of the neutral disk area. DdH find that many of the compact holes have OB associations inside and conclude that many such holes have been excavated by massive star formation activity. While DdH also allow for the possibility that holes are formed by the local collapse of H I to form GMCs, our CO survey shows that this is unlikely, since nearly all GMCs are outside the holes. In this section we consider evidence of a physical connection between GMCs and the H I shells.

There is, in fact, a tendency for the catalog GMCs to be spatially clustered around the DdH H I holes. We quantified this effect by measuring the fraction of GMCs located within annuli of width $`\mathrm{\Delta }R`$ bounding the holes, and comparing this to the fraction expected if the GMCs and holes are uncorrelated. We restricted the analysis to holes with diameters less than 400 pc because these are more likely to have been created by high mass star formation. The results are shown in Figure 23. We find a significant correlation between GMCs and H I hole edges for $`\mathrm{\Delta }R`$ between 60 and 300 pc, most significantly for $`\mathrm{\Delta }R`$ = 150 pc (binomial probability $`P_{rand}10^3`$). The large separation (150 pc) between hole edges and GMCs suggests, however, that expanding shells do not directly trigger formation of the GMCs. Figure 24 plots the GMCs (boundaries defined by the 2$`\sigma `$ clipping levels) and the 150 pc annuli (shaded) around the (compact) hole edges. This figure shows that there is no significant CO emission inside the DdH holes.

Since essentially all GMCs are on filaments but not all GMCs are at the edges of holes, we conclude that the filamentary structure is more fundamental to GMC formation than the action of expanding shells. These shells may create filamentary structure on small scales, especially in the presence of turbulence (e.g. Hartmann, Ballesteros-Paredes, & Bergin, 2001), but are unlikely to produce the largest filaments. The energies required to evacuate the large holes defining these filaments are $`10^{53}`$ ergs (see e.g. DdH). Such energies would leave bright stellar remnants at the centers of the large holes, which are not generally observed. In addition, these holes would be sheared apart on time scales of $`2\times 10^7`$ years by differential galactic rotation unless they were stabilized by some process. Thus some organizing principle other than stellar feedback seems to be responsible for much of the filamentary H I structure, and ultimately, for GMC formation in M33. One possibility is MHD instabilities in a self-gravitating, shearing disk (Elmegreen, 1987; Kim, Ostriker, & Stone, 2002).

Although nearly all GMCs are outside holes on filaments of H I, an exception is DdH 28 (Figure 24). The catalog GMCs in this hole may result from a high velocity cloud (HVC) colliding with the H I disk. This scenario for hole creation was first elaborated by Tenorio-Tagle (1981). GMCs 66 and 88 have velocities blue-shifted relative to local H I by 72 and 24 km s-1, respectively. The formation of GMCs 40, 51, 54, 123, and 133 at the edge of the hole – all blueshifted relative to local H I– may have been triggered by this collision as well. The combined kinetic energy of these clouds ($`1.4\times 10^{52}`$ ergs) could have been supplied by an atomic HVC. Alternatively, GMCs in DdH 28 may have formed on the front surface of neutral shells. That is, they may be similar to the Orion A and B molecular clouds, which appear to have formed from gas lifted 140 pc out of the plane of the Milky Way by stellar winds (Dame, Hartmann, & Thaddeus, 2001; Hartmann, Ballesteros-Paredes, & Bergin, 2001).

### 5.6. Spatial Comparison of GMC Distribution with H II Regions

Figure 25 compares the locations of GMCs with H$`\alpha `$ emission, showing that some GMCs follow the weak H$`\alpha `$ spiral arms of M33. Molecular clouds in the overlay are represented by filled circles with area proportional to mass; the underlying halftone is a CCD image from Cheng et al. (1996). Most of the H$`\alpha `$ flux comes from H II regions, which are cataloged in Boulesteix et al. (1974), Wyder, Hodge, & Skelton (1997), and Hodge, et al. (1999). As shown in §5.1, the radial profiles of CO flux and H II region flux have comparable scale lengths; however, not all H II regions are associated with GMCs. In fact, there are over 20 times as many H II regions as catalog GMCs, and many of the H II regions are far from detected clouds. Since our GMC catalog has a lower dynamic range than the H II region catalog, it’s possible that low mass GMCs are associated with most of these “orphan” H II regions; §4.1 argues that there could be $``$ 2000 molecular clouds in M33, comparable to the number of H II regions. Alternatively, H II regions may rapidly dissipate molecular clouds, accounting for the discrepancy in numbers (Hartmann, Ballesteros-Paredes, & Bergin, 2001).

We calculate the fraction of the H$`\alpha `$ flux that originates from H II regions that spatially overlap catalog GMCs. Using the positions and radii of H II regions from the combined catalogs (Boulesteix et al., 1974; Wyder, Hodge, & Skelton, 1997; Hodge, et al., 1999), we find that 40% of the H$`\alpha `$ flux originates from H II regions which are tangent to, or contain, catalog GMCs. If the remaining $`2.2\times 10^7M_{}`$ of molecular gas (§3.4) produced a proportional amount of H$`\alpha `$ flux, at least 80% of the emission would arise from H II regions associated with GMCs.

One may also ask what fraction of the GMCs are actively forming stars. Figure 26 shows the positional correlation of catalog GMCs with H II regions (circles). We counted the fraction of GMCs that have at least one H II region within distance $`\mathrm{\Delta }r`$. The thin curve is the fraction expected from random association (see Appendix B). There is a statistically significant clustering of GMCs and H II regions out to a separation of 150 pc. As many as 100 GMCs (67%) are apparently forming massive stars – that is, their centroid position is within 50 pc of some H II region. This is a similar fraction of association as found in a more limited study by Wilson & Scoville (1991). Of the remaining GMCs, we estimate that only a few contain obscured H$`\alpha `$ sources. From our correlation curve, there are only 7 H II regions coincident with a GMC position ($`\mathrm{\Delta }r10`$ pc). By symmetry, we expect that only $`7`$ catalog GMCs are obscuring signs of massive star formation.

Winds and ionizing radiation from massive stars destroy molecular clouds. Since at least 2/3 of our sources are associated with H II regions, we infer that GMCs spend less than 1/3 of their life in a quiescent phase prior to the onset of massive star formation. The fraction of time spent in the quiescent phase appears to be even smaller for the more massive clouds. Figure 26 shows that 85% of the higher mass GMCs are within 50 pc of H II regions versus 55% of the lower mass GMCs. A reciprocal correlation exists between luminous H II regions and GMCs: luminous H II regions (top 10%) are nearly twice as likely to have an associated catalog GMC than the population as a whole.

Between the cutoff of the molecular disk ($`R=4.5`$ kpc, §5.1) and the edge of the star forming disk ($`R=6.7`$ kpc), there are 7 cataloged GMCs and 600 H II regions. Comparing the total H$`\alpha `$ luminosity in this region to the total cataloged molecular mass, we find the ratio $`L_{\mathrm{H}\alpha }/M_{\mathrm{H2}}`$ for this region is 8 times larger than in the inner galaxy ($`R<4.5`$ kpc). Correcting for different sensitivities and survey coverage in the two regions can only account for a factor of two discrepancy. Based on the H$`\alpha `$ flux, we would expect to have cataloged at least 4 times as much molecular mass in the outer galaxy as we actually do. It is unlikely that a burst of star formation in this region has depleted the supply of GMCs. We conclude either that dispersal times for GMCs are significantly shorter in the outer galaxy, that the mass spectrum of the GMCs becomes even steeper there, or the metallicity significantly affects the CO-to-H2 conversion factor in this region.

## 6. Summary

In this paper we present the results of an unbiased interferometric survey of the star forming disk of M33 in CO(1$``$0); the observations were done using the BIMA array. The 50 pc linear resolution of our survey map is comparable to the size of most GMCs. We derive a catalog of GMCs for M33 complete down to 1.5 $`\times 10^5M_{\mathrm{}}`$. From simple statistics, we expect that no more than 15 of the 148 sources listed in our catalog are spurious. We estimate that approximately 60 clouds with a total mass of $`4.3\times 10^6M_{}`$ have been falsely rejected from the catalog. The interferometer data were compared with single dish fields observed at the UASO 12 m telescope to estimate the effects of spatial filtering and masking on the catalog. Our results generally agree with those of Wilson & Scoville (1990), but the greater spatial coverage and flux completeness of the BIMA survey leads, in some cases, to different physical inferences.

From the survey data, we conclude the following:

1. The total mass of GMCs in our catalog is $`2.3\times 10^7M_{}`$. From a comparison of the interferometer and single dish data in selected fields, we estimate that the total molecular mass of the galaxy is $`4.5\times 10^7M_{}`$, 2% of the atomic mass.

2. GMCs with masses above our completeness limit are described by a power-law distribution with $`dN/dMM^{2.6\pm 0.3}`$. The power-law index of GMCs in M33 is steeper than that in the Milky Way ($`dN/dMM^{1.6}`$). We infer that there exists a low mass cutoff or a change in the index to $`<`$ 2 between $`3\times 10^4M_{}`$ and $`10^5M_{}`$. The cutoff or change in slope implies that GMCs in M33 form with a characteristic mass of $`7\times 10^4M_{}`$.

3. The surface density of molecular gas decreases exponentially in radius, with a scale length of $`1.4\pm 0.1`$ kpc. The scale length agrees well with the surface brightness of H$`\alpha `$ emission ($`1.7\pm 0.2`$ kpc) implying a scaling between star formation and molecular gas surface density of SFR $`\mathrm{\Sigma }_{\mathrm{CO}}^{0.9\pm 0.1}`$. Based on the H$`\alpha `$ luminosity of M33, we estimate a star formation rate of 0.24 $`M_{}`$ yr-1, implying a molecular gas depletion time of $`1.9\times 10^8`$ yr. This very short time is not unreasonable given the large reservoir of atomic gas available at all radii.

4. The rotation curve of the galaxy derived from CO emission is in excellent agreement with that of the H I. We find no evidence for large scale radial motions of molecular gas in the galaxy.

5. Giant molecular clouds are preferentially found on bright H I filaments where the H I column is approximately $`10^{21}`$ cm-2. Regions with H I column densities $`<5\times 10^{20}`$ cm-2 are devoid of catalog GMCs.

6. We estimate a GMC lifetime of 10–20 Myr based on the close association between GMCs and H I filaments. Given the typical velocity difference of 8 km s-1, the clouds would drift away from the filaments on longer time scales.

7. Since the filamentary structure of the H I gas appears to be independent of the existence of nearby GMCs, we conclude that it is the H I filamentary structure that forms first. The fraction of gas that becomes molecular is then determined by some other parameter, such as hydrostatic pressure, that is inversely correlated with galactic radius.

8. At least 40% of the H$`\alpha `$ flux is associated with catalog GMCs. If the undetected molecular mass produces a proportional amount of H$`\alpha `$ flux, at least 80% of the emission would arise from H II regions associated with GMCs. At least 2/3 of catalog GMCs are within 50 pc of an H II region, so it appears that GMCs spend less than 1/3 of their lifetime in a quiescent phase prior to the onset of star formation. Clouds above the median mass spend less than 15% of their lifetime in the quiescent phase. H II regions above the 90th percentile in luminosity are nearly twice as likely to have a GMC within 50 pc as compared to the entire population.

9. Nearly all catalog GMCs are exterior to H I holes. Clearly, H I holes do not result from conversion of atomic to molecular gas. GMCs and compact H I holes ( $`d<400`$ pc) are clustered over separation scales of 60 — 300 pc. However, since many GMCs are not associated with holes, the holes probably play only a secondary role in GMC formation.

10. GMCs 66 & 88 are the only GMCs interior to an H I hole (DdH 28); they are highly blueshifted with respect to H I. GMCs 40, 51, 54, 123 & 133, located along the edge of DdH 28, are also blueshifted with respect to H I. Conceivably these clouds were formed by the impact of a high velocity cloud on the atomic disk.

11. Although the molecular disk shows a sharp decline in the number of GMCs beyond $`R4.5`$ kpc, the H$`\alpha `$ disk extends to $`R6.7`$ kpc. This suggests that, for the outer disk, either dispersal times of GMCs are significantly shorter than for the inner disk or that there exists a large population of GMCs below our detection threshold.

This work was partially supported by NSF grant AST-9981308 to the University of California. This research made extensive use of NASA’s Astrophysics Data System (ADS) and the NASA/IPAC Extragalactic Database (NED). ER’s work is supported in part by a NSF Graduate Fellowship.

## Appendix A Probability calculation

This appendix provides a brief description of the calculation used to rank the likelihood of candidate molecular clouds in the data cube.

As discussed in the text, the flux density $`S`$ at each pixel in the data cube is first normalized by the RMS noise $`\sigma `$, as a function of position and velocity, to produce a signal-to-noise, or significance, image $`y`$:

$$s(x,y,v)=\frac{S(x,y,v)}{\sigma (x,y)\sigma (v)}.$$
(A1)
Suppose that at some location in the data cube we have identified a spectral feature, occupying $`n`$ adjacent velocity channels, which we suspect may be CO emission from a molecular cloud. What is the likelihood that this is a false detection?

Assume that the signal-to-noise ratios of the $`n`$ channels have been sorted in ascending order such that $`s_1s_2\mathrm{}s_n`$. We wish to calculate the probability of randomly drawing $`n`$ values from a normalized Gaussian distribution (with zero mean and standard deviation 1), such that all $`n`$ values are $`s_1`$, at least $`(n1)`$ values are $`s_2`$, …, and at least 1 value is $`s_n`$ – that is, the probability of selecting $`n`$ values which are at least as unlikely as the actual spectrum.

Suppose that one draws a single value $`s`$ from the Gaussian distribution. There are $`(n+1)`$ mutually exclusive outcomes:

$`\begin{array}{cc}s<s_1,\hfill & \mathrm{with}\mathrm{probability}p_0=0.5\{1\mathrm{erfc}(s_1/\sqrt{2})\}\hfill \\ \mathrm{}\hfill & \\ ss_n,\hfill & \mathrm{with}\mathrm{probability}p_n=0.5\{\mathrm{erfc}(s_n/\sqrt{2})\}\hfill \end{array}`$

where

$$\mathrm{erfc}(x)\frac{2}{\sqrt{\pi }}_x^{\mathrm{}}e^{t^2}𝑑t.$$
(A2)
Then, if one draws a total of $`n`$ values from the distribution, the probability that $`k_0`$ have values $`s<s_1`$, $`k_1`$ have values in the interval $`s_1s<s_2`$, …, and that $`s_n`$ have values $`sy_n`$, is given by the multinomial distribution:

$$p(k_0,k_1,\mathrm{},k_n)=\frac{n!}{k_0!k_1!\mathrm{}k_n!}p_0^{k_0}p_1^{k_1}\mathrm{}p_n^{k_n}.$$
(A3)
The probability of a false detection is the sum of all the multinomial terms which fulfill the original requirement, namely that $`k_0=0`$, $`k_n1`$, $`(k_{n1}+k_n)2`$, and so on:

$$P_{false}=p(0,0,\mathrm{},n)+p(0,0,\mathrm{},1,n1)+\mathrm{}+p(0,1,\mathrm{},1).$$

The M33 data cube contains approximately 20500 independent positions, each of which corresponds to a 220-channel spectrum. For an $`N_{ch}`$-channel spectrum, there are $`(N_{ch}n+1)`$ ways of selecting $`n`$ adjacent channels. Thus, over the 94 channel range covering the velocity interval $`290<`$ $`V_{\mathrm{LSR}}`$$`<80`$ there are of order $`N_{trials}=20500(94n)2\times 10^6`$ ways of selecting $`n`$ adjacent channels from the data cube. If $`P_{false}1`$ and $`N_{trials}1`$, the probability of $`r`$ false detections is given by Poisson statistics:

$$P(r)=\frac{(NP)^r}{r!}e^{NP}.$$
(A4)
The probability of any given detection being real is $`P_{real}=P(r=0)=e^{NP}1NP`$ for $`NP1`$. If there are a given number of detections at this significance level, the probability of $`q`$ of them being false is simply $`P(q)`$.

As an example, suppose we detect 3 adjacent channels with $`s_1=s_2=s_3=2.5\sigma `$. In this case the probability of a false detection is simply the product of the probabilities of a false detection in each of the 3 channels: $`P_{false}=\{0.5\mathrm{erfc}(2.5/\sqrt{2})\}^3=2.39\times 10^7`$. Using the Poisson formula, one finds a 60% likelihood of all sources at this significance being real.

If the signal-to-noise ratios of the 3 channels are unequal, then the probability of a false detection is greater than the simple product of the 3 probabilities, essentially because one has more than one chance to draw the less likely values from the population of Gaussian deviates. For example, if $`s_1=2.0\sigma ,s_2=2.5\sigma `$, and $`s_3=3.4\sigma `$,

$$P_{false}=p(0,0,0,3)+p(0,0,1,2)+p(0,0,2,1)+p(0,1,0,2)+p(0,1,1,1)=2.39\times 10^7,$$
(A5)
whereas the simple product of the probabilities is $`4.76\times 10^8`$. Again, Poisson statistics calculates about a 60% likelihood that all sources at this level are real. Because of the non-Gaussian nature of the noise, we ultimately use this probability as a relative measure for ranking the significance of cloud candidates (§3.1).

## Appendix B GMC–H II Region Spatial Correlation Function

To understand in a statistical sense how massive star formation affects GMCs, it is a matter of interest to determine how many GMCs in our catalog are physically linked with H$`\alpha `$ emitting H II regions. The filling factors $`f`$ for H$`\alpha `$ and CO in the survey region are 2.4% and 0.41%, respectively, where both types of emission show strong positional correspondence with H I filaments. H II regions are so widespread that, to some degree, observed proximity to GMCs may be random. If H II regions and GMCs are all randomly distributed in a survey area $`A_{survey}`$, the average separation $`r`$ between a GMC and the nearest H II region obeys the relation

$$\frac{\pi r^2f_{\text{H}\alpha }}{A_{\mathrm{HII}}}1,$$
(B1)
where $`A_{\mathrm{HII}}`$, the average area of an optical H II region, is $`700\text{pc}^2`$. If we assume an H II region and GMC are closely associated if they are separated by no more than $`r_{}`$, say, then the expectation value for the number of close associations is

$$N_{exp}=N_{gmc}\left(\frac{r_{}}{r}\right)^2=\frac{\pi r_{}^2A_{survey}f_{\text{H}\alpha }f_{CO}}{A_{CO}A_{\mathrm{HII}}}$$
(B2)
Assuming that $`A_{CO}`$ has approximately the same area as a resolution element of our survey, $`A_{\mathrm{HII}}/A_{CO}0.4`$, and $`r_{}40`$ pc – the average projected separation of an H II region lying at the edge of a GMC – then $`N_{exp}34`$. A comparison of our GMC catalog to the H II region catalogs previously referenced yields 100 H II regions in close association with GMCs, about a factor of three in excess of what might be expected if GMC and H II positions were random.


format: text

text:
# 1 Spanning probability function for a fixed lattice size 𝐿 (solid line). For larger and larger sizes (dotted lines), this function approaches a step. By fixing some value 𝑟 at the vertical axis, one can find a sequence of values 𝑝_𝐿(𝑟) at the horizontal axis approaching the critical threshold 𝑝_𝑐, for increasing lattice sizes.

Are the Tails of Percolation Thresholds Gaussians ?

P.M.C. de Oliveira1,2, R.A. Nóbrega2 and D. Stauffer1,3.

1) Laboratoire de Physique et Mécanique des Milieux Hétérogènes

École Supérieure de Physique et de Chimie Industrielles

10, rue Vauquelin, 75231 Paris Cedex 05, France

2) Instituto de Física

Universidade Federal Fluminense

av. Litorânea s/n, Boa Viagem, Niterói, Brasil 24210-340

3) Institute for Theoretical Physics

Cologne University

D-50923 Köln, Euroland

e-mails: pmco@if.uff.br; rafaella@if.uff.br; stauffer@thp.uni-koeln.de.

PACS numbers: 02.70.-c, 05.10.Ln, 64.60.Ak, 05.70.Jk

Abstract: The probability distribution of percolation thresholds in finite lattices were first believed to follow a normal Gaussian behaviour. With increasing computer power and more efficient simulational techniques, this belief turned to a stretched exponential behaviour, instead. Here, based on a further improvement of Monte Carlo data, we show evidences that this question is not yet answered at all.

In reference , the percolation on a $`N`$-site square lattice is treated with high numerical accuracy. Indeed, the best known estimate for the critical threshold, $`p_c=0.59274621(13)`$, comes from this work. In order to study this kind of problems, the authors follow a very fruitful Monte Carlo approach which allows one to obtain continuous functions of $`p`$, the concentration of occupied sites, namely the canonical-like average

$$R(p)=\underset{n}{}C_N^np^n(1p)^{Nn}R_n,$$
(1)
of some quantity $`R`$. Here, $`R_n`$ is a uniform average over all configurations with just $`n`$ occupied sites, i.e. a microcanonical-like average. $`C_N^n`$ are simple binomial factors. By filling the initially empty lattice, site by site at random, and repeating this process many times, one is able to get the discrete set of microcanonical averages $`R_n`$ accumulated into an $`n`$-histogram, over the entire range, $`n`$ = 0, 1, 2 $`\mathrm{}`$ $`N`$. From this set of numbers, the determination of the continuous $`p`$-function $`R(p)`$ is straightforward.

In particular the authors of  fix attention at the horizontal wrapping probability $`R_L(p)`$ around a $`L\times L`$ torus, i.e. a square lattice with periodic boundary conditions on both directions. In the thermodynamic limit, this function approaches a step: $`R_{\mathrm{}}(p)=0`$ below the critical threshold $`p_c`$, and $`R_{\mathrm{}}(p)=1`$ above $`p_c`$. For finite sizes, $`R_L(p)`$ presents a sigmoid aspect similar to figure 1. A good approach to $`p_c`$ is to choose some fixed value $`r`$, and solve the equation $`R_L(p)=r`$, getting the root $`p`$ shown at the horizontal axis. Here, one can appreciate the advantage of knowing $`R_L(p)`$ as a continuous function of $`p`$. Keeping the same value $`r`$ and repeating this task for a series of increasing lattice sizes (dotted lines), one gets a series of roots $`p_{L_1}(r)`$, $`p_{L_2}(r)`$, $`p_{L_3}(r)`$ $`\mathrm{}`$ which converges to the desired threshold $`p_c`$ in the thermodynamic limit.

The above reasoning is valid no matter which is the fixed value for $`r`$ one chooses. However, for the very particular choice $`r^{}=0.521058290`$, a universal probability exactly known through conformal invariance arguments , the convergence becomes fast, i.e. the root $`p(L)`$ differs from $`p_c`$ as $`L^{21/\nu }=L^{2.75}`$, where $`\nu =4/3`$ is the correlation length critical exponent. The above quoted accurate value for $`p_c`$ was obtained in this way. For details, see  and references therein.

Reference  proposes the mathematical form

$$p_L(r)=p_c+\frac{1}{L^{1/\nu }}\left[A_0(r)+\frac{A_1(r)}{L}+\frac{A_2(r)}{L^2}+\mathrm{}\right],$$
(2)
for estimators $`p_L`$ obtained from quantities like $`R(p)`$. The option for the wrapping probabilities around the torus and the convenient choice of Pinson’s number $`r=r^{}`$ lead to vanishing values for the two first terms $`A_0(r^{})=A_1(r^{})=0`$, a lucky coincidence which accelerates very much the convergence.

The quantity $`R_L(p)`$ is obtained, as quoted before, by filling up the initially empty lattice site by site, at random. Clusters of neighbouring occupied sites grow. As soon as the horizontal wrapping along the torus is set, one books the corresponding value of $`n`$, the number of occupied sites so far, and stops the process. For that particular sample, the wrapping probability $`R_L`$ is a step function, i.e. $`R_L=0`$ below $`n`$ and $`R_L=1`$ above. The same routine is repeated many times, in order to have a probability distribution for $`n`$. The various step functions are then superimposed to get the microcanonical averages $`R_n`$ in equation (1), stored in an $`n`$-histogram. Finally, the continuous canonical average $`R_L(p)`$ can be calculated for any value of $`p`$.

Each process of filling-up the lattice (one sample) yields a single value $`n`$ for the statistics, i.e. just one more entry on the $`n`$-histogram. In , we decided to improve this point, by changing the definition from wrapping to spanning probability, figure 1. We fix two parallel horizontal lines on the $`L\times L`$ torus, separated by a distance of $`L/2`$, for instance lines $`i=1`$ and $`i=1+L/2`$. The measured quantity is now the probability of having these two lines connected by the same cluster of neighbouring occupied sites, instead of the wrapping probability along the whole torus. The advantage is that we can measure the same thing for lines $`i=2`$ and $`i=2+L/2`$, also for lines $`i=3`$ and $`i=3+L/2`$, and so on. Moreover, also vertical parallel lines can be included into this counting. At the end, from a single sample we store $`L`$ new entries into our $`n`$-histogram, instead of just one more entry. Note that this advantage even increases for larger and larger lattices.

Within the same computational effort, our approach allows the test of larger lattices. Because of that, we were able to confirm the validity of equation (2) with high precision, by sampling 27 different lattice sizes from $`L=18`$ up to $`L=1594`$, a 8-thousand factor in number of sites. On the other hand, our definition does not allow the chance of vanishing both $`A_0(r)`$ and $`A_1(r)`$ at once. We can have only $`A_0(\widehat{r})=0`$, for a particular universal probability $`\widehat{r}=0.984786(11)`$ numerically determined within the same work . Independently, Cardy  tried to determine it by conformal invariance arguments, however, in looking for configurations which link two parallel lines, he was forced to disregard configurations which also wrap along the other direction. As a result of using larger lattices but a slower convergence rate of $`L^{11/\nu }=L^{1.75}`$, we get the same figure $`p_c=0.59274621(33)`$  as in , but within a 3 times larger error bar.

The non-Gaussian behaviour of the finite-lattice-threshold distribution near the infinite-lattice critical point is already stablished . Here, we profit from the same simulational data in order to investigate the distribution tails, far from the critical point. Which is the mathematical form of the tail observed in figure 1, below the root $`p`$? One possible answer is a simple Gaussian form

$$R_L(p)\mathrm{exp}[K(pp_c)^2].$$
(3)
Another alternative is a stretched exponential

$$R_L(p)=\mathrm{exp}[C(p_cp)^\nu ],$$
(4)
where the strict equality (for large $`L`$ and far from $`p_c`$) is a consequence of the periodic boundary condition , which holds for our data. Profiting from this strict equality, one can test equation (4) by constructing a plot of $`\mathrm{ln}[\mathrm{ln}(R)]`$ against $`\mathrm{ln}(p_cp)`$. This was done in , and we repeat the same for our data, in figure 2. Note that our range for $`\mathrm{ln}(p_cp)`$ (up to $`3.4`$ for $`L10^3`$) is larger than in reference  (up to $`4.6`$ for the same size). This means that we are testing more deeply the distribution tails, thanks to our trick of sampling $`L`$ new entries for each run. Even so, the conclusion in favour of either equation (3) or (4) is by no means obvious. Note a further difficulty in what concerns equation (3), because the leading multiplicative constant in front of the exponential is not necessarily 1.

Another, perhaps better way to address the same question is by plotting $`\mathrm{ln}(R)`$ twice, against $`(pp_c)^2`$ and $`(p_cp)^{4/3}`$. Figure 3 shows the result for our data. Note that our range for $`\mathrm{ln}(R)`$ (down to $`16`$) doubles the one presented in . The would-be Gaussian case (up) presents clear positive curvatures, whereas the would-be stretched exponential (down) presents negative curvatures although not so pronounced. The exponent $`4/3`$ seems to fit better, but one cannot extract a clear conclusion from these data.

Still more undefined is the situation of the other tails on the right of figure 1, above $`p_c`$. In this case (not shown), our accuracy limit for $`\mathrm{ln}(1R)`$ (down to $`16`$) is reached much closer to $`p_c`$ than the case shown in figures 2 and 3, below $`p_c`$.

Concluding, we present new Monte Carlo data concerning the distribution probability of percolation thresholds on a finite square lattice. We address the question of the mathematical form of the distribution tails, equation (3) against (4). Even considering that our statistics is over than 1000 times larger than previous works, no definitive conclusion can be extracted from our data, in what concerns the asymptotic tail exponent.

Aknowledgements: This work is partially supported by Brazilian agencies FAPERJ and CNPq.


text:
# Constraints on intergalactic dust from quasar colours

## 1 Introduction

Scattering and absorption of light by a homogeneous distribution of intergalactic large grains giving a nearly achromatic cross-section – so called grey dust – has been proposed  as an alternative to the dark energy explanation for the observed faintness of Type Ia supernovae (SNIa) at $`z0.5`$

The presence of a significant component of dark energy has been inferred from complementary cosmological tests, e.g., from recent WMAP background radiation measurements  in combination with limits on the Hubble parameter from the HST Key Project  and large scale structure measurements from the 2dFGRS survey . However, in Ref.  it was demonstrated how the recent CMB anisotropy observations can be explained also by an Einstein-de Sitter cosmology with no dark energy if one allows for a low Hubble parameter, $`H_050`$ km/s/Mpc, the value currently favoured by gravitational lensing time delay measurements . It is further argued that conflict with large scale structure measurements can be avoided by suppressing cluster scale mass fluctuations assuming a small but non-negligible hot dark matter component.

Thus, probing the redshift-distance relation through Type Ia SN observations remains essential for establishing and exploring the dark energy of the universe. Ruling out alternative explanations to the dark energy interpretation and minimizing systematic effects on the measured distances to high-$`z`$ supernovae is still crucial. In this note, we discuss the possibility to constrain the effects of intergalactic grey dust.

It has been argued that the amount of intergalactic dust required to make the SNIa results compatible with an Einstein-de Sitter universe is disfavored by the far-infrared background measured by the DIRBE/FIRAS instruments . However, the derived bounds rely on the assumption that optical and UV photons are absorbed and re-emitted rather than just scattered by the dust grains.

It is expected that large dust grains will also cause some differential extinction for high-$`z`$ sources (although less than for smaller dust grains), allowing the possibility to distinguish between extinction and a cosmological origin of the dimming. In Ref. , it was shown that with 1 % relative spectrophotometric accuracy or broadband photometry in the wavelength interval 0.7–1.5 $`\mu `$m for high-$`z`$ SNIa observations, one could reduce the uncertainty from grey dust extinction down to $`\mathrm{\Delta }m=0.02`$ magnitudes (the target error for the proposed SNAP SN survey ).

In this paper, the differential extinction for high-$`z`$ quasars (QSOs) caused by extragalactic dust along the line of sight is used to put an upper limit on the possible effects on SNIa observations at $`z=1`$ of $`\mathrm{\Delta }m=0.2`$ (comparable to current SNIa survey errors) for a wide range of models for evolution of dust density and extinction properties.

## 2 Extinction by dust at cosmological distances

For a given emission redshift $`z_e`$, the attenuation $`\mathrm{\Delta }m_\mathrm{d}`$ at observed wavelength $`\lambda _o`$ due to dust can be written (see also Ref. )

$$\mathrm{\Delta }m_\mathrm{d}(z_e,\lambda _o)=\frac{2.5}{\mathrm{ln}10}_0^{z_e}a[\lambda _o/(1+z),R_V]h(z)/D_V(z)𝑑z,$$
(1)
where $`a(\lambda ,R_V)`$ is the dimensionless wavelength-dependent attenuation coefficient where we use the parameterization given in  and $`D_V(z)[\sigma \rho _{\mathrm{dust}}(z)]^1`$ is the interaction length for scattering or absorption by dust particles in the observed $`V`$-band ($`\lambda =0.55\mu `$m). Here, $`\rho _{\mathrm{dust}}(z)`$ is the physical dust density at redshift $`z`$, $`\sigma `$ the interaction cross-section and the proportionality constant will be related to the mass-distribution of the dust particles (since $`D_V(z)=[\sigma n_{\mathrm{dust}}(z)]^1`$ where $`n_{\mathrm{dust}}(z)`$ is the dust number density). The cosmology-dependent function $`h(z)`$ is given by

$$h(z)=\frac{1}{H_0\left(1+z\right)\sqrt{(1+z)^2(1+\mathrm{\Omega }_\mathrm{M}z)z(2+z)\mathrm{\Omega }_\mathrm{\Lambda }}}.$$
(2)
Note that we have set $`c=1`$ and $`h(z)`$ and $`D_V(z)`$ should be given using the same set of units, e.g., Mpc.

Thus, our parameters are the reddening parameter $`R_V`$ and the interaction length at zero redshift $`D_{0V}(\sigma \rho _{\mathrm{dust}}^{})^1`$.

The Monte-Carlo simulation program SNOC  was used to perform the integral in Eq. (1) numerically. Note that the model is valid also for a patchy dust distribution, as long as the scale of inhomogeneities is small enough, i.e., $`1/\sqrt{N}1`$ where $`N`$ is the number of dust clouds intersected by the light-ray.

Perfectly grey dust (i.e., with wavelength-independent absorption) can obviously mimic the effects of any dark energy component if one allows for arbitrary fine-tuning of the dust distribution. We have investigated whether any realistic dust model related to astrophysical sources are consistent with the available data. We note first that the dust density has to be related to astrophysical sources, such as star formation, and second, that physically reasonable dust models generally implies a wavelength-dependence in the absorption and scattering properties.

As in Ref. , two dust distributions are considered: $`\rho _{\mathrm{dust}}(1+z)^\alpha `$, where

$`\alpha (z)=\{\begin{array}{cc}3\mathrm{for}\mathrm{all}z,\hfill & \mathrm{Model}\mathrm{A}\\ 0\mathrm{for}z>0.5(3\mathrm{for}\mathrm{lower}z).\hfill & \mathrm{Model}\mathrm{B}\end{array}`$ (5)
Thus, Model A has a constant comoving dust density whereas Model B with $`\rho _{\mathrm{dust}}(z>0.5)=\rho _{\mathrm{dust}}(z=0.5)`$ implies dust creation at $`z>0.5`$. Models A and B are considered to be extreme in opposite directions, i.e. they should capture the full range of physical models leading to the production of intergalactic dust.

In our simulations we explore values of $`R_V`$ in the interval ranging from 1 to 9 and the interaction scale-length at zero redshift $`D_{0V}`$ in the range $`10\mathrm{10\hspace{0.17em}000}`$ Gpc. In a $`[\mathrm{\Omega }_\mathrm{M}=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7]`$-cosmology, $`z=[0.5,1.0,1.5]`$ correspond to proper distances of 1.55, 2.37 and 2.84 Gpc, respectively. Note that since the dust density varies with redshift according to Eq. (5), the interaction length will vary correspondingly.

## 3 Spectral homogeneity of QSOs

QSOs have been found to have remarkably homogeneous spectral properties in spite of their wide range in brightness. The Sloan Digital Sky Survey collaboration (SDSS, ) has produced empirical spectral templates based on a couple of thousand objects in the redshift range $`0.044<z<4.789`$ showing a spectrum-to-spectrum $`1\sigma `$ difference of approximately 20 % in the rest system wavelength range 1200 to 8500 Å . HST QSO data was used by  to extend the spectral template down to 300 Å in order to obtain full wavelength coverage for all broadband filters (see Sec. 4) over the entire redshift interval used in this paper.

## 4 The QSO sample

The SDSS and 2dF  Galaxy Redshift Survey will eventually include about 105 QSOs up to redshift $`z5`$. We have used a data sample of 3814 QSOs from the SDSS Early Data Release (EDR) . These have redshifts between $`z=0.15`$ and $`5.03`$ distributed according to Fig. 1.

The optical magnitudes are measured through five broadband filters ($`ugriz`$) with limiting point-source magnitudes of 22.3, 22.6, 22.7, 22.4, and 20.5, respectively, for a signal-to-noise ratio of 5:1.

The quoted photometric calibration errors are 3 % in $`g,r`$ and $`i`$, and 5 % in $`u`$ and $`z`$ . The instrumental response functions are given in Ref.  where we use the second column (respt) corresponding to 1.3 airmasses at APO for a point source.

We exclude QSOs that are classified as extended to avoid contamination from the host galaxy. This cut leaves 3692 of the original 3814 QSOs. Applying the 5:1 signal-to-noise ratio limiting point-source magnitude cut leaves 3556 of the point-like QSOs. In order to avoid confusion with Seyfert galaxies, we also exclude objects at $`z<0.5`$ . To avoid putting excessive weight in the analysis on a possible dust distribution at very high redshifts (and also to minimize a possible bias from flux limitations) we have only considered QSOs with redshift $`z<2`$ in our primary analysis, leaving 2740 point-like QSOs at $`0.5<z<2`$. An alternative method to constrain intergalactic dust at $`z>2`$ by studying the thermal history of the intergalactic medium is presented in Ref. .

## 5 Method and results

We divide the QSO sample into redshift bins of size $`\mathrm{\Delta }z=0.1`$. In each redshift bin, we calculate the mean colour, e.g., $`(ug)_{\mathrm{obs}}`$ as a measurement error weighted average. Objects that have colours that are more than $`2\sigma `$ off from the mean value are rejected to avoid contamination from object with anomalous colours, e.g., QSOs that have suffered severe extinction in the host galaxy.

For each set of parameter values $`R_V=[1,9],D_{0V}=[10,\mathrm{5\hspace{0.17em}000}]`$ Gpc and redshift, we simulate the attenuation due to grey dust and add to the median QSO spectrum obtained in Ref. .

For a pair of broadband filters $`X`$ and $`Y`$, we define

$$\mathrm{\Delta }(XY)_{z,R_V,D_{0V}}=(XY)_z(𝒳𝒴)_{z,R_V,D_{0V}}$$
(6)
where $`(XY)_z`$ is the observed colour of a QSO at redshift $`z`$ and $`(𝒳𝒴)_{z,R_V,D_{0V}}`$ the simulated colour for a specific set of dust parameter values, i.e, the colour from the template spectrum corrected for differential extinction and convolved with the SDSS filter functions. We then calculate the $`\chi ^2`$ as a function of $`R_V`$ and $`D_{0V}`$ as

$$\chi ^2[R_V,D_{0V}]=\underset{i,j=1}{\overset{N}{}}\mathrm{\Delta }(XY)_{z_i,R_V,D_{0V}}(V(XY)^1)_{i,j}\mathrm{\Delta }(XY)_{z_j,R_V,D_{0V}}$$
(7)
where the sum is over $`N`$ redshift bins and $`V(XY)_{i,j}`$ is the covariance matrix defined as

$$V(XY)_{i,j}=\sigma (X)_{\mathrm{sys}}^2+\sigma (Y)_{\mathrm{sys}}^2+\delta _{ij}\left[\sigma (X)_{\mathrm{tem}}^2+\sigma (Y)_{\mathrm{tem}}^2+\sigma (XY)_{\mathrm{obs}}^2\right]$$
(8)
Here, $`\sigma _{\mathrm{sys}}`$ denote the systematic errors associated with each broadband filter (see Sec. 4), $`\sigma _{\mathrm{tem}}`$ the template error (Sec. 6.2) and $`\sigma _{\mathrm{obs}}`$ the error in the mean observed colour in each redshift bin (after the cut of outliers described above). The scatter in each redshift bin is approximately 0.2 mag so with the order of 100 QSOs in each redshift bin we have $`\sigma _{\mathrm{obs}}0.02`$ mag.

In Fig. 2, the colour-redshift evolution for a couple of different dust scenarios are compared. The upper left panel shows $`\mathrm{\Delta }(ui)`$ vs redshift for the Model A best fit parameters $`R_V=9`$ and $`D_{0V}=100`$ Gpc and the lower left panel shows $`\mathrm{\Delta }(ui)`$ vs redshift for the Model A best fit parameters $`R_V=2`$ and $`D_{0V}=1000`$ Gpc. These parameters yield $`\chi ^2/\mathrm{dof}0.5`$, indicating that systematic errors may have been overestimated (see Sec. 6). These plots are to be compared with the case of $`R_V=4`$ and $`D_{0V}=10`$ Gpc that give high $`\chi ^2`$-values (i.e., low probabilities) for the same dust distribution model (right panels).

We have also investigated the dimming effect on high-$`z`$ sources at different wavelengths, in particular the restframe $`B`$-band flux of SNIa. In Fig. 3, we have combined the results from this investigation with the results from the $`\chi ^2`$-test. In the left panels, the results for $`\mathrm{\Delta }(ui)`$ assuming the Model A (upper panel) and Model B (lower panel) dust distribution are shown. The right panels show the corresponding results for $`\mathrm{\Delta }(gz)`$. The red, dark orange, light orange, and yellow regions indicate 68, 90, 95 and 99 % confidence level allowed regions from the $`\chi ^2`$-test. The four contour levels labeled \[0.2,0.1,0.05,0.02\] correspond to the amount of attenuation caused by grey dust when integrating over the restframe $`B`$-band of a set of SNIa at a redshift of $`z=1`$.

In Fig. 4, we show combined results from colours $`\mathrm{\Delta }(ui)`$ and $`\mathrm{\Delta }(gz)`$. Allowing for any value of $`R_V`$ within the parameter range, we are able to rule out $`\mathrm{\Delta }m>0.2`$ at 99 % confidence level for SNIa at $`z=1`$. Assuming intergalactic dust properties more similar to galactic dust ($`R_V4`$) we are able to do almost a factor of ten better ($`\mathrm{\Delta }m0.03`$). Note also that all our results are consistent with a no dust hypothesis since the $`\chi ^2`$-values flattens out when $`D_{0V}\mathrm{}`$. The small scale features in Fig. 3 and 4 are due to low resolution in the parameter grid.

Redoing the analysis combining $`\mathrm{\Delta }(uz)`$ and $`\mathrm{\Delta }(gi)`$ yields very similar results.

## 6 Systematic effects

### 6.1 Selection effects

The selection algorithm for the SDSS quasar sample varied throughout the time that the EDR observations were obtained, so the objects in this catalog were not found via a uniform set of selection criteria . Richards et al  analysed the limitations on the colour range of QSOs and found that for $`z=2`$, the SDSS survey was sensitive to reddened QSOs with $`E(BV)0.54`$. We conclude that selection effects do not significantly affect our ability to compute the average colour of the QSO sample out to redshift $`z=2`$.

### 6.2 Spectral evolution of QSOs

The possible spectral evolution is limited by the nearly Gaussian 20 % variations found for the whole sample. As we only consider less than half of the survey’s redshift range, systematic colour trends with redshift would not exceed $``$ 10 %. In our analysis, we use a systematic uncertainty on the magnitude derived from the median spectrum vs redshift of 10 % in $`u`$ and $`z`$ and 5 % in $`g,r`$ and $`i`$. In general, unless the spectral evolution trends exactly cancel the differential extinction from dust, our limits on grey dust are conservative in the sense that the reddening that we attribute to possible extinction by dust may stem from a redder population of QSOs. Of course, a possible problem is the fact that the template spectrum we use is derived from spectra that themselves might be affected by dust absorption and extinction. Such an effect would not cause any major change in differential evolution of the QSO colours vs redshift. However, it should be noted that it could cause a close to constant offset in the colour-redshift diagram that might bias the outcome of the $`\chi ^2`$-analysis. An iterative method where this is corrected for by shifting the colour of the template to match the colour of the current best fit model at the mean redshift of the QSO sample ($`z1.4`$) yields results similar to Fig. 4 with the exception that for Model B, the upper right corner (corresponding to high values of $`R_V`$ and $`D_{0V}`$) is also excluded. However, the main conclusion that $`\mathrm{\Delta }m>0.2`$ is excluded for SNIa at $`z=1`$ remains unchanged. Improved statistics and systematic errors in the QSO surveys will improve the control of systematic effects from the template spectrum, see Sec. 7.

### 6.3 Host galaxy dust

Internal reddening may come from depletion by material associated with the torus of gas and dust that is thought to surround the central QSO engine . Richards et al.  estimate that 6 % of the SDSS quasar sample shows signs of reddening by dust, probably in the host environment. The computation of average colours for each redshift bin takes this into account by an iterative rejection mechanism where QSOs more than $`2\sigma `$ away from the average colour at that redshift are rejected, see Sec. 5. Thus, the tails of the distribution are suppressed. Remaining reddened QSOs would tend to make our limits more conservative, as the residual host galaxy dust would be counted as possible intergalactic dust.

### 6.4 Evolution of dust properties

Although we have assumed a constant value of $`R_V`$ for the redshift range $`0.5<z<2.0`$, our derived limits show only moderate dependence on the actual value of the reddening parameter. A variation of $`R_V`$ within the considered range would thus not lead to a significant difference in the observed attenuation of high-$`z`$ SNIa.

## 7 Future prospects

The Sloan Digital Sky Survey has recently made available the Data Release One (DR1) with spectroscopic parameters for 17 700 medium-$`z`$ ($`z<2.3`$) and 980 high-$`z`$ QSOs ($`z>2.3`$), including reprocessed data from the Early Data Release. A catalog of DR1 QSOs is being prepared  and we plan to redo the analysis with the full data set when the catalog is completed. The increased number of QSOs will of course help in decreasing the statistical error. Hopefully, the added data will also be helpful in order to improve our understanding of different systematic effects. For example, it should be possible to constrain the possible evolution of QSO colours and minimize the systematic error from the spectral template by studying the homogeneity of QSO spectra over a range of redshifts.

## 8 Summary and conclusions

Studying the redshift evolution of colours of QSOs from the Sloan Digital Sky Survey Early Data Release, we constrain the effect from extinction by homogeneous intergalactic dust on high-$`z`$ observations. We compare observed spectra from 2740 QSOs at $`0.5<z<2`$ with simulated spectra obtained by adding the attenuation due to different grey dust scenarios to a reference median QSO spectrum. For each dust model, we also compute the related extinction effects on high-$`z`$ SNIa.

For a wide range of models for the dust distribution and extinction properties, we are able to rule out $`\mathrm{\Delta }m0.2`$ at 99 % confidence level for SNIa at $`z=1`$. This is comparable to current SNIa survey errors and allows us to rule out the “grey” dust scenario as an alternative to the cosmological constant explanation for the observed faintness of high-$`z`$ SNIa. It should be noted that these results are based on the specific parameterization of the wavelength-dependent dust attenuation given in Ref.  and that the validity of this specific parameterization for very high values of $`R_V`$ remains untested.

To constrain the dark energy properties (e.g., the equation of state parameter) with future SN surveys such as SNAP, yet better constrains on dimming by grey dust are required. Assuming intergalactic dust properties similar to galactic dust ($`R_V4`$) the present analysis rules out $`\mathrm{\Delta }m0.03`$. However, for a more general dust scenario, larger statistical samples and better control of the systematics of the QSO sample are required. Hopefully this will be within reach with future data sets from, e.g., the Sloan Digital Sky Survey.

## Acknowledgments

The authors would like to thank Anthony Aguirre and the anonymous referee for careful reading and valuable suggestions improving the quality of the manuscript. They would also like to thank Mamoru Doi and Gordon Richards for help regarding the SDSS response functions. A.G. is a Royal Swedish Academy Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the US Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The Participating Institutions are the University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, the Johns Hopkins University, the Max Planck Institute for Astronomy (MPIA), the Max Planck Institute for Astrophysics (MPA), New Mexico State University, Princeton University, the United States Naval Observatory, and the University of Washington.