EMERGING MASSIVE STAR CLUSTERS REVEALED: HIGH-RESOLUTION IMAGING OF NGC 4449 FROM THE RADIO TO THE ULTRAVIOLET

Radio observations of NGC 4449 were obtained with the NRAO VLA from 2001 April to 2002 November. Observations using the C-band (6 cm, 5 GHz) were carried out with the A-array, observations using the X-band (3.6 cm, 8 GHz) were carried out with both the A- and C-array, and observations using the K-band (1.3 cm, 23 GHz) were carried out with the C-array. The sources 3C48 and 3C286 were used as flux density calibrators, and we estimate that the absolute flux density calibration is ≲5%, based on the scatter in the VLA Flux Density Calibrator database. The high-frequency observations at the K-band utilized fast-switching to a nearby phase calibrator, displaced by ∼5°, in order to mitigate the effect of the atmosphere.

These data were reduced and calibrated using the Astronomical Image Processing System (AIPS). In order to better match the largest spatial scale to which the observations at each wavelength is sensitive, antenna spacings of <20 kλ were not used. Additionally, in order to achieve synthesized beams of approximately the same size at each wavelength, the three data sets were imaged with different weighting of the uv-coverage, ranging from purely natural weighting at 6 cm to slightly uniform weighting at 1.3 cm. The data were convolved to identical beams of 13 × 13 in the imaging process to better facilitate comparison between wavelengths. Despite our best efforts, it should be noted that it is not possible to precisely match the uv-coverage at each wavelength, and there are slight variations in sensitivity to different spatial scales. After imaging the data, a primary beam correction was applied at each wavelength. The resulting parameters are listed in Table 1.

We identify 39 compact radio sources toward NGC 4449 which are shown in Figure 1 and labeled in order of increasing right ascension (R.A.). Source identification is based on a minimum 3σ detection (locally) in the 3.6 cm image, which has the highest sensitivity of the three radio images. Figure 2 shows an archival HST F550M (∼V-band) image in yellow overlaid with our VLA 3.6 cm image in blue.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A wealth of data on NGC 4449 is contained in the HST⁴ archive and we have obtained broadband and narrowband images of the galaxy spanning ultraviolet to near-IR wavelengths. These observations were taken over a period of several years with the Advanced Camera for Surveys (ACS), the Wide Field and Planetary Camera 2 (WFPC2), and the Near Infrared Camera and Multiobject Spectrometer (NICMOS). The HST observations have spatial resolutions ranging from ∼01 to 02. A summary of the HST observations used in this project is given in Table 2.

The ACS Wide Field Channel (WFC) has a field of view of 202''× 202'' and a plate scale of ∼005 pixel⁻¹. The ACS pipeline employs the MultiDrizzle (Koekemoer et al. 2002) software to produce cosmic ray cleaned, combined images and we have obtained these data on NGC 4449 through the F435W, F550M, F814W, F658N, and F660N filters. These observations provide coverage of all of the locations of the radio sources in NGC 4449 identified in Figure 1.

WFPC2 is composed of four CCDs: three Wide Field (WF) chips and the Planetary Camera (PC). The plate scales of the PC and WF chips are ∼005 pixel⁻¹ and ∼01 pixel⁻¹, respectively. The WF chips provide an L-shaped field of view of 150''× 150'' and the PC fills in another 35''× 35''. All but two of the radio sources (6 and 29) in NGC 4449 have WFPC2 coverage through the F170W and F336W filters. Two exposures were taken through each of these filters and we combined the calibrated images, chip by chip, and rejected cosmic rays.

A snapshot NICMOS observation through the F160W filter was obtained with the NIC 3 Camera and we have obtained this calibrated image from the archive. This image, with a field of view of 51''× 51'' and a plate scale of ∼02 pixel⁻¹, only covers seven radio sources (15, 16, 17, 19, 23, 24, and 25) in the nuclear region of the galaxy.

We register all of the HST images to match the astrometry of the VLA images, which is accurate to ∼01. Coincident sources (typically 5–10) are found in pairs of images and the mean pixel offsets are found between the two images. We begin by registering the F658N (Hα) image to the 3.6 cm radio continuum image since both observations trace H ii regions. The other HST images can then be registered with respect to this or subsequent registered images. The HST images are aligned with respect to the VLA images to a precision of ∼02.

We have obtained IR images (3.6–24 μm) of NGC 4449 from the Spitzer Space Telescope archive. We primarily use these data in this paper to qualitatively investigate the IR properties of NGC 4449 since these images suffer from significantly lower resolution (∼3–6'') than the VLA and HST images. However, despite their limitations, the IR data are informative, providing the general IR morphology of NGC 4449 and highlighting regions associated with warm dust.

The high-quality data used in this study of embedded massive clusters allows for a detailed inspection of the starburst morphology of NGC 4449. Figure 3 shows a selection of images of the galaxy at multiple wavelengths. "Close-up" three-color images (3.6 cm, Hα and V-band) are shown in Figure 4 with the sources labeled according to Figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The similarities between the 3.6 cm and Hα observations are immediately apparent in Figure 3. This likeness is not surprising since both observations trace dense ionized gas. The 3.6 cm and 24 μm images are also strikingly similar, indicating that the presence of warm dust is almost exclusively associated with the star-forming regions in the galaxy.

