THE SPITZER SURVEY OF STELLAR STRUCTURE IN GALAXIES (S⁴G): PRECISE STELLAR MASS DISTRIBUTIONS FROM AUTOMATED DUST CORRECTION AT 3.6 μm

For nearby galaxies, the light detected in the IRAC $3.6\;\mu {\rm{m}}$ filter arises mainly from two components: the photospheres of (old) stars and from dust emission. Old stars dominate the stellar flux, as their atmospheres are cold and their blackbody curves peak close to $3.6\;\mu {\rm{m}}$; additionally, due to CO absorption at $4.5\;\mu {\rm{m}}$, old K and M giants exhibit blue [3.6]–[4.5] colors (Pahre et al. 2004; Willner et al. 2004). Younger (hotter) stars are not expected to contribute significantly to the observed stellar emission. However, due to their strong UV fluxes, these younger stars can heat their surrounding dust, which, in turn, re-radiates at longer wavelengths and can also account for a significant fraction of the light at $3.6\;\mu {\rm{m}}$ (Meidt et al. 2012b). Models show that the dust emission arises from the PAH bands (specifically, the $3.3\;\mu {\rm{m}}$ feature) and the continuum radiation from very small grains (e.g., Li & Draine 2001), and this (hot) dust emission becomes more prominent near the sites of star formation, where the presence of young stars leads to stronger radiation fields.

According to the measurements of Flagey et al. (2006), PAHs should have colors [3.6]–[4.5] ∼ 0.3, including both the $3.3\;\mu {\rm{m}}$ PAH feature and the underlying PAH continuum detected, e.g., in Milky Way reflection nebulae. In the absence of the PAH continuum, the [3.6]–[4.5] color is bluer (as low as [3.6]–[4.5] ∼ − 0.1), although this scenario will likely be uncommon here, given unavoidable mixing at our resolution. For the continuum dust component, we expect [3.6]–[4.5] ∼ 1.0, adopting the power-law approximation ${f}_{\nu }\propto {\nu }^{-\alpha }$ to the Wien side of the dust spectral energy distribution (SED) with $\alpha =2$, which Blain et al. (2003) suggest better accounts for the observed dust SED shape than a modified blackbody in the near-IR. As shown in Figure 1, the diffuse mixture of PAH and continuum dust emission will exhibit resulting colors in the range 0.2 < [3.6]–[4.5] < 0.7.

Zoom In Zoom Out Reset image size

The relative flux contributions of the PAH and continuum at $3.6\;\mu {\rm{m}}$ can be estimated from the ratio of the $6.2\;\mu {\rm{m}}$ PAH feature to the $6.2\;\mu {\rm{m}}$ continuum measured by Peeters et al. (2002). For normal spiral galaxies, this ratio is between 0.72–1.16 (on average ∼1) and can be as low as 0.1 in QSOs and the nuclei of Seyfert galaxies. The ratio of the 3.3–6.2 $\mu {\rm{m}}$ PAH features measured in H ii regions, star-forming regions and in planetary nebulae falls in the range 0.15–0.25 (Hony et al. 2011). With the Blain et al. power-law dust continuum, this implies that the ratio of the $3.3\;\mu {\rm{m}}$ PAH emission to the underlying continuum is ∼0.4–0.9, or as low as 0.1 in galaxy nuclei or other regions dominated by hot dust.

Figure 1 shows the resulting colors for the diffuse dust at $3.6\;\mu {\rm{m}}$ adopting this range of PAH flux fractions for a wide range of possible PAH colors and assuming a fixed continuum color. As found by Meidt et al. (2012b), the primary non-stellar emission detectable with ICA is in the form of this "diffuse dust" component, i.e., the mixture of PAH and the dust heated by the ambient interstellar radiation field, away from star-forming regions. As considered later in Section 4.3 and in Appendix B, some galaxies also exhibit additional emission from hot dust, isolated in star-forming H ii regions. This dominant component in H ii regions, which can be thought of as the far end of the spectrum exhibited by the "diffuse dust," should have [3.6]–[4.5] colors closer to 1.

In the implementation of our S⁴G pipeline we make use of the fact that regions containing hot dust (and negligible PAH) have colors that are distinguishable from the nominal diffuse dust component. We also make use of the fact that the colors of the dust and the old stellar population are very different. The colors for old stellar populations lie in the range $-0.2\lt [3.6]-[4.5]{| }_{\mathrm{stars}}\lt 0$, according to the observed colors of giant stars (Meidt et al. 2014). This is consistent with the observed colors of early-type galaxies, in the absence of significant dust emission (Peletier et al. 2012, and Norris et al. 2015). The non-stellar $[3.6]-[4.5]$ color, on the other hand, always appears positive in these bands ($[3.6]-[4.5]{| }_{\mathrm{dust}}\gt 0$), although it can span a wider range, as seen in Figure 1.

Table 1 summarizes the colors expected for the different sources of emission at 3.6 μm. These sources combine together to produce the observed [3.6]–[4.5] colors plotted in Figure 2. There we show representative colors for young, relatively metal-poor, late-type galaxies and their old, metal-rich early-type counterparts (with bluer [3.6]–[4.5] colors),²² adopting the dust contamination fractions estimated in Appendix A. The latter consistently exhibit [3.6]–[4.5] < 0, whereas the addition of even a little bit of dust emission will move the [3.6]–[4.5] colors of late-type galaxies above zero.

Zoom In Zoom Out Reset image size

