Spectral Energy Distribution and Radio Halo of NGC 253 at Low Radio Frequencies

At low radio frequencies the intensity of the synchrotron radiation may become sufficiently high (optically thick regime) for reabsorption, termed SSA, to take over. The process may be important, or even dominant, for compact sources (Slysh 1990; Chevalier 1998). We model the synchrotron radio spectra that may turn over due to self-absorption at low radio frequencies as (e.g., Tingay & de Kool 2003)

$$\begin{eqnarray}&&{S}_{\nu }={S}_{\tau =1}{\left(\displaystyle \frac{\nu }{{\nu }_{\tau =1}}\right)}^{-\alpha }\left(\displaystyle \frac{1-{e}^{-\tau (\nu )}}{\tau (\nu )}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\tau (\nu )={(\nu /{\nu }_{\tau =1})}^{-(\alpha +2.5)},\end{eqnarray} \tag{ 4 }$$

where ${\nu }_{\tau =1}$ is a frequency at which the optical depth (τ) reaches unity.

The self-absorbed bremsstrahlung (i.e., free–free absorbed) radio spectrum can be expressed as (e.g., McDonald et al. 2002)

$$\begin{eqnarray}&&{S}_{\nu }={S}_{\tau =1}{\left(\displaystyle \frac{\nu }{{\nu }_{\tau =1}}\right)}^{2}(1-{e}^{-{\tau }_{\mathrm{ff}}(\nu )}),\end{eqnarray} \tag{ 5 }$$

where the opacity coefficient is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}(\nu )={(\nu /{\nu }_{\tau =1})}^{-2.1}.\end{eqnarray} \tag{ 6 }$$

As discussed in Section 1, radio emission from NGC 253 is a mixture of synchrotron emission from cosmic rays from SNRs and thermal emission from H ii regions. The FFA is expected to start dominating at low radio frequencies, where the intensity of the electrons in the ionized gas becomes high (optically thick regime). For NGC 253 it is a natural assumption that the thermal plasma coexists with the synchrotron-emitting electrons, hence the radio spectrum can be modeled as a synchrotron power law with an internal free–free absorbing screen (SFA; Tingay & de Kool 2003),

$$\begin{eqnarray}&&{S}_{\nu }={S}_{0}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-\alpha }\left(\displaystyle \frac{1-{e}^{-{\tau }_{\mathrm{ff}}(\nu )}}{{\tau }_{\mathrm{ff}}(\nu )}\right).\end{eqnarray} \tag{ 7 }$$

Radio continuum spectra of NGC 253 between 76 MHz and 10.7 GHz are plotted in Figure 4 and the flux density measurements are tabulated in Table 4. Background radio sources (Section 4.1) located within the diffuse emission of NGC 253 were subtracted from the total flux density measurements. In the construction of the radio spectrum we used archival data provided that the measurements were of angular resolution comparable to GLEAM or were sensitive to low-brightness emission on angular scales of at least 05, the total flux density was integrated over the diffuse emission of NGC 253 and not fitted by Gaussian components, and the absolute flux density scale and the uncertainties of the measurements were quoted. We do not use measurements at an angular resolution of >20 arcmin because of confusion of NGC 253 with nearby sources.

Zoom In Zoom Out Reset image size

We find the best fitting model to be a second-order polynomial with S₀ = 7.30 ± 0.04 Jy, α = 0.56 ± 0.01, and a curvature c₁ = −0.12 ± 0.01 at a reference frequency of 1 GHz (χ² = 140, with degrees of freedom: dof = 45; Figure 4), which is significantly preferred to a simple power law (Δln(Z) = 100.5 ± 0.3).

The angular resolution of the MWA data is too low to resolve the central starburst region of NGC 253; the highest angular resolution achieved is 102 arcsec at 227 MHz (GLEAM) and 138 arcsec at 169 MHz (EoR), which is over three times the size of the starburst region (Condon et al. 1982; Antonucci & Ulvestad 1988; Carilli 1996). We construct radio spectra of the central starburst region using data from the literature and new measurements from the TGSS ADR1 survey (Figure 4, Table 5). We limit the measurements to those that are at an angular resolution comparable to the size of the central starburst region (approximately 20–30 arcsec).

