A New, Deep JVLA Radio Survey of M33

The majority of the emission in our radio maps is in the form of point or discrete sources.⁸ Hence, we have constructed a catalog of these sources to permit comparisons between these new radio data and catalogs at other wavelengths.

The radio source catalog was created by averaging the data in the 1.4 and 5 GHz bandpasses to produce a single detection image. For the 1.4 GHz data, the entire bandpass was used, and 14 of the 16 125 MHz channels from the 5 GHz observations were used. (Channels 1 and 2 of the 5 GHz data were contaminated by RFI.) The 5 and 1.4 GHz bands are equally weighted in the combined image, resulting in source lists that are not a strong function of the source spectral index. This equal weighting is not optimal for a typical extragalactic source (which has a power-law spectrum ν^α with α ≈ −0.7) but is a reasonably "spectrum-neutral" strategy for the source detection, which is more appropriate for our study.

The complexity of the radio emission causes a significant spatial variation in the background level with position. We developed a multiresolution (median pyramid) image processing algorithm, both for use in source detection and to subtract the background from the images to make the source detection threshold more uniform. The multiresolution median pyramid images were also used for both source detection and flux density measurements for the sources. We describe this new and somewhat complex algorithm in detail in the Appendix.

An island map (segmentation map) technique is used to indicate image regions where sources were detected. Figure 4 shows an enlargement near the complex nuclear region as an example of the results of this process. Here we provide a brief summary of the algorithm used to compute the islands and the positions listed in our catalog:

• 1.  
  The detection image is separated into a stack of images of equal size that have structures on scales ranging from compact to very extended. The sum of this stack is equal to the original image.
• 2.  
  At each stack level, the local rms noise is estimated, and pixels that exceed that noise by a factor of 2 are included in a segmentation map indicating potential source positions. Contiguous pixels are grouped into "islands" for different sources.
• 3.  
  The lowest resolution image from the stack is treated as a variable sky level and is not used for source detection or in flux calculations.
• 4.  
  Overlapping islands at different levels are merged and/or split as necessary to group them with islands at other levels. When this is complete, each pixel of the image is assigned to a single island in one or more levels of the stack or is not assigned to any source.

Zoom In Zoom Out Reset image size

As seen in Figure 4, the resulting islands can be irregular in shape and sometimes include what to one's eye might appear as two or more sources (or a point source with adjacent more diffuse emission). However, in complex regions, the islands usually fit together snugly and thus capture all of the radio flux from the region.

Each island corresponds to a single source in the catalog. The fluxes for sources are calculated by computing similar multiresolution image stacks for each frequency band image. Pixels that are assigned to an island at a given resolution level are included in the summed flux for the source. Note that the island size may change at each resolution level, which effectively leads to an adaptive weighting of the image pixels that depends on the source profile. There is a PSF-dependent correction factor (close to unity) for flux that spills outside the island. The island shapes are identical for all the different frequency bands (as are the PSFs after the beam sizes have been matched), which leads to accurate and reliable spectral indices. See the Appendix for further details.

The moments of the mean detection-image flux densities in each island are used to determine the centroid position and a representative elliptical morphology (FWHM major and minor axis and position angle) for each source, as listed in the catalog. The size quoted in the catalog has not been corrected for the round 59 FWHM beam size. For an unresolved source, the major and minor axes will be approximately equal to 59. Note that the size is not constrained to be larger than this value; sources with sizes smaller than the beamwidth (due to noise fluctuations) can have integrated flux densities smaller than their peak flux densities. In such cases, the peak flux density may be a more reliable estimate of the source brightness.

The integrated flux densities for each of four bands (two from 1 to 2 GHz and two from 4 to 6 GHz) are fitted with a power law to determine the spectral index α for each source, F_ν = F_int(ν/ν_p)^α. The noise in the band fluxes is used to determine the errors in F_int and α. The pivot frequency ν_p is chosen as the frequency where the covariance between F_int and α is zero, which also is the frequency of the optimal signal-to-noise ratio ${F}_{\nu }/\sigma ({F}_{\nu })$. The spectral index fit is constrained to −3 ≤ α ≤ 3; sources with spectral indices at these limits have σ(α) = 0.

The accuracy of the spectral indices degrades in the outer region where fewer frequency bands are available (Figure 5). In the center of the galaxy, sources are measured across the full 1–6 GHz bandpass. The frequency-dependent field of view reduces the sky area covered at higher frequencies, so outer regions of the galaxy are covered by fewer frequency bands (Figure 1). The uncertainty in the spectral index is determined by both the signal-to-noise ratio of the flux density measurements and by the frequency lever-arm of the data. The small region with three bands has only slightly worse noise than the four-band data, but the two-band data (with measurements only in the 1–2 GHz bandpass) have spectral index errors that are six times larger. As a result, only the brightest sources with two-band detections have accurately measured spectral indices. (Obviously, the outermost region with only a single band detected has no spectral index information at all.) Fortunately, the majority of the sources are found in the inner region of the galaxy where there is full-frequency coverage.

Zoom In Zoom Out Reset image size

In addition to the island centroids, the catalog also includes the peak flux density and the position of the peak for each island. These values are computed from the detection image, just as for the flux-weighted positions. The peak is simply the brightest pixel in the island after background subtraction. The position of the peak (which may differ from the flux-weighted mean in asymmetrical sources) is determined in the sharpest channel of the multiresolution image stack where the source is detected. This means it represents the position of the peak for the most compact source component. We found that this choice gives more accurate peak positions in crowded regions.

The results are summarized in Table 2,⁹ which includes all sources detected at a signal-to-noise ratio F_int/σ(F_int) ≥ 4. It also includes seven sources that fall below that detection threshold but that have associations in external X-ray or H ii region catalogs; those sources are flagged as subthreshold objects. A full description of the columns in the table is found in Table 3.

Our final catalog contains 2875 sources with 1.4 GHz flux densities ranging from about 2.2 μJy to 0.17 Jy. (The faintest 1.4 GHz fluxes are for objects with inverted spectra that are brighter at 5 GHz.) Of these, 670 have spectral indices determined with an accuracy σ(α) ≤ 0.15. Essentially all sources with 5 GHz data (2/3 of the catalog) and fluxes in excess of 150 μJy have well-measured spectral indices. The typical noise in the integrated flux density is <5 μJy, although there are a few crowded regions where the noise is significantly higher (Figure 6).

Zoom In Zoom Out Reset image size

With the sensitivity and resolution of these new data, some of the larger optical emission structures like giant H ii regions or H ii emission complexes incorporate numerous individual radio emission peaks. Figure 7 shows two examples, IC 133 and NGC 592, where the color coding and labeling are similar to those in Figure 3. IC 133 shows a bright, relatively high-excitation (strong [O iii]-emitting) bi-lobed structure in the optical with five radio peaks on portions of the main nebula. The radio contours extend to the lower left, incorporating another eight radio peaks that are listed in the catalog. Several optical SNRs are indicated nearby but are not directly related to IC 133. In contrast, NGC 592 has a very bright core of optical emission with a number of fainter loops of surrounding emission that are mapped reasonably well by the radio contours, including a faint arm of emission to the west. Several SNRs are associated with this region as well. Several radio peaks correspond with the bright core region, while the extended radio contours incorporate a dozen or more additional radio peaks.

Zoom In Zoom Out Reset image size

