NEAR-INFRARED DETECTION OF A SUPER-THIN DISK IN NGC 891

All observations were obtained with the WIYN³ 3.5 m telescope at Kitt Peak National Observatory during an observing run in 2011 October. We used the WIYN High-resolution InfraRed Camera (WHIRC; Meixner et al. 2010) NIR imager, and took data in the broadband J, H, and K_s filters. WHIRC has a 33 field-of-view and 01 pixel scale, which allows it to fully exploit the excellent seeing at WIYN (median seeing ∼07).

A log of the observations is presented in Table 1. Three regions of NGC 891 (the center, along the major axis NE of the center, and along the major axis to the SW) were observed to provide coverage on both sides of the galaxy out to ∼2.5 and 14 K-band thin disk scale-lengths and scale-heights, respectively (10 and 4.7 kpc; Xilouris et al. 1999). Enough overlap was left between these pointings to allow for combined mosaics to be produced. Multiple exposures were taken in each band at each pointing in a dither pattern (with offsets between dithers ∼15'' for J and H and 25'' for K_s) and averaged into a mosaic; the seeing estimates provided in Table 1 are for the image mosaics and come from Gaussian fits to foreground stars in the image using the IRAF⁴ task IMEXAM. In addition to the on-source data, we also intersperse dithers of a sky field to aid in sky background subtraction; the total exposure time of the sky dithers is of the same order as the total on-source exposure time. All observations were taken under photometric conditions. Mean seeing conditions were 08, 067, and 063 FWHM in J, H, and K_s, respectively. Full details of the data reduction can be found in Appendix A, and a false-color image of our final mosaics is shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Log of Observations

Filter	Date	Region	Exposure Time (s)		Seeing FWHM
	(UT)		Individual	Total
J ...	2011 Oct 18	SW	80	640	06
	2011 Oct 19	Center	80	640	09
	2011 Oct 19	NE	80	640	09
H ...	2011 Oct 18	SW	80	640	08
	2011 Oct 19	Center	80	640	07
	2011 Oct 19	NE	80	640	05
Ks ...	2011 Oct 18	SW	40	1080	09
	2011 Oct 19	Center	40	1080	05
	2011 Oct 19	NE	40	1080	05

Download table as:  ASCIITypeset image

We downloaded "post-basic" calibrated data of NGC 891 from the Spitzer archive⁵ in the 3.6 and 4.5 μm IRAC channels. The images were inspected for defects. Very low level surface gradients were found across the images and surface-subtracted using the iterative method described in Schechtman-Rook & Hess (2012) and Martinsson et al. (2013). The resolution of the IRAC data is ∼1.8 arcsec⁻¹ (FWHM), twice as large as the worst quality JHK_s data. At these wavelengths, however, the SB distribution of NGC 891 is fairly smooth and slowly varying. Therefore we interpolate the IRAC data to match the WHIRC mosaics, assuming that oversampling the IRAC images does not contribute significant additional uncertainties to our computed attenuation correction and will not significantly degrade the resolution of the mosaics.

Without knowing the 3D distribution of both dust and stars it is impossible to perfectly attenuation-correct images of NGC 891. However, we can compute a statistically accurate correction by modeling the dependence of color on the attenuation at multiple points in a simulated galaxy, which produces the average attenuation as well as the statistical uncertainty due to the varying distributions of stars and dust in different sight-lines. This correction is then broadly applicable for real galaxies with similar physical parameters (e.g., dust clumpiness) and orientations.

We base our models on the NGC 891 models of Bianchi (2008), with a bulge, two stellar disks (an old, thin disk and young super-thin disk), and a dust disk. We use the same physical scales and unattenuated stellar SEDs for all components as in Bianchi (2008), however instead of individual molecular cloud-sized clumps we choose to first place half of the total dust mass in fractal clumps. We choose not to model spiral structure: the rearrangement of dust and starlight might change the distribution of attenuation across the galaxy but should not effect the color resulting from a given amount of attenuation. Furthermore, the addition of spiral structure has a negligible effect on the integrated galaxy SED (Popescu et al. 2011); therefore ignoring spirality will not significantly alter the parameters of a given RT model. Both the percentage of dust in the fractal component as well as the details of the construction of the fractal grid follow from the median values for these parameters found in Schechtman-Rook et al.'s (2012a) modeling of NGC 891—briefly, roughly 50% of the total dust mass is in clumps. In our models here the clumpy component is divided up into 130 fractal clumps. We choose to use the stellar emission parameters from Bianchi (2008) instead of this more recent work because Bianchi (2008) fit the entire SED of NGC 891 whereas Schechtman-Rook et al. (2012a) were most concerned about the clumpy nature of the dust in absorption and therefore only explored models without dust emission. The young super-thin disk component is placed proportionately to the dust density in each grid cell. In this way the super-thin disk is not only clumpy, but the clumpy stellar emission is co-spatial with the dust clumps. In order to fit the integrated SED both in the UV and the MIR, we find that we must also set a minimum dust density threshold below which there is zero emission from young stars. This threshold mimics the effect of having a gas surface density threshold for star formation (Martin & Kennicutt 2001).

The unattenuated stellar SEDs (again, following Bianchi 2008) are set to that of an Sb stellar SED for a galaxy with no dust from Fioc & Rocca-Volmerange (1997). The light from stars with M < 3 M_☉ is used as the input for the bulge and thin disk, while light from stars with M > 3 M_☉ is assigned to the super-thin disk. The dust is included based on the prescription of Draine & Li (2007) with the mass fraction of PAHs q_PAH = 4.6%, with a size distribution to approximate R_V = 3.1 from Weingartner & Draine (2001). For distribution in the model galaxy, the dust is split into three components: large grains, with grain size a > 200 Å, ultra-small (PAH dominated) grains, with a < 20 Å, and small grains, with sizes between those of the large and ultra-small grains. These three components make up 80.63%, 5.86%, and 13.51% of the total dust mass in a given voxel (volume pixel), respectively. This combination of physical parameters we designate as Model A, with the full set of parameters used to create the model given in Table 2.

As seen in Table 2, the super-thin disk of the model has a scale-length double the size of the thin disk's scale-length but a scale-height half as large as for the thin disk. The thin disk has an axial ratio (h_R/h_z) = 10, common for spiral galaxies (Bershady et al. 2010b), while the super-thin disk has a much larger (h_R/h_z) = 40. However, due to the minimum dust density threshold and the addition of clumpiness, the axial ratio of the super-thin disk has less physical meaning than for the thin disk: at large distances from the center of the galaxy and/or the midplane the super-thin disk will be essentially truncated. The physical parameters of the thin disk and dust component are based on the results of Xilouris et al. (1999). Popescu et al. (2011) take a similar approach, basing dust and thin disk parameters on Xilouris et al. (1999), but they choose a super-thin disk scale-height based on that found in the MW. We also note that Popescu et al. (2011) have an extra component of diffuse dust with the same scale parameters as the super-thin disk.

For the model used here we mapped a 100 × 100 × 100 cell Cartesian grid to a 40 × 40 × 5 kpc volume, where the 5 kpc dimension represents vertical height above and below the disk.⁶ We chose to employ a rectangular grid rather than a cubical one to properly sample the vertical extent of the dust grid, which (as shown in Table 2) has a scale-height in our model of only 200 pc. In our previous work (Schechtman-Rook et al. 2012a), we used a RT code that only tracked scattering and absorption of photons emitted at a specific wavelength. The RT is performed by the HYPERION software package (Robitaille 2011), which computes the full SED of a model. HYPERION uses an iterative scheme to compute the dust temperature, and we use 100 million photon packets per iteration (an average of 100 photon packets per grid cell) for both for the temperature calculation iterations and to produce SEDs and images. We compare the output SED of the model to data in Figure 2. While the model has difficulty reproducing the MIR features associated with PAH emission and photo-dissociation regions, Bianchi (2008) indicate that this is likely due to a resolution effect in our models. Regardless of the reason for this discrepancy we are most concerned about the SED between 1 <λ (μm) < 5, where the model appears to fit the data well.

Zoom In Zoom Out Reset image size

To correct for the effects of dust on star-light we ran the model (described in Section 3.1) twice: once with and once without any dust. We produced images in J, H, K_s, 3.6 μm, and 4.5 μm for both model runs at the inclination of NGC 891 (897; Xilouris et al. 1999). The attenuation ${A}_{\lambda }^{\rm e}$ at any given pixel is then simply

$$\begin{equation} {A}_{\lambda }^{\rm e}=-2.5\log \left(\frac{F_{\lambda }}{F_{\lambda ,{\rm nd}}} \right), \end{equation} \tag{ 1 }$$

where F_λ is the flux of the model with dust and F_{λ, nd} is the flux of the model without dust. Note that we use "attenuation" here broadly to include the threefold effects of dust: absorption, scattering, and emission. The latter is of particular importance moving into the MIR.