Table 1.  Main Sources at $3.6\;\mu {\rm{m}}$ and Their [3.6]–[4.5] Colors

Source	Typical [3.6]–[4.5] Range
Old starsa	$-0.2-0$
Diffuse dust	$\sim 0.2-0.7$
PAH emissionb	∼0.3
Dust in H ii regionsc	∼1.0

Notes. Main sources identified with ICA, and typical values of their colors.

^aWillner et al. (2004), Pahre et al. (2004), Peletier et al. (2012), Meidt et al. (2014), Norris et al. ^bDominated by PAH emission. Approximate estimation based on Flagey et al. (2006). ^cRepresentative number based on Blain et al. (2003); it can either correspond to hot dust, e.g., in H ii regions, or, in the most extreme case, to hot dust heated to large temperatures near an AGN.

Download table as:  ASCIITypeset image

Here, we introduce the ICA method and the nomenclature that we will use. The method has already been described and validated by Meidt et al. (2012a), so we will only provide a short summary here.

Similar to Principal Component Analysis (PCA), ICA is a means of blind source separation which we use to extract measurements of the flux and the wavelength-dependent scaling of individual components from linear combinations of the input data. But in contrast to PCA, ICA maximizes the statistical independence of the sources rather than requiring that the sources are orthogonal (with zero covariance). We use the fastICA realization of the method developed by Hyvärinen (1999) and Hyvärinen & Oja (2000), which achieves statistically independent solutions by maximizing the non-Gaussianity of the source distributions.

In practice, given N input images (e.g., at N different wavelengths), ICA will identify at most N underlying distinct sources that contribute to the flux in each of the N images, under the assumption that each image is a linear combination of the sources. Because S⁴G provides images of each galaxy in two channels, ICA can identify two distinct components through solution of the following equation:

$$\begin{eqnarray}&&{\boldsymbol{x}}={\boldsymbol{A}}\cdot {\boldsymbol{s}},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{x}}$ is the $2\times P$ measurement set, with as many columns as pixels P in the analysis region, ${\boldsymbol{s}}$ is the $2\times P$ source solution set, and ${A}_{i,j}$ is an invertible 2 × 2 matrix of "mixing coefficients" determined simultaneously with ${\boldsymbol{s}}$.

As first considered by Meidt et al. (2012a), ICA provides a way to distinguish between the old stellar population that dominates the light in the IRAC bands and additional emission present to varying degrees (depending on the nature of the source), without a priori knowledge of the relative proportions of the sources or their colors between 3.6 and 4.5 μm. In six prototypical star-forming disk galaxies, Meidt et al. (2012a) showed that ICA can identify and remove the combined PAH and continuum dust emission (tracing star formation, as observed at longer wavelengths) and localized emission from dusty asymptotic giant branch (AGB) and red supergiant (RSG) stars. The ICA correction leaves a cleaned, smooth map of the old stellar light consistent with expectations for an old, dust-free stellar population.

Our pipeline implementation of ICA includes an estimate of the uncertainties on both components by running the ICA sequence 48 times, based on 48 perturbations to the mixing matrix. As an initial guess for the mixing coefficients we adopt the expected color range of stars and dust ($-0.2\lt [3.6]-[4.5]{| }_{\mathrm{stars}}\lt 0$ and $0\lt [3.6]-[4.5]{| }_{\mathrm{dust}}\lt 1.5$), but we find that ICA solutions quickly depart from the initial seeds and converge to the final solution. The perturbations to the mixing matrix correspond to small steps in dust and stellar color in a count-controlled nested loop. The primary loop implements a change in the stellar color by 0.04 mag, during which the secondary iterates the dust color in steps of 0.3 mag (at fixed stellar color).

For each of the two original input images we obtain two source images, which we refer to as s1 and s2, providing a total of four images (which we could in principle name $s{1}_{3.6}$, $s{2}_{3.6\mu {\rm{m}}}$ and $s{1}_{4.5}$, $s{2}_{4.5\mu {\rm{m}}}$). At either wavelength it will always hold that $(s1+s2)=\{\mathrm{original}\ \mathrm{image}\}$. However, since the two source images at $4.5\;\mu {\rm{m}}$ are scaled replicas of the two source images at $3.6\;\mu {\rm{m}}$, with scalings set by the colors of the solutions ($[3.6]-[4.5]{| }_{s1}$ and $[3.6]-[4.5]{| }_{s2}$), throughout the paper we will, by default, refer to the $3.6\;\mu {\rm{m}}$ solutions as s1 and s2.

With this notation, ICA typically separates the emission arising from old stars (which we identify with s1) from dust (s2), as is illustrated by Figure 3. The secondary emission that arises from dust can have a number of physical origins, covering different color ranges (see Table 1), but is always redder than old stars in these bands, which allows ICA to perform the separation into two distinct components.

Zoom In Zoom Out Reset image size

Because our implementation of ICA seeks to maximize the non-Gaussianity of the source distributions (Hyvärinen & Oja 2000), ICA is sensitive to extreme outliers in color, even if these cover only a small region. This means that a few pixels with very different color from the dominant color can bias the whole solution (not only the area they cover), yielding an unrepresentative separation of the original image into two components. Consequently, properly masking any foreground stars or background galaxies is particularly important.

We have developed a specific masking strategy that builds on masks already available from S⁴G Pipeline 2, but includes further corrections to ensure that bright regions belonging to the galaxy are not masked in a first attempt. Pipeline 2 masks were developed primarily for galaxy photometric decompositions (Pipeline 4), and that is the reason why this subtle yet important modification is necessary to make the masks applicable for our purposes here.