The radio catalog described above was constructed without any reference to the positions of known optical SNRs in M33 (although many of the remnants were detected). However, we have also constructed a separate catalog of radio properties at the positions of known optical SNRs by integrating the radio images over the elliptical region determined for each optical SNR, as measured previously by either L10 or LL14. (When a SNR appeared in both lists, we somewhat parochially adopted the regions of L10.) This method allows us to establish either measured flux densities or appropriate upper limits for the radio emission from all SNR candidates, whether or not a radio source was independently detected at the SNR position. We will refer to this version of the SNR radio data as the "forced photometry" SNR extraction catalog.

The calculations of flux densities, spectral indices, and flux-weighted positions proceeds much as described above for the radio catalog. The island map for the forced photometry catalog is determined from the elliptical regions in the optical SNR catalog rather than from the radio map itself. Because the island map makes no reference to the radio morphology, we do not have multiresolution islands but instead sum the radio flux density over the entire region using the single background-subtracted image in each radio band. The catalog does not include the radio peak information (which was not found to be useful) but does include the input positions for the regions along with the radio-flux-weighted centroid positions.

Associations between SNRs and objects in our master radio catalog are determined by the overlap between the radio island map and the SNR island map. Many SNRs have unambiguous matches with cataloged radio sources, but there are also ambiguous cases where a single SNR island overlaps several radio islands and vice versa. The information on the overlapping sources and a flag that captures information on the ambiguity of the association are also included in the table.

Of the 217 SNRs and SNR candidates in the merged list of L18, 216 are in the region covered by our radio images (see Figure 1). The only object falling completely outside the radio region is LL14-195. There are 188 SNRs in the central region covered by both the 1.4 and 5 GHz frequency bands plus two more that have data from 1.4 GHz and the low-frequency 5 GHz bandpass only. Of the 26 sources that have only 1.4 GHz data, only one (L10-003) has a reasonably accurate spectral index; we have little spectral index information in the outer region of the galaxy.

We find that 155 of the 216 SNRs and candidates are detected above 3σ at radio wavelengths, a detection rate of 72%. That detection rate rises to 76% (145 of 190) in the central region where we have full-frequency coverage. The remainder have radio upper limits from the forced photometry exercise. (Some of the sources below the 3σ threshold are also likely to be detected: there are, for example, 12 2σ detections discussed further below.) Of the 155 3σ radio-detected SNRs, 126 came from the list of 137 candidates in L10, while 29 came from the 79 SNRs that were first suggested as candidates by LL14. The fact that a much lower percentage (37% versus 92%) of the extended list of SNR candidates identified by LL14 were detected in the radio is undoubtedly related to the fact that many of the LL14 candidates are optically fainter than those in the L10 list.

An important consideration in setting the detection threshold at 3σ is the expected rate of chance coincidences between the SNR islands and the radio map. The radio images are sufficiently crowded that randomly placed ellipses will sometimes fall on radio sources. This is particularly true in the center of the galaxy and in the spiral arms. We estimated the rate of chance radio detections by shifting the SNR islands by ±1 arcmin in R.A. and decl. The size of the shift was chosen to be larger than the largest island's major axis. The shift is large enough to move off the local object while still being small enough to be sampling the same environment (spiral arm, galaxy center) as the true SNR position. The shifted calculations also use the same distribution of region sizes as the real catalog, which has a significant influence on the background rate. Eight shifts were done, including shifts only in decl., only in R.A., and in both R.A. and decl. For every shift, forced photometry was performed on all the regions using exactly the same approach as for the real positions.

A comparison between the measured flux densities for the real SNR regions and a typical example of the shifted SNR region fluxes is shown in Figure 8. While positive detections are much rarer at the shifted positions, they are common enough that they cannot be ignored.

Zoom In Zoom Out Reset image size

All eight shifted catalogs for each object class were combined to compute the distribution of the signal-to-noise ratio, F_int/σ(F_int), for random positions on the sky. With eight simulations of the random detections, the simulated noise distribution is reasonably well determined. Given the cumulative distribution of the number of sources exceeding a given signal-to-noise ratio, we then compute for every source in the real catalog the probability that a random source of equal or greater flux would be detected.

Finally we compare the cumulative distribution of the actual number of sources detected as a function of the detection threshold with the cumulative distribution of the random detection probability above that threshold (Figure 9). That gives an estimate of the number of false detections that are included as a function of the detection threshold.

Zoom In Zoom Out Reset image size

We conclude that a detection threshold of 3σ leads to acceptable contamination by false detections. Of the 141 SNRs above 4σ, ∼5 are likely to be due to false matches. Adding the SNRs between 3σ and 4σ increases the number detected by 14 while adding only two additional false matches. Even detections between 2σ and 3σ are likely to be fairly reliable, with ∼75% of the 12 detections in that range expected to be true detections. (These coincidence rates are lower than intuition might suggest because the high true detection rate for SNRs means that there have been relatively few "rolls of the dice" at empty sky positions.)

The basic characteristics of the SNRs as well as results of our forced photometry extraction are presented in Table 4,¹⁰ where the nonradio data are from L18. The [S ii]:Hα column in this table indicates whether optical spectra exist and whether these showed a [S ii]:Hα ratio of 0.4 or greater. The X-ray column indicates whether an object was detected at 3σ or greater with Chandra by L10 or with XMM-Newton by Garofali et al. (2017). A full description of the columns in the table is found in Table 5. Our forced photometry radio catalog of SNRs contains 178 SNRs that have been optically confirmed by the [S ii]:Hα criterion; there are another 27 objects that have spectra but show lower [S ii]:Hα ratios, and 12 for which there are not yet any optical spectra. A total of 112 of the objects are also detected in X-rays.

The SNR forced photometry table also identifies radio sources from the master catalog (Table 2) that have detection islands that overlap the SNR regions. In cases where more than one radio source overlaps the SNR, all the matches are listed. Multiple matches are sorted by the number of overlapping pixels, so the radio source sharing the largest number of pixels with the SNR region is listed first. Similar information is included in the radio catalog for cases when a radio source overlaps one or more SNRs, and the same information is available for radio matches to the forced photometry X-ray catalog (described below). The Sflag (Table 2) and Rflag (Tables 4 and 6) columns are bit flags with values that help determine whether a match is reliable. For the Sflag column, for example, the bit values b indicate that

• 1.  
  b = 1: the SNR has a radio match,
• 2.  
  b = 2: the match is unambiguous, so only one radio source matches the SNR, and
• 3.  
  b = 8: the match is mutually good, meaning that the best radio source for this SNR also has this SNR as its own best match, where "best" is determined by the number of overlapping pixels between the islands.

The most reliable matches will have flag bit 8 set. All of those will have bit 1 set as well, and most of them will also have bit 2 set. Flag values of 9 and above indicate good matches, with the best having flag values of 11. The distributions of the values of Sflag (in the radio catalog) and Rflag (in the SNR and X-ray catalogs) are shown in Table 7.

Of the 155 SNRs detected at 3σ in the forced photometry catalog, 134 are associated with a source in the master radio catalog. Of these, there are 19 sources in the SNR catalog associated with two (or three) radio sources. There are 21 SNRs detected as faint sources in the forced photometry catalog that are not found in the master catalog. Because the master radio catalog has a 4σ detection limit (compared with 3σ for the forced photometry), it is not surprising that there are sources below the radio catalog limit.