We find the best fitting model to be a second-order polynomial with S₀ = 2.28 ± 0.02 Jy, α = 0.20 ± 0.01, and a curvature c₁ = −0.24 ± 0.01 at a reference frequency of 1 GHz (χ² = 12.8, dof = 13; Figure 4). We further attempt to model the spectral turnover, and we find SFA (Equation (7)) to be the best fitting model with ${S}_{\tau =1}=4.43\pm 0.14$ Jy, α = 0.43 ± 0.01, and ${\nu }_{\tau =1}=238\pm 15\,\mathrm{MHz}$ (χ² = 42.9, dof = 13). Based on the Bayesian evidence the SFA model (synchrotron plasma absorbed by an internal free–free absorbing screen) is preferred to the pure SSA model (${\rm{\Delta }}\mathrm{ln}(Z)=8.3\pm 0.3).$

Using the total flux density images at 200 MHz (GLEAM), 169 MHz (MWA/EoR), and 1.46 GHz (Carilli et al. 1992), we created maps of the spectral index distribution and corresponding maps of the uncertainty and signal-to-noise ratio (Figure 5). To create the spectral index maps we convolved the 1.46 GHz image with the resolution of the GLEAM image (for the 200 MHz–1.46 GHz spectral index map), and the MWA/EoR image (169 MHz–1.46 GHz spectral index map). The spectral index maps suggest an apparent variation in α across the disk and halo of NGC 253. The central regions of the galaxy are dominated by a flat component (α = 0.31–0.34), coinciding with the central starburst. The gradual steepening along the minor axis found by Heesen et al. (2009a) is seen only on the northern side of the galaxy in our maps. Further, the region extending southwest (SW) from the central starburst region seems to be flatter (α ≲ 0.53) than the other parts of the galaxy outside the central starburst (α ∼ 0.60–0.65).

Zoom In Zoom Out Reset image size

To verify the significance of the variation in spectral index across the galaxy, we used the T–T method (Turtle et al. 1962). The method allows one to estimate a spectral index within defined regions of a source between two frequencies. We define nine regions within NGC 253 (Figure 5). Due to the low angular resolution of our observations and to avoid oversampling, only one data point is associated with each region. The results are shown in Table 6 and Figure 6.

Zoom In Zoom Out Reset image size

Table 6.  Radio Spectral Indices (α) Measured for Selected Regions of NGC 253 (Figure 5)

Region	${\alpha }_{1.4\mathrm{GHz}}^{169\mathrm{MHz}}$	${S}_{169\mathrm{MHz}}$	${S}_{1.4\,\mathrm{GHz}}^{169\mathrm{res}}$	${\alpha }_{1.4\mathrm{GHz}}^{200\mathrm{MHz}}$	${S}_{200\mathrm{MHz}}$	${S}_{1.4\,\mathrm{GHz}}^{200\mathrm{res}}$
		(mJy)	(mJy)		(mJy)	(mJy)
1	0.54 ± 0.06	242 ± 26	75 ± 10	0.57 ± 0.06	233 ± 25	75 ± 10
2	0.63 ± 0.04	586 ± 59	149 ± 13	0.68 ± 0.03	581 ± 51	149 ± 13
3	0.57 ± 0.02	883 ± 89	257 ± 11	0.62 ± 0.02	884 ± 73	257 ± 11
4	0.67 ± 0.04	415 ± 42	97 ± 9	0.72 ± 0.04	412 ± 36	97 ± 9
5	0.31 ± <0.01	6558 ± 656	3341 ± 14	0.34 ± <0.01	6589 ± 528	3341 ± 14
6	0.60 ± 0.07	224 ± 24	61 ± 10	0.70 ± 0.07	246 ± 25	61 ± 10
7	0.58 ± 0.03	488 ± 49	139 ± 10	0.63 ± 0.03	487 ± 42	139 ± 10
8	0.50 ± 0.02	616 ± 62	211 ± 9	0.51 ± 0.02	578 ± 49	211 ± 9
9	0.53 ± 0.06	216 ± 23	69 ± 9	0.51 ± 0.06	192 ± 21	69 ± 9