Specifically, we remove from the Pipeline 2 masks any contiguous regions that are small (maximum area in pixels of $8000\;{\mathrm{Mpc}}^{2}/{d}^{2}$) and have colors $-0.3\lt [3.6]-[4.5]{| }_{\mathrm{orig}}\lt 0.75$, as they tend to be bright regions intrinsic to the target galaxy (e.g., H ii regions). The distance-dependent cut in size corresponds to a radius of ∼150 pc, covering well even the largest H ii regions (see Whitmore et al. 2011), whereas the color criterion matches roughly the range spanned by emission from old stars and dust, leaving out saturated field stars and observational artifacts, which usually exhibit colors outside this range.

Therefore, the set of pixels to which we apply ICA (the ICA solution area) is defined to contain all emission out to the edge of the galaxy and avoid external sources such as field stars and background galaxies. The area is centered on the photometric center defined in Pipeline 3 and extends out to the 25.5 mag arcsec⁻² isophote (Muñoz-Mateos et al. 2015, in preparation), excluding the masked regions. At this early stage, we also correct for PSF effects by convolving both images with the IRAC 4.5–3.6 μm PSF kernel from Aniano et al. (2011).

The expected color for an old stellar population of ages $\tau \;\sim$ 2–12 Gyr is $-0.2\lt [3.6]-[4.5]\lt 0$ (see references in Table 1). Some of the galaxies in S⁴G have original [3.6]–[4.5] global colors in that range, which implies that they are already compatible with an old stellar population. As shown in Figure 2, an originally blue color is indicative of little to no dust emission; we therefore do not apply ICA to these galaxies. We have found that these cases can be best identified by calculating the weighted mean of the original color for the pixels in our ICA solution area; the weights are chosen as the inverse of the variances (${w}_{i}=1/{\sigma }_{i}^{2}$, where ${\sigma }_{i}$ is the original color error, as described in Section 4.7). By adopting the weighted mean we avoid the bias of low signal-to-noise regions, from which color information is less reliable. A total of 376 galaxies in S⁴G (16% of the sample) have colors originally compatible with old stars, and are therefore excluded from the further ICA analysis.

Inspection of the 2D color maps confirms very little contamination from dust in galaxies where the original negative color is consistent with that of old stars, although dust may not be entirely absent. Still, we prefer not to run ICA in these cases, since our testing suggests that more uncertainty in the old stellar light map can be introduced under these circumstances than if the non-stellar emission is simply retained. In particular, we have found that when emission from dust becomes negligible, the separation starts to be dominated by spatial fluctuations in the original color due to noise. As shown in Appendix A, a low fractional dust contamination leads to the largest errors on the stellar flux obtained from the ICA separation.

As described in Section 3.1, the primary non-stellar emission detectable with ICA is in the form of a "diffuse dust" component, i.e., a mixture of PAH and dust continuum. Some galaxies also exhibit secondary non-stellar emission from hot dust isolated in star-forming H ii regions or near an active nucleus. Even in these cases, it is most common that the diffuse dust, which is spatially more extensive than the dust in H ii regions, dominates the ICA solution. But as the hot dust regions are assigned a color that is unrepresentative (and, in particular, less red than their true color), the flux in this secondary dust component can be overestimated (see Section 4.6 below).

In other cases (e.g., a high number of H ii regions, or very pervasive star formation), ICA will favor the hot dust. While this may provide the best description of the non-stellar emission in some galaxies (with genuinely more hot than diffuse dust), often the emission from the diffuse dust component is still present, but underestimated by ICA (since it is assigned a color redder than its intrinsic value). As we are more interested in isolating dust emission in the disk (tracing star formation, heating from hot evolved stars, etc.), we prefer to completely avoid these very red, hot dust-dominated regions by performing a second iteration of ICA with the corresponding areas masked.

This same iterative process is useful in general for reducing the number of sources of emission present in the solution area. Extremely bright nuclear point-like sources have a very similar effect on ICA solutions as the bright field stars and background galaxies described in Section 4.1. In some cases, they initially completely dominate the secondary component identified with ICA. Ignoring the central source helps ICA identify other types of secondary sources of emission, as can be seen in Figure 5.

Zoom In Zoom Out Reset image size

The relative contributions of the additional sources determines the degree to which ICA identifies a realistic secondary source. When, for example, the second component is dominated by the diffuse dust emission, and hot dust appears in only small localized H ii regions, ICA successfully describes the dominant source (in this case diffuse dust; see Meidt et al. 2012a, and Appendix B). This is close to the optimum decomposition and is characteristic of most star-forming disk galaxies. However, when the contributions of the two sources become comparable (e.g., when they cover similarly sized areas in the disk), ICA finds a compromise between the two. In some cases, this compromise is an acceptable outcome, but it is often possible to obtain a noticeably improved solution for one of the components by running ICA again, now with pixels containing the other source masked. This is confirmed by the tests described in Appendix A, and, for the interested reader, the details of this empirically optimized strategy are described below (Section 4.3.1).

