INFRARED SPECTRAL ENERGY DISTRIBUTIONS OF SEYFERT GALAXIES: SPITZER SPACE TELESCOPE OBSERVATIONS OF THE 12 μm SAMPLE OF ACTIVE GALAXIES

IRAC observations were centered on the NED coordinates of the sample galaxies based on catalog names listed in Rush et al. (1993). A beam splitter and supporting optics centers the target on two detectors simultaneously (Fazio et al. 2004), either at 3.6 and 5.8 μm or 4.5 and 8.0 μm, and so each observation consisted of two pointings at a common orientation of the focal plane relative to sky. These observations were taken as snapshots with no attempt to mosaic or dither. To guard against saturation, we used the high-dynamic range (HDR) mode which provides 0.6 s and 12 s integrations at each pointing.

The data were initially processed and calibrated by the IRAC basic calibration data (BCD) pipeline, version S14.0. The pipeline performs basic processing tasks, including bias and dark current subtraction; response linearization for pixels near saturation; flat-fielding based on observations of high-zodiacal background regions; saturation and cosmic ray flagging, the latter indicated by signal detections more compact than the point-spread function (PSF); and finally flux calibration based on observations of standard stars.¹² The nominal photometric stability for compact sources is better than 3% in all detectors (Reach et al. 2005). Corrections for extended sources¹³ are accurate to ∼10%, but the contribution of extended emission and the resulting correction is small for most of the sample galaxies. Color corrections are typically much greater, especially in the 8.0 μm channel where in-band PAH emission can result in factors of 2 or greater corrections.

The data were further processed for remaining artifacts, including cosmic rays and detector artifacts not corrected by the BCD pipeline. The steps performed for artifact removal and photometric extraction are detailed below in Sections 3.1.1–3.1.4. Figure 3 illustrates the effect of our artifact removal techniques.

Zoom In Zoom Out Reset image size

The IRAC photometry is listed in Table 2. Detailed below, the photometry is presented in the IRAC magnitude system with zero-point flux densities 280.9, 179.7, 115.0, and 64.13 Jy for the 3.6, 4.5, 5.8, and 8.0 μm channels, respectively (Reach et al. 2005). The photometry includes corrections for extended emission and color corrections.

3.1.1. Cosmic Ray and Bandwidth-effect Mitigation

3.1.2. Bias Artifacts

3.1.3. Saturation and Distortion

3.1.4. Photometric Extraction

IRAC observations of 16 sample galaxies were obtained through public release to the Spitzer archive. Data processing largely followed the techniques described above for our observations, except that, where our observations comprise snapshots, the archival observations all employed dithering or mosaicking techniques. The archival data were initially mosaicked using standard procedures with MOPEX, including background matching based on overlaps.

Artifacts, particularly banding, seriously affected the matching of mosaic overlap regions and so we developed a modified mosaicking procedure building on the algorithm of Regan & Gruendl (1995). We modified their least-squares approach to employ least absolute deviations as a more robust measure of the quality of overlap matching, and we further masked bright source and artifact regions prior to overlap matching. The final images were assembled using a median stack of astrometrically aligned images, background subtracted, and finally sub-imaged to roughly 5' square to match our survey observations.

Residual bias artifacts were removed by computing column or row biases where sufficient background was available for bias determination. Cosmic rays were largely mitigated by the median stack, but residual bad pixels were identified and flagged interactively, as described in Section 3.1.1.

We observed the sample galaxies using the Spitzer IRS in spectral mapping mode and using the Short-Low (SL) and Long-Low (LL) modules. The modules cover a wavelength range of ∼5–36 μm with resolution R = λ/Δλ ranging from 64 to 128. The integration times were 6 s per slit pointing.

The spectral maps were constructed to span >20'' in the cross-slit direction and centered on the target source. The observations were stepped perpendicular to the slit by half slit-width spacings. For the SL data, the mapping involved 13 observations stepped by 18 perpendicular to the slit, and, for the LL data, 5 observations stepped by 525 perpendicular to the slit. The resulting spectral cubes span roughly 252 × 546 × (5.3–14.2 μm) for the SL data and 291 × 151''× (14.2–36 μm) for the LL data.

The raw data were processed through the Spitzer BCD pipeline, version S15.3.0. The pipeline handled primary processing tasks including identification of saturated pixels; detection of cosmic ray hits; correction for "droop," in which charge stored in an individual pixel is affected by the total flux received by the detector array; dark current subtraction; and flat-fielding and response linearization. Details are provided in the IRS Data Handbook.¹⁶