The circles in Figure 3 are centered on two example radio sources: one with and one without counterparts in the other wavebands. Source 30 (to the NE) has both optical and IR counterparts, suggesting this stellar cluster is in the process of emerging from its birth material, all the while ionizing the surrounding gas and heating the dust. Source 28, on the other hand, does not have detectable counterparts in any of the other wavebands. We will return to the radio sources without optical or IR counterparts in Section 5.4.

In order to study the emergence of young massive star clusters, we first use radio continuum imaging to identify potential clusters surrounded by ultra-dense H ii regions (see Figure 5 for some examples). We then measure the ultraviolet to radio flux densities of the radio-selected sources (in the VLA and HST images) using a new photometry code that we developed (SURPHOT), allowing for consistent irregular apertures and background annuli across multiple wavebands. SURPHOT is first used to detect a specified contour level in a selected region of the 3.6 cm image which is used as a reference. A source is selected and the contour may be used for the aperture in all of the images. Alternatively, an ellipse, a circle, or a free-form aperture can be used if appropriate (e.g., nearby yet distinct sources that are contained within the same contour). A background annulus is constructed by radially expanding the chosen aperture by specified factors (typically 1.75× and 2.5×). The pixel coordinates of the aperture and background annulus are transformed between images using their FITS World Coordinate Systems (WCS) which relate image pixels to sky coordinates. In this way, flux densities are computed in all (11) of the VLA and HST images using identical apertures and background annuli. (A more detailed description of SURPHOT is provided in the appendix.) The resulting flux densities are used in conjunction with stellar population synthesis models to derive physical properties such as extinctions, ages, and masses.

Zoom In Zoom Out Reset image size

Our choice to use radio-selected apertures is both scientifically and practically motivated. The focus of this work is on extremely young massive clusters and centimeter radio emission is a good indicator of youth. In addition, the resolution of the VLA images (13) is larger than the resolutions of the HST images (∼01–02). We note that using a larger aperture than the intrinsic size of a stellar cluster (in the HST data) allows for the possibility of contamination from an underlying or surrounding stellar population within the apertures. However, this effect is not likely to be very significant since the stellar light is dominated by the clusters of interest and we perform background subtraction in the photometry.

Thermal radio sources are recognized as having flat or positive radio spectral energy distributions (SEDs; S_ν ∼ ν^α,  α ≳ 0), indicative of free–free emission from dense ionized gas associated with young massive star clusters. Non-thermal sources have decreasing radio SEDs characteristic of supernova remnants and (background) active galactic nuclei (AGNs). Here we are concerned with the thermal sources in NGC 4449, which must be disentangled from the non-thermal sources by calculating the spectral indices, α^6.0 cm_3.6 cm and α^3.6 cm_1.3 cm, using the measured flux densities.

We use the 3.6 cm radio image as our reference for both source identification and as the reference image in SURPHOT. The 1.3 cm image is somewhat less likely to suffer from non-thermal contamination or self-absorption but the sensitivity is lower than the 3.6 cm image by a factor of ∼2 at best (the image quality is degraded farther from the center of the 1.3 cm image). Source selection is based on a 3σ detection threshold, where σ is the local rms noise (see Figure 1). We estimate the minimum uncertainty in the measured radio flux densities as this local rms noise. In many cases, however, the dominant source of uncertainty comes from the choice of aperture. For faint sources, estimates of the background level can also lead to large uncertainties. We therefore make multiple flux measurements for each source, varying these parameters. The fluctuations in these measurements are accounted for in our final reported values of flux densities and spectral indices, which are given in Table 3 along with the coordinates of the 39 radio sources identified in NGC 4449. In the cases of weak or non-detections in the 6.0 and 1.3 cm images, 3σ detection limits are given for the flux densities. Corresponding limits are given for the spectral indices. Sources are designated as either (likely) thermal, non-thermal, or mixed/uncertain based on their spectral indices as described above. Of the 39 radio sources, 13 are thermal, 22 are uncertain or composed of mixed sources, and 4 are non-thermal.

Of the four definitive non-thermal sources, three appear to be associated with the galaxy (Sources 13, 14, and 39) while Source 6 is most likely a background AGN with no optical counterpart. Source 14 is a well-studied, young (60–200 years) luminous supernova remnant (SNR) (J1228+441) (see Lacey et al. 2007, and references therein). Lacey et al. (2007) have given a modern summary of the 30-year history of the radio light curve of J1228+441. Source 13 may also be an SNR, although it appears to be at least partially embedded in an H ii region. Source 39 is quite intriguing. The radio emission is partially associated with an obvious point-like optical counterpart (probably an SNR) within a larger bi-conical-shaped H ii region (see Figure 4). The non-thermal sources within NGC 4449 are interesting in their own right, but they are excluded from further analysis in this work.