As can be seen Figure 2 (and will be discussed again in Section 3.3), there is a very strong, narrow PAH feature which lies within the IRAC 3.6 μm bandpass. The 4.5 μm filter, on the other hand, has no discrete dust components, although it is somewhat affected by continuum dust emission. We find that dust emission contributes 12% and 9% of the total integrated 3.6 μm and 4.5 μm model flux, respectively. Local values, along sight-lines involving dust clumps, are likely higher: Meidt et al. (2012) find ranges of 5%–13% contamination at 3.6 μm in the integrated light of Spitzer Survey of Stellar Structure in Galaxies (S⁴G) targets but local amounts around 20% in star-forming regions. Therefore, based on the reduced contamination from dust emission at 4.5 μm as well as the fact that it provides a larger wavelength baseline for other attenuation effects (absorption and scattering), we expect that the 4.5 μm data will be more useful for estimating the NIR attenuation than the 3.6 μm data.

With these considerations in mind we inspected the relationship between color and K_s-band attenuation of the model (Figure 3). An ideal attenuation correction metric would have both a tight and steep relationship between attenuation and color—i.e., a small error in color leads to a small error in the estimated attenuation. The J−K_s and H−K_s corrections both have very small slopes at ${A}_{K_{\rm s}}^{\rm e}\gtrsim 0.75$ mag, making them poor choices for use as attenuation correctors. All NIR-IRAC colors (NIR−4.5 μm colors are shown in Figure 3, but the trends are very similar for NIR−3.6 μm colors) have similar slopes at ${A}_{K_{\rm s}}^{\rm e}\gtrsim 0.75$ of ∼1 mag in color per magnitude in attenuation, which indicates that either the attenuation in the 4.5 μm band is slowly varying or dust emission is, on average, offsetting the light losses due to attenuation at that wavelength.

Zoom In Zoom Out Reset image size

While the NIR-IRAC colors have similar responses to increasing attenuation, the K_s−4.5 μm (also K_s−3.6 μm) color clearly has the least scatter in the relationship. This is likely due to two factors: the K_s-band has less dust attenuation than either J or H, and that K_s is significantly redward of the Wein peak for even the coolest stars, which makes the K_s−IRAC color insensitive to the effects of different stellar populations. This feature is especially critical because it means that the assumptions we made about the parameters of the super-thin disk in Section 3.1 will not bias our results. Indeed, the attenuation corrections we produce use data at all positions in the model galaxy, and show a single trend applicable across the entire disk. For these reasons and those given concerning MIR dust contamination earlier in this section, we choose to use the K_s−4.5 μm color to correct the attenuation. We take the 1σ error on the attenuation correction in ${A}_{K_{\rm s}}^{\rm e}$ and propagate that with the errors on the data to compute the final errors for the corrected SB images.

To test the applicability of Model A we overlaid color–color profiles of the model on the data and checked for discrepancies. We compared the colors as a function of both R and z to isolate any trends with position. A sample set of K_s−3.6 μm versus J−K_s color–color distributions is shown in Figure 4; K_s−4.5 μm versus J−K_s distributions are shown in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At large distances from the midplane (high z values), both the model and the data trend toward K_s−3.6 μm and K_s−4.5 μm colors of zero and J−K_s colors of ∼1 mag. The K_s−IRAC color is a result of even the coolest stars being on the Rayleigh–Jeans tail in those spectral regions—the slopes of all stars are self-similar at these wavelengths, so Vega colors go to zero. We note that the large J−K_s color implies the presence of very cool stars—a J−K_s color of 1 mag is as red as the coldest dwarf in Bessell & Brett (1988). If the NIR light at these heights is dominated by giants, the J−K_s color is roughly equivalent to a M2III star or an O-rich asymptotic giant branch (AGB) with T_eff  ∼ 3500 K (Lançon & Mouhcine 2002).

It is also evident that at small height Model A fails to reproduce the colors of the data; model colors are increasingly too red in K_s−3.6 and K_s−4.5 when z ≲ 0.5 kpc. Properly modeling the colors at z < 0.5 kpc is crucial, as the midplane is where attenuation from dust is the strongest and any signal from a super-thin disk component would be most apparent. While there is better overlap when using the K_s−4.5 μm color (top panel of Figure 5), Model A still produces redder K_s−4.5 μm colors than the data.

The culprit for the discrepancy in the colors appears to be an excess of PAH line emission in the model, as seen in Figure 2: there is a weak emission line at ∼3.5 μm which falls within the IRAC 3.6 μm bandpass, while the wings of several lines contribute to the flux in the 4.5 μm band. There are (at least) two possible reasons for this discrepancy. The first is that there are fewer PAHs relative to other grain sizes near the midplane of NGC 891. Reducing the amount of PAHs decreases the total 3.6/4.5 μm flux, increasing the K_s−3.6/4.5 μm color without affecting the J−K_s color. Lowered PAH levels near the midplane could be the result of PAH grain destruction in high density areas, similar to what is seen in the MW (e.g., Povich et al. 2007). On the other hand there is also empirical evidence that, contrary to what is frequently assumed in dust models (including the model used here), the extinction law flattens in the MIR (Indebetouw et al. 2005; Zasowski et al. 2009). Flattening the extinction law for the IRAC bands would also have the effect of reducing the K_s−3.6/4.5 μm color at higher attenuations while leaving the NIR colors unchanged.

A full investigation of why a MW dust model has difficulty producing proper K_s−3.6/4.5 μm colors is beyond the scope of this work, given that we are merely concerned with computing a proper attenuation correction rather than conclusively determining the root cause of the attenuation. For simplicity, we choose to adjust the PAH levels in high density regions instead of changing the entire extinction law. Therefore we ran a second RT model with two prescriptions for the fraction of different sized dust grains. In this new model (called Model B) voxels with dust densities ⩽3.4 × 10⁻²⁷ g cm⁻³ are weighted 15% for ultra-small grains, 10% for small grains, and 75% for large grains, while high-density voxels use the following breakdown of grain sizes: 1% ultra-small grains, 35% small grains, and 64% large grains. All other quantities are held constant save the bolometric luminosity of the thin disk, which is decreased by ∼7%. These changes were made empirically by eye to better reproduce the infrared colors of the data while still accurately fitting the integrated SED. The SED of Model B is shown alongside that of Model A in Figure 2, and color–color K_s−3.6/4.5 μm versus J−K_s profiles are shown in the bottom panel of Figures 4 and 5.

Model B's SED has only 9% and 6% of the total 3.6 μm and 4.5 μm flux contributed by dust emission, generally better replicates the color distribution of NGC 891 than Model A. There are still some minor discrepancies between the data and the model in the MIR, however, which are shown in Figure 6. The model 3.6 μm–4.5 μm color is systematically bluer than the data. Near the midplane this is largely caused by the mismatch in the 3.6 μm flux as discussed earlier, but at large heights some of the discrepancy arises from the model being slightly too faint at 4.5 μm (visible in the bottom panels of both models in Figure 5). This mismatch is ≲0.08 mag, which implies that the y-intercepts of the attenuation correction formulae in Appendix C should be shifted by approximately that amount. The systematic error resulting from this mismatch will therefore be at the ≲8% level, which should only impact flux measurements and (as will be seen later) does not have a significant effect on any of our conclusions. While the 3.6 μm model flux does match the data at large heights, due to the aforementioned significant mismatch in dust-obscured regions we still consider the 4.5 μm flux as more reliable for this purpose, systematic error notwithstanding. We will investigate ways to improve our input SEDs for future publications, but for now we will focus on demonstrating the rest of our methodology on NGC 891.

Zoom In Zoom Out Reset image size

Functionally, adjusting the dust distribution decreases the reddening in the K_s−IRAC μm color while leaving J−K_s essentially unchanged. In addition to the systematic error in intrinsic MIR color, to have some leverage on our systematic uncertainty caused by our imperfect knowledge of the dust distribution we generally perform subsequent analyses with both models, using only Model B when we do not expect to see different responses to analysis from the two models and the savings in computation time is significant.

The exact behavior of the attenuation as a function of color will be different in each bandpass, so we fit ${A}_{\lambda }^{\rm e}$ separately for J, H, and K_s. Fitting the attenuation from the RT models is non-trivial, largely because at low ${A}_{\lambda }^{\rm e}$ the slope of the attenuation curve changes significantly while at large attenuations the function becomes linear. This problem is compounded by the dearth of model data at very high attenuations. These factors lead to poor fits of single polynomial functions, especially at large ${A}_{\lambda }^{\rm e}$. To deal with this issue we employ piecewise fits of the attenuation curve. Details of our method and the formulae for the resulting fits are shown in Appendix C.