Note. Radio spectral index maps are used at 169 MHz (MWA/EoR0; S_{169 MHz}), 200 MHz (GLEAM; S_{200 MHz}), and 1.465 GHz convolved to 169 MHz (${S}_{1.4\,{\rm{G}}{\rm{H}}{\rm{z}}}^{169{\rm{r}}{\rm{e}}{\rm{s}}}$) and separately to 200 MHz (${S}_{1.4\,{\rm{G}}{\rm{H}}{\rm{z}}}^{200{\rm{r}}{\rm{e}}{\rm{s}}}$) angular resolution.

Download table as:  ASCIITypeset image

We find that between 200 MHz and 1.465 GHz the apparently flat regions (Regions 8 and 9) are statistically different from the other regions within NGC 253 apart from Region 1 (Figure 6, Table 6). This spectral flattening is not present in the map of the radio spectral index distribution between 330 MHz and 1.46 GHz of Heesen et al. (2009a), even though the same high-frequency map is used. This clearly indicates that the flattening occurs at <300 MHz. There is also a slight flattening of the spectral index in the northeast (NE) region perpendicular to the major axis (Region 1), although we find this flattening to be statistically different only from Regions 2 (eastern northwest (NW) halo), 4 (radio spur), and 5 (including central starburst) in the ${\alpha }_{1.4\,\mathrm{GHz}}^{200\mathrm{MHz}}$ map. All regions flatten further at 169 MHz, reducing the differences between spectral indices of the regions.

We measured the observed, projected maximum vertical extent of the disk and halo of NGC 253 at 169 MHz as a function of the distance from the nuclear region along the major axis (Figure 7, Table 7). The extent is measured perpendicular to the major axis (at PA = −38°, z-direction) in steps of 132 arcsec (2.4 kpc) separately for the north (filled circles) and south (empty circles) sides of the disk and halo as divided by the major axis (PA = 52°). We find the projected radio halo to extend up to 4.75 kpc above the optical edge (B band including 90% of total light, Lauberts & Valentijn 1989) and 6.3 kpc above the infrared edge (total K_s band, Jarrett et al. 2003) of the galaxy, reaching up to a total 7.9 kpc in the z-direction. This is consistent with previous radio measurements at higher radio frequencies (Reynolds & Harnett 1983; Hummel et al. 1984; Carilli et al. 1992, Figure 1; but see Heesen et al. 2009a), as well as with broadband X-ray observations of the galaxy's extended extraplanar emission (Pietsch et al. 2000). Heesen et al. (2009a, 2009b) attributed the decrease in scale height they measured and modeled to the increased synchrotron losses in the central regions, where the magnetic field is highest.

Zoom In Zoom Out Reset image size

The shape of the radio halo in our MWA radio images (Figure 2) resembles the "horn-like" or "X-shaped" structure seen at gigahertz radio frequencies (Heesen et al. 2009a), in H i (Boomsma et al. 2005; Lucero et al. 2015), X-rays (Fabbiano 1988; Pietsch et al. 2000; Bauer et al. 2008), Hα (G. Meurer, private communication; Figure 2), UV (Hoopes et al. 2005), and far-IR (Kaneda et al. 2009).