Sky subtraction used off-source orders. In SL and LL observations, the source is centered in the first- or second order separately, with the "off-order" observing the sky at a position offset parallel to the slit: 79'' away for SL and 192'' away for LL. Sky frames were constructed using median combinations of the off-source data and subtracted from the on-source data of matching order. No detectable contamination of the sky frame appeared in our observations based on inspection of the sky frame data. In addition, a portion of the first-order slit spectrum appears on the second-order observation and is used as an additional check of spectral features that appear near the first-/second-order spectral boundary.

Data cubes were constructed by first registering the individual, single-slit spectra onto a uniform grid (slit position versus wavelength). The IRS slits do not align with the detector grid, and the detector pixels undersample the spatial resolution. Interpolating the undersampled slit onto a uniform grid produces resampling noise (Allington-Smith et al. 1989), which significantly impacts the spectra of unresolved sources. We instead used the pixel re-gridding algorithm described by Smith et al. (2007a) that minimizes resampling noise. Pixels from the original image are effectively placed atop the new, registered pixel grid. Uniform surface brightness is assumed across the original pixel, and that original signal is weighted and distributed based on the fractional overlap with pixels on the new grid. After registering the single-slit frames, the data were re-gridded by bilinear interpolation, one wavelength plane at a time, into the final data cube. Flux uncertainty images were similarly processed, with modifications to accommodate variance propagation, to produce an uncertainty cube.

The spectra were extracted from synthetic, 20'' diameter circular apertures centered on the brightest, compact IR source nearest the target coordinates. Fractional pixels at the edge of the aperture were accounted for by assuming uniform surface brightness across the pixel and weighting by the area of intersection between the pixel and aperture mask.

The cube-extracted spectra were optimally weighted based on a three-dimensional modification of Horne's (1986) two-dimensional slit extraction algorithm. Successful application of this algorithm requires estimation of spatial profiles P_xλ, the probability that a detected photon falls in a given pixel x rather than some other pixel on the spectral map at wavelength λ. Optimal extraction requires that the fractional uncertainty of the estimator for P_xλ is less than the fractional uncertainty of the original spectral image. Generation of P_xλ therefore requires some smoothing along the wavelength axis. To help preserve real variations of P_xλ with wavelength, we employed Savitzky-Golay (1964) polynomial smoothing. In this smoothing scheme, P_xλ is calculated by a polynomial fit to the spectrum at pixel x over the fitting window δλ either centered on λ or limited by the ends of the spectrum. We performed trial-and-error smoothing experiments on CGCG381-051, which shows a relatively weak continuum at short wavelengths but high EQW PAH features, to decide on the polynomial order and window smoothing parameters; ultimately we selected quadratic polynomials and 5-pixel smoothing windows to balance improved signal to noise and the tracking of real variations of P_xλ with wavelength.

The effect of optimal extraction is illustrated in Figures 4 and 5, which compare the optimally extracted and non-weighted extraction of the IRS SL spectrum of NGC 3079 and F01475-0740, and Table 3, which compares the measurements derived from these extractions (see Section 4.1 for a description of the measurement technique). NGC 3079 presents a challenging case because it is edge-on and shows bright, extended PAH emission. The spatial profile is therefore a complex function of wavelength, and coarse smoothing in the wavelength direction could potentially affect the measurement of the PAH fluxes and EQWs. We find however that the fractional difference between the optimally weighted spectrum and the unweighted spectrum is typically only a few percent, comparable to the statistical uncertainties in the unweighted spectrum. Furthermore, the measured fluxes and EQWs agree to within the measurement uncertainties; in fact, the measurement uncertainties are dominated by the systematics of the measurement technique.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared to NGC 3079, F01475-0740 is a compact source at relatively low signal to noise. Figure 5 clearly illustrates the advantage of optimum weighting. The continuum shape and PAH spectral features are preserved, but the formal statistical uncertainties are reduced by a factor of 2–3 over the SL spectral range. Again, the fluxes and EQWs of lines agree to within the measurement uncertainties, but the optimally weighted spectrum produced significant (>3σ) detections of PAH 6.2 μm, Ne v 14.3 μm, and S iii 18.7 μm that are too faint for the unweighted spectrum.

Fringing at the 5%–10% level was apparent in the LL observations, particularly those of sample objects with a strong point-source contribution. The spectra undersample the fringing, and so the fringe pattern could not be removed using a conventional filtering in Fourier space. We employed instead the technique described by Kester et al. (2003), which involves fitting sinusoids in wavenumber to line-free sections of the spectrum. Successive fringe components are added to the model until the next model fringe amplitude falls below the noise level of the spectrum. The effect of this technique is illustrated in Figure 6.