The attenuation fits are plotted against the model data in Figure 7. In all cases it is apparent that the model trends are well fit. In Figure 7 we also compare the fits between Models A and B. At low attenuations Model B has much more complex behavior, although for all three bands the two models are within ∼0.1 mag of each other in ${A}_{\lambda }^{\rm e}$. At large attenuations the attenuation relations have very similar slopes but are offset from each other by about 0.3 mag in J, H, and K_s. Model B has more extinction at a given K_s−4.5 μm color than Model A, for K_s−4.5 μm ≳ 0.5 mag. The difference in behavior between the models at large values of K_s−4.5 μm arises because the depletion of PAH grains in Model B acts to lower the amount of line/continuum emission from dust at ∼5–6 μm, the tail of which is blue enough to affect the 4.5 μm IRAC band. The attenuation in the K_s-band, however, is not strongly affected by changes in PAH densities and therefore remains essentially constant. The result is that Model B has bluer K_s−4.5 μm colors for a given ${A}_{\lambda }^{\rm e}$, seen in Figure 7. This also indicates that Model A overpredicts the 4.5 μm SB near the midplane by ∼0.3 mag arcsec⁻².

Zoom In Zoom Out Reset image size

Assuming a MW extinction law and a foreground screen of dust, Aoki et al. (1991) correct their K-band emission by the formula

$$\begin{equation} K_{c}=K_{r} - \frac{(H-K)-0.25}{0.56}, \end{equation} \tag{ 2 }$$

where K_r is the uncorrected K-band SB and K_c is the corrected K-band SB. As the reddest filter and therefore the data least likely to be affected by errors in the attenuation correction, we follow Aoki et al. (1991) and focus on our K_s-band data. We show the uncorrected K_s data along with corrected K_s profiles in Figure 8 on both sides of NGC 891 at 3, 5, 7, and 9 kpc in radius. We plot the vertical profiles at these positions because they show data at high signal-to-noise ratio (S/N) while avoiding bulge contamination, and are representative of the profiles at all radii. While the K and K_s-bands are not identical, we use the Aoki et al. (1991) attenuation correction on our data to obtain a schematic understanding of the differences between it and our RT methodology.

Zoom In Zoom Out Reset image size

We find that generally both dust models in this work have larger peak attenuation than the Aoki et al. (1991) correction. Model A is on average brighter by ∼0.16 mag, while Model B is brighter by ∼0.33 mag. This discrepancy is likely due to the saturation of the H−K_s color at large optical depths (Matthews & Wood 2001 and see Figure 3), which we are less sensitive to since we are using K_s−4.5 μm color and models that accurately represent this effect. Additionally, while the attenuation corrected vertical profiles of both Model A and Model B peak at a relatively consistent location at all radii, the peak of the Aoki et al. (1991) dust corrected data has much more variation (especially visible at 5 and 9 kpc). Our two RT models are very similar to each other, but (as predicted by Figure 7) the model with a PAH-depleted grain size distribution in high-density regions (Model B) leads to larger inferred ${A}_{\lambda }^{\rm e}$ for a given K_s−4.5 μm color than in Model A, resulting in a brighter corrected $\mu _{{K_{\rm s}}}$.

We show NGC 891's JHK_s radial SB profile along the major axis (i.e., z = 0) uncorrected for dust attenuation in Figure 10. The vertical bin size is 09 or ∼40 pc. The radial profiles with dust are very similar to those shown in Figure 9 of Aoki et al. (1991), except that our data shows a distinct bump in all bands at ∼3 kpc on the SW side while they only see that feature in H and K. The reason for this discrepancy is unclear, but it indicates that this 3 kpc feature is a stellar overdensity, not a dusty spiral arm as hypothesized by Aoki et al. (1991). We also present uncorrected J−K_s and H−K_s major axis color profiles in the left panel of Figure 11. The J−K_s profile seems to show a weak trend toward bluer colors with increasing radius, but there is significant scatter. The H−K_s profile shows a similarly weak trend. No color trend is visible as a function of radius when using attenuation corrected profiles (right panel of Figure 11). There is also an asymmetry between the NE and SW colors of NGC 891 around R ∼ 3 kpc, in the same position with the apparent stellar overdensity discussed previously in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our major axis radial profiles also show the northeast side of NGC 891, which was not shown in Aoki et al. (1991) due to noise issues in the H band on their detector. Comparing the radial profiles on both sides of the bulge shows the significant amount of stochastic structure in the midplane. In addition to the bumps noted by Aoki et al. (1991) the southwest side has another light enhancement at ∼7.5 kpc and a marked dearth of light at ∼5.5 kpc. The northwest side appears to have fewer of these strong features, however, although it certainly has its fair share of smaller bumps and wiggles. All of these features point to the strong observable effect of clumpy dust and starlight, even at these long wavelengths.

We compare the K_s-band major axis radial profile from Figure 10 to its dust-corrected counterpart (using Model B) in the left panel of Figure 12. The dust-corrected profile is uniformly brighter than the fiducial profile, which is to be expected. Several local SB peaks seen in the observed profile are not present in the dust-corrected profile (at R ∼ 3 kpc, for instance, on the SW side) while others have been enhanced. The apparent star-forming regions mentioned in the previous section all correspond to local maxima in the dust-corrected profile; of particular note is point "B," which does not stand out on the original profile. Also interesting is the "trough" on the SW side of the dust-corrected profile between 4 and 6.5 kpc. While the original profile has a local minima in the same region, it is much smaller in radial width. It seems likely that striking appearance of this feature is due to a combination of a real intensity drop-off at R ∼ 4 kpc and the enhancement of the underlying smooth profile at ∼7 kpc (point "C").

Zoom In Zoom Out Reset image size

Just from looking at the attenuation-corrected radial SB and color profiles (Figures 11 and 12) it is difficult to discern whether the asymmetries in the flux and color profiles are due to stellar density enhancements, intrinsic color variations, or both. Figure 13 shows the magnitude and color differentials between the SW and NE sides of the major axis profiles, where ΔX is defined as the magnitude difference between corresponding ∼40 pc² cells at the same radii but on opposite sides of NGC 891's midplane. Figure 13 clearly indicates that asymmetries in SB are correlated with color—brighter regions of NGC 891's midplane are also bluer. The mean attenuation-corrected J−K_s and H−K_s colors for ΔJ and ΔH = 0 are ∼0.77 and ∼0.12 mag, respectively, roughly the color of a K2III star. There are minor (∼0.03 mag) differences in these colors between the two sides, but are well within the dispersions of the data.

Zoom In Zoom Out Reset image size

To approximate the effect of star formation on these Δ metrics we also compute how the colors would change if the excess midplane light between the two sides of NGC 891 is due to a hot star. We find that the ΔJ−K_s index is well fit by an O-type star, while the trend for ΔH−K_s versus ΔH is underpredicted for any stellar type. A full investigation of this minor discrepancy is beyond the scope of this work, but the general trend is clear. Rather than merely being overdensities of old stars, asymmetries in the radial major-axis profile of NGC 891 are being caused by the addition of young stellar populations to the underlying stellar density. In some cases these asymmetries are likely due to individual regions of enhanced star-formation, but the broad nature of several features are at much larger scales than would be expected for localized star-forming regions. For instance, the SB is relatively constant on the NE side of NGC 891 between 3 and 6 kpc, while the aforementioned SW "trough" is ∼2 kpc wide. The simplest explanation for these features is that we are viewing projected spiral structure, which NGC 891 is believed to possess (Kamphuis et al. 2007).

To verify the validity of our dust correction we compare IRAC 3.6 and 4.5 μm major axis profiles to the dust-corrected K_s profile in the middle and right panels of Figure 12. The IRAC data in general shows the same non-axisymmetric features as the dust-corrected K_s-band profile, with a few local exceptions (e.g., at ∼3 and ∼5.5 kpc on the SW side). The IRAC profiles are also globally fainter by up to 0.5 mag arcsec⁻², whereas for starlight the unattenuated K_s−3.6/4.5 μm color should be zero (in our Vega photometric system) due to all filters being in the Rayleigh–Jeans regime. These discrepancies are entirely due to dust attenuation in the 3.6 and 4.5 μm bands: we computed dust corrections for these filters in the same way as for the NIR, which when applied to the IRAC major axis profiles produces profiles virtually identical to the K_s-band data (Figure 14).

Zoom In Zoom Out Reset image size

A series of K_s-band vertical profiles as several radii are shown in Figure 15. Gray points show the profiles including the effects of dust, while black points have been attenuation corrected. The midplane attenuation is roughly constant at ∼1.5 mag arcsec⁻² out to R = 8.5 kpc. The amount of attenuation quickly decreases with z, however, and by z ∼ 0.5 kpc the effects of dust on the K_s-band vertical profiles are almost negligible. At R ∼ 9 kpc there is a truncation in the central SB of the dust corrected profiles but not of the uncorrected ones. This scenario can only result from a cutoff in both the light near the midplane as well as the dust attenuation. This feature illustrates the importance of properly accounting for dust attenuation, as the uncorrected profile appears to be smoothly decaying across this threshold.