In some cases where the forced photometry does not yield a 3σ detection, radio emission is still apparent at the optical remnant location. These cases arise from the incommensurability of the flux density distribution in the optical and radio bands, the presence of a nearby bright radio source, or radio emission just below the detection threshold. Examples of these cases are shown in Figure 10. In some cases (L10-028, L10-030, LL14-038, and L10-131), only the brightest part of the optical shell has coincident radio emission, and the overall source falls below threshold (Figure 10(a)). In other cases (L10-058, LL14-096, L10-132, and LL14-142) a bright, nearby (and apparently unrelated) source raises the local background sufficiently that the radio counterpart falls below threshold (Figure 10(b)). In yet other cases (LL14-125, LL14-128, L10-048, L10-080, L10-092, and LL14-058), a clear radio counterpart is apparent but falls just below the threshold (Figure 10(c)); note that the division between the latter two categories is somewhat arbitrary. In summary, over one-quarter of the remnants not listed here as having radio emission do have modest levels of radio flux associated with them that a deeper and/or higher-resolution survey would definitely have detected.

Zoom In Zoom Out Reset image size

The list of X-ray source positions from the ChASeM33 catalog of Tüllmann et al. (2011) was also used to construct a forced photometry catalog with radio flux densities. The procedure followed was very similar to the SNR forced-photometry catalog described above. The resulting table (shown as a sample in Table 6) has a 3σ threshold for radio detections (but includes the radio measurement regardless of the signal-to-noise ratio). It includes information from the X-ray catalog including the diameter and count rates in various energy bands. The Rflag column is defined exactly as for the SNR forced photometry table, with values of 8 or greater indicating confident (or at least unconfused) associations. The format is similar to Table 4; a full description of the columns is found in Table 5.

Of the 662 X-ray sources, 319 are detected in our radio survey. All of these X-ray sources fall in the region that has both 1.4 and 5 GHz radio data, and 102 of the detections have radio spectral indices determined to an accuracy σ(α) < 0.15.

Many of the radio sources in the catalog are associated with star-forming H ii regions, with optically thin thermal bremsstrahlung radiation having a spectral index α ∼ 0. In order to determine which sources in our radio catalog are associated with H ii regions, we have computed the Hα surface brightness over the radio island regions for each source in the catalog using the Hα image obtained by one of us (P.F.W.) using the Burrell Schmidt camera at Kitt Peak in Arizona (see L10). The tables for both the radio catalog and the forced photometry catalogs include the total Hα flux, F(Hα), and the average flux surface brightness (per square arcsec), Σ(Hα).

To determine an appropriate Σ(Hα) threshold to identify H ii regions, we have used the Lin et al. (2017) catalog with spectroscopic observations of 413 star-forming regions in M33. This catalog is not a complete compilation of H ii regions in M33 because its sample was shaped by the constraints of the fiber-fed spectrograph used for these observations. Despite that limitation, it represents a good subset of M33 ionized regions for our analysis.

Figure 11 compares the Hα surface brightnesses for the Lin et al. (2017) fiber positions with the background distribution computed from a set of shifted regions. Five different shifts of 1 arcmin were used to determine the background distribution. We conclude that a surface brightness threshold Σ(Hα) > 3 ×10⁻¹⁶ erg cm⁻² s⁻¹ arcsec⁻² separates H ii regions from background sources, generating a reasonably reliable and complete sample. The background rate is not negligible (about 10% of the background positions exceed this threshold), but some contamination is unavoidable in a galaxy this crowded with star formation.

Zoom In Zoom Out Reset image size

Of the list of 217 SNRs and SNR candidates compiled by L18, 216 are in the region covered by our radio survey and 155 are detected above 3σ via forced photometry, including 84 that have radio spectral indices with an uncertainty of less than 0.15. There are 122 of the detected objects where spectroscopy has confirmed [S ii]:Hα flux ratios exceeding 0.4, the usual optical criterion that an emission nebula is a SNR, and 92 of these have X-ray detections.

The derived radio luminosities at 1.4 GHz and spectral indices of these objects are shown as a function of SNR diameter in Figure 15. Neither the luminosity nor the spectral indices show a significant correlation with diameter. Radio luminosities for the entire sample and the various subsamples show very large dispersions at all diameters. Smaller diameter SNRs are not significantly brighter or fainter than those at large diameter. Similarly, radio spectral index does not appear to evolve with diameter. In contrast, a correlation between spectral index and diameter has been reported for the LMC (Bozzetto et al. 2017); we discuss this difference further in Section 4.3.

Zoom In Zoom Out Reset image size

Comparisons between the 1.4 GHz luminosity and the Hα and X-ray luminosities are shown in Figure 16. Although there is considerable scatter, there appears to be a trend between Hα and radio luminosity. More luminous SNRs in Hα appear to be more luminous at radio wavelengths. Whether this trend is physical is hard to determine because the scatter is large. The SNRs in M33, like those in essentially all other galaxies, were initially identified through optical interference filter imaging and are, thus, surface brightness limited. Larger SNRs tend to have higher Hα luminosities as a result. By contrast, the X-ray sample is luminosity limited. This accounts for the fact that there are almost no SNRs with L_x less than 4 × 10³⁵ erg s⁻¹, the approximate limit of the Chandra survey. That said, it is interesting that the objects with highest X-ray luminosities also tend to have the highest radio luminosities.

Zoom In Zoom Out Reset image size

3.1.1. Comparison to Previous M33 Radio Detections of SNRs

As noted in the Introduction, the last detailed radio study of SNRs in M33 was carried out by G99, who reported the radio detection of 53 of 98 optically identified SNR candidates that were known at the time. Of these 53 objects, 52 are in the region covered by our radio survey.¹¹ All of these are detected at 3σ or greater in our forced photometry catalog. As shown in Figure 17, there is a fairly good correlation between the flux densities we measure at 1.4 GHz (by constructing a weighted average of our two 1.4 GHz bands) and those measured by G99. The G99 flux densities tend to be higher, and in a few cases considerably higher. Many of the objects that are higher were identified by L10 as likely to be associated H ii emission rather than emission from the SNR (or the SNR only), but the median of the ratio of the flux densities they measured is only 50% higher than those we measure, suggesting that we and they are identifying the same radio sources as SNRs, but with differing amounts of contamination by thermal emission.

Zoom In Zoom Out Reset image size

The Tüllmann et al. (2011) catalog of X-ray sources in M33 contains 662 X-ray sources, the vast majority of which are due to background sources (galaxies and AGNs). We detect radio emission at the position of 319 of these sources at 3σ or greater via forced photometry. Of these, 55 are positionally coincident with SNRs and 25 with H ii regions observed by Lin et al. (2017). As shown in the X-ray hardness ratio diagram in the left panel of Figure 18, most of the X-ray sources we detect lie in the region of the diagram where background sources and X-ray binaries are expected to fall. As expected, the sources identified with SNRs mostly have soft X-ray spectra. The X-ray sources associated with H ii regions are not concentrated in a particular region of the diagram, suggesting a mix of source types (X-ray binaries in M33, unrecognized SNRs, and background objects). The right-hand panel of the figure shows X-ray hardness ratios as function of the radio spectral index. The SNRs form a well-defined region in the diagram; this suggests that one might be able to identify SNRs directly from such a diagram with sufficiently sensitive X-ray and radio observations, even without a prior optical identification.

Zoom In Zoom Out Reset image size