We use the 3.6 cm radio map as the reference image in SURPHOT (see the appendix) to obtain photometry of the radio sources in the HST data. HST system magnitudes with Vega zeropoints are reported in Table 4. Error estimates include contributions from Poisson noise and uncertainties from varying the apertures (see Section 3.2).

In order to estimate the physical properties of the radio-selected sources in NGC 4449, we compare our data to the latest STARBURST99 population synthesis models (Version 5.1) of Leitherer et al. (1999). We have run simulations at two metallicities, Z = 0.004 and Z = 0.008, as appropriate for NGC 4449 (Lequeux et al. 1979). We adopt an instantaneous burst of 10⁴ M_☉ with a Kroupa IMF, the Geneva evolutionary tracks with high mass loss, and the Pauldrach/Hillier atmospheres. Ages, masses, and extinctions are estimated using two distinct methods: (1) comparing the nebular emission from the H ii regions to the model output (Section 4.2) and (2) fitting model SEDs to the broadband HST flux densities (Section 4.3).

4.2.1. Ages

We estimate ages of the radio sources in NGC 4449 by comparing their measured Hα equivalent widths with the model predictions as a function of age. For an instantaneous burst of star formation, the equivalent width of Hα emission is a reliable age indicator for young clusters since it measures the ratio of the ionizing flux (from massive stars) to the total continuum flux (e.g., Alonso-Herrero et al. 1996). This ratio is strongly dependent on age for stellar populations ≲10 Myr, so long as the most massive stars have begun to die, thereby decreasing the ionizing flux. In addition, if the extinction is the same to the stars and the gas, Hα equivalent width will be independent of extinction.

We calculate Hα equivalent widths, W_Hα, using Equations (1) and (2):

$$\begin{equation} W_{{\rm H}\alpha } = {F_{{\rm H}\alpha } \over f_{\rm cont}^{({\rm H}\alpha)}}, \end{equation} \tag{ 1 }$$

where F_Hα is the total Hα flux (erg s⁻¹ cm⁻²) and f^(Hα)_cont is the continuum flux density (erg s⁻¹ cm⁻² Å⁻¹) at 6563 Å. The F658N filter contains both Hα and [N ii]λ6584 emission so the contribution from the [N ii] line must be subtracted. F_Hα is given by

$$\begin{eqnarray} F_{{\rm H}\alpha } &=& \big[\big(f_{\rm F658N} - f_{\rm cont}^{({\rm H}\alpha)}\big) \Delta \lambda _{\rm F658N}\big]\nonumber\\ && - \big[\big(f_{\rm F660N} - f_{\rm cont}^{({\rm N}\,{\rm II})}\big) \Delta \lambda _{\rm F660N}\big], \end{eqnarray} \tag{ 2 }$$

where the first term in square brackets is the total Hα + [N ii]λ6584 flux and the second term in square brackets is the total [N ii]λ6584 flux. Differences in throughput of the filters at the wavelengths of Hα and [N ii] are accounted for in Equation (2). The measured flux densities and widths of the filters are given by f_F658N, f_F660N, and Δλ_F658N, Δλ_F660N, respectively. The flux densities of the continuum values, f^(Hα)_cont and f^(N II)_cont, are found by interpolating between the F550M and F814W data.

We compare the calculated Hα equivalent widths, W_Hα, to the models in order to estimate ages for the radio sources in NGC 4449. The thermal radio sources in NGC 4449 have Hα equivalent widths between ∼340 and ∼2460 Å, implying they have ages ≲5 Myr. However, it should be noted that the STARBURST99 models reflect a fully sampled stellar IMF. For clusters with masses ≲10⁶ M_☉, it becomes increasingly unlikely that a given cluster will contain extremely massive stars. As a result, lower mass clusters will tend to produce less ionizing flux for the same age than would be predicted by simply scaling the STARBURST99 predictions. In addition, the possibility of contamination from an underlying or surrounding stellar population to the continuum measurements as a result of using radio-selected apertures (see Section 3.1) could artificially decrease the measured Hα equivalent widths. Therefore, the radio-detected clusters may be even younger than derived here (Table 5).

Table 5. Properties of the Thermal Radio Sources

Source	Log WHα (Å)	Age (Myr)	QLyc (1049 s−1)	Mass (103 M☉)	AV (mag)
2	>3.10	<3.7	17(3)	<7	4.3(0.3)
3	3.00(0.28)	(1.2)4.0(0.8)	40(4)	21(2)	0.3(0.1)
4	2.68(0.07)	(0.2)4.9(0.1)	39(5)	40(5)	0.2(0.1)
5	2.73(0.16)	(0.4)4.8(0.3)	19(4)	18(3)	1.5(0.3)
7	3.12(0.22)	(0.6)3.3(1.0)	18(3)	6(1)	1.0(0.2)
11	2.53(0.50)	(2.1)5.2(2.2)	33(4)	44(5)	1.5(0.4)
18	3.26(0.07)	(0.2)2.9(0.1)	21(2)	6(1)	0.1(0.1)
22	2.79(0.03)	(0.1)4.6(0.1)	23(2)	19(2)	0.7(0.1)
23	3.21(0.08)	(0.1)2.9(0.7)	68(3)	19(1)	1.2(0.1)
27	3.13(0.14)	(0.8)3.6(0.5)	42(3)	16(1)	0.8(0.1)
29	3.20(0.58)	(2.9)2.9(2.1)	23(4)	6(1)	1.0(0.3)
30	3.39(0.09)	(1.0)2.6(0.2)	27(2)	6(1)	1.2(0.1)
33	3.24(0.48)	(2.9)2.9(1.8)	18(3)	5(1)	1.2(0.2)