Zoom In Zoom Out Reset image size

4.2.1. Model Profiles

The intensity distribution of an edge-on spiral with a radially exponential disk luminosity density distribution at all heights can be described, in the absence of attenuation, by the analytic function

$$\begin{equation} I(R,z) = I_{0}\frac{|R|}{h_{R}}K_{1}\left(\frac{|R|}{h_{R}}\right) f(z), \end{equation} \tag{ 3 }$$

where I₀ is the central flux, K₁ is the modified Bessel function of the first order, and f(z) is a function governing the vertical light distribution (van der Kruit & Searle 1981a).

A simple analytic description of the vertical density profile of a spiral galaxy has been proposed by van der Kruit (1988) to be of the form:

$$\begin{equation} f(z) = 2^{-2/n}\mathrm{sech}^{2/n}\left(n|z|/(2h_{z})\right), \end{equation} \tag{ 4 }$$

where n varies between 1 and ∞. The limiting behavior of this formula produces a sech² distribution

$$\begin{equation} f(z) = \frac{\mathrm{sech}^{2}(|z|/(2h_{z}))}{4} \end{equation} \tag{ 5 }$$

when n = 1 and the exponential distribution

$$\begin{equation} f(z) = \mathrm{exp}(-|z|/h_z) \end{equation} \tag{ 6 }$$

when n = ∞. For convenience we include the factor of four in the sech² profile with I₀ and use

$$\begin{equation} f(z) = \mathrm{sech}^{2}(|z|/(2h_{z})) \end{equation} \tag{ 7 }$$

instead.

The first attempts to fit the detailed vertical SB profiles of edge-on spirals began by treating the disk as a locally isothermal self-gravitating sheet (van der Kruit & Searle 1981a), which results naturally in the sech² model (Spitzer 1942). As a result this model is generally referred to as the "isothermal" distribution. Note, however, that this strict physical connection is lost when considering multiple disks, rotating disks with radial gradients, and/or disks embedded in a dark matter halo. Nonetheless, the sech² profile is not a bad descriptor of many observed disk vertical light profiles. More recent work has found that a steeper profile appears to fit some data better—usually an exponential profile (Wainscoat et al. 1989; Aoki et al. 1991; Xilouris et al. 1999), although sometimes the intermediate n = 2 "sech" profile is employed (Barnaby & Thronson 1992; Rice et al. 1996). A theoretical framework which produces vertical SB distribution much closer to that of an exponential from physical conditions in spiral galaxies (by considering the gas in addition to the stars) has been developed partly in direct response to these experimental findings (Banerjee & Jog 2007). Comerón et al. (2011b) also find that the vertical density profiles for dynamically consistent models of multiple self-gravitating disks are steeper than a sech² function. Discussed further in Section 5.2, their models are well characterized by the intermediate sech vertical profile.

In general, however, there is currently no concrete consensus for which type of profile (or superposition of profiles) better represents the underlying distribution of stars, due largely to the fact that both the exponential and sech² profiles have significant deviations from each other only at small heights, exactly the place where dust obscuration is most prevalent. In this work, since we are able to accurately remove the deleterious effects of dust on the data, we will use both exponential and sech² profiles in our vertical fits and see if we are able to distinguish one as clearly better in a statistical sense. Additionally, since NGC 891's thin and thick disks are well known (e.g., Morrison et al. 1997; Ibata et al. 2009) and there is a clear super-exponentiality in the dust corrected profiles near the midplane in Figures 8, 9, and 15, we expect our fits to require three separate vertical components in order to constrain the data at all z.

4.2.2. Seeing Convolution

Normally, observational studies of edge-on spirals have not required a detailed treatment of the effects of seeing on the observed profiles. This is a result of the combination of large pixel-sizes on detectors and the censoring of regions near the midplane to avoid dust effects (e.g., de Grijs & van der Kruit 1996; de Grijs et al. 1997). However, in this work the seeing correction is of significant importance, especially in distinguishing between sech² and exponential fits, where the seeing serves to smooth out the exponential profile and make it appear more like a sech².

We make the simplifying assumption that the seeing profile can be approximated by a 2D Gaussian point-spread-function, as in Pritchet & Kline (1981). Then, the convolution is

$$\begin{eqnarray} I_{{\rm conv}}(R,z) &=& \frac{1}{2\pi \sigma ^2}\int _{-\infty }^{\infty }I(r^{\prime },z^{\prime })\, \mathrm{exp}\left(\frac{-(R-r^{\prime })^{2}}{2\sigma ^{2}}\right) \, \nonumber\\ && \times \mathrm{exp}\left(\frac{-(z-z^{\prime })^{2}}{2\sigma ^{2}}\right)\,dr^{\prime }dz^{\prime }, \end{eqnarray} \tag{ 8 }$$

where σ is the FWHM of the seeing. Since all the profiles we are considering here are of the form of Equation (3) we can perform the convolution separately in R and z:

$$\begin{eqnarray} I_{{\rm conv}}(R,z) = &\frac{I_{0}}{2\pi \sigma ^2}\int _{-\infty }^{\infty }\frac{|r^{\prime }|}{h_{R}}K_{1}\left(\frac{|r^{\prime }|}{h_{R}}\right)\,\mathrm{exp}\left(\! \frac{-(R-r^{\prime })^{2}}{2\sigma ^{2}} \!\right)\,dr^{\prime } \nonumber \\ & \times \int _{-\infty }^{\infty }f(z)\,\mathrm{exp}\left(\frac{-(z-z^{\prime })^{2}}{2\sigma ^{2}}\right)\,dz^{\prime }. \end{eqnarray} \tag{ 9 }$$

For both exponential and sech² profiles we used the "convolve" function in the NumPy Python library. For every value to be convolved the SB profile was sampled with a spacing of 0.1 times the seeing (FWHM = 09) over a range of values equal to six times the seeing in either direction. The Gaussian kernel was sampled at the same spacing between −3 and +3 times the seeing. This scheme ensured that the convolution did not miss any significant amounts of flux due to edge effects from either the convolution kernel or the function, while minimizing the computation time required to perform the integral.

When the relevant scale-length or -height is much larger than the seeing, the effect of the convolution on the original profile is negligible. This is especially true at large distances from the center of the profile. Since the radial scale-length is much larger than the seeing we use the unconvolved radial intensity profile in our fitting analyses to minimize computation time.

4.2.3. 2D Fitting

With the SB profiles corrected for the attenuation and convolutions for seeing in place, we now fit the intrinsic distribution of starlight. While one-dimensional (1D) fits have been successfully employed to compare single-component SB models to data, Morrison et al. (1997) make a compelling argument, based on degeneracies between free parameters, for the necessity of fully 2D fits when multiple stellar distributions overlap. We therefore opt to model our data exclusively with image-based fits. We choose to model the K_s-band, as it is where both the dust correction and stellar population-mismatch errors will be smallest (due to the decreasing attenuation with increasing wavelength and all stellar fluxes being on the Rayleigh–Jeans tail).

Before performing our fits, however, we first mask the image based on several criteria. All pixels with flux from foreground or background contaminants are masked, as well as pixels with errors >1 mag (while large-error pixels can still contribute to a χ² fit in aggregate, gains from leaving these pixels in will be minimal and removing them allows for computational speedup in the fitting process). Contaminant masking was done by hand and whenever possible we choose to be conservative. We also mask out all pixels with R < 3 kpc to remove any pixels dominated by "bulge" light. We refer to the bulge light in quotations because, as we shall see in Section 4.2.5, the central light distribution is complex and likely non-axisymmetric. For exactly these reasons we wish to focus first on the outer (in radius) disks.

We use the Python lmfit⁷ extension to the Levenberg–Marquardt nonlinear least-squares minimization routine in scipy.optimize.leastsq (which is itself based on the MINPACK-1 library). Although we report SB values in mag arcsec⁻², all fitting is performed on flux values to preserve error estimates. To avoid false minima, we ran the fit 100 times for one- and two-component models and 200 times for three-component models. The initial guesses of all free parameters for each run were chosen randomly within upper and lower limits based on what we considered to be reasonable ranges for spiral galaxies (Table 3). The model with the lowest χ² was selected and used as the input guess to 100 bootstrapping iterations, finding, as Morrison et al. (1994) did, that standard error estimates from the nonlinear fit were often too small. Note that only random uncertainties are shown in Table 3; the ≲8% systematic error on the central SBs from input stellar SED mismatch (Section 3.3) is not shown.

Table 3. Two-dimensional Fits