A Venn diagram showing the number of objects detected in the various wavelength bands is shown in Figure 19. To be counted as a confirmed optical SNR, we require not only that a nebula was suggested by either L10 or LL14 as a candidate on the basis of interference filter imagery, but also that a spectrum exists that confirms the [S ii]:Hα ratio is at least 0.4, the conventional criterion for declaring that an emission nebula is a valid optical SNR candidate. There are 98 emission nebulae satisfying this criterion that also show up as positive detections in our radio forced photometry catalog and that have been detected at 3σ or greater with either Chandra or XMM-Newton. Additionally, there are 31 optical SNRs with spectral confirmations but without detected X-ray or radio counterparts, 41 that have a radio counterpart but no X-ray, and eight objects with optical–X-ray detections but no radio counterpart.

Zoom In Zoom Out Reset image size

The small numbers in the nonoptical parts of the diagram deserve specific discussion. The location of an individual object in the Venn diagram depends on many factors, not the least of which is our application of the optical spectral confirmation criterion above and beyond the imaging candidate status. There are optical SNR candidates that were identified from imaging but have no optical spectra, and there are optical candidates for which the spectra obtained showed [S ii]:Hα ratios somewhat below 0.4. The objects that lie in the X-ray (3), radio (13), and radio–X-ray (3) portions of the Venn diagram align with some of these objects.

For the 13 "radio-only" objects, five do not have spectra, two were too contaminated for accurate spectral measurements, and six have spectra with derived [S ii]:Hα ratios below 0.4. Three of these objects had [S ii]:Hα ratios between 0.36 and 0.39, and the other three were near 0.25. However, for fainter objects or objects with nearby H ii contamination, even the spectrally derived ratios can have significant errors bars. If the background Hα is even slightly undersubtracted in the spectra, the derived [S ii]:Hα ratio can be artificially lowered. The presence of associated radio emission for these eight objects significantly enhances the probability that these nebulae are SNRs.

Of the three "X-ray-only" objects, one has no optical spectrum, one is too contaminated to make a good measurement, and one has a spectrum with a [S ii]:Hα ratio of 0.27, significantly below the threshold (LL14-008). The presence of an X-ray source at the optical position strengthens these as possible SNRs as well.

The three objects shown as having both radio and X-ray emission would already seem to be strong candidates. All three have optical spectra with derived [S ii]:Hα ratios of 0.37 (L10-035), 0.33 (L10-040), and 0.37 (LL14-005). Again, slight Hα background subtraction problems could have caused these to fall below our normal threshold ratio, so these are quite likely SNRs.

There are 20 emission nebulae that, other than their original selection as nebulae with enhanced [S ii]:Hα ratios from imaging, have no additional supporting information to suggest that they are SNRs and, thus, do not appear in this diagram. This does not necessarily mean that they are not SNRs since some objects have not been observed spectroscopically, and one other was simply outside of the region surveyed with the JVLA. Fourteen of the 20 have spectra that show ratios below the normal threshold, so without further corroboration the majority of these are probably not SNRs.

With an X-ray luminosity of ∼10³⁹ erg s⁻¹, M33 X-8, the source coincident with the the nucleus of M33, is one of the brightest persistent X-ray sources in the Local Group (Long et al. 1981). The nuclear star cluster is compact but otherwise unprepossessing optically; its spectrum can be modeled in terms of two bursts of star formation, one with an age of 40 Myr involving stars with a total mass of 9000 M_☉, and another of age 1 Gyr with a mass of 76,000 M_☉ (Long et al. 2002). If a black hole (BH) exists at the center of the cluster, it must have a mass <1500 M_☉ (Gebhardt et al. 2001). It has been suggested (although not confirmed) that the nuclear source has a period of 109 days (Dubus et al. 1997). Taken together, these observations imply that the source is not an AGN but is instead an ∼10 M_☉ X-ray binary radiating near the Eddington limit (see, e.g., Foschini et al. 2004; Middleton et al. 2011).

As shown in Figure 20, a nonthermal radio source, W19-0683, is coincident with the nucleus. The spectral index is −0.39 ± 0.11, and the flux at 1.4 GHz is 0.61 ± 0.04 mJy. The radio source is extended, with a ratio of peak to integrated flux density of 0.35 (compared with unity for a point source), and a nominal major axis of 18'', or 70 pc. In the forced photometry catalog, where source T11-318 corresponds to the nuclear source, the 1.4 GHz flux is similar, 0.60 ± 0.04 mJy, with a slightly steeper spectral index of −0.68 ± 0.17. (The spectral indices in the main catalog and the X-ray catalog are consistent, differing by 1.4σ.) The nuclear source was also observed by G99, who reported a consistent flux density (0.6 ± 0.2 mJy at 1.4 GHz) and spectral index (−0.8 ± 0.3). For comparison, SgrA* in our Galaxy, which arises from the accretion disk of a 4 × 10⁶ M_☉ BH, is a 0.7 Jy source at 1.4 GHz, corresponding to ∼0.07 mJy at the distance of M33, and has a spectral index of −0.3 (Duschl & Lesch 1994).

Zoom In Zoom Out Reset image size

Whether M33 X-8 is a very low mass BH at the nucleus of M33 or simply, as is generally thought, a BH binary, the existence of radio emission is not surprising, given that BH systems on all mass scales produce radio emission (see, e.g., Plotkin et al. 2012). The 1.4 GHz radio luminosity of M33 X-8 is about 6.7 × 10³² erg s⁻¹, less than 10⁻⁶ the X-ray luminosity. Exactly how to fit M33 X-8 into the zoo of BH binaries remains unclear. It has properties that resemble those of a microquasar trapped in the so-called high soft state, where disk accretion is near the Eddington rate (Long et al. 2002; Foschini et al. 2004). Most microquasars, however, exhibit major outburst events, and M33 X-8 has not been seen to vary appreciably. During the rise to X-ray maximum, microquasars show active radio jet emission (with flat spectral indices) and hard X-ray spectra. At outburst maximum, however, the X-ray spectra transition to something similar to that of M33 X-8, and the radio spectral indices steepen. The idea is that active jet generation has ceased (or ebbed at least), and the residual emission arises from the interaction with the surrounding medium (Fender et al. 2004).

Because of its X-ray luminosity and the fact that the luminosity is persistent, M33 X-8 has also been characterized as an ultraluminous X-ray source (ULX). It is at the low luminosity end of sources classified as ULXs that also have X-ray spectra fit in terms of steady-state accretion onto a massive BH. Other ULXs in this state observed at radio wavelengths have spectral indices characteristic of optically thin synchrotron emission. The radio emission for a number of the nearest ULXs (and microquasars) has been resolved and appears to arise from a jet interaction with the surrounding ISM (Pakull & Mirioni 2003). The typical sizes of the resolved bubbles around these sources are 100–300 pc. At the resolution of our observations (about 20 pc), the bulk of the radio emission in M33 X-8 appears pointlike, although, as mentioned above, the source is extended at the ∼70 pc level. Thus, source extent does not rule out the jet-interaction scenario in M33 X-8. The prototypical high-mass X-ray binary Cygnus X-1 is surrounded by a 5 pc diameter radio ring that appears to have been inflated by a jet from the central source (Gallo et al. 2005). Clearly, sensitive, higher-resolution radio observations are needed to clarify the situation for M33 X-8.