Notes. Hα equivalent widths, ages, ionizing fluxes, masses, and extinctions of the thermal radio sources in NGC 4449. These values were derived from the nebular emission from the H ii regions (Section 4.2).

Download table as:  ASCIITypeset image

4.2.2. Ionizing Fluxes

The ionizing flux of a starburst region, Q_Lyc, can be estimated from its radio luminosity density following Condon (1992):

$$\begin{eqnarray} \left({Q_{\rm Lyc} \over {\rm s^{-1}}}\right) &\gtrsim& 6.3\times 10^{52} \left({T_e \over 10^4\ {\rm K}}\right)^{-0.45} \left({\nu \over {\rm GHz}}\right)^{0.1}\nonumber\\ && \times\, \left({L_{\nu, \rm thermal} \over 10^{27}\ {\rm erg\ s^{-1}\ Hz^{-1}}}\right). \end{eqnarray} \tag{ 3 }$$

Equation (3) is applicable to thermal radio emission but provides only a lower limit since some of the ionizing photons may be absorbed by dust and the radio emission may be optically thick. Studies of Lyman continuum extinction in H ii regions find that the fraction of ionizing photons absorbed by dust can be quite large (Inoue 2001; Dopita et al. 2003), leading to derived ionizing fluxes (and therefore masses) that can be artificially low by up to a factor of 10. Nonetheless, we adopt an electron temperature of T_e = 10⁴ K and use the 3.6 cm radio luminosities to calculate Q_Lyc using Equation (3) for the radio sources in NGC 4449. The 1.3 cm radio luminosities are somewhat less likely to be contaminated by non-thermal emission or be optically thick, but the 1.3 cm image does not have the requisite sensitivity. The derived Q_Lyc of the thermal radio sources range from ∼20–70 ×10⁴⁹ s⁻¹, or, equivalently (Vacca et al. 1996), ∼20–70 O7.5 V stars.

4.2.3. Masses

4.2.4. Extinctions

Global extinction of the ionized nebular gas is estimated by comparing the measured Hα flux to the predicted Hα flux (with no extinction) derived from the thermal radio flux density, S_ν,thermal. The radio emission does not suffer from extinction per se, although again we note that the UV photons can suffer significant absorption before ionizing the gas, which is not accounted for in Equation (3). We follow Condon (1992) to obtain the predicted Hα flux:

$$\begin{eqnarray} \left({F_{{\rm H}\alpha,{\rm predicted}} \over {\rm erg\ s^{-1}\ cm^{-2}}}\right) &\sim& 0.8 \times 10^{-12} \left({T_e \over 10^4\ {\rm K}}\right)^{-0.59} \left({\nu \over {\rm GHz}}\right)^{0.1}\nonumber\\ && \times\, \left({S_{\nu, \rm thermal} \over {\rm mJy}}\right). \end{eqnarray} \tag{ 4 }$$

The extinction of Hα is then given by

$$\begin{equation} A_{{\rm H}\alpha } = -2.5 \times {\rm log}\left({F_{{\rm H}\alpha,{\rm observed}} \over F_{{\rm H}\alpha,{\rm predicted}}} \right). \end{equation} \tag{ 5 }$$

A 30 Doradus extinction curve (Fitzpatrick 1985; Fitzpatrick & Massa 1990; Misselt et al. 1999) is used to convert A_Hα to A_V, the extinction at 5500 Å. The thermal radio sources have measured nebular extinctions of A_V ≲ 1.5 with the exception of Source 2, which is highly obscured and has A_V ≈ 4.3. Measured extinctions, Hα equivalent widths, ages, ionizing fluxes, and masses of the thermal radio sources in NGC 4449 are listed in Table 5.

In addition to using the nebular emission, we also estimate the physical properties of the radio-selected sources in NGC 4449 using the stellar emission by comparing a grid of model SEDs to the broadband HST flux densities. The grid consists of models with ages up to 10 Myr in increments of Δt = 0.1 Myr and extinctions in increments of ΔA_V = 0.1. We typically run the grid from A_V = 0–5 mag which yields a total of 5,151 model SEDs to compare to each observed SED. Galactic foreground extinction is also accounted for (NGC 4449 has a Galactic reddening of E(B − V) = 0.019 (Schlegel et al. 1998)).