Parametera	Allowed Range		One Disk				Two Disks				Three Disks
	Min	Max	Exp	Δb	sech2	Δb	exp	Δb	sech2	Δb	Exp	Δb	sech2	Δb
μ0, 1	25	10	14.00 ± 0.018	0.04	14.52 ± 0.020	0.13	13.26 ± 0.040	−0.58	13.53 ± 0.026	−0.66	13.08 ± 0.040	−0.94	13.41 ± 0.023	−0.67
hR, 1	0.1	10	3.55 ± 0.027	0.12	3.67 ± 0.025	0.11	2.42 ± 0.045	−0.89	2.54 ± 0.032	−0.70	1.96 ± 0.038	0.13	2.32 ± 0.022	−0.37
hz, 1	0.01	2.5	0.28 ± 0.002	0.00	0.19 ± 0.002	0.01	0.12 ± 0.004	−0.10	0.07 ± 0.001	−0.02	0.08 ± 0.004	0.02	0.06 ± 0.001	0.00
Ltot, 1c	⋅⋅⋅	⋅⋅⋅	1.37 × 1011(100)	⋅⋅⋅	1.19 × 1011(100)	⋅⋅⋅	7.92 × 1010(47)	⋅⋅⋅	7.57 × 1010(48)	⋅⋅⋅	5.05 × 1010(29)	⋅⋅⋅	6.62 × 1010(38)	⋅⋅⋅
μ0, 2	25	10	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	15.29 ± 0.101	−1.91	15.68 ± 0.067	0.07	14.52 ± 0.096	0.47	15.32 ± 0.067	0.35
hR, 2	0.1	10	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	4.44 ± 0.241	0.63	4.30 ± 0.178	0.89	3.87 ± 0.118	0.46	4.11 ± 0.114	0.51
hz, 2	0.01	2.5	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	0.48 ± 0.016	−0.77	0.33 ± 0.007	−0.02	0.29 ± 0.010	0.04	0.25 ± 0.006	0.03
Ltot, 2c	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	8.97 × 1010(53)	⋅⋅⋅	8.33 × 1010(52)	⋅⋅⋅	9.60 × 1010(55)	⋅⋅⋅	8.41 × 1010(48)	⋅⋅⋅
μ0, 3	25	10	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	17.81 ± 0.015	0.04	18.70 ± 0.015	0.05
hR, 3	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	4.80d	⋅⋅⋅	4.80d	⋅⋅⋅
hz, 3	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	1.44d	⋅⋅⋅	1.44d	⋅⋅⋅
Ltot, 3c	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	2.85 × 1010(16)	⋅⋅⋅	2.52 × 1010(14)	⋅⋅⋅
Δz	−0.1	0.1	0.05 ± 0.001	0.00	0.05 ± 0.001	0.00	0.05 ± 0.001	0.00	0.05 ± 0.001	0.00	0.05 ± 0.001	0.00	0.05 ± 0.001	0.00
Reduced χ2	⋅⋅⋅	⋅⋅⋅	7.32	1.09	11.7	1.70	4.21	−0.19	4.61	−0.58	3.77	−0.47	3.65	−0.55

Notes. ^aμ values in mag arcsec⁻¹, h_R, h_z, and Δz values in kpc. ^bReported values in "exp" and "sech2" columns are for Model B, and Δ = Parameter_Model B − Parameter_Model A. ^cTotal luminosities are given in solar units where L_{☉, K} = 3.909 × 10¹¹ W Hz⁻¹ (from Binney & Merrifield 1998). Values in parenthesis are the percentage of the total disk luminosity in that component. ^dFixed values from Ibata et al. (2009).

Download table as:  ASCIITypeset image

We fit models with one-, two-, and three-components for both attenuation corrections, motivated by the presence of three disk components in the MW: the young, star-forming disk, the thin disk, and the thick disk (van Dokkum et al. 1994; although Bovy et al. 2012 indicate that perhaps the boundaries between the disks are not so well-defined). It is important to note that we do not force the physical parameters of the super-thin or thin disks to resemble those seen in NGC 891 and/or the MW by prior studies by way of limiting our parameter space (although we do restrict the parameter space for the thick disk; see next paragraph for details). On the contrary our externally imposed constraints are extremely permissive and allow for fitted values that can either agree or disagree with the literature.

For models with one or two components each component had three free parameters: μ₀, h_R, and h_z. Additionally, a single vertical offset for each model was fit simultaneously, to remove the effects of any small errors from our choice of galaxy center. Because any such errors in R will be much smaller than h_R (and therefore have minimal effect on the fits) to avoid additional computational overhead we do not allow a radial offset in the fits. Therefore, for a model with 1 (2) vertical component(s), there are 4 (7) free parameters. Preliminary three-component fits were unable to consistently select physically plausible fits due to the low S/N at large heights; however, visual inspection of the two-component residuals showed an excess of light at large heights above the midplane, especially when fitting with sech² profiles. This arises from the shallower inner profile of the sech² function which drives the fits (in a χ² sense) to a smaller h_z. To reduce our parameter space (and therefore the degeneracies) we restrict the scale-length and -height of the third (most extended) component to values given by Ibata et al. (2009), who used HST to obtain very precise scale parameters with resolved star-counts of NGC 891's thick disk. We do, however, leave the flux normalization for this component a free parameter, to allow for band-pass differences and an incomplete understanding of the integrated SED from the star-count analysis. Therefore our three-component fits involve eight free parameters. The results from all fits are shown in Table 3.

Generally there is very good agreement between the fits from the two attenuation correction models. The one-component fits are nearly identical, especially for the exponential case. The main distinction of note for the two- and three-component models is that the central SBs of the innermost (super-thin) disk fits are smaller for Model A than for Model B, which is to be expected given the increased amount of attenuation correction for a given color in Model B (see Figures 7 and 8). Since none of the differences in fitted parameters between Model A and B would substantially affect our initial analysis (although the differences between the models are more relevant for our discussion of the total disk luminosities; see Section 5.1.2), for simplicity we will focus the discussion on the fits from Model B because it appears to provide a better fit to the colors of the data (cf. Figures 4 and 5 and Section 3.3).

Several vertical slices of the data, overlaid with the relevant 1D portions of the fits, are shown in Figures 16–18. The single component fits (Figure 16) are clearly too simple, failing to match the peak intensity at the midplane as well as at high latitude. Figure 17 illustrates the dramatic improvement between the one- and two-component fits. Whereas the single component fit had difficulty fitting the data at any height, the two-component sech² fit does a reasonable job within ±0.5 kpc, while the exponentials appear to accurately fit almost the entire profile. The difference in reduced χ² between using sech² functions versus exponentials is relatively small compared to the statistical improvement in adding the second component to the fits, however, since the regions of largest discrepancy in the two-component fits are where the noise is largest. This discrepancy arises because a single exponential (even with the seeing correction) is more peaked than a sech² and is therefore better able to fit the (most statistically relevant) midplane regions, leaving the second exponential to fit the outer regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While adding a third component decreases the reduced χ² of the exponential fits only by a small amount, the improvement in the sech² fits is more significant, both in the χ² sense and upon visual inspection of Figure 18. The three-component exponential and sech² fits are virtually indistinguishable from each other, both in a visual and statistical sense, and both are clearly able to fit the vertical profile at all heights. A representative plot of vertical SB fits with the contributions of the individual disks is shown in Figure 19. First, the visual improvement resulting from adding the third component is clearly noticeable, especially at heights above ∼0.75 kpc. Additionally, while the total SB profiles of the three-component exponential and sech² fits are very similar, there are noticeable differences in the form of the individual components. These differences are largest in the relative μ₀ between the exponentials and the sech²s, mostly to compensate for the flatter central profile of the sech² profiles.

Zoom In Zoom Out Reset image size

Even with the three-component models, we find reduced χ² > 1. This could be due to an underestimation of our errors, most likely in the attenuation correction. However, given the relative simplicity of our model compared to the complex nature of a galactic system (e.g., the model is perfectly axisymmetric), the fact that we do not recover a reduced χ² ∼ 1 is not surprising. Our 2D fits do not take into account the potential effect of warps, star-forming regions, spiral arms, or any other non-axisymmetric structure. Despite these limitations, however, the reduced χ² is still effective at distinguishing between the relative goodness-of-fit of models, as shown by the improvement from one-component fits to three-component fits both in reduced χ² and upon visual inspection.

4.2.4. Disk Asymmetries

To investigate the impact of non-axisymmetric structures on our model fits to NGC 891, we constructed model − (dust-corrected) data residual SB images for both three-component models (top two panels of Figure 20). These images have been masked (bad values set to zero) in the same way as for the 2D fitting (except for the inner 3 kpc, the boundaries of which are marked by dashed lines but not masked out), and smoothed by a σ = 27 Gaussian to bring out larger-scale detail. Additionally, the K_s-band bulge determined by Xilouris et al. (1999) has been added to the model image before producing the residual. White contours show the boundaries of regions with errors <0.2 and <0.5 mag pixel⁻¹ (unsmoothed)—while we fit all points with errors <1 mag pixel⁻¹, we weight each pixel by its error and as a result the fit is most sensitive to pixels with small uncertainties. Therefore we focus our analysis of the quality of the fit on the regions of the image with smaller uncertainties.