In addition to detecting a nuclear radio source, radio sources were also detected at the positions of several other known X-ray binaries in M33. Among these is a source with a 1.4 GHz flux density (in the forced photometry catalog) of 0.63 ± 0.02 mJy at the position of M33 X-7, the well-studied eclipsing BH binary, composed of a 16 M_☉ BH orbiting a 70 M_☉ O-star every 3.45 days (Orosz et al. 2007). However, the spectral index of this source is 0.08 ± 0.02, and the Hα surface brightness at this position is high, about 3.4 × 10⁻¹⁵ erg cm⁻² s⁻¹ arcsec⁻², suggesting that most, if not all, of the radio emission arises from H ii regions along the line of sight. A similar situation pertains to M33 X-4, another bright X-ray source thought to be an X-ray binary (Grimm et al. 2007), where the forced photometry 1.4 GHz flux density is 0.17 ± 0.01 mJy, but the apparent spectral index is +0.3 ± 0.1. This source is also located along a line of sight with considerable Hα flux, 1.1 ×10⁻¹⁵ erg cm⁻² s⁻¹ arcsec⁻². Higher-resolution observations are required to determine whether either of these sources has any intrinsic radio emission.

Our large sample of SNRs with both radio and X-ray luminosities enables some useful comparisons to models for SNR radio emission. Sarbadhicary et al. (2017) describes a relatively simple analytical model for the radio luminosities of SNRs as a function of their radius, evolutionary phase, ISM density, explosion energy, and so forth. According to Equation (A12) in Sarbadhicary et al. (2017), the radio luminosity density L_ν at 1.4 GHz is given by:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{1.4} & = & 2.2\times {10}^{24}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}{\left(\displaystyle \frac{R}{10\mathrm{pc}}\right)}^{3}\\ & & \times \,\left(\displaystyle \frac{{\epsilon }_{e}}{{10}^{-2}}\right){\left(\displaystyle \frac{{\epsilon }_{b}^{u}}{{10}^{-2}}\right)}^{0.8}{\left(\displaystyle \frac{{v}_{s}}{500\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{3.6}\end{array}\end{eqnarray} \tag{ 1 }$$

where R is the shock radius, v_s is the shock velocity, and _e and ${\epsilon }_{b}^{u}$ are the fractions of the postshock energy density $\rho {v}_{s}^{2}$ converted into relativistic electrons and the amplified upstream magnetic field in a diffuse acceleration model, respectively.

If this model were correct, and if all other factors except the shock velocity were the same, L_1.4 would be a strong function of v_s. Indeed, a factor of 2 difference in shock velocity would correspond to a factor of 12 difference in L_1.4. The radio luminosity also increases in proportion to the SNR volume, L_1.4 ∼ R³. As a result, during the SNR free-expansion phase, when the shock velocity is constant, the radio luminosity increases rapidly as R³ ∼ t³. The radio luminosity peaks at the beginning of the Sedov phase, declines through the adiabatic expansion phase as the shock slows, and finally decreases rapidly in the radiative phase as v_s drops off and the SNR expansion slows.

Sarbadhicary et al. (2017) assume that _e is constant with time, citing Soderberg et al. (2005) and Chevalier & Fransson (2006). However, ${\epsilon }_{b}^{u}$ is not a constant in this equation but depends on the magnetic field, shock velocity, and other parameters. Following the discussion in their appendix, the definition of ${\epsilon }_{b}^{u}$ is:

$$\begin{eqnarray}&&{\epsilon }_{b}^{u}=\displaystyle \frac{{B}^{2}/8\pi }{\rho {v}_{s}^{2}}.\end{eqnarray} \tag{ 2 }$$

Citing Bell (2004), Sarbadhicary et al. (2017) argue that at high Alfvén Mach number (M_A), magnetic amplification saturates at a value given by:

$$\begin{eqnarray}&&\displaystyle \frac{{B}^{2}}{8\pi }\sim \displaystyle \frac{1}{2}\displaystyle \frac{{v}_{s}}{c}{\xi }_{\mathrm{cr}}\rho {v}_{s}^{2}\end{eqnarray} \tag{ 3 }$$

where ξ_cr is the fraction of shock energy that goes into relativistic particle acceleration (including electrons and ions). They then assume ξ_cr is a constant value of 0.1 for high Alfvén Mach numbers but declines for smaller M_A:

$$\begin{eqnarray}{\xi }_{\mathrm{cr}}=\left\{\begin{array}{ll}0.1 & \mathrm{when}\ {M}_{{\rm{A}}}\gt 30,\\ {10}^{-2}\left(0.15\,{M}_{{\rm{A}}}+6\right) & \mathrm{when}\ {M}_{{\rm{A}}}\lt 30.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

Here M_A < 30 is an approximate limit—ξ_cr increases in proportion to M_A until it saturates at a value ξ_cr = 0.1.

The final equation that Sarbadhicary et al. (2017) use for ${\epsilon }_{b}^{u}$ is

$$\begin{eqnarray}&&{\epsilon }_{b}^{u}=\displaystyle \frac{{\xi }_{\mathrm{cr}}}{2}\left(\displaystyle \frac{{v}_{s}}{c}+\displaystyle \frac{1}{{M}_{{\rm{A}}}}\right).\end{eqnarray} \tag{ 5 }$$

That makes a transition from values that are proportional to the shock velocity v_s (at high M_A values) to values that are proportional to 1/M_A, meaning ${\epsilon }_{b}^{u}$ scales as 1/v_s, at lower shock velocities. Coincidentally, Figure A1 in Sarbadhicary et al. (2017) shows that the transition between these scalings occurs close to the beginning of the Sedov phase of evolution. In the example they show, the value of ${\epsilon }_{b}^{u}$ varies by only a factor of 3 from the initial SN explosion to an age of 2 × 10⁴ yr; given the modest dependence of L_1.4 on ${\epsilon }_{b}^{u}$, this implies that ${\epsilon }_{b}^{u}$ contributes only a factor of ∼2 to the variation in L_1.4 over a SNR lifetime. That is why the treatment of ${\epsilon }_{b}^{u}$ as a constant in Equation (1) is approximately correct.

The ${\epsilon }_{b}^{u}$ term also introduces dependencies on the density n_H and magnetic field B₀ in the ISM. Also, the magnetic field itself depends on the ISM density because it is coupled to the ISM state. Equation (A4) of Sarbadhicary et al. (2017) can be rewritten to give the ISM magnetic field as

$$\begin{eqnarray}&&{B}_{0}=9\,\mu {\rm{G}}{n}_{{\rm{H}}}^{0.47},\end{eqnarray} \tag{ 6 }$$

where n_H is the ISM hydrogen number density.¹²

To see what this model for the radio emission means, it is useful to recall the Sedov equations:

$$\begin{eqnarray}&&R=5.2\,\mathrm{pc}{\left(\displaystyle \frac{{E}_{51}}{{n}_{{\rm{H}}}}\right)}^{1/5}{\left(\displaystyle \frac{t}{1000\mathrm{yr}}\right)}^{2/5}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{v}_{s}=2000\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{E}_{51}}{{n}_{{\rm{H}}}}\right)}^{1/5}{\left(\displaystyle \frac{t}{1000\mathrm{yr}}\right)}^{-3/5}\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{2}{5}\displaystyle \frac{R}{t}\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&=\ 3900\,\mathrm{km}\,{{\rm{s}}}^{-1}\displaystyle \frac{R/10\,\mathrm{pc}}{t/1000\,\mathrm{yr}}.\end{eqnarray} \tag{ 8c }$$

In the Sedov phase, the radio luminosity density L_1.4 from Equation (1) can be expressed as a function of t:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{1.4} & = & 4.8\times {10}^{25}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}{\left(\displaystyle \frac{t}{1000\mathrm{yr}}\right)}^{-0.96}\\ & & \times \,\left(\displaystyle \frac{{\epsilon }_{e}}{{10}^{-2}}\right){\left(\displaystyle \frac{{\epsilon }_{b}^{u}}{{10}^{-2}}\right)}^{0.8}{\left(\displaystyle \frac{{E}_{51}}{{n}_{{\rm{H}}}}\right)}^{1.32}.\end{array}\end{eqnarray} \tag{ 9 }$$