Each model SED, with a given age and extinction, is convolved with the HST total throughput curves corresponding to the broadband filters listed in Table 2 before it is compared to the observations. Once we have obtained the model flux densities for each broadband filter, we compute a goodness-of-fit statistic, σ, given by the standard deviation of the difference between the log of the observed flux densities and the log of the model flux densities, weighted by the measurement errors of the observed flux densities. In the case where we only have an upper limit for an observed flux density, we use this limit to exclude models but do not include this point in the fit. The model that minimizes σ is chosen as the best fit. A mass estimate is obtained by scaling the input model mass by the weighted mean offset between the observed and best-fit model flux densities (again done in the log).

We have performed SED fitting on the radio-selected sources in NGC 4449 with detections in the HST data. The most heavily obscured sources are not fit since we only have upper limits for all of the flux densities. The strongest Balmer lines, Hα and Hβ, are not contained in any of the broadband filters used here, although the F435W filter does include Hγ and Hδ. Originally, we included up to five filters in the fitting process: F170W (∼UV), F336W (∼U), F435W (∼B), F550M (∼V), and F814W (∼I). However, we found that even the best-fitting models did not visibly fit the data in many cases because of an apparent excess in the I-band. To check there were no calibration problems in the F814W data, we obtained WFPC2 F814W images of NGC 4449 from the archive. The results were consistent with our findings in the ACS data.

To investigate further, we ran our SED fitting routine 12 times on each possible source: including and excluding the F814W data point, using models with metallicities of Z = 0.004 and Z = 0.008, and applying three different extinction curves: the Galactic extinction curve (Cardelli et al. 1989), the starburst obscuration curve (Calzetti et al. 2000), and a 30 Doradus extinction curve adopted from Misselt et al. (1999, Table 3) and Fitzpatrick (1985, Table 6), using the parameterization given by Fitzpatrick & Massa (1990). For each of the 12 runs, we have computed the mean σ (goodness-of-fit statistic) of all the sources and its standard deviation (Table 6). These values measure how well a family of models (metallicity, extinction curve, including/excluding I-band data) fit the data set as a whole. Table 6 shows that excluding the I-band data and using the 30 Doradus extinction curve fits the data significantly (∼20–50%) better than the other families of models, but that there is no significant difference between the fits using either metallicity.

It is also evident in Table 6 that excluding the I-band data significantly improves the fits, regardless of which extinction curve is used. A better fit is expected if using fewer data points (i.e., excluding the I-band); however, we ran Monte Carlo simulations which showed that including one more data point in the fits should only increase σ by ∼2–5%, not 33–91% as seen in Table 6. Data were simulated by adding noise (representative of measurement errors) to model SEDs and then convolving the SEDs with the appropriate throughput curves. The simulated data were compared to the original SEDs and the goodness-of-fit statistic, σ, was calculated with and without the I-band data.

It should be emphasized that the extinction of the stellar continuum cannot account for the observed SEDs when including the I-band data. Figure 6 shows a plot of B − V (F435W–F550M) versus V − I (F550M–F814W) illustrating that the reddest clusters in V − I are the bluest in B − V, contrary to the expected trend if extinction were responsible for the red V − I values. The model evolutionary track and a reddening vector (using the 30 Doradus extinction curve) are also shown for comparison. This anomalous color–color diagram, along with our results from SED fitting, leaves little doubt that the I-band excess is real.

Zoom In Zoom Out Reset image size

On average, the I-band flux densities are ∼60% higher (0.5 mag brighter) than the predicted values from the best-fitting model SEDs. In the most extreme cases, the measured I-band flux densities exceed the model values by more than a factor of 2.5. An even greater excess is present in the F160W (H) filter (also not included in the fits) with the data brighter than the best-fitting SED by an average of 1.6 mag, or ∼4× the flux. We will discuss possible origins of the I- and H-band excesses in Section 5.1.

Figure 7 shows the observed flux densities and best-fitting model SEDs for the thermal sources that have detections in the HST data. The models have a metallicity of Z = 0.004, use the 30 Doradus extinction curve, and exclude the I- and H-band data in the fits. The grids of σ are also shown to give a sense of the uncertainty in the estimates of ages and extinctions. A green cross indicates the best-fitting model and, for comparison, an orange cross indicates the values obtained from the nebular emission as described in the previous section. Overall, there is good agreement between the physical properties derived from the nebular emission and SED fitting.

Zoom In Zoom Out Reset image size

By combining results from the nebular and stellar emission, a self-consistent picture relating the observed properties of young emerging massive star clusters to model predictions is beginning to develop. It is noted that caution must be applied when interpreting the physical properties of observed clusters derived from population synthesis models, since stochastic fluctuations can bias the results. However, we do not believe that this is significantly affecting the results presented in this paper since the cluster masses are ∼10⁴ M_☉ and we are not solely relying on broadband colors (Lançon & Mouhcine 2000; Cerviño & Luridiana 2004). In this section, we combine the results from the nebular (Section 4.2) and stellar (Section 4.3) emission and explore various trends in the data.