Zoom In Zoom Out Reset image size

Several asymmetries are clearly present in Figure 20. Based on a visual inspection of the data before resampling, at least three of them (labeled in Figure 20 and also shown in Figures 1 and 9) are likely to be regions of enhanced local star formation in NGC 891. Oosterloo et al. (2007) find a mild H i warp in NGC 891 (position angle = 15). There are some slight indications of this warp in the data, chiefly in the excess of light over the model above the midplane at radii larger than region "C" on the SW side. On extremely close inspection of the non-dust-corrected mosaics (see Figure 1, although it is difficult to distinguish in that image stretch) it is also possible to discern a similar offset in the dust, but due to the clumpiness of the dust and the faintness of this feature we are unable to make a stronger claim about whether or not it is real.

4.2.5. Fitting the Central Region of NGC 891

Beyond the specific asymmetries called out above, it is clear that these models do a poor job of fitting NGC 891's K_s-band light within 3 kpc. In addition to the large amount of light above and below the midplane from the galaxy's bar that is not reproduced in our model there is also clearly excess light in NGC 891's nucleus—even with the addition of the Xilouris et al. (1999) bulge to the model. In contrast, just outside of the nucleus the model greatly overpredicts the flux near the midplane. The simplest explanation for this feature is that at R ∼ 3 kpc there is a transition between disk and bar dominated regions, such that the three-component disk is truncated within this radius and the bar does not extend beyond this radius. There is ample evidence for such truncations in other galaxies (e.g., Anderson et al. 2004). To the extent that bars arise from dynamical instabilities in the inner disk, then this is not a truncation per se, but a redistribution of the disk from an axisymmetric, dynamically cold stellar mass distribution to a dynamically hotter distribution with more structure.

This transition is more apparent as a function of radius in vertical profiles (Figure 21). While the three-component model well fits the data outside of 3 kpc, the data clearly has a dearth of flux inside of 3 kpc right up until the very central regions of the galaxy. At small radii, however, there is a clear excess of flux at larger heights. We suggest this excess is due to a bar because previous studies have found evidence for such a feature in NGC 891 (e.g., Garcia-Burillo & Guelin 1995). Further, there are clearly kinks in the vertical light profile (most evident at R = 1.5 and z = 0.8 kpc) caused by the "X" shape of the projected bar (Combes & Sanders 1981; Mihos et al. 1995; Bureau & Freeman 1999; Bureau et al. 2006). Other than these small kinks, however, the bar appears to be very close to having a pure exponential vertical light distribution.

Zoom In Zoom Out Reset image size

Understanding the galaxy's behavior near the center is critical to estimating the total luminosity of each disk component and to fully characterize their radial distribution. Therefore we perform a second set of 2D fits on the data attenuation-corrected by Model B. While using the same basic scheme as in Section 4.2.3, this time we leave the central region of NGC 891 unblanked. This necessitates several new model parameters. We approximate the bar in edge-on projection as an outer-truncated disk. This appears to be a good approximation to the light distribution of the bar seen edge-on, and the functional form of a truncated disk simplifies our fitting procedure. Our bar model therefore has three free parameters: a central SB (μ_{0, bar}), scale-length (h_{R, bar}), and scale-height (h_{z, bar}).

Since there is an excess of light at small R and z that appears morphologically disconnected from the bar emission (see the top two panels of Figure 20), it is necessary to introduce an additional component. We first attempted to fit a classical R^1/4-law bulge to NGC 891, but found that this form made for extremely poor fits. We found instead that the central component of NGC 891 was well-fit by a nuclear disk with the same scale-height as the super-thin disk. We therefore fixed the scale-height of this component h_{z, nuc} = h_{z, ST} and allowed the central SB μ_{0, nuc} and scale-length h_{R, nuc} to vary.

Finally, we add a truncation radius R_trunc, which functions as an inner cutoff for all three main disks and an outer cutoff for the bar and nuclear disk. For simplicity this truncation is abrupt: at all points inside (outside) of R_trunc the main disks (bar and nuclear disk) have zero flux. In reality there will be a transition length instead of a discontinuous truncation, but we do not find that the data warrants this complication in the model. The scale-heights of the three main disks should be well determined by our original fitting scheme, and tests indicated that there were minimal changes in central SBs and scale-lengths of the thin and thick disks when the central region of NGC 891 was uncensored. Therefore we fixed the parameters of the thin and thick components at their best-fitting values in Table 3, but we allowed μ_{0, ST} and h_{R, ST} to vary. We perform these new fits only for sech² disks.

The parameters of the best-fitting model are listed in Table 4 and compared to the data and untruncated three-component sech² disk model (from Table 3) in Figure 21. The residual image from this model is shown in the bottom panel of Figure 20. The outstanding residuals that remain are associated with the "X" shape pattern attributed to the bar. It would be very useful to have a simple analytic parameterization of this light distribution, but we are not aware of one in the literature. While still not perfect, adding a disk truncation, a simple representation for the bar, and a nuclear disk dramatically improves our ability to fit the inner portion of NGC 891. The central SB and scale-length of the super-thin disk change minimally from the original fits. The bar has a scale-height ∼1.5 times that of the thin disk, but substantially smaller than the thick disk, which agrees with a visual inspection of our images. The reduced χ² of this fit is somewhat larger than for the fits with only the three disk components; this is due to the unmodeled "X" shape bar pattern.

Table 4. Two-dimensional sech² Fits, with Bar and Nuclear Disk

Parametera	Allowed Range		sech2
	Min	Max
μ0, ST	25	10	13.39 ± 0.02
hR, ST	0.1	10	2.30 ± 0.023
hz, STb	⋅⋅⋅	⋅⋅⋅	0.06
Ltot, STc	⋅⋅⋅	⋅⋅⋅	4.04 × 1010(25)
μ0, Tb	⋅⋅⋅	⋅⋅⋅	15.32
hR, Tb	⋅⋅⋅	⋅⋅⋅	4.11
hz, Tb	⋅⋅⋅	⋅⋅⋅	0.25
Ltot, Tc	⋅⋅⋅	⋅⋅⋅	6.86 × 1010(42)
μ0, Thb	⋅⋅⋅	⋅⋅⋅	18.71
hR, Thb	⋅⋅⋅	⋅⋅⋅	4.80
hz, Thb	⋅⋅⋅	⋅⋅⋅	1.44
Ltot, Thc	⋅⋅⋅	⋅⋅⋅	2.13 × 1010(13)
μ0, bar	25	10	15.42 ± 0.01
hR, bar	0.1	10	1.28 ± 0.009
hz, bar	0.01	2.5	0.38 ± 0.001
Ltot, barc	⋅⋅⋅	⋅⋅⋅	2.49 × 1010(15)
μ0, nuc	25	10	13.23 ± 0.04
hR, nuc	0.1	10	0.25 ± 0.007
hz, nucb	⋅⋅⋅	⋅⋅⋅	0.06
Ltot, nucc	⋅⋅⋅	⋅⋅⋅	8.32 × 109(5)
Rtrunc	1.5	4.5	3.09 ± 0.014
Δzb	⋅⋅⋅	⋅⋅⋅	0.05
Reduced χ2	⋅⋅⋅	⋅⋅⋅	4.55

Notes. ^aμ values in mag arcsec⁻¹, h_R, h_z, Δz, and R_trunc values in kpc. ^bValues fixed using best-fitting values for the three-component sech² model. ^cTotal luminosities are in units of L_{☉, K} = 3.909 × 10¹¹ W Hz⁻¹ (from Binney & Merrifield 1998). Values in parenthesis are the percentage of the total disk luminosity in that component.

Download table as:  ASCIITypeset image

5.1.1. Radial and Vertical Disk Scale-lengths

5.1.2. Disk Luminosities

For our single component fits (as well as for the thin disk component of our three-component exponential fits) we compute comparable scale-lengths and scale-heights but ∼1 mag arcsec⁻² brighter central SBs than found by Xilouris et al. (1999). This result implies that our single component models predict ∼50% higher total K_s-band luminosities for NGC 891. Our images do go deeper than those of Xilouris et al. (1999), which would boost the luminosity of our model. Our attenuation correction, which takes into account the effects of clumpy dust, may also contribute to this discrepancy. We note that our total luminosities are in much better agreement with more recent works (e.g., Popescu et al. 2011), which will be discussed in more detail later in this section.

For exponential and sech² disks the total flux can be computed by integrating over R and z. Including the effect of an inner disk truncation gives