Zoom In Zoom Out Reset image size

The flux scale was calibrated against archival staring and spectral mapping observations of the stars HR7341 (SL) and HR6606 (LL). The staring mode spectra were extracted using SPICE,¹⁷ and the spectral mapping observations were processed into cubes as described above for our sample galaxies. The flux calibration curve was determined by the ratio of the flux-calibrated staring mode spectra to the uncalibrated, cube-extracted spectra and smoothed by a polynomial fit.

No correction was produced for resolved, extended emission. The flux calibration includes an aperture loss correction factor appropriate for compact sources, but extended sources do not suffer as much aperture loss. Therefore, wavelength bands containing extended emission relative to the PRF will be calibrated systematically high. Discussed below in Section 3.6, the systematic error introduced is of order 10%, particularly near the 8 μm PAH features.

We note that these IRS observations have been independently and differently processed and recently presented in Wu et al. (2009). The main differences in data processing are (1) Wu et al. used the CUBISM tool and its calibrations for spectral cube construction (Smith et al. 2007a); (2) it is not clear that Wu et al. extracted the spectra optimally; (3) a rectangular aperture, which covers a similar solid angle to our extraction aperture, was commonly used (204 × 153; larger for SINGS galaxies); (4) there appears to have been no attempt to de-fringe the long wavelength (LL) data; (5) for an unknown number of sources, Wu et al. had to scale the SL data to match up with the flux calibration of the LL data. We address these differences in turn, excepting (3), the minor difference of extraction aperture geometry.

We found that (1) CUBISM forced the data onto a new grid in which the original pointing center was shifted, commonly resulting in nuclei that were centered between pixels rather than on a pixel. This shift apparently introduced a systematic error in the resampling of the surface brightness, because we noticed artifacts and a few percent reduced signal in the spectrum extracted from cubes produced by CUBISM. Our cubing technique was designed to require minimal resampling of the observational grid and mitigated these artifacts.

We also found that (2) non-optimal extractions decreased the signal to noise of the extracted spectra of fainter sources by factors of ∼50% or more (see Figure 5 and the discussion above). Although it varies source to source, (4) fringing can produce ∼10% artifacts between 14 and 36 μm (cf. their extraction of MRK 1239 and M-6-30-15, among others). Finally, (5) we did not need to apply any additional scaling to the spectra extracted from SL cubes; rather our calibrations of the SL module show excellent agreement with both the IRAC measurements and LL measurements, even though the calibrations for each respective Spitzer module were based on different calibration stars; this result suggests that the calibration against spectrally mapped stars, as presented here, was more robust for calibration and systematic corrections, accepting the ∼10% calibration uncertainty for sources dominated by extended emission.

Archival data of our sample galaxies, such as obtained by the SINGS collaboration, were processed in the same manner as for our program data. Spectral mapping data were available in LL for all sources, but there are however a few sources for which only SL staring mode observations are available: NGC 526A, NGC 4941, NGC 3227, IC 5063, NGC 7172, and NGC 7314. The staring mode spectra were extracted using SPICE and subjectively scaled to achieve the best match to align with LL measurements and, where possible, IRAC 5.8 and 8.0 μm measurements. Not surprisingly, there result disagreements between the IRAC measurements and scaled, IRS staring observations owing to the contribution of the host galaxy in the 20'' aperture.

The sample galaxies were observed using the 70 μm low-resolution spectrometer (SED mode) of the MIPS instrument. The spectrometer provides spectral resolution R ∼ 15–25 between λ55 and 95 μm. The slit dimensions are 20'' × 120''.

Each observation consisted of three pairs of on-source and off-source measurements, where the off-position was located 1'–3' from the on-source position. Integration times were 3 or 10 s depending on the IRAS 60 μm flux density of the source. The observations include measurements of built-in calibration light sources, called stimulators, for flux calibration.

We used post-BCD products of the S14.4.0 pipeline; the pipeline processing includes flux calibration, background subtraction, co-addition of the on-source pointings, and uncertainty images. Details are provided in the MIPS data handbook.¹⁸

The MIPS spectra were extracted by summing across a 29''-wide synthetic aperture centered on the target source (three-column extraction). The spectra were further corrected for aperture losses assuming a point-source model; the present extractions include no corrections for spatially extended emission. The photometric accuracy is expected to be 10% for compact sources and ∼15% for extended sources (Lu et al. 2008).