Figure 8 shows the ratio of radio flux densities (3.6 cm) to the optical (V-band) flux densities plotted against Hα equivalent width. The ionizing flux, and therefore Hα equivalent width, are expected to decrease as a cluster ages. However, the ionizing flux is also a function of a cluster's mass whereas Hα equivalent width is roughly mass independent. Normalizing the radio flux densities (which traces the ionizing flux) by the optical flux densities removes the mass dependence. It is clear that the mass normalized radio flux densities and Hα equivalent widths do indeed decrease with cluster age. A model evolutionary track from 1–6 Myr is also shown in Figure 8. The model flux density at 3.6 cm was obtained from inserting the predicted ionizing flux from starburst99 into Equation (3). The vertical offset between the data and the model evolutionary track can be attributed to extinction in the V-band (top of Figure 8). The data match the model track extremely well after correcting for the extinctions derived from the nebular emission (bottom of Figure 8). This result suggests that the extinction to the stars and the gas is approximately the same, which is significant since the stellar light and nebular gas may have different geometries (i.e. the nebular emission is typically more extended than the stellar emission).

Zoom In Zoom Out Reset image size

Figure 9 also indicates that the extinction to the stars and the gas is roughly the same in most cases. The majority of the thermal radio sources have consistent extinction estimates for the stars and the gas, although approximately one-third of the sources have nebular extinctions that are ∼1 mag higher than their derived stellar extinctions. Based on the results from Figure 8, these differences may simply be due to uncertainties in the extinction estimates from SED fitting, which are largely dependent on the ultraviolet (and lowest-quality) data. The SEDs also suffer a degeneracy between age and extinction.

Zoom In Zoom Out Reset image size

However, it is worth mentioning here that discrepancies in extinction estimates of extragalactic H ii regions have been recognized since the 1970s. In particular, early Westerbork data for H ii regions in M33, M101, and M51 (Israel & van der Kruit 1974; Israel et al. 1975; van der Kruit 1977) showed that comparing radio continuum emission to Hα implied about 1.2 mag excess extinction as compared to the Balmer decrement procedure. These results combined with other determinations were summarized by Israel & Kennicutt (1980), who discuss two possible origins for the "missing dust": (1) the dust may be internal and mixed with the ionized gas, (2) over large sizes of many hundreds of parsecs, a non-uniform distribution of gas and dust will cause the Balmer decrement to be mainly weighted toward regions with lower extinction, while the radio-derived absorption will be determined mainly by the more absorbed regions. It is unclear whether the discrepant extinctions derived for four of the thermal radio sources in this work are a result of these types of effects, or simply due to uncertainties in the SED-derived extinctions.

The age estimates from SED fitting also have relatively high uncertainties (see Figure 7) and tend to be concentrated around ∼3 Myr. These results are probably due to an insensitivity of the model SEDs, which do not change significantly at these young ages. Hα equivalent width, however, is a more sensitive age indicator at young ages since ionizing flux is a strong function of time once the most massive stars begin to die (although these age estimates may be systematically high for reasons discussed in Section 4.2.1). A comparison of the age estimates from SED fitting and Hα equivalent widths is shown in Figure 10.

Zoom In Zoom Out Reset image size

We find two interesting correlations with Hα equivalent width, and hence cluster age. First, the I-band excess (Section 4.3) is largest in the youngest clusters and decreases with age. We define the I-band excess in magnitudes as ΔI = −2.5 × log(f^(I)_observed/f^(I)_model), where f^(I)_observed is the observed flux density in the F814W filter and f^(I)_model is the flux density of the best-fitting model SED. Figure 11 shows a plot of ΔI versus age for the thermal radio sources and the best linear fit is given by

$$\begin{equation} \left({\Delta I \over {\rm mag}}\right) = 0.4 \left({t \over {\rm Myr}}\right) - 2.1. \end{equation} \tag{ 6 }$$

Setting ΔI equal to zero implies that whatever is causing the I-band excess in these young clusters is unlikely to affect clusters older than ∼5 Myr. We will discuss possible origins of the I-band excess in Section 5.1.

We also observe a puzzling trend that seems to indicate cluster mass is positively correlated with age for ages ≲5 Myr. Figure 12 shows a plot of the extinction-corrected M_V versus Hα equivalent width. The model evolutionary tracks for various cluster masses are also shown. A strong anti-correlation between M_V and Hα equivalent width is evident, implying more massive clusters tend to be older than less massive clusters. The same overall trend is observed without correcting M_V for extinction, while a positive correlation is found between M_V and Hα total flux (as expected). An observational bias could possibly explain why optically faint clusters with low Hα equivalent widths are not detected, since these clusters may be less likely to show up in the radio and make it into our sample. However, it is curious that we do not detect clusters that are relatively bright in V with high Hα equivalent widths. It is not clear what is causing this trend in the data, but we explore various possibilities in Section 5.2.