For this relationship and assuming ${\epsilon }_{b}^{u}$ is approximately constant, then according to Sarbadhicary et al. (2017), we should expect L_1.4 ∝ t^−0.96. A SNR in the Sedov phase should be brightest when it enters the Sedov phase and should decline linearly with time (ignoring the variation in ${\epsilon }_{b}^{u}$).

The radio luminosity can also be expressed as a function of R:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{1.4} & = & 1.0\times {10}^{25}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}{\left(\displaystyle \frac{R}{10\mathrm{pc}}\right)}^{-2.4}\\ & & \times \,\left(\displaystyle \frac{{\epsilon }_{e}}{{10}^{-2}}\right){\left(\displaystyle \frac{{\epsilon }_{b}^{u}}{{10}^{-2}}\right)}^{0.8}{\left(\displaystyle \frac{{E}_{51}}{{n}_{{\rm{H}}}}\right)}^{-9/5}.\end{array}\end{eqnarray} \tag{ 10 }$$

Thus, we expect small-diameter SNRs to be significantly brighter, all other things being the same (L_1.4 ∝ R^−2.4). Clearly, we do not see our observations reflecting this dependence (see Figure 15).

On the other hand, this relationship does allow for significant variations in the radio luminosity at a specific diameter if the ratio E₅₁/n_H varies significantly. A SNR with an explosion energy of 10⁵¹ erg expanding into an ISM with density 1 cm⁻³ would be 63 (4000) times brighter than one expanding into an ISM with density 0.1 (0.001) cm⁻³, assuming both are in the Sedov phase. (See also Asvarov 2006, who reached a similar conclusion about radio luminosity variations associated with variations in E₅₁/n_H from models made with slightly different assumptions.)

If the variations in radio luminosity are primarily due to variations in n_H, then we would expect that to have consequences at X-ray wavelengths, where emission arises from the same shock that powers the radio emission. Simply put, SNRs expanding into a dense ISM are brighter at the same diameter in soft X-rays than those expanding into a tenuous one, so we would expect radio-bright SNRs to also be X-ray bright and vice versa. Berkhuijsen (1986) argued that a correlation of this type does exist using a collection of Galactic and extragalactic SNRs; here we discuss the specific case of M33.

The X-ray luminosity in the Sedov phase has a complex dependence on the physics involved (e.g., nonequilibrium ionization, electron conduction, etc.). However, White & Long (1991) give a simple approximation based on the observation that the X-ray emissivity function Λ_x(T) is constant to within 50% over a wide range of temperature (Hamilton et al. 1983). From Equation (21) in White & Long (1991):