Note that the MIPS aperture covers a solid angle ∼1.8× larger than the IRAC and IRS apertures. Modeling SEDs of sources with extended (>20'') far-infrared emission will therefore require an (unknown) aperture correction that is potentially much greater than the reported calibration uncertainty. Candidates for MIPS aperture corrections will show a jump in the MIPS flux relative to a suitable extrapolation of the IRS spectrum, although a real spectral peak in the 40–50 μm range might similarly result in an observed MIPS SED that appears too blue even in the absence of aperture effects.

The IRS spectra overlap with the 5.8 and 8.0 μm IRAC channels. To compare the relative spectrophotometry, we interpolated the IRS spectra to find the flux densities at the effective wavelengths of the color-corrected IRAC data, 5.731 μm and 7.872 μm, respectively. The average relative difference of the overlapping flux densities, F_ν(IRS)/F_ν(IRAC) − 1, are 4 ± 1% at 5.8 μm and 6 ± 1% at 8.0 μm. The frequency distributions are shown in Figure 7. For comparison, the average standard scores, (F_ν[IRS] − F_ν[IRAC])/σ, are 0.00 ± 0.03 and 0.18 ± 0.05.

Zoom In Zoom Out Reset image size

These results indicate that the IRS flux densities are, on average, systematically higher than the IRAC measurements by a few percent. However, the excess at 5.8 μm is dominated by statistical uncertainties, as the average standard score is consistent with zero.

Strong, extended PAH emission however contributes significantly to the 8.0 μm IRAC channel. This spectral feature is difficult to color correct accurately owing to numerical integration errors that arise at sharp spectral features. In addition, PAH emission is often spatially resolved in this sample, and the present IRS calibration includes no correction for extended emission. The systematic offset of 6% may result from the overcorrection of aperture losses for extended PAH sources. The IRAC photometry of extended sources is moreover accurate to only ∼10%.

Figure 8 shows the 8 μm IRS/IRAC flux ratio as a function of the fractional contribution of the central point source to the 20'' nuclear aperture. The fainter sources provide appreciable scatter, but the brighter sources reveal a trend in which point-source-dominated objects show better agreement, but more extended objects present ∼10% IRS excesses. This trend supports the interpretation that residual extended source calibration uncertainties are the primary cause for the IRS–IRAC discrepancies at 8 μm.

Zoom In Zoom Out Reset image size

In addition to continuum radiation from dust grains and associated Sil emission and absorption bands, the infrared spectrum of Seyfert galaxies includes diagnostics such as fine-structure lines tracing a range of ionization states, H₂ lines, and PAH emission. We used the spectrum fitting tool PAHFIT (Smith et al. 2007b), which is tailored to low-resolution IRS spectroscopy. The fits to the SEDs of NGC 4151 and NGC 7213 are provided in Figure 9 as examples of the decomposition. The results are summarized in the following tables: integrated PAH fluxes are given in Table 4; PAH EQWs in Table 5; H₂ line fluxes (mostly upper limits) in Table 6 and EQWs in Table 7; ionic fine-structure line fluxes in Table 8, and their EQWs in Table 9. The line measurements are not corrected for model extinction, but for completeness we list the best-fit model dust opacity (normalized to 10 μm) in Table 10.

Zoom In Zoom Out Reset image size

As provided, PAHFIT is best suited for nearly normal or star-forming galaxies; its model includes thermal continuum radiation from dust grains, PAH features, fine-structure lines from lower ionization state species, H₂, and Sil absorption, whether by assumed mixed or foreground screen extinction. We added two components to the PAHFIT model to fit the Seyfert SEDs: (1) fine-structure lines from high-ionization states, such as [Ne v] and [Ne vi]; and (2), to fit silicate emission features, a simple model for warm dust clouds that are optically thin at infrared wavelengths.

By default, PAHFIT models extinction based on the dust opacity law of Kemper et al. (2004). Sirocky et al. (2008) found that the cold dust model of Ossenkopf et al. (1992) better matches the high 18 μm Sil/10 μm Sil absorption found in active ultraluminous infrared galaxies. From inspection, the present spectra similarly show relatively strong 18 μm features, whether in emission or absorption, and so we further modified PAHFIT to use the cold dust model of Ossenkopf et al. (1992).

The thin, warm dust model assumes clouds with simple, slab geometry and opacity at 10 μm, τ₁₀ < 1. These warm clouds are further assumed to be partially covered by cold, absorbing dust clouds. The model spectrum is then