In our automatic pipeline, the second ICA iteration is determined to be effective when the second ICA s2 is bluer than the first ICA s2 (i.e., the reddest sources have been removed from the dust map). This is the case for 66% of the objects, in which $[3.6]-[4.5]{| }_{s2}$ is effectively reduced after the second iteration and ICA₂ is chosen. The improvement in the solution for these cases is also associated with a decrease in the color uncertainties (in 79% of those objects $[3.6]-[4.5]{| }_{s2}^{{\mathrm{ICA}}_{2}}$ has a smaller uncertainty than $[3.6]-[4.5]{| }_{s2}^{{\mathrm{ICA}}_{1}}$). Even in these cases, the fractional change to the total old stellar flux (and, thus, total mass) is not dramatic, typically on the order of ∼5%, but of course this can locally become more significant.

We have confirmed that running a third iteration of ICA would introduce a change which is negligible, smaller than the uncertainty in 98% of the cases. Therefore, for simplicity, we have only implemented two iterations in the final pipeline.

4.3.1. Technical Details on Second Iteration, ICA₂

ICA provides a solution for each pixel in the solution area using information from all other pixels in the area. But not all pixels contain a genuine secondary source; consequently, some low-level noise is systematically introduced throughout the analysis area. These low-level pixels can either be positive or negative, since ICA does not impose non-negative solutions, and they are clearly artificial, as they fluctuate nearly uniformly around zero. To prevent arbitrary removal (or addition) of flux from the stellar map, we conservatively impose a minimum flux of s2 above which the emission is arguably genuine. We base our threshold on the map of propagated uncertainties, which we are able to compute following the recipe described and justified below.

4.4.1. Technical Details on Thresholding

We define the threshold relative to the map of propagated uncertainties, with the additional simple assumption that the noise randomly introduced by ICA is symmetric about zero. In principle, there is no mechanism within ICA that should asymmetrically bias the noise toward positive values, and this is confirmed by the histograms of the flux distribution of s2 over areas where no dust emission is expected (e.g., in the outermost part of the galaxy). For a randomly selected sample of galaxies, these histograms are indeed very symmetric distributions centered around zero (some are shown in Figure 6 for illustrative purposes).

Zoom In Zoom Out Reset image size

In light of the symmetric distribution of noise around zero, we adopt a variable threshold, which we set to M times the propagated map of uncertainties, where M is defined by the 95th percentile in the noise distribution measured from the negative pixels in s2 divided by the average propagated noise:

$$\begin{eqnarray}\mathrm{flux}\;\mathrm{threshold} & = & \left(-\displaystyle \frac{95\mathrm{th}\;\mathrm{percentile}\;\mathrm{of}\;\mathrm{negatives}}{\mathrm{average}\;\mathrm{propagated}\;\mathrm{noise}}\right)\\ & & \times \mathrm{map}\;\mathrm{of}\;\mathrm{propagated}\;\mathrm{noise}.\end{eqnarray} \tag{ 2 }$$

As we assume that the distribution of negative pixels is representative of the global noise distribution, such that a symmetrical cut in the 95th percentile toward positive values removes all negative noise and also 95% of the artificially introduced positive noise. What is left is primarily the significant emission in s2.

The fluxes at these threshold locations are returned back to the map s1, which, as a result, once again becomes identical to the original image for those positions. We note that some genuinely identified low-level dust emission may be returned to the stellar map. (At low flux values, where the distributions of genuine dust emission and noise overlap, there is no way to uniquely determine their true contributions.) In effect, we modify the picture of the old stellar light map supplied by the original image by removing flux only if we are confident that it corresponds to dust emission. This prevents the introduction of additional structure by artificially removing flux consistent with noise.

Because ICA optimally separates two sources in two images, the second iteration, for which the number of sources has been reduced, supplies in principle our best solution. In practice, this translates into a dust color which becomes less red after the second iteration, as a combination of hot dust, evolved star regions and possibly red nuclei have been masked, to converge on a color representative of the spatially dominant diffuse dust emission.

However, when the second iteration returns solutions that are unphysical (i.e., the ICA colors of one or both sources fall outside the range of colors expected for old stars or dust emission, or the identified dust color turns redder after the second iteration), then we take this as an indication that the solution has been biased by the masking of noisy pixels rather than by a true additional secondary source. The solution from the first iteration is chosen in this case (34% of the cases).

Independently of the choice of solution, a posteriori processing is necessary, as the light from old stars will be over- or under-estimated in pixels containing dust properties that deviate from the dominant ones. There are two cases when we choose to adjust a posteriori the information in map s1. In the first case, when pixels contain dust with redder colors than the ICA $[3.6]-[4.5]{| }_{s2}$ color and these are unmasked, the solution overestimates the true secondary flux in these pixels. This manifests itself in oversubtractions on the stellar flux map, which we find to be an acceptable price to pay for the high quality with which a larger fraction of the dust emission can be described with ICA. In these cases we opt to linearly interpolate over the regions of oversubtraction, which are typically small, effectively filling them in with information about the old stellar light from neighboring pixels.

A second case arises when those red regions have been masked in the second iteration (i.e., regions with the reddest dust emission). Since we do not want to leave that extra flux in the stellar map, we also linearly interpolate over all those masked regions (since they are typically of the order of a few pixels). The map of dust emission is then adjusted in this case, by taking the difference between the original image and the adjusted old stellar map.

Naturally, the colors of the sources are also recalculated after these slight flux modifications.

Several uncertainties combine to set the level of accuracy that we can obtain in our maps of the distribution of old stellar light, due to both systematic and measurement errors. For clarity, Table 2 summarizes the typical values of the different uncertainties involved.

Table 2.  Uncertainties Involved in the Pipeline^a