$$\begin{equation} F_{{\rm tot,exp}} = 2\pi h_{z}f_{\nu ,0}10^{-\mu _{0}/2.5}e^{-R_{{\rm trunc}}/h_{R}}(h_{R}+R_{{\rm trunc}}) \end{equation} \tag{ 10 }$$

$$\begin{equation} F_{{\rm tot,sech}^{2}} = 4\pi h_{z}f_{\nu ,0}10^{-\mu _{0}/2.5}e^{-R_{{\rm trunc}}/h_{R}}(h_{R}+R_{{\rm trunc}}) \end{equation} \tag{ 11 }$$

for the super-thin, thin, and thick disks, where R_trunc is the truncation radius and f_{ν, 0} is the flux zero-point. In this work we use the Two Micron All Sky Survey (2MASS) value of 667 Jy for the K_s-band zero-point (Cohen et al. 2003). The equations for the total flux of the bar and nuclear disk are similar:

$$\begin{equation} F_{{\rm tot,exp}} = 2\pi h_{z}f_{\nu ,0}10^{-\mu _{0}/2.5}(h_{R}-e^{-R_{{\rm trunc}}/h_{R}}(h_{R}+R_{{\rm trunc}})) \end{equation} \tag{ 12 }$$

$$\begin{equation} F_{{\rm tot,sech}^{2}} = 4\pi h_{z}f_{\nu ,0}10^{-\mu _{0}/2.5}(h_{R}-e^{-R_{{\rm trunc}}/h_{R}}(h_{R}+R_{{\rm trunc}})). \end{equation} \tag{ 13 }$$

We show the total luminosities of each component for both three-component models as well as the model which adds a bar and nuclear disk (the "truncated" model) in Tables 3 and 4.

For the truncated model we find that the bar contributes roughly the same total luminosity as the thick disk, despite being confined to the central regions of NGC 891. The nuclear disk is by far the faintest, with a total luminosity only 40% of the second least-luminous component (the thick disk). The total K_s-band luminosity of the truncated model is slightly (∼5%) smaller than the three-component model. The fact that the K_s-band light lost by truncating the disks of the three-component model is almost entirely made up for by the light from the bar and nuclear disk is tantalizing evidence that disk instabilities have merely rearranged the orbits of existing stars in NGC 891's inner regions (e.g., Combes & Sanders 1981; Combes et al. 1990; Raha et al. 1991).

The super-thin disk has a total K_s-band luminosity of 5.05 × 10¹⁰ (6.62 × 10¹⁰) L_{☉, K} for the exponential (sech²) three disk model. The truncated super-thin disk has a total luminosity ∼62% of its non-truncated counterpart. Assuming the nuclear disk is an inner continuation of the super-thin disk and adding its luminosity to that of the truncated super-thin disk increases this ratio to ∼74%. Assuming that the "young" stellar populations used as inputs to RT models correspond exactly to our super-thin disk, the truncated disk luminosity is about an order of magnitude larger than 5.36 × 10⁹ L_{☉, K}, the value assumed by Popescu et al. (2011) for modeling the integrated SED of NGC 891, and roughly four times larger than assumed for our RT modeling. We again note that due to the nature of our attenuation correction Model B tends to have more luminous central SBs, but even the truncated total super-thin disk luminosity of the Model A fits are still ≳2 times as large than the luminosity used in Popescu et al. (2011). These discrepancies are all significantly larger than either our random or systematic errors, which are of order ∼10%.

However, when the total flux (bulge, young disk, and old disk) from the model used by Popescu et al. (2011) is compared to the total flux of our truncated model, they agree to within 10%. It therefore seems likely that the input SEDs used by Popescu et al. (2011) need to be adjusted, with the young stellar SED getting boosted in the NIR at the expense of flux from the old stellar SED. Adjusting the bolometric luminosities of the two components will not solve this problem—boosting the luminosity of the young component at all wavelengths will cause an overprediction of the UV flux (or an overprediction of the FIR light if the extra UV emission is suppressed by dust). It is plausible that Popescu et al.'s (2011) young stellar SED underestimates the impact of cool, evolved stars, especially thermally pulsating AGB stars which are known to be difficult to model and are very bright in the NIR (Maraston 2005). If this is the case it would explain why Popescu et al. (2011) are able to fit the UV and FIR flux of NGC 891 while underpredicting the NIR flux of the young stellar disk component of the galaxy. Depending on the spectral type of these cool stars the peak of the old component's SED will shift, but it is not possible with our current data to predict exactly how these input SEDs will change.