Overall, the mass estimates derived from the nebular emission and SED fitting are consistent within a factor of ∼3. However, a slight trend is evident in Figure 13 suggesting that mass estimates from the nebular emission tend to be higher than those from SED fitting for the lower-mass clusters, while the mass estimates from SED fitting tend to be higher than those derived from the nebular emission for the higher-mass clusters. We will return to this in Section 5.2

We have demonstrated that the I-band fluxes of the radio-selected clusters in NGC 4449 exhibit an excess compared to model SEDs (Figure 7). In the B − V versus V − I color–color diagram, the reddest clusters in V − I are simultaneously the bluest in B − V (Figure 6). The I-band excess is also positively correlated with Hα equivalent width, implying the excess decreases with age (Figure 11). Based on these observations, we do not expect the origin of this I-band excess to affect clusters older than ∼5 Myr (Equation (6)).

Zoom In Zoom Out Reset image size

Possible origins of the I-band excess are discussed below, but first we note that other studies have also revealed the presence of an I-band excess in a small number of YMCs. Whitmore & Schweizer (1995) discovered about a dozen extremely red objects in the SE quadrant of the merging Antennae galaxies. This region also contains the strongest radio continuum sources in the system as well as strong mid-IR emission and a few of these sources are probably associated with the red objects. Whitmore & Schweizer (1995) interpreted these red objects as extremely young, embedded clusters that have not yet cleared out their natal dust cocoons.

In their study of the circumnuclear starburst galaxies NGC 5248 and NGC 1512, Maoz et al. (2001) found a number of YMCs with excess emission in the H-band, and one of these sources had a detectable I-band excess. Maoz et al. (2001) were able to fit the SED of this source by adding a 2000 K blackbody to their model and determined that if the emission were due to dust heated by stars, it must be circumstellar.

Possible origins of the observed I-band excess in the young radio-detected clusters of NGC 4449 are considered here.

5.1.1. Red Supergiants

5.1.2. Emission Lines

5.1.3. Thermal Emission from Hot Dust

5.1.4. Optically Veiled Photospheric Continuum from Embedded Stars

5.1.5. "Extended Red Emission"

Figure 12 illustrates that while M_V and Hα flux are positively correlated (as expected), M_V and Hα equivalent width are anti-correlated. In other words, fainter clusters have relatively more Hα emission with respect to the stellar continuum than brighter clusters. This can be interpreted as a mass–age correlation, since M_V generally traces cluster mass (modulo extinction) and Hα equivalent width traces age. The two most trivial, and unsatisfying, interpretations of this trend are that (1) the low-mass clusters formed more recently than the high-mass clusters, and (2) smaller clusters coalesce and make larger clusters over time. However, we note that there is no observed correlation between cluster size (measured by aperture size) and Hα equivalent width.

Zoom In Zoom Out Reset image size

A third interpretation of the apparent mass–age correlation for clusters in this work (with ages ≲5 Myr), although purely speculative, is the following. It is possible that the observed correlation is actually the imprint of an evolutionary sequence, where the youngest clusters merely appear less massive than their older counterparts. First, relative to more evolved clusters, younger clusters are likely to contain a higher fraction of dust enshrouded, optically veiled stars. These heavily embedded stars would not contribute to the integrated V-band flux, yet the young optically visible stars within the clusters would still produce large Hα equivalent widths revealing their young ages. Therefore, younger clusters may tend to appear less massive. In addition, if the younger clusters have a majority of their mass in the optically veiled sub-population of stars, it would not be surprising that the masses inferred from the nebular emission (i.e., radio and Hα equivalent width, which are mostly unaffected by extinction) are higher than those inferred from the stellar emission (which only samples the optically visible "emerged" stars) as is apparent in Figure 13.

Zoom In Zoom Out Reset image size

A second factor that could contribute to the apparent mass–age correlation is stellar winds. More massive clusters will have disproportionately strong stellar winds, which could preferentially clear away the gas. A lower gas fraction will result in a lower equivalent width, and therefore a larger inferred age. As a result, more massive clusters would tend to appear older by virtue of their dearth of ionized gas. A deficiency of ionized gas could also explain why the more massive clusters tend to have higher mass estimates from SED fitting than those derived from the nebular emission (Figure 13).

NGC 4449 is thought to have undergone some type of interaction in its past (Hunter 1998; Theis & Kohle 2001), which has likely influenced the global star-formation history of the galaxy over the last ∼1 Gyr (e.g. Annibali et al. 2007). Local processes such as stellar winds and supernovae (perhaps the result of a previous burst), however, seem to be playing a significant role in triggering the current bursts of star formation, resulting in the formation of massive clusters detected in the radio.