Uncertainty	Typical Value (mag)
Error of original global color	0.001
ICA bootstrapping error	0.074
Propagated photometric error $[3.6]-[4.5]{| }_{s1}$	0.028
2nd iteration change on $[3.6]-[4.5]{| }_{s1}$	0.036

Note.

^aTypical values understood as the median values of the distributions.

Download table as:  ASCIITypeset image

The 3.6 and $4.5\;\mu {\rm{m}}$ images on which we base our separation have photometric uncertainties, which propagate through the ICA method into the final maps we produce. The original photometric uncertainties are quantified via the sigma maps (${\sigma }_{1}$, ${\sigma }_{2}$) that we obtain from the S⁴G weight maps and, according to Equation (1), they propagate in to the following uncertainties:

$$\begin{eqnarray}{\sigma }_{11}^{2} & = & \displaystyle \frac{{\mathrm{ZP}}^{2}}{{10}^{0.8({cs}1)}-{10}^{0.8({cs}2)}}\left[({\sigma }_{1}^{2}+{\rm{\Delta }}{\mathrm{sky}}_{1}^{2})\right.\\ & & +\left.\displaystyle \frac{{10}^{0.8({cs}2)}}{{\mathrm{ZP}}^{2}}({\sigma }_{2}^{2}+{\rm{\Delta }}{\mathrm{sky}}_{2}^{2})\right],\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}{\sigma }_{12}^{2} & = & \displaystyle \frac{{\mathrm{ZP}}^{2}}{{10}^{0.8({cs}2)}-{10}^{0.8({cs}1)}}\left[({\sigma }_{1}^{2}+{\rm{\Delta }}{\mathrm{sky}}_{1}^{2})\right.\\ & & +\left.\displaystyle \frac{{10}^{0.8({cs}1)}}{{\mathrm{ZP}}^{2}}({\sigma }_{2}^{2}+{\rm{\Delta }}{\mathrm{sky}}_{2}^{2})\right],\end{eqnarray} \tag{ 4 }$$

where $\mathrm{ZP}=280.9/179.7$, ${cs}1=[3.6]-[4.5]{| }_{s1}$, ${cs}2=[3.6]-[4.5]{| }_{s2}$, and ${\rm{\Delta }}{\mathrm{sky}}_{1}$ and ${\rm{\Delta }}{\mathrm{sky}}_{2}$ refer to the uncertainties in the determination of the sky for a given galaxy in each of the bands. In particular, there are two sky uncertainties computed within S⁴G Pipeline 3: one reflects Poisson noise, while the other refers to the large-scale background errors (rms calculated within and among different sky boxes; Muñoz-Mateos et al. 2013, see also Muñoz-Mateos et al. 2009 for a detailed discussion). Our global estimate of the sky error is calculated as ${\rm{\Delta }}\mathrm{sky}=\sqrt{{\rm{\Delta }}{\mathrm{sky}}_{\mathrm{Poisson}}^{2}+{\rm{\Delta }}{\mathrm{sky}}_{\mathrm{large}\;\mathrm{scale}}^{2}}$, which has a typical (median) value of 0.013 MJy/sr for both bands. Here, we neglect any uncertainties associated to the photometric zero points.

There is also uncertainty intrinsic to the ICA method, namely the reliability or uniqueness of a solution identified in any given measurement set. Our algorithm includes a quantification of the uniqueness of the solution by performing ICA on each measurement set N times, each time with a different initial seed matrix of mixing coefficients (bootstrapping). Perturbations to the initial seed are fixed for all galaxies in the sample and represent N possible realistic mixtures of old stars and dust emission. Optimally, ICA quickly converges to its final solution independent of the initial seed (Meidt et al. 2012a). The range of final mixing coefficients sets the uncertainty on the final average solution, and this ICA color uncertainty corresponds typically to values of order ∼0.07 mag. The ICA tests presented in the appendices confirm the meaningfulness of these uncertainties.

The s1 color change associated with the second iteration is typically of the order of ∼0.04 mag. Interestingly, in 88% of the cases, this color change is well constrained by the initial ICA bootstrap uncertainty (i.e., the change occurs within the s1 color error bar). Taking this into account, and considering the typical values listed in Table 2, we can state that the ICA boostrap error is the dominant source of uncertainty, and provides a good estimate of how accurate our solutions are (see Appendix B).

To account for the varying quality of solutions, flags are provided along with the final products. Flagging has been independently performed by two of the authors (Sharon Meidt, Miguel Querejeta) and an external classifier (Emer Brady), according to well-defined criteria.

• 1.  
  Is the physical distribution of the identified dust in accordance with the signatures that appear in the color map?
• 2.  
  Are there significant oversubtractions in the map of old stellar flux?
• 3.  
  Are there any artifacts that prevent ICA from obtaining a correct solution?

Three quality flags have been established: 1—excellent, 2—acceptable, 3—bad, and the statistical mode is chosen as the final classification (cases in which all three classifiers disagree have been revised and individually discussed). Depending on the specific application of the mass maps, only those galaxies classified as either 1 or both 1 and 2 will be suitable; also, some of the solutions classified as 3—bad correspond to galaxies in which the input data was contaminated by artifacts (e.g., mux bleeds or saturated PSFs), and a personalized case-by-case treatment of the masks can potentially improve the quality of those solutions (but we have not done it to avoid introducing a subjective component, and for homogeneity of the pipeline results).