The radio halo was investigated in detail by Heesen et al. (2009b) who, through modeling of the large-scale magnetic field, attributed its origin to disk wind, confirming previous suggestions (Carilli et al. 1992; Beck et al. 1994). This is also in line with the Hα and optical analyses of the inner starburst-driven superwind (Westmoquette et al. 2011). In our MWA maps both the NE and NW halo regions are pronounced. The extended soft X-ray emission (<1 keV) of the halo, detected in the northwestern direction from the disk of NGC 253, is interpreted as bubbles of hot low-density gas (Figure 2; Pietsch et al. 2000; Strickland et al. 2002). Heesen et al. (2009b) postulate that the large-scale magnetic field of the halo follows the walls of these bubbles, where it may be compressed, producing the X-shaped synchrotron radio halo as well as heating up pre-existing cold gas to X-ray energies. The northern halo can also be easily seen, in projection, in our Figure 2 where we overlay MWA/EoR intensity contours on the XMM-Newton image of soft X-ray emission.

The southeast (SE) region of the extended halo, the "spur" (Carilli et al. 1992), is contaminated by a background source. We modeled the background source as unresolved (at MWA angular resolution) and subtracted it from the 169 MHz EoR and deep 200 MHz GLEAM images as described in Section 4.1. The residual emission, which we consider intrinsic to the "spur," is shown in Figure 2 overplotted on the Hα and X-ray images. Although slightly offset (as expected, see Figure 19 in Heesen et al. 2009b), the feature broadly coincides with the extended outflows at both frequencies as clearly seen in our figure. The spur has been previously interpreted as originating from the active star formation region in the NE end of the bar (Waller et al. 1988; Carilli et al. 1992). We confirm the association of the radio "spur" with an extended Hα outflow, which can be clearly seen in our Figure 2 (G. Meurer, private communication; Kennicutt et al. 2003). In our radio spectral map between 200 MHz and 1.46 GHz the spur displays the steepest spectral index within the galaxy (Figure 5), gradually steepening from α ∼ 0.7 to α ∼ 0.9 outward from the disk of the galaxy (though it still contains background source no.1 with α = 0.57), which is indicative of ageing unabsorbed synchrotron plasma.

The broadband spectrum of the total radio emission from NGC 253 is steep, although flattening at megahertz radio frequencies (Figure 4). The total radio emission originates from SNRs, H ii regions (predominantly the central starburst region; e.g., Ulvestad & Antonucci 1997; Ulvestad 2000; but see Waller et al. 1988; Hoopes et al. 1996) and electrons (cosmic rays) freely spiralling in the large-scale magnetic field (radio halo; Heesen et al. 2009a). The steep spectrum is understood to be of a synchrotron origin. The flattening of the spectrum, however, may be due to a number of reasons, including some degree of absorption of the synchrotron emission (Equation (3)) and a low-energy cutoff of the electron population, where in general it is assumed the spectrum of the electron energy (E) can be described as a power law $N(E)\propto {E}^{-p}$ with the index p related to the radio spectral index as $\alpha =(p-1)/2$. The flattening of the spectra of starburst galaxies at low radio frequency is not unusual and it has been observed previously (e.g., Marvil et al. 2015).

We attempt to separate the contributions of the extended and central starburst emission to the total flux densities to see whether the spectral flattening can be attributed mostly to the absorption occurring in the central starburst region. To do this we simultaneously fitted the spectra of the total flux density and the central starburst region, assuming an underlying two- or three-component spectrum composed of the central starburst region (component C) and extended emission (component E, consisting of one or two components). As an example of our fitting method, in the case where component E is modeled as a power law and component C as free–free absorbed synchrotron emission, the model equation (${S}_{\nu }^{\mathrm{mod}})$ takes the form