King et al. (1990) report a B-band thin disk integrated luminosity of 6.7 × 10⁹ L_☉, which (even though it's in the optical) should be minimally affected by dust because it is computed away from the midplane of NGC 891. Using our fitted thin-disk K_s-band integrated luminosity, the thin disk has a B−K_s color of ∼4.5 mag (Vega), similar to a K4III (Johnson 1966) star. Using the values from what Popescu et al. (2011) call the old disk, which is most similar to our thin disk, gives a B−K_s color of ∼5.2 mag, close to a K5III or M0III. Based on a galactic evolution model Just et al. (1996) find I−K colors for the thin disks of edge-on galaxies to be ∼1.2 mag, similar to the colors of a K3III star. While these models were not designed specifically for NGC 891, they are another indication that the NIR luminosity in Popescu et al. (2011) needs to be redistributed between the young and old stars.

5.2.1. The Luminosity Profile and Total Light

5.2.2. Component Scales and Scale Ratios

Of value to our investigation would be corroborating evidence of a super-thin disk in NGC 891 in another dataset, ideally in the MIR which has been claimed to be ideal for studies of starlight with minimal attenuation (Sheth et al. 2010). Toward that end we analyzed the archival IRAC 4.5 μm data looking for a super-thin component consistent with our NIR fits. While the IRAC data has resolution only half that of the lowest resolution of our WHIRC data, we note that our K_s-band fits indicate that NGC 891's super-thin disk has a central SB ∼1.5–2 mag arcsec⁻² brighter than its thin disk, and ∼4 mag arcsec⁻² brighter than its thick disk. At 4 kpc in radius and at the approximate resolution of the IRAC data in height (18 ∼ 80 pc), the super-thin disk SB should be comparable to the fitted thin disk SB (∼0.2 mag fainter for the exponential fit, and ∼0.9 mag brighter for the sech² fit, assuming no color gradients between K_s and 4.5 μm). Therefore it is not unreasonable to expect to see at least some evidence of the super-thin disk in the IRAC data, even though it may never be the dominant component in the ranges we can probe.

We computed 2D fits to the IRAC 4.5 μm data, using the same methodology as our fits of the K_s-band WHIRC data, with a three-component disk model. The exponential disk model had difficulty converging on a physically plausible fit, due to degeneracies introduced by the lowered resolution. The sech² model, however, produced a good fit (Table 6 and Figure 22) with parameters very similar to that for the K_s-band data. The scale-heights of the super-thin and thin disk in the fits to the 4.5 μm data are identical to the K_s-band values.

Zoom In Zoom Out Reset image size

In hindsight it is perhaps unsurprising that we were able to identify NGC 891's super-thin disk at 4.5 μm: Comerón et al. (2011b) set out to decompose the galaxies in the S⁴G sample into thin and thick disk components, but several (e.g., IC 2058, NGC 1495, and NGC 5470) have fits that are closer to what would be expected for a super-thin+thin disk. This suggests not only that the IRAC bands are capable of identifying super-thin disks in (very) nearby galaxies, but that such disks may be common in the local universe. Attempting a three-component decomposition of these S⁴G galaxies would be a logical next step in determining if they are truly super-thin+thin+thick disk systems. Indeed, Comerón et al. (2011a) perform a three-component (1D) analysis on S⁴G data of NGC 4013, and find a system that appears to have a MW-like super-thin disk but extended (in height) thin and thick disks. Comerón et al. (2011a) do not identify the narrowest disk as a super-thin disk but rather term the components thin, thick, and "extended"; this difference is semantic until there is a clear astrophysical distinction between the components (e.g., age).

Although the lower-resolution IRAC data appears good enough to identify super-thin disks in at least some nearby edge-on galaxies there are still several advantages to our methodology. First and foremost is the increased resolution in the NIR, which allows super-thin disks to be constrained at much larger distances. This is especially important for the study of MW-sized galaxies, which are relatively rare in the local universe. The increased resolution is helpful even in nearby galaxies: lower spatial resolution increases the likelihood of degeneracies in the fits, which is already known to be an important issue in performing edge-on disk decompositions (Morrison et al. 1997). The primary motivation for undertaking vertical disk decompositions in the MIR is the limited dust attenuation at those wavelengths and the availability of deep imaging data from Spitzer. However, our RT models indicate there is as much as 0.5 mag of attenuation at 4.5 μm (Figures 12 and 14). While most of the 4.5 μm light experiences very little attenuation given the concentration of dust near the midplane, a full, accurate accounting of the MIR light will still require an attenuation correction.

Further, by using K_s-band data we ensure that we exclusively fit the light from the stellar disk(s) and avoid contamination from dust, which Meidt et al. (2012) find to be ⩾5% of total integrated light and up to 20% of the light locally from star-forming regions in spiral galaxies. While MIR dust emission does affect our attenuation correction in the K_s-band by skewing the K_s−4.5 μm color, this is a second-order effect—a correction to a correction in the K_s-band itself. The ability of our model to reproduce the data at a variety of radii and heights (Figure 5) indicates that it has been properly accounted for in our modeling. In contrast, dust emission is both a first and second order correction at 3.6 and 4.5 μm.

Finally, there is additional value in making the attenuation correction in the NIR as it opens up this important spectral region for future imaging and spectroscopic efforts. For instance, there are several molecular band-heads between 1 and 2.3 μm which may be able to separate different types of cold giant stars (Lançon et al. 2007), allowing for age and composition gradients to be determined in the stellar disks.

Before (and during) the data reduction process we inspected the individual images for artifacts and defects. Three images from both the SW and NE K_s-band pointings contained elevated signal levels over a large region near the left-middle of the detector and were removed from the reduction. Additionally, during a few exposures of the central K_s-band and SW H-band the sky background appears to undergo a rapid change, faster than our sky-source duty cycle. We therefore also removed four SW H-band exposures and five central K_s-band exposures from the reduction pipeline.

Data reduction generally followed the procedures described in the WHIRC Quick Guide to Data Reduction manual,⁸ although we modified several steps to optimize the reduction for very extended targets like NGC 891. We detail all reduction steps here for clarity, going into depth for the non-standard portions of our reduction.

First, excess data on the images were trimmed off. Then, images taken with Fowler-4 sampling (J-band) were divided by four to remove the additive nature of the sampling. We then applied a linearity correction to all images with the IRAF task IRLINCOR using coefficients provided in the data reduction manual. The coordinate matrix of the image headers was then adjusted to take into account the rotation of the telescope during the observations.

We took 20 dome flat exposures for each filter: 10 with the dome lamps on and 10 with the lamps off. The 10 lamp-on flat images were scaled to have the same mean flux in the central quarter of the images and then averaged together to produce a combined lamp-on flat; this process was then repeated with the lamp-off flats. The combined lamp-off flat was then subtracted from the combined lamp-on flat to remove any signal not coming from the flat lamps (e.g., dark current, scattered light, or thermal emission in the dome).

Before reducing the flat-fields any further, a first-order correction must be applied because WHIRC has a "pupil ghost," likely due to internal reflections in the imager's optics. It appears in the central 80'' of the detector. this feature is most prominent in the K_s-band, with an amplitude about 15% of the background, but it is also present (albeit at much lower levels) in J and H as well. The ghost was modeled and removed from the flat-fields (it is removed from the object frames naturally during the sky-subtraction process) using the representative pupil template provided on the WHIRC Web site⁹ and the IRAF MSCRED.RMPUPIL task. While pupil templates can be constructed from J and K_s flats, for our purposes the provided template was more than adequate, with errors <0.3% on the corrected images. We note that the flats we used are from an earlier night in our run (UT date 2011 October 16); we found night-to-night deviations in the flat-fields of <1% over the entire image area. All object and sky images were flattened before proceeding with further data processing.

Finally, the central quarter of the pupil-ghost-removed flat was normalized to an average value of 1. Small and negative pixel values (caused by bad pixels and later masked out) were replaced by unity, to prevent erroneously large and/or negative flux values from appearing in the flat-fielded science images.

With extended sources like NGC 891, sky subtraction can be very difficult. To combat this problem, we interleaved observations of a relatively empty sky field ∼10' away (with similar dither amplitudes between sky images as for the data) with our observations of NGC 891. For each set of on-source dithers (four for J and H, nine for K_s) we combine sky images taken before (two for J and H, five for K_s) and after (again, two for J and H and five for K_s) to produce a sky frame suitable for use on the data images. Each sky frame was scaled to have the same median value to avoid biasing the statistics. The sky frames were then median combined iteratively pixel-by-pixel with 2σ clipping until the median converged or all the pixels in the individual sky frames had been clipped out.

After all pixels have been median combined, we used an iterative relaxation algorithm (described in Appendix B) to interpolate sky values for sky image pixels where all the constituent pixels were clipped out. The relaxation procedure provides accurate interpolation and converges very quickly when the pixels needing interpolation are relatively few and widely spaced; we deliberately chose to obtain dedicated sky images and ensured that dither amplitudes were larger than the stellar point-source-functions (PSFs) in order to make this process as efficient and accurate as possible. For a given set of on-source images the appropriate combined sky frame was scaled to have the same median value as the data frame and then subtracted. The resulting sky-subtracted image was then field flattened by the normalized flat-field image.

Prior to registration, the sky-subtracted data images were cleaned of (most) chip defects using the IRAF task FIXPIX (a linear interpolator) and a pixel mask created from the ratio of two lamp-on flats with different exposure times. The data images were then distortion corrected using the IRAF task GEOTRAN and distortion files downloaded from the WHIRC Web site.

WHIRC has a very small dark current (∼0.2 DN s⁻¹ pixel⁻¹) compared to the sky background in most filters. For comparison, WHIRC J-band has an average sky brightness of ∼6 DN s⁻¹ pixel⁻¹. (Dark current is not an issue even in principle for dome flats since we utilize "on-off" differences.) Additionally, the sky subtraction procedure functions to remove most of the (already small) contribution to the image flux from the dark current. However, the dark current exhibits a fixed pattern structure. Further, the correct procedure for dark current removal is subtraction before field-flattening. Nonetheless, due to our limited set of dark frames, we were unable to demonstrate an improvement in image flatness with dark subtraction. Therefore a separate dark subtraction has minimal impact on the reduced images. For future studies focusing on low SB extended structure we do recommend a more thorough investigation of dark-current subtraction.

When compared to a super-sky frame made up of an average of a given night's sky images, the dome flats appear to have small (on the order of 1%) illumination errors. However, the illumination appears to vary between nights and potentially even during nights, making correction difficult. We tested an illumination correction on a small subset of our data and found it made very little difference in the reduced images. Due to the limited improvement and variability of the illumination we do not employ an illumination correction, but again we recommend a more careful exploration of this issue for detailed low SB studies.

Because of WHIRC's relatively small field of view and the extremely extended nature of NGC 891 on the chip, fully automatic registration procedures (e.g., SCAMP; Bertin 2006) were found to be unreliable. Therefore we performed relative image registration between individual exposures by hand, first computing source catalogs for each image using Source Extractor (Bertin & Arnouts 1996), then hand-identifying individual sources that overlapped between images. Based on visual comparisons between individual registered images and the PSF of the final mosaics, our registration scheme is accurate to ≲1 pixel.

With the images all registered in relation to each other, we created full J, H, and K_s mosaics using the Swarp software package (Bertin et al. 2002). Because of the extended nature of NGC 891 compared to the WHIRC field-of-view the estimate of the average flux in the individual data frames used to subtract the combined sky frame is biased, and therefore the background of the sky-subtracted images is oversubtracted. To combat this, we used the high S/N image mosaics to hand-produce very conservative masks of NGC 891 (and foreground stars), essentially masking out all pixels except those very close to the edge of the mosaic where galaxy contamination is minimized. This global mask was then converted into a pixel mask for each individual frame. We then re-ran the pipeline, starting from the sky-subtraction step, this time normalizing the sky images to the source images using only unmasked pixels. The relative World Coordinate System (WCS) offsets were held constant from the first iteration of the pipeline. This two-step process effectively removes the oversubtraction bias from the extended nature of NGC 891 in our images.

The next step was to register the image mosaics in an absolute astrometric system. First we hand-constructed catalogs of the centroids of stars visible on both the WHIRC K_s-band image and the IRAC 4.5 μm image (see Section 2.2 for details on the IRAC images) using Ds9 (Joye & Mandel 2003). Using the IRAC image as an astrometric reference is useful not only because we will be using K_s−IRAC colors in our analysis and therefore want good registration between these datasets but also because both the IRAC and WHIRC data go much deeper than standard NIR calibration catalogs. We then used the IRAF task CCMAP to perform the registration. To preserve the excellent relative agreement between bands the K_s-band absolute calibration was applied manually to the J and H bands.

To flux calibrate the images, stars from the 2MASS PSC with reliable photometry were matched to Source Extractor integrated fluxes on the mosaics. We then fit the relationship between log(DN) and Vega magnitude, obtaining fits with standard errors of 0.02 (0.03, 0.02) mag in the J (H, K_s) bands.