It should be noted that, within the group of more than 1200 galaxies that constitute our final science sample (Section 6), very few have the bad quality flag, "3" (only 3.5%). This is, to a great extent, because we have conservatively excluded from the analysis all galaxies with low S/N, which systematically lead to solutions of poorer quality.

Using the IRAC zero magnitude flux density at 3.6 and $4.5\mu {\rm{m}}$, ${\mathrm{ZP}}_{3.6\mu {\rm{m}}}=280.9$ Jy and ${\mathrm{ZP}}_{4.5\mu {\rm{m}}}=179.7$ Jy (Reach et al. 2005), and the corresponding absolute magnitude of the Sun, ${M}_{\odot }^{3.6}=3.24$ mag (Oh et al. 2008), we obtain the following relationship for the $3.6\mu {\rm{m}}$ IRAC band:

$$\begin{eqnarray}&&1\;\mathrm{MJy}/\mathrm{sr}=704.04\;{L}_{\odot }/{\mathrm{pc}}^{2}.\end{eqnarray} \tag{ 5 }$$

Starting from an ICA-corrected flux surface density measurement ${S}_{3.6\mu {\rm{m}}}$ in Spitzer units of MJy/sr, the stellar mass contained by a pixel (0.75'') can be obtained as

$$\begin{eqnarray}&&\displaystyle \frac{M}{{M}_{\odot }}=9308.23\times \left(\displaystyle \frac{{S}_{3.6\mu {\rm{m}}}}{\mathrm{MJy}/\mathrm{sr}}\right)\times {\left(\displaystyle \frac{D}{\mathrm{Mpc}}\right)}^{2}\times \left(\displaystyle \frac{M/{L}_{3.6\mu {\rm{m}}}}{{M}_{\odot }/{L}_{\odot }}\right).\end{eqnarray} \tag{ 6 }$$

For a detailed discussion on how to choose the optimal $3.6\;\mu {\rm{m}}$ M/L, the reader is referred to Meidt et al. (2014) and Norris et al.. Here we adopt a single M/L = 0.6 (assuming a Chabrier IMF), which according to both sets of authors can convert the $3.6\;\mu {\rm{m}}$ old stellar flux (with dust removed) into stellar mass with an accuracy of ∼0.1 dex. Given that the dependence on age and metallicity of the M/L at $3.6\;\mu {\rm{m}}$ is relatively small for old stellar populations, Meidt et al. (2014) advocate for this constant M/L and its uncertainty assuming a universal age–metallicity relation, together with the constraint on metallicity (and thus age) provided by the [3.6]–[4.5] color. Following an independent, empirical strategy, Norris et al. argue for a comparable value, without invoking such an argument. A M/L of 0.6 is also found to be representative in stellar population synthesis models, extended to the wavelength range of 2.5–5 μm using empirical stellar spectra by Röck et al.

To give users the opportunity to choose their preferred distance and M/L, we release the map of the old stellar flux (map s1), along with the dust map (s2). To allow for the choice of a spatially varying M/L, we also provide a color map for the old stars. This color map conserves the original [3.6]–[4.5] color in all pixels without non-stellar emission (i.e., where s2 = 0 after thresholding), and $[3.6]-[4.5]{| }_{s1}$ elsewhere. Additionally, the quality flags, colors and integrated fluxes of each component are made public in a table format.

Our recommended post-processing strategy is explained in Section 4.6, which includes interpolation over masked areas. However, a myriad of interpolation techniques exist, and some users may even prefer to leave all the original flux unchanged for those regions. Therefore, we also provide the masks used, making it possible for different strategies to be applied. In a final step, for aesthetical purposes, dust maps have been smoothed in the areas of significant dust emission using a Gaussian kernel of $\sigma =2\;\mathrm{px}$ and conserving total flux between s1 and s2, but unsmoothed maps are available upon request. The release of Pipeline 5 products takes place on the NASA/IPAC Infrared Science Archive (IRSA).

Here we explore how the stellar flux fraction (s1/total) changes as a function of Hubble type (Figure 8). We can see that the observed $3.6\;\mu {\rm{m}}$ flux becomes a poorer tracer of stellar mass as the Hubble type increases, because dust emission becomes more significant. This contamination from dust increases slightly toward late-type galaxies, as shown by the running medians, and seems to decrease again for the latest Hubble types.

The reason for this final decline in contamination is probably due to the fact that very late-type galaxies tend to have less dust as they have lower metallicity. For standard spiral galaxies (Sa-Sc), the stellar component contributes ∼70%–80% of the flux (i.e., as much as 20%–30% is coming from dust), with both upper and lower outliers. Moreover, these are only global estimates, and these values can be significantly higher in individual star-forming regions, as shown by Meidt et al. (2012a). To account for completeness, we include as lower limits the galaxies to which ICA was not applied due to original blue colors (compatible with old stars, which implies little dust emission): they are shown as gray arrows marking the upper limit of 15% contamination that is estimated in Appendix A (and the arrow thickness is proportional to the number of galaxies per bin).

Zoom In Zoom Out Reset image size