$$\begin{eqnarray}\begin{array}{rcl}{L}_{x} & = & 3.0\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{R}{10\mathrm{pc}}\right)}^{3}{n}_{{\rm{H}}}^{2}\\ & & \times \,\left(\displaystyle \frac{{{\rm{\Lambda }}}_{x}}{1\times {10}^{-22}\,\mathrm{erg}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where we set the contribution of evaporating clouds to zero and are using an ISM density n = ρ/m_H = 0.75 n_H. Combining Equations (10) and (11), we can calculate the ratio of the radio luminosity L_R = νL_ν to the X-ray luminosity:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{L}_{R}}{{L}_{x}} & = & 4.6\times {10}^{-4}{\left(\displaystyle \frac{R}{10\mathrm{pc}}\right)}^{-5.4}{n}_{{\rm{H}}}^{-1/5}{E}_{51}^{-9/5}\\ & & \times \,{\left(\displaystyle \frac{{{\rm{\Lambda }}}_{x}}{{10}^{-22}}\right)}^{-1}\left(\displaystyle \frac{{\epsilon }_{e}}{{10}^{-2}}\right){\left(\displaystyle \frac{{\epsilon }_{b}^{u}}{{10}^{-2}}\right)}^{0.8}.\end{array}\end{eqnarray} \tag{ 12 }$$

Note that this decreases rapidly with radius but depends very weakly on the ISM density. It predicts that large-diameter SNRs should have much lower L_R/L_x ratios, regardless of the local ISM density (which may vary widely among remnants).

This appears to be in sharp conflict with the observations (Figure 21). There is a wide range of radio-to-X-ray ratios at all SNR diameters, and there is no indication of a marked decline for larger SNRs. The dashed lines show the model prediction from Equation (12). We have included in the model calculation the additional dependence on B and n_H that comes through the ${\epsilon }_{b}^{u}$ factor by using Equation (5) rather than assuming ${\epsilon }_{b}^{u}$ is constant, but even large changes in the ISM density and the B field scaling in Equation (6) have little effect on the general trend. Moreover, Elwood et al. (2019) found that the scatter in n_H for SNRs in M31 and M33 is relatively small ($\mathrm{log}{n}_{{\rm{H}}}\sim -1.35\pm 0.7$), which makes explanations that rely on a wide variation in the ISM density problematic. The predicted L_R/L_x ratios from the Sarbadhicary et al. (2017) model are too low at D ∼ 30 pc by one to two orders of magnitude compared with the observations, even using extreme parameter values.

Zoom In Zoom Out Reset image size

We have considered several possible explanations for the discrepancy in the L_R/L_x ratio predictions compared with the observations. Our calculation assumes that the SNR is in the Sedov phase. That is likely to be a good assumption because the radio and X-ray luminosities both decline dramatically at the onset of the late radiative phase (Asvarov 2006), while the early free-expansion phase is short (∼10² yr) and can be relevant only for small-diameter remnants (Elwood et al. 2019).

Another possible explanation is that the forced photometry SNR radio fluxes are contaminated by emission from nearby H ii regions, leading to elevated values of L_R/L_x. To test this hypothesis, Figure 21 shows (via the symbol colors) the spectral indices of the SNRs. If H ii contamination were important, all the large-diameter SNRs would have the thermal spectral indices typical of H ii regions, α ∼ 0, rather than the nonthermal α ∼ −0.5 indices observed for most SNRs. While there is an indication that some of the largest diameter SNRs have thermal contamination, there are also large-diameter SNRs with nonthermal spectra. Care is required to account for the overestimated radio luminosities for SNRs embedded in H ii regions, but that effect cannot explain the large L_R/L_x ratios in M33.

As shown in Figure 12, the luminosity function for the SNRs we detect can be represented as a power law, as was pointed out by Chomiuk & Wilcots (2009) for M33 (and a number of other galaxies) on the basis of earlier data. Chomiuk & Wilcots (2009) sought to interpret this fact in terms of a model that was similar to that of Sarbadhicary et al. (2017). While it may well be true that the radio luminosity function of SNRs can be represented as a power law, our analysis makes it clear that a physical connection between the observations and the theory is lacking.

Although we have chosen to focus on the specific model for radio emission in SNRs described by Sarbadhicary et al. (2017), this model embodies the same elements as most other models of radio emission in SNRs. We conclude that the L_R/L_x ratio is a strong diagnostic for testing models of radio emission from SNRs in the presence of ISM density variations. Apparently current models still have significant shortcomings when confronted by the radio and X-ray observations.

4.3.1. Large Magellanic Cloud

4.3.2. M31

4.3.3. M83

4.3.4. NGC 7793

4.3.5. Efficacy of SNR Identification Using Radio Emission

Most extragalactic SNRs and SNR candidates have been identified optically according to [S ii]:Hα emission line ratios. However, some attempts have been made to identify SNRs in other galaxies on the basis of their radio properties, usually coupled with more limited optical data. For example, Lacey & Duric (2001) identified 35 radio-selected SNRs/SNR candidates in NGC 6946 on the basis of their nonthermal radio spectral indices and associated Hα emission. These SNRs were concentrated in the spiral arms of NGC 6946. They compared this sample to one obtained by Matonick & Fesen (1997) optically and found very little overlap between the two samples. The Matonick & Fesen (1997) sample was distributed in both arm and interarm regions of the galaxy, and so Lacey & Duric (2001) argued that the fact that the samples were largely disjointed was a selection effect. Optical SNRs are simply difficult to pick out against the light of an H ii region. Most SNe occur in star-forming regions of a galaxy and hence most SNRs should exist in the spiral arms. Star-forming regions are also the regions of a galaxy with the most Hα emission, and therefore the fact that optical SNRs are not concentrated there is prima facie evidence that they are being missed in optical surveys.

The general argument made by Lacey & Duric (2001) that there are portions of a galaxy where optical SNR detection is difficult has to be correct at some level, but it raises a number of questions, especially as many of the spectral index measurements in radio maps are not precisely determined and because, with current instrumentation, it has been difficult to obtain any confirming evidence that a radio SNR candidate in an external galaxy based on radio criteria alone is indeed a SNR. On the theoretical side, these questions include how important it is that after the first SN explodes in a star cluster, subsequent SNe will explode into the region evacuated by the first SN. On the observational side, one of the questions is at what Hα surface brightness will what fraction of SNRs actually be missed, and how many background interlopers will be counted as SNRs. Unfortunately, there has very little systematic work to investigate any of these problems.

We can use the observations of M33 to gain some insight into this. There are 284 nonthermal sources (with a spectral index of less than −0.3) in our survey that are brighter than 0.1 mJy and are in the region covered by our 5 GHz data. Essentially all of those sources have well-determined spectral indices. Of these, 28 lie at positions where the Hα surface brightness exceeds 3 × 10⁻¹⁶ erg s⁻¹ arcsec⁻². Of these sources, 14 are known SNRs. However, most of the sources that are not known to be SNRs are quite faint. As shown in Figure 22, the number of interlopers drops quite rapidly as a function of the radio source brightness. At 200 μJy, the number of interlopers and known SNRs is about the same, and at 600 μJy essentially all of the interlopers have disappeared. At 200 μJy, the typical spectral index error is 0.05 (1σ) in our survey, so confusion with H ii regions is not a problem; nearly all of the non-SNRs are background sources. Our results are not very sensitive to the specific Hα surface brightness limit chosen, and the apparent distribution of fluxes of background sources does not depend on the distance to the foreground galaxy. This indicates that, as long as a radio SNR candidate is relatively bright (≳0.5 mJy), radio SNR identifications should be fairly reliable, as Lacey & Duric (2001) suggested on the basis of their location.

Zoom In Zoom Out Reset image size

The JVLA's combination of high spatial resolution and wide frequency coverage provides an excellent tool for surveys of M33 and other nearby northern galaxies. Future observations that are higher resolution might enable algorithms that separate SNR emission from H ii region emission using the morphology and spectral index of the radio flux. Higher angular resolution would also allow one to separate background objects from SNRs.

Improved algorithms for processing the very-wide-band JVLA data could also enable more information to be extracted from our current data. Although we have put considerable effort into development of new methods for analyzing multiresolution radio images, we have deliberately degraded the resolution to generate images that are matched across the full-frequency bandpass. One can imagine methods that might avoid this step.

Elsewhere, the Low-Frequency Array (LOFAR) radio telescope is now being utilized for low-frequency (0.15 GHz) surveys of the northern sky with a resolution of 6'' (Shimwell et al. 2019). While these data alone will not produce accurate spectral indices from its 0.12 to 0.17 GHz bandpass, it should be readily combined with our similar-resolution, higher-frequency JVLA observations. There are plans to do LOFAR observations of nearby galaxies including M33 to an rms sensitivity of 15 μJy at 0.15 GHz (van Haarlem et al. 2013); for a SNR with a spectral index −0.5, that corresponds to an equivalent rms ∼5 μJy at 1.4 GHz, which is also well-matched to our survey.

The Square Kilometer Array (SKA) pathfinder telescopes (the Australian SKA Pathfinder (ASKAP), MeerKAT) are located in the southern hemisphere and so are not directly relevant to studies of M33, but they should contribute interesting observations of other nearby galaxies. According to their specifications, they will generate very deep images with noise comparable to or better than that of our JVLA data (Johnston et al. 2008; Norris et al. 2011, 2013). Their lower resolution of 10''–15'' will make studies in the vicinity of H ii regions more challenging, and their initial configurations have a more limited frequency range than the JVLA. Nonetheless, with a flood of high-quality radio data on the horizon, we can expect dramatic advances in our understand of radio emission from SNRs over the next decade.

The median pyramid calculation begins with a median filter over a circular region of radius R = 8 pixels (8'' for the full-resolution image). The median-filtered image is the "smooth" channel, and the difference between the original image and the median-filtered image is the "sharp" channel. We also estimate the local rms noise as a function of position in the sharp channel by computing a second median filter, operating on the absolute value of the sharp channel and using an annular region with an inner radius of R pixels and an outer radius of 2R pixels. The median absolute value is multiplied by 1.4826, a factor that converts the value to an equivalent Gaussian sigma when the noise is purely Gaussian. (This robust noise estimation approach is insensitive to outliers in the image.)

To create additional levels in the pyramid, the process is then repeated using the smooth image, with the radius of the filtering region expanded by a factor of 2. This is iterated as many times as desired, with the filtering radius doubling at each new level. The choice of the number of levels determines the largest scale structure that remains in the image. For M33, we compute L = 4 levels of the pyramid.

The median-filtered smooth channel has only faint residuals remaining at the positions of sources that are small compared with the filter area. We further reduce the contamination by sources in the smooth channel by doing a second filtering pass after removing pixels that are found to be larger than twice the local rms estimate. Those pixels are replaced by interpolating values from nearby pixels (using yet another median filter). Then the smooth channel is computed again. This single iteration of source rejection leads to very good removal of sources from the smooth channel.

We also improve the algorithm's performance, which would be slow for higher levels of the pyramid when the filtering radius becomes large, by reducing the size of the smoothed image at each level by a factor of 2 in each dimension. We simply bin the smoothed image by averaging blocks of 2 × 2 pixels. To compute the sharp image, we then re-expand the smoothed image back to the original size using bilinear interpolation. The result is a stack of smoothed images of decreasing size. With this modification, the iterated smoothing uses a median filter region that is the same size (in pixels) as used for the initial iteration. The effective radius is increased by a factor of 2 because the smoothed image was shrunk by that factor. The result is that the algorithm is actually faster for the higher levels of the pyramid.

When the median pyramid is complete, we are left with a stack of L = 4 images (L − 1 sharp images and the remaining smooth background image) along with L − 1 local rms estimates for each sharp image. There are several adjustable parameters in the algorithm. The filter radius R should be slightly larger than the PSF FWHM so that the first sharp images includes the point sources. The number of levels L, combined with R, determines the scale of the largest sources in the lowest resolution sharp image. Larger structures are in the final smooth image and are treated as background.

A single background-subtracted image can be constructed by summing the L − 1 sharp channels. However, it is more effective to do source detection directly on the pyramid levels. At each level, we use the FellWalker clump-finding algorithm (Berry 2015) to identify islands of source emission. FellWalker identifies distinct peaks in the image and assigns nearby pixels to the nearest peak that is reachable by going uphill. It has parameters that prevent noise bumps from being included, that determine when two neighboring peaks have a significant valley between them, and that merge islands that are consistent with being noisy features of the same peak. The output of FellWalker is an island (or segmentation) map that divides the image into empty regions and contiguous regions at the locations of detected sources. We morphologically filter the FellWalker maps to eliminate islands that are smaller than expected given the resolution at that pyramid level.

There is one FellWalker island map for each level of the pyramid. We combine the independent island maps by merging islands that are detected in more than one level (which is a frequent occurrence). When necessary, islands that have multiple high-resolution peaks with overlapping low-resolution emission are split among the peaks using rules that attempt to maintain the morphological characteristics of the existing islands.

The merging process starts by creating a merged island map M from the highest resolution image. The island numbers at each higher level are then updated using these rules:

• 1.  
  Islands that do not overlap an existing island in M get a new number larger than any existing island number.
• 2.  
  Islands that overlap a single existing island in M are assigned the same number as that island. This is the most common case (e.g., when a core-halo source is detected at several levels of the pyramid).
• 3.  
  Islands that overlap more than one existing island in M are split and are assigned numbers chosen from among those islands:

  • (a)  
    Pixels already assigned in M to an island number retain their values.
  • (b)  
    Unassigned pixels are assigned to the closest island using a distance metric that weights the Euclidean distance to the island center by the elliptical moments of the island shape in M.
  • (c)  
    The extended islands are checked to ensure that they are contiguous. If noncontiguous regions are found, those pixels are reassigned to the second closest island.
• 4.  
  After renumbering is complete, M is updated by filling blank pixels with the corresponding values in the new level's island map.

There are other algorithms that could be used for island splitting (case 3), but our approach is simple and has been found to produce reasonable results. The effect of our choice is to assign a pixel that is near two different islands to the island that reaches most closely toward that pixel, either because it is large or because it is elliptical and is extended in the direction of the pixel.

There are several advantages to this source detection approach. The detection of sources at each level is simplified because more extended background structures have already been subtracted. Sources that are extended with low surface brightnesses are much more easily recognized in the lower resolution levels of the pyramid. The island merging process provides a clear and simple approach for combining the source lists at different resolution levels. Finally, having a multiresolution set of detection maps enables the computation of higher signal-to-noise ratio flux density estimates for sources with complex morphologies that include sharp components embedded in broad halos. That algorithm is described below.

It is possible to combine the multiresolution detection maps into a single map that assigns every pixel of the image to a particular source (or to the empty background). We use that map, for example, to display the locations of detected sources. In some applications it might be sufficient to sum the fluxes within islands in order to determine source brightnesses. In our images, however, that leads to relatively high noise for core-halo sources, which have a sharp component surrounded by a lower surface brightness halo that can spread across many pixels.

To reduce the noise, we instead use the segmentation maps at each level of the image stack along with the image from that level of the stack. For each level, we sum only the pixels that are detected in the island at that level. For a sharp core embedded in a broad halo, the sharp channel island ("level 0") only includes the few pixels needed to cover the PSF, but the broad channel island ("level 1") with the halo includes many more pixels. The flux for the object is the sum of the few pixels from level 0 plus the many pixels from level 1. It does not include all the level 0 pixels that are covered by the level 1 island. That decreases the noise in the sum by not adding in a bunch of empty pixels that contribute only noise. Moreover, the level 1 image is significantly smoothed, reducing the rms noise in that level. As a result, the pixels that are included over a large region for the smooth halo contribute significantly less noise to the sum.

Note also that sources that have only a sharp component (e.g., point sources) typically have islands that have no components at any of the lower resolutions. Omitting those pyramid levels from the sum also reduces the noise for compact sources.

The end result is that this approach effectively creates an adaptive set of weights for the image pixels that adjusts itself depending on the actual distribution of flux in the object.

Figure 23 shows a comparison between the simple (single level) island flux densities and our multiresolution island flux densities (left panel) and a comparison of the fractional flux density errors σ(F)/F (right panel). Sources are color-coded using the island flux. There is a bias in the measured flux densities, with the multiresolution values being systematically lower for fainter sources. This is not unexpected because the multiresolution islands are more tightly delineated, and some flux spills outside the summed regions. On the other hand, the right panel shows that the signal-to-noise ratio in the measurements is about 20% better in the multiresolution measurements. The fractional flux density errors are systematically lower for both bright and faint sources using the multiresolution algorithm.

Zoom In Zoom Out Reset image size

Figure 24 shows the effect of the two algorithms on the spectral indices. The flux densities are computed using exactly the same islands for the (resolution-matched) 1.4 and 5 GHz images. Each band is divided into two frequency channels, and a power law is fitted to the four channels using a constrained algorithm that requires the spectral index α (${F}_{\nu }\propto {\nu }^{-\alpha }$) to lie in the range −3 ≤ α ≤ 3. The variances and covariances of the fitted parameters are derived from the noise in the individual bands. Figure 24(a) shows the spectral index from the simple island fluxes (x-axis) versus the index from the multiresolution fluxes (y-axis). The sample of sources was restricted to objects with σ(α) ≤ 0.1 in the multiresolution measurement. As in Figure 23, the points are color-coded by flux. Clearly there is excellent agreement, particularly for the brighter sources. There is no significant bias between the two spectral index estimates.

Zoom In Zoom Out Reset image size

Figure 24(b) compares the noise in the spectral index for the two algorithms. Here the advantage of the multiresolution islands can be seen clearly, with the median noise being about 30% lower compared with the simple islands.

Our conclusion is that the multiresolution islands lead to significantly higher signal-to-noise ratio measurements of the flux densities (although with some bias toward lower values) and also to more accurate estimates of the spectral indices. The algorithm is robust, giving reliable results even in very crowded regions of the image. We find that it has significant advantages compared to the other algorithms that we tried.

We foresee several possible areas for further improvement. It might be effective to use a continuous weighting scheme for pixels rather than the simple in-or-out approach used here. Pixels near the edges of islands could get a lower weight. Pixels that have some but not much flux could be weighted down. A better weighting algorithm would probably reduce the bias and might ideally be able match Gaussian fits in signal-to-noise ratio while retaining the ability to robustly determine source flux densities for complex objects that are poorly measured using Gaussian models.

Another idea would be to perform multiresolution Gaussian fits at each level of the image stack. That might have a couple of benefits: (1) the fitting process would converge more easily for the sharp channel because sources are better separated, and (2) the resulting image could model the actual flux distribution much more accurately using a sum of multiple Gaussians, while avoiding the difficulties of simultaneously fitting overlapping Gaussians.