Figure 14 illustrates that the formation of a number of radio sources appears to have been triggered by an expanding superbubble. Taken alone, the fact that seven sources (23, 24, 25, 28, 30, 37, and 39) lie along the same circular arc (radius = 44'' = 832 kpc) could be viewed as a coincidence, albeit an unlikely one. However, previous observations of gas and dust in NGC 4449 suggest otherwise. H i, CO, and cold dust appear to be depressed or almost absent in this region of the galaxy (e.g., Figure 1 in Böttner et al. 2003) while soft diffuse X-ray emission fills it in (Figure 5 in della Ceca et al. 1997 and Figure 10 in Summers et al. 2003). In addition, kinematic studies of NGC 4449 are consistent with the stellar component of the galaxy observed face-on (Hunter et al. 2002, 2005), which would account for the sources lying on a circular arc in projection. All of this evidence supports the view of an expanding superbubble triggering several sites of star formation in NGC 4449.

Zoom In Zoom Out Reset image size

Further evidence of triggered star formation is apparent within NGC 4449. An Hα shell with radius ∼250 pc is observed in the SW region of the galaxy (see Figures 3 and 4) surrounding numerous optically visible stars and clusters, whose stellar winds and supernovae likely produced the shell, while six radio sources (1, 2, 3, 4, 5, and 7) are found along the periphery of this shell. It is possible that a collision between this shell and the larger superbubble is responsible for compressing the gas in the southern region of the galaxy and triggering the formation of sources 18, 20, and 21.

Another group of radio sources (31–36) is found along a smaller Hα loop (radius ∼40 pc) in the NE region of the galaxy, again containing optically visible stellar clusters that could have triggered the current local burst of star formation.

A close inspection of Figures 3 and 4 reveals that at least a few radio sources in NGC 4449 do not have optical counterparts (e.g. sources 20, 24, and 28). Sources 24 and 28 also do not appear to have any detectable 24 μm emission associated with them. It is unclear whether this is the case for the other radio sources without optical counterparts due to confusion in the 24 μm image. In addition, the "pure" radio sources have only been detected at 3.6 cm and we cannot determine whether they are thermal. We do believe these sources are in the galaxy and not background objects because their locations are associated with star-forming regions. With the current data, however, we can only speculate about the nature of these sources.

If the sources are thermal, the lack of Hα could be explained by heavy extinction from dust. For Sources 24 and 28, the extinctions would have to be at least ∼7 visual magnitudes (Section 4.2.4). Again, if we assume that these sources are thermal, they are probably younger than the radio-detected clusters with optical counterparts (≲2.5 Myr). They would also have relatively low ionizing fluxes (the equivalent of ∼6–15 O7.5 V stars) and masses. These young ages and low masses would be consistent with the observed correlation between derived masses and ages for the thermal radio sources with optical counterparts. However, if the visibly obscured sources are heavily embedded clusters, we would expect them to produce thermal emission from heated dust and be observable in the mid-IR which does not seem to be the case (at least for Sources 24 and 28). It may be more likely that these "pure" radio sources are supernova remnants, since non-thermal emission is anti-correlated with bright mid-IR emission (e.g. Brogan et al. 2006).

We have developed a new multi-wavelength photometry code (SURPHOT) to measure flux densities within identical irregular apertures in multiple images. SURPHOT is written in the Interactive Data Language (IDL) and the basic algorithm is as follows: the user supplies a reference image and selects an aperture type, including "contour" or "free-form" (in addition to a circle or an ellipse). If a contour aperture is chosen, SURPHOT will detect a specified contour level to define a region of interest. Alternatively, the user can choose to draw a free-form aperture by clicking around a region of interest. Once an aperture is defined in the reference image, a background annulus is assigned by expanding the aperture out radially by two specified factors, matching the overall shape of the aperture. The pixel coordinates of vertices in the aperture and background annulus are transformed into R.A. and declination (decl.) using astrometry keywords in the FITS header of the reference image. The R.A. and decl. values can then be transformed into pixel coordinates of any other image using its own FITS header astrometry keywords. In this way, SURPHOT uses the same aperture and background annulus to compute flux densities in any number of user supplied images, provided they are registered to match the astrometric reference frame of the reference image.

For a given source, four background levels and their errors are estimated in each image from which the user can choose mode, mean, resistant mean, and low mean. The mode and resistant mean are calculated using the MMM and RESISTANT_MEAN routines, respectively, which can be found in the IDL Astronomy User's Library.⁵ The low mean is found by using MEANCLIP (also found in the IDL Astronomy User's Library) to compute a 3σ clipped mean of the 25% lowest valued pixels in the background annulus. In general, the mode is our preferred estimate of the background level, however, the other estimates can provide a consistency check or an alternative value if appropriate. The estimated background value (per pixel) is multiplied by the number of pixels in the aperture and then subtracted from the sum of the pixels in the aperture to obtain the total background-subtracted source counts. In the case of the VLA images, the image units of Jy beam⁻¹ are first converted to Jy pix⁻¹, where the conversion factor is given by the pixel area of the 2D Gaussian beam. The background-subtracted flux density (in Jy) is then computed as described above.

Error estimates include contributions from three sources: Poisson noise from the source counts, Poisson noise from the background counts, and the error in measuring the background level. These error terms are calculated in units of electrons, added in quadrature and converted back to counts using the gain. A 1σ detection limit is found in the same way, but only includes contributions from the Poisson noise due to the background counts and the error in measuring the background level.