In an attempt to investigate the main driver behind the varying flux fraction due to old stars at $3.6\;\mu {\rm{m}}$, in Figure 9 we plot that same quantity (${F}_{3.6,\mathrm{stars}}/{F}_{\mathrm{total}}$) as a function of the specific star formation rates (SSFRs), and find a declining correlation. We calculate SSFRs as the SFRs divided by our stellar mass, which is derived from the $3.6\;\mu {\rm{m}}$ flux corrected for dust with a M/L = 0.6, and integrated within the 25.5 mag arcsec⁻² isophote. The SFRs come from IRAS 60 and $100\;\mu {\rm{m}}$ fluxes (the weighted average of the values reported in NED), following the recipe from Larsen & Richtler (2000). Therefore, the subsample shown in Figure 9 is all S⁴G galaxies with ICA solutions compatible with old stars (within their uncertainties) and with availability of SFRs, which makes a total of 819 galaxies. Galaxies with higher SSFRs have on average higher fractional contamination from dust, which makes the fractional flux due to stars drop (see running median in Figure 9).

Zoom In Zoom Out Reset image size

Finally, there is no obvious trend of dust emission fraction with SFR alone: the SFR relative to the total mass of the galaxy (SSFR) is a better first-order indicator of how much flux is arising from dust and not from old stars at $3.6\;\mu {\rm{m}}$.

We have successfully isolated for old stellar emission the $3.6\;\mu {\rm{m}}$ images of more than 1600 galaxies in S⁴G, which provides us with an unprecedentedly large sample of stellar mass estimates for nearby galaxies. This statistically powerful tool allows us to produce an empirical calibration of the M/L and its scatter as a function of observed original color of the galaxy.

Our empirical effective M/L as a function of [3.6]–[4.5] global color is shown in Figure 10. This is the M/L required to convert the uncorrected 3.6 micron flux—including both stellar and non-stellar emission—into stellar mass, calculated by dividing the stellar mass implied by the s1 map for each galaxy by the total original luminosity in the analysis area. The stellar mass is calculated by multiplying the total luminosity in the s1 map by a constant mass to light ratio $0.6{M}_{\odot }/{L}_{\odot }$ (Meidt et al. 2014), and divided by the total original luminosity within our area of analysis. It is important to note here that the abscissa is the original integrated color that would be available, for instance, to an observer that cannot spatially resolve a galaxy; this is different from the weighted mean color (based on the pixel-to-pixel photometric uncertainty) that we used to discriminate which galaxies to apply ICA to in the first step of our pipeline.

Zoom In Zoom Out Reset image size

Overplotted in Figure 10 as a thick, continuous line is the regression to the data points in the range $-0.1\lt [3.6]-[4.5]\lt 0.15$, which provides the calibration

$$\begin{eqnarray}&&\mathrm{log}(M/L)=-0.339(\pm 0.057)\times ([3.6]-[4.5])\\ &&-0.336(\pm 0.002).\end{eqnarray} \tag{ 7 }$$

We can also express this in the following equivalent form, which provides stellar mass in terms of the measured 3.6 and 4.5 $\mu {\rm{m}}$ flux densities (expressed in Jy) and the distance (in Mpc):

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{*}}{{M}_{\odot }}={10}^{8.35}\times {\left(\displaystyle \frac{{F}_{3.6}}{\mathrm{Jy}}\right)}^{1.85}\times {\left(\displaystyle \frac{{F}_{4.5}}{\mathrm{Jy}}\right)}^{-0.85}\times {\left(\displaystyle \frac{D}{\mathrm{Mpc}}\right)}^{2}.\end{eqnarray} \tag{ 8 }$$

The decreasing trend in Figure 10 suggests that color is indeed a good first-order indicator of the level of contamination and, therefore, of the M/L for a given galaxy. The impact of contamination on the M/L is demonstrated by the color-coding based on SSFRs: since high SSFRs are typically associated with redder galaxies (see color coding in Figure 9), this provides a physical connection to the trend found in Figure 10. Galaxies in which star formation is relatively more prominent (higher SSFRs) have more dust emission in these bands, leading to redder colors and thus also requiring a lower M/L to convert the light into mass.

While Figure 10 suggests that the color can be used to constrain the effective M/L, there is still considerable scatter at fixed [3.6]–[4.5] color. We therefore recommend adopting an uncertainty of∼0.2 dex on our empirical relation, which also accounts for the 0.1 dex uncertainty associated with our adopted stellar M/L (see Section 5.1). (We note, however, that a method that can correctly remove dust emission, such as the one that we have presented here, should be the preferred approach for esimating stellar masses.)

Our relation between the effective M/L and [3.6]–[4.5] color, calibrated using our optimal stellar mass estimates, compares well with other relations in the literature. As part of the WISE Enhanced Resolution Galaxy Atlas, Jarrett et al. (2013) also found a declining correlation between W1 and W2 color (the 3.4 and $4.6\;\mu {\rm{m}}$ WISE bands) and M/L at $3.6\;\mu {\rm{m}}$, based on their observational stellar mass estimates of a few (∼15) galaxies. Figure 11 shows that the slope that we find is considerably shallower, based on a sample nearly 100 times larger. In this plot we show ${\rm{\Delta }}{\text{}}M/L$ in order to highlight the IMF-independent slopes of each relation and avoid differences in vertical offset due to choice of the IMF.

Zoom In Zoom Out Reset image size

We also find good agreement with the stellar M/L estimated by Eskew et al. (2012). Using a different method, based on spatially resolved star formation histories (SFHs) in the LMC, they calibrated a conversion between 3.6 and $4.5\;\mu {\rm{m}}$ fluxes and stellar mass. This can be translated into an equivalent color-dependent M/L ($\mathrm{log}{\text{}}M/L=-0.74$ $([3.6]-[4.5])-0.236$), which is overplotted on Figure 11, also normalized so that log(M/L) = 0 at [3.6]–[4.5] = 0.