$$\begin{eqnarray} F_{\nu } &=& (1 - C_f) B_{\nu }(T_W)(1 - e^{-\tau _W(\nu)}) + C_f B_{\nu }(T_W)\nonumber\\ &&\times\, (1-e^{-\tau _W(\nu)}) e^{-\tau _C(\nu)}, \end{eqnarray} \tag{ 1 }$$

where C_f is the covering fraction of cold, foreground clouds, modeled independent of the galaxy extinction; τ_W is the opacity through the warm clouds; τ_C is the opacity through the cold, foreground clouds, and B_ν(T_w) is the source function, for which we adopt a scaled Planck spectrum at (fitted) temperature T_W for simplicity. Note that the cold dust opacities described by Equation (1) are taken to be independent of the global PAHFIT dust opacity model; the opacity values listed in Table 10 are determined by the global dust opacity fit. The intention of including this additional model component is to provide a realistic continuum baseline for emission line fits and Sil strength determination rather than an interpretable, radiative transfer model, and details of the opacity law or more realistic source functions are beyond the scope of the present spectral decomposition.

We further modified PAHFIT better to accommodate the IRAC broadband and MIPS-SED measurements. Both instruments provide much lower sampling density in wavelength than IRS data, and consequently they receive lower weight in the χ² minimization procedure. We therefore employed the sampling weight correction described by Marshall et al. (2007), which effectively re-weights individual data points based on the sampling density local to that data point; regions of low sampling density receive increased weight. The weighting is normalized so that each data point carries, on average, unity sampling weight. We also introduced as a fitted parameter an aperture correction factor for the MIPS-SED wavelength range that boosts the model SED by a factor up to 1.81, which is the areal aperture ratio of the MIPS-SED extraction to the nominal 20'' diameter extraction aperture. Table 10 includes the best-fit aperture corrections.

There are a few caveats to the results of this decomposition. These low spectral resolution data are not best suited for measuring line fluxes, particularly in spectrally crowded regions. PAHFIT takes a conservative approach, severely restricting the centroid wavelengths and widths of the fitted lines; for example, the fine-structure lines are assumed to be unresolved and their widths are fixed at the instrumental resolution. Even so, the fine-structure lines near 35 μm are crowded, and furthermore the IRS measurements are very noisy at that region of the spectrum, with sensitivity ∼10× poorer than at 25 μm, depending on background ([http://ssc.spitzer.caltech.edu/documents/som/]Spitzer Observer's Manual). These limitations are reflected in the uncertainties and upper limits reported in Tables 4–9.

Extracted line fluxes were however affected by occasional, residually poor fits to the local continuum or PAH features. For example, if the PAHFIT model placed the continuum too low locally to an emission line, PAHFIT would grow the emission line to meet the data and therefore produce a line flux greater than observed. Similarly, the PAHFIT model might produce a local continuum that is too high and the line flux is reported too low. A good illustration of this problem is the too-high continuum model surrounding the [Ne v] λ14.3 μm line of NGC 4151 (Figure 9). To compensate for locally poor continuum models, the mean and rms of the model-subtracted spectrum was evaluated in spectral regions surrounding each line. Each fitted line peak was accordingly adjusted by subtracting the local, mean residual continuum. Line fluxes were similarly adjusted; line widths are unaffected as they were held fixed to the instrumental resolution during the fit.

There is a weak PAH feature near 14.3 μm that potentially contaminates the [Ne v] λ14.3 μm fine-structure line (Sturm et al. 2000). Deblending these features uses the fact that the PAH feature is somewhat broader (FWHM ∼0.4 μm) than the unresolved (FWHM <0.1 μm) [Ne v] line (Smith et al. 2007b). In the worst-case scenario, the code may fail to fit a weak but present PAH 14.3 feature resulting in an artificially enhanced [Ne v] line strength. This effect is at least partially ameliorated by the residuals analysis and is reflected in the large uncertainties of Table 8. The [Ne v] 14/24 ratio provides however a good check for contamination. This ratio is a density diagnostic with lower limit ∼0.9 (e.g., Alexander et al. 1999). We demonstrate in a companion paper (Baum et al. 2010) that for all of the sources where there is a [Ne v] 14.3 μm detection, the line ratio is consistent with the low density limit (cf. Sturm et al. 2002); contamination from PAH 14.3 would push the ratio to a (forbidden) value below the low density limit. We conclude that the PAH 14.3 μm feature does not significantly contaminate the [Ne v] λ14.3 μm line strengths in this study or that any contamination is within the uncertainties of the line strength.