$$\begin{eqnarray}&&{S}_{\nu }^{\mathrm{mod}}={({S}_{\nu }^{{\rm{E}}}+{S}_{\nu }^{{\rm{C}}})}^{\mathrm{tot}}+{({S}_{\nu }^{{\rm{C}}})}^{\mathrm{cor}},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{S}_{\nu }^{{\rm{E}}}={S}_{0}^{\mathrm{ext}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-{\alpha }^{\mathrm{ext}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{S}_{\nu }^{{\rm{C}}}={S}_{0}^{\mathrm{cor}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-{\alpha }^{\mathrm{cor}}}\left(\displaystyle \frac{1-{e}^{-{\tau }_{\mathrm{ff}}(\nu )}}{{\tau }_{\mathrm{ff}}(\nu )}\right).\end{eqnarray} \tag{ 14 }$$

τ_ff(ν) is given by Equation (6) and the indices indicate: tot—total, cor—core, C—component C, E—component E, ext—extended. This model equation is then compared to our observed data (${S}_{\nu }^{\mathrm{obs}}$), where ${S}_{\nu }^{\mathrm{obs}}={S}_{\nu }^{\mathrm{obs},\mathrm{tot}}+{S}_{\nu }^{\mathrm{obs},\mathrm{cor}}$. In the two-component model, we fit component E with either a simple power law or a curved spectrum (second-degree polynomial; Equation (1)). In the three-component model, we fit component E with a combination of a power law and either a curved spectrum, a synchrotron self-absorbed component (Equation (3)), or a synchrotron free–free absorbed component (Equation (7)). Component C is modeled with either a second-degree polynomial, self-absorbed synchrotron emission, or synchrotron power-law emission with a free–free absorbing screen.

We find the best fitting model to be the three-component model, with (1) component E modeled as a combination of a simple power law and a second-order polynomial, with a total flux density S₀ = 5.08 ± 0.50 Jy and α = 0.71 ± 0.01 at reference frequency 1 GHz, and 34% of S₀ becoming absorbed at low frequencies as described by a second-order polynomial with c₁ = −0.76 ± 0.17, and (2) component C modeled as a synchrotron plasma with an internal free–free absorbing screen, with ${S}_{\tau =1}=4.34\pm 0.11$ Jy, ${\nu }_{\tau =1}=231\pm 14\,\mathrm{MHz}$, α_SSA = 0.41 ± 0.01. The three-component model (χ² = 173, dof = 58) is favored over any two-component model, even if component E is modeled as a second-order polynomial (${\rm{\Delta }}\mathrm{ln}(Z)\gt 4.3\pm 0.3).$

Preference of the three-component model indicates that flattening of the spectrum of component E at the lower frequencies is non-negligible. In principle, this flattening could be attributed to SSA caused by shock reacceleration of the halo/disk plasma, an external free–free absorbing screen, or an intrinsic low-energy cutoff of the electron distribution. Thermal FFA can be largely excluded based on the limited evidence for high thermal content in the halo of NGC 253, especially in the SW region (see Section 5.2.2 for details). Although our Bayesian inference tests indicate that the flattening caused by the SSA is moderately favored over the low-energy cutoff in the electron distribution, which could be inferred from the second-order polynomial fit (${\rm{\Delta }}\mathrm{ln}(Z)=2.5\pm 0.3$), we find that the former fit is associated with very high uncertainties. We also find that any model invoking multiple internal components that we tested is strongly favored over an external free–free absorbing screen. Our results do not change in the absence of the data points between 300 and 600 MHz that may seem unusually high, which further strengthens our result that the radio-emitting plasma in the disk and halo of NGC 253 is composed of at least two spectral components that behave differently. This result is also in line with our findings on the variation of the radio spectral index across the galaxy (discussed further in the next section).

Furthermore, another important result of our spectral modeling is that the central starburst region is best modeled by the SFA model. Although, in principle, the second-order polynomial is statistically favored, the SFA model is more realistic. Curved radio spectra can be explained by a low-energy cutoff of the electron population, SSA, SFA, or FFA models. As we have shown, the SSA model is statistically ruled out (Section 4.2.2). Given the overwhelming evidence of a significant thermal component in the central starburst (e.g., Figure 2; Ulvestad & Antonucci 1997; Keto et al. 1999; Kepley et al. 2011) coexisting with synchrotron plasma, the synchrotron FFA model is more likely than the low-energy cutoff in the electron population. The plasma becomes optically thick around a frequency of 230 MHz. Given this result, and under a simplified assumption that a uniform optical depth holds across the region, we estimate a typical emission measure (Oster 1961; Spitzer 1978) of the absorbing gas toward the central starburst region to be very high, of the order 4 × 10⁵ pc cm⁻⁶ (assuming electron temperature T_e = 4000 K; see discussion in Carilli 1996).

We now consider the origin of the variation in α across the disk and halo of NGC 253. The southern flattening occurs beyond the SW spiral arm, in the halo region. In Figure 8 we overplot the ${\alpha }_{1.4\,\mathrm{GHz}}^{200\,\mathrm{MHz}}$ and extraplanar H i contours on the XMM-Newton soft X-ray image of NGC 253. The diffuse X-ray emission indicates ionized hot gas. As pointed out by Lucero et al. (2015), the neutral cold H i gas seems to surround the X-ray-emitting regions. The spatial variations of the radio spectral index seem to follow the distribution of the X-ray emission, with the steepening of α occurring in the regions of intense soft X-ray emission (radio spur and NW halo) and the flattening around the voids of diffuse X-ray plasma (western SE halo, eastern NW halo). This distribution seems to also match the extraplanar H i emission, especially in the western SE halo region. In Hα we detect faint diffuse emission in the NE halo and the southern "spur" (Figure 2), with the line fluxes measured down to 3 × 10⁻¹⁸ erg cm⁻² s⁻¹ arcsec⁻². In these regions the spectral index seems to steepen (Figure 8(B)). There is, however, almost no Hα emission present in the western SE halo, while most of the flattening of component E in the modeling of Section 5.2.1 is due to this region (Figure 8(C), Regions 8 and 9 in Figure 5).

Zoom In Zoom Out Reset image size

It has previously been suggested that the halo gas originates from both galactic "fountains" from the star-forming disk and a galactic superwind (Pietsch et al. 2000; Bauer et al. 2008). This strong superwind may be pushing and collimating the neutral cold gas in the halo (Lucero et al. 2015). In the case of strong collimation shocks one may observe flattening of radio spectra due to SSA in transverse shocks. Our results seem to favor such a scenario. It is also worth noting that the spectral flattening of the SW halo corresponds to an extended loop, or arch, seen in optical images (blue filter; Beck et al. 1982). However, if the SW halo is predominantly diffuse and of low density, the flattening may be rather due to an intrinsic low-energy cutoff of the electron distribution.

Another important note is that, based on our modeling of the broadband radio spectrum, the SW halo region cannot be fully responsible for the total spectral flattening. Radio emission that becomes absorbed at lower frequencies constitutes more than 30% of the total extended radio emission at 1 GHz (1.73 ± 0.36 Jy), while the SW region is only 0.77 Jy at that frequency, which means that the flattening must also be occurring, although to a lesser degree, in other regions across the disk and halo.

Although the SW flattening of the radio spectral index is most likely of a synchrotron origin, an external free–free absorbing screen was also previously suggested. The foreground absorption model was favored by Bauer et al. (2008) based on their modeling of X-ray data and apparent differences in radio and X-ray halo morphologies. As proved later (e.g., Heesen et al. 2009a, 2009b; this publication, Section 5.1), deep continuum and polarization radio observations at both gigahertz and megahertz frequencies reveal the horn-like structure of the radio halo, which directly resembles the X-ray diffuse emission. Based on the assumptions of equipartition, Heesen et al. (2009a, 2011) find the magnetic field within the halo to be very high, 7–12 μG, reaching as much as 160 ± 20 μG in the central regions and 46 ± 10 μG in the starburst outflow. The magnetic field in the central regions is strong enough for synchrotron emission to contribute a few per cent to the total X-ray emission (Lacki & Thompson 2013). As discussed in the previous section, we also find that an external free–free absorbing screen is not a statistically preferred model. These new findings support models in which the total X-ray emission may indeed come from a combination of thermal and synchrotron plasma rather than multi-temperature pure thermal plasma with an externally caused absorption (see Bauer et al. 2008).