It is remarkable that, in spite of the very different method used, their declining slope is very similar to the fit to our data. Finally, it is also worth mentioning that McGaugh & Schombert (2014) recently point to a constant $M/{L}_{3.6\;\mu {\rm{m}}}=0.47$, with a scatter of 0.1 dex, which is also compatible with our result.

We conduct our tests by creating two-dimensional realistic distributions of the emission from old stars and dust. As input models we utilize the spatial distributions of the components s1 (old stellar light) and s2 (dust emission) identified with ICA for four galaxies that cover a range in morphology. For these specific galaxies, the s1 and s2 components were shown to be representative of old stars and dust in Meidt et al. (2012b). Based on these, our model input maps at 3.6 and $4.5\;\mu {\rm{m}}$ are generated in the following way: (a) the color of old stars is fixed at [3.6]–[4.5] = −0.12; (b) the color of dust is varied between [3.6]–[4.5] = 0 and [3.6]–[4.5] = 0.6; (c) the fractional contribution of dust emission to the total flux is varied from 5% to 40%.

In each case, we verified that the resulting colors of each of the combined s1+s2 models are consistent with the observed colors in the S⁴G sample, which cover the range $-0.05\lt [3.6]-[4.5]\lt 0.20$ (if we exclude the 5% tail at the upper and lower ends).

We run a total of 60 tests per galaxy, covering the total dust color range $0.2\leqslant [3.6]-[4.5]\leqslant 0.7$ in steps of 0.1, and varying the fractional contribution of dust by 10%, 20%, 30%, and 40%. We also explore the dependence on the dust spatial distribution by repeating each battery of experiments on four different galaxy models; for reference, models 1 through 4 are based on NGC 2976, NGC 3184, NGC 4321, and NGC 5194, respectively.

A.2.1. Relative Color Difference Between s1 and s2

A.2.2. Fractional Contribution of Dust Emission

The difference between the stellar and dust colors imposed on each model and the s1 and s2 colors found by ICA are shown in Figure 12, as a function of the dust color of the model. To improve the visibility of the plotted data points, we introduce a small artificial random scatter around each discrete dust color value (a random dithering within ±0.05). The symbol size of each point is proportional to the dust fraction, while the colors denote the original galaxy on which the models are based (highlighting the effect of two-dimensional dust distribution). Most striking is the clear trend for increased scatter in the ICA-determined dust colors (and stellar colors, although less so) when the "intrinsic" dust colors are redder. ICA underestimates the dust color more severely for low contamination fractions, and, conversely, it tends to overestimate the stellar color (but not as much) at the highest contamination fractions.

Zoom In Zoom Out Reset image size

However, as shown in Figure 13, this large uncertainty in dust and stellar colors does not necessarily imply larger uncertainties in the estimation of fluxes. Plotted there is the difference between the stellar flux imposed on the model and the total flux in the s1 map retrieved by ICA for all possible combinations of dust colors and contamination fractions in our model grids. Flux uncertainties in both components are sensitive to the relative stellar and dust fluxes (i.e., the level of dust contamination) but this tends to be balanced by an additional dependence on the actual difference between the stellar and dust colors. As we are ultimately interested in the error in the corrected stellar fluxes found by ICA, Figure 13 suggests that ICA can perform quite well. We note that many of the models exhibit original (weighted) mean [3.6]–[4.5] colors that are negative, so that we would not have applied ICA to them (see Section 4.2). Here we see that it is reasonable to exclude such galaxies from the analysis given that their low contamination fraction and/or very low dust colors (both of which can produce original colors $[3.6]-[4.5]\lt 0$) lead to the highest uncertainties.

Zoom In Zoom Out Reset image size

Figure 13 also shows that most models lead to final uncertainties in flux within 10%, with some outliers covering a band up to 20% of uncertainty. However, these are mainly associated with the highest contamination fractions in our models (40%), and such extreme contamination from dust can be only very rarely expected (it is the upper limit for star-forming spirals, as argued in Section A.2.2).

Additionally, these tests confirm that the ICA uncertainties determined through the bootstrap method are reliable measures of the errors in color. This bootstrap uncertainty is obtained by running ICA based on 48 different seeds (different starting points in s1 and s2 colors). According to our models, in 94.6% of cases the difference in stellar color was well confined by the bootstrap uncertainty, confirming its meaningfulness.

To test the impact of a third component on ICA solutions assuming two components (and to confirm the usefulness of our second ICA iteration), we have included an additional component with an intrinsically redder color representative of hot dust in H ii regions, i.e., [3.6]–[4.5] = 1. For the spatial flux distribution, we use truncated De Vaucouleurs profiles (de Vaucouleurs 1948), and place them randomly in the regions of the galaxy model where the more extended dust emission is significant (higher than the mean dust flux). We vary the total number of regions (10, 20, 30) and their relative brightness (5%, 10%, 20% of total dust flux). We also explore a range of colors, starting from the color of the spatially extended dust itself, up to a maximum color of 0.6.

When run through our pipeline, we can indeed see that our simulated H ii regions get effectively masked, and the second iteration leads to a dust color closer to the nominal spatially extended "diffuse dust" in 94.7% of the experiments (see Figure 14). Moreover, the bootstrap uncertainty is effectively reduced after the second iteration (from ∼0.8 to ∼0.2 mag), suggesting that our measure of errors is sensitive to the problem posed by a third component of intrinsically different color, which is relieved by the iterative ICA approach that we have constructed.

Zoom In Zoom Out Reset image size