© ESO, 2014

1. Introduction

Dust and gas are crucial ingredients for star formation in the Universe. Dust grains provide a surface on which molecular hydrogen (H₂) gas forms from atomic hydrogen (Hi) gas via catalytic reactions (Gould & Salpeter 1963; Hollenbach & Salpeter 1971; Cazaux & Tielens 2002). Regions of cooling and fragmentation of molecular gas in giant molecular clouds (GMCs) are the dominant formation sites of stars in galaxies (see Fukui & Kawamura 2010, and references therein). These stars produce the heavy elements via nucleosynthesis in their cores. Upon the death of a star, metals are expelled into the interstellar medium (ISM) and either mix with the gas phase or condense to form dust grains in the metal-rich cooling gas (see e.g., Nozawa & Kozasa 2013), such as in the expanding ejecta of novae and supernovae (e.g., Kozasa et al. 1991; Nozawa et al. 2003; Barlow et al. 2010), the stellar winds from cool, dense atmospheres of asymptotic giant branch stars (e.g., Ferrarotti & Gail 2006; Zhukovska & Henning 2013), and the common-envelope in close binary systems (e.g., Lü et al. 2013). Growth, modification and destruction processes may affect the overall abundance of interstellar dust (e.g., Clemens et al. 2010; Mattsson et al. 2012), whilst surviving grains have the potential to aid future episodes of star formation. Thus, the evolutionary processes governing gas, stars, metals and dust are mutually related, and constraining the relationships between these four components is indispensable for understanding galaxy formation and evolution.

Detailed studies of the interaction between gas and stars have been technologically feasible for a long time. Global star formation rates (SFRs) of galaxies are observed to be well correlated with their atomic gas content (e.g., Schmidt 1959; Kennicutt 1998). There is a power-law correlation between the surface densities of gas (Σ_gas) and SFR (Σ_SFR), with numerous studies attempting to constrain the power-law index via global quantities from galaxy samples and resolved observations of individual systems (see Kennicutt & Evans 2012, for a review). Spatially resolved data has uncovered that Σ_SFR is more strongly correlated with the surface density of molecular hydrogen (Σ_H2) than atomic hydrogen (Σ_Hi) gas (Wong & Blitz 2002). More recent studies have begun to probe this relationship at sub-kpc scales, approaching the typical size of GMCs and the intrinsic physical scale of star formation, for larger samples of galaxies (e.g., Bigiel et al. 2008; Momose et al. 2013). Strong observational evidence also shows a correlation between the stellar mass and the gas-phase metallicity (e.g., Lequeux et al. 1979; Tremonti et al. 2004). Recent efforts are attempting to observationally constrain the relationships between stellar mass, metal content and Hi gas content (e.g., Hughes et al. 2013; Bothwell et al. 2013).

Studying the dust content has proved a difficult endeavour mainly due to technological limitations. Dust absorbs up to 50% of the optical and ultraviolet (UV) photons from stars. The dust can then be traced via the re-emission of this energy at mid-infrared (MIR), far-infrared (FIR) and submillimetre (submm) wavelengths. Previous space missions − necessary to avoid atmospheric absorption in these bands − such as IRAS, ISO, Spitzer, and AKARI (respectively: Neugebauer et al. 1984; Kessler et al. 1996; Werner et al. 2004; Murakami et al. 2007) typically operated at wavelengths λ ~ 240 μm, potentially missing a great part of the spectral energy distribution (SED) associated with a cold dust component. Fortunately, the success of the Herschel Space Observatory (Pilbratt et al. 2010) satellite has substantially improved this situation. With two of its instruments, the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) and the Spectral and Photometric Imaging REceiver (SPIRE, Griffin et al. 2010), Herschel observations are capable of constraining the FIR SED, and thus detect emission at λ> 250 μm originating from a cold dust component (e.g., Bendo et al. 2010a), via imaging in six wavebands from 70 to 500 μm with significantly higher sensitivity and angular resolution than previous FIR/submm experiments.

Recent observational studies have used Herschel data to investigate the correlations between the dust component and other galaxy properties (see e.g., Smith et al. 2012a; Bourne et al. 2013; Rémy-Ruyer et al. 2013). For example, scaling relations have been found between the total masses of dust, H₂ and Hi in late-type Virgo cluster galaxies (Corbelli et al. 2012), between the dust-to-stellar mass ratio and NUV − r colour (Cortese et al. 2012), and between the dust temperature and SFR per unit dust mass (Clemens et al. 2013). A growing number of studies are using the unprecedented angular resolution of the Herschel observations to perform pixel-by-pixel fitting of the FIR/submm SED and map the main dust properties of nearby galaxies (e.g., Smith et al. 2010; Galametz et al. 2012; Boquien et al. 2013; Draine et al. 2014; Tabatabaei et al. 2014). Observations obtained by the Herschel Exploitation of Local Galaxy Andromeda (Fritz et al. 2012) enabled highly-detailed maps of the dust properties of M31 to be produced via such pixel-by-pixel SED fitting (Smith et al. 2012b), and allow for an in-depth study of Andromeda’s ISM (Ford et al. 2013; Viaene et al. 2014; see also Draine et al. 2014). Similarly, Foyle et al. (2012) and Mentuch Cooper et al. (2012) have examined the gas and dust using Herschel observations to carry out a pixel-by-pixel analysis of the inner regions of M83 and M51, respectively. Neither galaxy shows strong variations in the dust-to-gas mass ratio, except for a small gradient in M51 (Mentuch Cooper et al. 2012), whereas many other galaxy studies have found strong gas-to-dust ratio gradients (e.g., Muñoz-Mateos et al. 2009; Bendo et al. 2010b; Magrini et al. 2011). Although much effort has been made to quantify the radial dust distribution in these face-on spiral galaxies, we also require information on the vertical | z | structure for a complete picture of the three dimensional properties of dust and gas in galaxies.

Edge-on spiral galaxies with inclinations at or near 90° offer the best opportunity to study the radial and vertical structure of the various galaxy components. The stellar component has been shown to be sensitive to the merging history of a galaxy, whereas the ISM may probe gravitational instability (see e.g., Dalcanton et al. 2004; Yoachim & Dalcanton 2006). Dust lanes feature prominently in optical images of massive edge-on spirals, enabling the study of dust in both emission and absorption. Furthermore, the increased surface brightness arising from the line-of-sight projection allows for both the vertical and radial dust distribution to be studied in detail at greater distances from the centre. Edge-ons therefore provide an important contribution to our knowledge of the dust structure. Fortunately, we have one such prototypical example of an almost perfect edge-on spiral right in our neighbourhood: NGC 891. Located at a distance of ~9.6 Mpc (e.g., Strickland et al. 2004) with an inclination angle of >89° (e.g., Xilouris et al. 1998), this bright, non-interacting SA(s)b galaxy (de Vaucouleurs et al. 1976) presents similarities to our own Milky Way in Hubble type, rotational velocity and optical luminosity (see Table 1). These properties make NGC 891 an ideal target for studying the stars, gas and dust in the disc, and so it has already been extensively observed at a range of wavelengths.

The stellar content in the galaxy, traced in the optical wavelengths (e.g., van der Kruit & Searle 1981; Ibata et al. 2009) and at 3.6 μm (Kamphuis et al. 2007), are found to be fairly symmetric along the galactic plane. However, the Hi gas shows a slight north-south asymmetry in both disc and halo thickness (Swaters et al. 1997), and extends up to at least 22 kpc above the galactic plane (Oosterloo et al. 2007). In contrast, molecular hydrogen gas observed via CO emission is predominantly found in a thin disc around 400 pc thick (Scoville et al. 1993), and also seen in the halo up to 1.4 kpc above the plane (Garcia-Burillo et al. 1992). The Hα emission, a tracer of ionized gas linked to star formation, is more prominent and extended in the | z | direction on the northern side of the galaxy than the southern (Dettmar 1990; Rand et al. 1990). Whether the north-south asymmetry arises due to higher SFR in the north (Rossa et al. 2004) or increased attenuation by dust in spiral arms along the line-of-sight through the disc (Kamphuis et al. 2007; Schechtman-Rook et al. 2012) remains an open issue, which may be resolved by studying the cold dust in the galaxy.

Evidence for cold dust in NGC 891 first came from IRAS observations at FIR wavelengths (Wainscoat et al. 1987), and later at millimetre wavelengths using the IRAM 30 m telescope (Guelin et al. 1993). The 1.3 mm dust emission was noted to correlate strongly with the H₂ gas traced by CO emission, but poorly with the Hi emission. Later observations at λ 450 and 850 μm with the Submillimeter Common-User Bolometer Array (SCUBA) at the James Clerk Maxwell Telescope (JCMT) found a large amount of cold dust (~15 K) in the disc (Alton et al. 1998). The submm emission was confirmed to correlate spatially with the molecular gas in the disc (Alton et al. 1998; Israel et al. 1999). The combination of the SCUBA observations, IRAS and ISO data, with maps from the French balloon-borne PRONAOS experiment (Serra et al. 2002) traced the complete FIR/submm SED of NGC 891 for the first time (Dupac et al. 2003a), finding a dust temperature of ~18−24K and estimating a global gas-to-dust mass ratio of ~240. Popescu et al. (2004) present deep maps of dust emission at 170 and 200 μm obtained with the ISOPHOT instrument. However, the low resolution of these datasets, typically ≥1′, meant a pixel-by-pixel analysis of the FIR/submm SED to map the cold dust distribution was not previously feasible.

New, high resolution (≤) Herschel PACS/SPIRE observations of NGC 891 have been obtained as part of the Herschel Guaranteed Time Key Project the Very Nearby Galaxies Survey (P. I.: C. D. Wilson), which aims to study the gas and dust in the ISM of a diverse sample of 13 nearby galaxies using Herschel. These observations now make it possible to trace the FIR/submm SED in detail and map the resolved dust properties in the galaxy disc. In this paper, we investigate the connection between dust and gas in NGC 891 in order to test for asymmetries in the distributions of the dust properties and the gas-to-dust ratio. We use these new Herschel PACS/SPIRE observations combined with a JCMT SCUBA 850 μm image to trace the FIR/submm SED. For the first time, we fit one-component modified blackbody models to the SEDs of each pixel to map the dust mass and temperature in the galaxy. Hi 21 cm line emission and ¹²CO(J = 1−0) emission maps, combined with a new JCMT HARP-B ¹²CO(J = 3−2) map (P. I.: T. Parkin), trace the atomic and molecular hydrogen gas, enabling the study of the gas-to-dust ratios of these components. Our paper is structured as follows. In the next section, we present the new and existing data used in our analysis. In Sect. 3.1, we present our detection of extraplanar dust. In Sect. 4, we describe our SED fitting methodology and the resulting integrated and resolved dust properties. Section 5 describes the relationships we find between the dust and gas, and Sects. 6 and 7 present our discussion and conclusions, respectively.

2. Observations

We have assembled together a multi-wavelength dataset capable of tracing the key ISM components in NGC 891 (see Fig. 1 and Table 2). Our new VNGS Herschel PACS photometric observations and a new JCMT HARP-B ¹²CO(J = 3−2) emission map are presented here for the first time.

Fig. 1 Images of NGC 891 at various wavelengths, centred on the galaxy’s coordinates from NED: , (J2000.0). Each image is presented in its native resolution and pixel size, with the open symbols at the bottom left of each panel indicating the beam size at each wavelength. The image characteristics are described in detail in the text (see Sect. 2). The solid white ellipse denotes the aperture for integrated photometry and the dashed white circles represent the apertures to reconstruct the analysis of Dupac et al. (2003b), as described in Sect. 4.4. North is up and east is to the left.

2.1. Mid-infrared

Spitzer Infrared Array Camery (IRAC; Fazio et al. 2004) 3.6 μm data were obtained in astronomical observation requests 3631872, 3632128, 3632384, and 3633152. Individual corrected basic calibration data frames were processed with version 18.25.0 of the IRAC pipeline and remosaicked using the standard IRAC pipeline within the MOsaicker and Point source EXtractor (Makovoz & Khan 2005). The final images have pixel scales of , and the images are oriented with north up and east to the left. The PSFs have FWHM of and calibration uncertainties of 3% (IRAC Instrument and Instrument Support Teams, 2013, IRAC Instrument Handbook, Version 2.0.3, JPL, Pasadena).

The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) imaged NGC 891 as part of an all-sky survey in four bands: 3.4, 4.6, 12, and 22 μm. For the purposes of our analysis, we only use the latter two MIR bands. The angular resolution is and at 12 and 22 μm, and the astrometric precision for high signal-to-noise (S/N) sources is better than . Although deconvolution techniques may enhance the WISE image resolution, the mosaics with standard beams are sufficient for our analysis.

Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) 24 μm data were reprocessed by Bendo et al. (2012b) using the MIPS Data Analysis Tools (Gordon et al. 2005) along with additional processing steps. The final mosaic has a pixel scale of , the PSF FWHM is 6″, and the calibration uncertainty is 4% (Engelbracht et al. 2007).

2.2. Far-infrared

VNGS Herschel PACS observations at 70 and 160 μm comprised four cross-scans each, taken at a 20″ s^-1 scan rate. We also include an additional PACS observation at 100 μm, taken as part of The Herschel EDGE-on galaxy Survey (HEDGES; P.I.: E. Murphy), which consists of only two cross-scans at a 20″ s^-1 scan speed (Obs IDs: 1342261791, 1342261793) but with a smaller separation between individual scan legs. The Herschel interactive processing environment (HIPE, v.11.0.1; Ott 2010) with PACS Calibration (v.32) was first used to bring the raw Level-0 data to Level-1, which means flagging of pixels, flat-field correction, flux conversion into units of Janskys, and assigning sky coordinates to each detector array pixel. Scanamorphos (v.22.0, Roussel 2013) was then used to remove the 1/f noise, drifts and glitches by exploiting the redundancy of the observations in each sky pixel to estimate the optimal baseline correction for each detector pixel. This process produced the final images. The maps have pixel sizes of , and for the 70, 100 and 160 μm maps respectively. These pixel sizes correspond to one quarter of the point spread function (PSF) full width at half maximum (FWHM) for the scan speed used in our observations (PACS Observers’ Manual 2013). We apply a colour correction factor, taken from the ICC PACS Photometer Passbands and Colour Correction Factors document¹, adopting values relative to a modified black body with temperature of 20 K and emissivity index β = 2, which are both appropriate for the SED of NGC 891 (see Sect. 4).

2.3. Submillimetre

VNGS Herschel SPIRE photometric observations were taken at 250, 350, and 500 μm. The galaxy was observed in large map mode, covering an area of 20′× 20′ with two cross-scans at a 30″ s^-1 scan rate. The SPIRE data were processed up to Level-1 using a script based on the standard pipeline, and was run in HIPE continuous integration build number 11.0.1200. There are two differences to the standard pipeline: firstly, we use the latest calibration due for inclusion in HIPE 12.0 (Bendo et al. 2013), and secondly we do not apply the default temperature drift correction and de-striping. Instead we use a custom method (BriGAdE; M. W. L. Smith, in prep.) to remove the temperature drift and bring all the bolometers to the same level. SPIRE instrument characteristics, e.g. PSF, beam size, calibration uncertainty, etc., are taken from the SPIRE Observers’ Manual (v.2.4, 2011). Artefacts introduced due to extended bright sources were avoided by masking the galaxy prior to removing the median baselines from each bolometer timeline. The HIPE mapmaking procedure was then used to produce the maps. The main beam size has a FWHM of , and at 250, 350 and 500 μm maps, respectively, and the resulting images have pixel sizes of 6, 8 and 12″. We note that in addition to the use of more up-to-date versions of the pipeline and calibration tree compared to the maps previously presented in Bianchi & Xilouris (2011), we adopt new beam areas of 446, 794 and 1679 arcsec² at 250, 350 and 500 μm, respectively². We also apply a calibration factor to account for extended sources and updated colour corrections based on an index of α_s = 4 (SPIRE Data Reduction Guide, v.2.3, 2013), which is appropriate for β = 1.8 (see Sect. 4.3).

We use the JCMT SCUBA observations at 850 μm presented in Alton et al. (1998). The dedicated SCUBA software package SURF (Jenness & Lightfoot 1998), was used to clean, flat-field, and calibrate the image according to atmospheric attenuation. The measured half-power beam width of the central beam is at 850 μm, with sidelobes estimated as 20% (from the same observing configuration, see Bianchi et al. 1998). Alton et al. estimate the error in the calibration to be less than 15% with a 1σ noise level of 3.5 mJy/ beam at 850 μm. The original data from Alton et al. (1998) has since been recalibrated and presented in Haas et al. (2002) and Whaley et al. (2009). Residual noise seen at the SCUBA field edges is avoided in our analysis (see Sect. 4.4). We apply a conservative correction of 5% to the flux density to account for contamination due to the ¹²CO(J = 3−2) line (see e.g., Israel et al. 1999).

2.4. CO and Hi line data

We obtained new ¹²CO(J = 3−2) observations, which were carried out on 19 December 2009 and 25 January 2010 using the HARP-B receiver mounted on the JCMT, as part of project M09BC05 (P. I.: T. Parkin). We used quarter-step spacing and the basket weave pointing technique (see Appendix of Warren et al. 2010) to produce a raster map of the galaxy 3.25′ × 8.75′ in size, at a resolution of 14″. The total integration time required for this observation was 2180 s. The sky background was quantified by position-switch chopping of the secondary mirror with a chop throw of 60″. We used the Auto-Correlation Spectrometer Imaging System (ACSIS) for the backend receiver, set to a bandwidth of 1000 MHz with 2048 channels. This configuration produces a spectral resolution of 0.43 km s^-1. We processed the raw data with Starlink³ (Currie et al. 2008). These data were flagged for bad pixels, combined into a data cube with a pixel-size of , and then baseline subtracted by fitting a third order polynomial to the line-free regions of the spectrum. Finally, the cube was rebinned in the spectral dimension to a resolution of 20 km s^-1. The average root mean-square (rms) noise in the spectrum then produced a signal-to-noise cube. The integrated intensity map was created by running the Starlink routine “find-clump” on the signal-to-noise cube to identify regions of emission with a S/N> 3. This cube was collapsed along the spectral axis to produce a two-dimensional signal-to-noise integrated intensity map, which was subsequently multiplied by the noise to produce the final integrated intensity map. The integrated intensity measurement uncertainty is a function of the channel width, the rms noise, and the number of channels in the spectral line and baseline. The calibration uncertainty is estimated to be 10% based on a comparison of spectra taken of calibrators during observations and their standard spectra. For more details on the data reduction process and uncertainty measurements, see Warren et al. (2010) and Parkin et al. (2012).

A high resolution map of the ¹²CO(J = 1−0) emission is taken from Scoville et al. (1993). The map was made using the Owens Valley Millimeter Array (OVMA) to observe nine fields along the semimajor axis, covering ± from the centre. The synthesized beam from the different telescope configurations ranged from to with a mean of . The spectral resolution is based on 32 channel filter banks with 5 MHz resolution and the continuum calibration was performed with a broad-band channel with an effective bandwidth of 375 MHz. Observations of 3C 84 and Uranus provided the phase and absolute flux calibrations. The nine fields were cleaned and mosaicked. The final map has a rms noise of 20 mJy beam^-1.

Finally, neutral hydrogen 21 cm line emission data necessary for quantifying the Hi content was available from Oosterloo et al. (2007), taken as pathfinder observations for the Hydrogen Accretion in LOcal GAlaxieS survey (HALOGAS; Heald et al. 2010, 2011) with the Westerbork Synthesis Radio Telescope (WSRT). Standard array configurations for the WSRT were used to complete 20 × 12 h exposures, with a total bandwidth of 10 MHz using 1024 channels. The data was reduced using the MIRIAD package (Sault et al. 1995). We use their intermediate resolution data, which has a beam FWHM of 30″ (see Fig. 1. in Oosterloo et al. 2007). The map is the deepest Hi observation available for NGC 891.

3. FIR/submm morphology

Before proceeding with our main analysis, we first take a brief diversion to exploit one of the main advantages of studying an edge-on galaxy like NGC 891 − the ability to study the vertical and radial distribution of the FIR/submm emission.

3.1. Detection of extraplanar dust

The presence of a spatially-resolved vertical profile in the FIR/submm emission of an edge-on galaxy may be related to dust ejected from the galactic plane, i.e. extraplanar dust. Determining the prevalence of extraplanar dust and the relationship to other galaxy properties may shed light on the mechanisms responsible for transporting ISM material from the galactic plane. Deep optical to near-infrared (NIR) imaging has revealed significant amounts of extraplanar dust (e.g., Howk & Savage 1999; Thompson et al. 2004) as well as PAHs and small dust grains (e.g., Whaley et al. 2009) in several edge-on spirals, including NGC 891. Most recently, extraplanar dust was detected at the 5σ level in the small edge-on spiral NGC 4244 (Holwerda et al. 2012a) and one (NGC 4013) of the seven edge-on spirals comprising the HERschel Observations of Edge-on Spirals (HEROES, Verstappen et al. 2013).

We extracted vertical profiles of NGC 891 from the PACS 100, 160 and SPIRE 250 μm images. We selected an area surrounding the galaxy and summed all the pixel values along the major axis to generate the profiles. Following Verstappen et al. (2013), we model these profiles with an exponential function appropriate for an exactly edge-on, double-exponential disc, given by (1)where h_z is the dust scaleheight. We first convolve the vertical profile model with the Herschel beams at each corresponding wavelength, using the circularised PSF images⁴described in Aniano et al. (2011). In order to obtain one-dimensional beams, the two-dimensional PSFs are averaged along one direction in the same manner as we obtain the vertical profiles. The optimal value of h_z that reconciles the observed and model profiles was found by using a χ² minimisation technique. Uncertainties on the value of h_z were also derived from the χ² probability distribution. We estimate the uncertainty in the profiles due to the galaxy’s position angle (PA), which varies in the literature from 22 to 23°, by comparing the resultant h_z values and profiles when systematically varying the assumed position angle from 21 to 24° in 0.1° increments. For position angles moving out of the range of 22.6 < PA < 23.1°, the profiles demonstrate increased broadening and an expected increase in the corresponding h_z value. The profile widths and scaleheights tend to reach minimum values within the 22.6 < PA < 23.1° range. Combining these results, the recent result from PA fitting performed by Bianchi & Xilouris (2011), and a visual inspection of the data, we decide to adopt PA = 22.9° for our analysis.

The resulting vertical profiles are shown in Fig. 2. We observe a slight asymmetry in the profiles towards the NW side of the disc (represented by positive vertical offsets in Fig. 2), which are also seen in the submm profiles obtained by Bianchi & Xilouris (2011, see their Fig. 4). For NGC 891, we derive FIR scaleheights of 0.24, 0.43 and 1.40 kpc at 100, 160 and 250 μm, respectively. The scaleheight at 100 μm is consistent with scaleheights found in previous studies using radiative transfer modelling of optical and NIR images. For example, Xilouris et al. (1998) found a mean dust scaleheight of 0.26 kpc with a range of 0.22 to 0.32 kpc from models based on B to K band images. Similar results were recently found in the modelling by Schechtman-Rook et al. (2012). In such a comparison of scaleheights derived from different methods, we should not forget that (i) the FIR profiles have a smaller resolution (e.g. ~7″ at 100 μm) and are much more dominated by the PSF FWHM than compared to optical images used for the radiative transfer modelling (~1″); and (ii) our implicit assumption of a perfectly edge-on disc likely differs from the inclination adopted by the radiative transfer models. However, NGC 891’s near perfect edge-on inclination (89.7°, e.g., Xilouris et al. 1998) means the latter issue is likely negligible.

Fig. 2 Vertical profiles at 100, 160 and 250 μm, adopting a position angle of 22.9°, are plotted as blue circles connected by a blue line. We fit with an exponential profile (solid red line), which has been convolved with the Herschel beam at the corresponding wavelength (dashed black line).

Given that the profiles are not dominated by the telescope beam, as evident in Fig. 2, and since the deconvolved scaleheight value we derive from the profile fitting is not consistent with zero at the 5σ level, we conclude that our vertical FIR profiles are spatially resolved. Such a conclusion confirms previous findings of extraplanar dust in NGC 891 inferred from optical-NIR imaging (see e.g., Howk & Savage 1999). However, whether this material has been expelled from the galactic plane or infalling from the intergalactic medium (e.g., Howk 2009) is still not understood.

3.2. Break in PACS radial profiles

A break in the radial profiles in the SPIRE bands, found by Bianchi & Xilouris (2011), was shown to occur in all bands at ~12 kpc from the centre, with the profiles sharply declining at larger radii with a reported radial scalelength of about 1 kpc. A similar break was also detected in the optical images of Xilouris et al. (1999). We now investigate whether we observe the same behaviour in the radial profiles of the PACS bands. We produce horizontal profiles of the FIR/submm emission from the PACS and SPIRE data via a similar method to that described for extracting the vertical profiles. The flux was summed along pixels parallel to the minor axis and divided by the number of pixels used for each strip, thus yielding an average surface brightness profile. The resulting profiles are presented in Fig. 3. The SPIRE 250 μm profile matches that found by Bianchi & Xilouris (2011). Our profiles from the PACS bands broadly display similar features as the SPIRE profiles: a central peak with two clear secondary peaks within radial distances of approximately ±5 kpc, a sharp decrease in surface brightness between ±11 and 13 kpc, and a flattening of the profile at radial distances greater than ±13 kpc. Thus, we also observe the same break in the PACS bands as the analysis of the SPIRE profiles performed by Bianchi & Xilouris (2011).

The most interesting observation to note is the correlation between the profiles of the IRAC 3.6 μm emission and the MIPS 24 μm emission, which respectively trace the total stellar population and star formation. Within the break radius at ±12 kpc (dashed line in Fig. 3) the FIR profiles appear to correlate best with the 24 μm emission. At radial distances outside the break radius, the FIR profiles on the SW edge of the disc (i.e. negative radial offsets) appear to flatten out to reflect the 3.6 μm emission, whereas the profiles on the NE edge (i.e. positive radial offsets) are only traceable to ~85% of the optical size. This appears to hint that the dust properties change across the break radius, possibly as the dominant source of the dust heating may change from star-forming regions to the stellar population. We return to the topic of the source of dust heating in more detail in the following section.

Fig. 3 Radial profiles at 100, 160 and 250 μm are compared with the IRAC 3.6 μm and MIPS 24 μm emission. The solid black line represents the optical diameter adopting D25 = 37.7 kpc and the two dashed black lines demark the break radius found at ±12 kpc from Bianchi & Xilouris (2011). The alternative y-axis scale refers only to the dust temperature profile obtained from modelling the SED (see Sect. 4.3), represented by the solid red line.

4. Dust properties of NGC 891

We now present our methodology to determine the global and spatially-resolved dust properties of NGC 891.

4.1. Image preparation

In our analysis, we fit simple models to both the integrated SED of the galaxy and to the individual SEDs of each pixel, so that we may derive the global dust properties and map their spatial distributions. However, the various images described above were taken using several different instruments and detectors, with differing PSFs of varying shapes and widths. To properly compare the data to each other and perform a pixel-by-pixel comparison, the physical scale of each image should be consistent. We therefore require convolution kernels that will transform all the images into a common PSF, so we can generate image cubes in which each pixel of each image corresponds to the same sky region. Fortunately, common-resolution convolution kernels have already been made publicly available ⁴for the instrumental PSFs of several space- and ground-based telescopes, as well as general purpose Gaussian PSFs (Aniano et al. 2011). We use the dedicated kernels for PACS, SPIRE and WISE images, which were generated with the latest available PSF characterisation for each instrument, and use the appropriate Gaussian kernels for the remaining images.

All images were first converted into surface brightness units (i.e. MJy sr^-1) and convolved to the resolution of the 500 μm image, since this band has a PSF with the largest FWHM (), using the IDL task convolve_image.pro (Aniano et al. 2011). The images were regridded to the pixel size of the 500 μm map using the Montage software package. We note that since the pixel size () is smaller than the 500 μm beam size (), adjacent pixels are not independent. Errors on each pixel were calculated by summing the flux calibration uncertainty, instrumental noise and sky background measurement in quadrature. For the pixels covering the galaxy, the flux errors are dominated by the calibration uncertainty. We use calibration uncertainties of 5% for PACS (PACS Observers’ Manual, 2011), 7% for SPIRE (SPIRE Observers’ Manual, 2011) and 25% for SCUBA (Haas et al. 2002). The sky background error was measured using the standard deviation of the flux from ten circular background apertures, which were carefully placed around the galaxy at a suitable distance to avoid diffuse galaxy emission. For each image, the mean sky flux estimated from these background apertures was subtracted.

4.2. Source of dust heating

Before fitting the SED data for NGC 891, we assessed how the dust was related to its heating sources. Several authors (e.g., Bendo et al. 2010a, 2012a; Boquien et al. 2011) have compared dust colour variations to emission from either star forming regions or the total stellar populations and have demonstrated emission seen at ≥250 μm by Herschel may originate from a different thermal component than the emission seen at ≤160 μm. We perform a similar analysis on the NGC 891 data by comparing the 100/160, 160/250, 250/350, and 350/500 μm flux ratios to 3.6 μm emission (a tracer of the total stellar population that is relatively unaffected by dust extinction) and 24 μm emission (a tracer of star formation that is also unaffected by dust extinction e.g., Calzetti et al. 2007). The galaxy is viewed edge-on, which may cause some issues with properly relating dust emission to heating sources, as emission from both star forming regions and diffuse emissions is integrated along the line of sight. On the other hand, the distribution of starlight and star forming regions is distinctly different in this galaxy. Emission from star formation peaks in the centre and in a couple of star forming regions in the disc, while 3.6 μm emission is relatively flat within the centre (because of the presence of the galaxy’s bulge) and then smoothly decreases with radius (see Fig. 1). Because of this, associating colour variations with dust heating sources is relatively straightforward.

The results are shown in Fig. 4, with data along the major axis shown in blue and data off-axis shown in grey. The off-axis pixels tend to exhibit scatter, potentially because of either low signal-to-noise in the data or issues with matching the extended structure in the PSFs following the convolution process described above, so we focus on the on-axis pixels. The 100/160, 160/250, and 250/350 μm ratios all correlate better with the 24 μm emission than with the 3.6 μm emission. In particular, the relations between the 100/160, 160/250, and 250/350 μm ratios and the 3.6 μm emission flatten out at high surface brightnesses, which indicates that the colours are independent of the galaxy’s stellar surface brightness. This suggests that the 100−250 μm emission originate from a component of dust heated by star forming regions. The results from the 350/500 μm ratios are more ambiguous. It is possible that the dust seen in this wavelength range is heated by a combination of light from star forming regions and diffuse light from the total stellar population. However, it is also possible that the dust emitting at ≥350 μm is still primarily heated by star forming regions and that the colour variations are difficult to detect because the emission originates from the Rayleigh-Jeans side of the SED.

These results suggest that, following integration along the line-of-sight, emission from dust in NGC 891 at >250μm heated by the diffuse interstellar radiation field is either weak or absent in the middle of the plane. This would contradict the results from other spiral galaxies (e.g., Bendo et al. 2010a, 2012a; Boquien et al. 2011; Groves et al. 2012; Smith et al. 2012b), which all found significant amounts of emission from dust heated by the total stellar population. However, a scenario where dust heating in NGC 891 is dominated by star-forming regions rather than the total stellar population could explain why the models of Bianchi & Xilouris (2011) required a stellar disc with radial scalelength of 5.7 kpc, more similar to that of young stars (traced via B-band imaging), instead of a 4 kpc radial scalelength derived from the old stellar population (emitting in NIR bands), to successfully reproduce the major axis profiles in the SPIRE bands. Furthermore, the results for NGC 891 would be consistent with those found for the dwarf irregular galaxy NGC 6822 (Galametz et al. 2010). Bendo et al. (in prep.) are performing a follow-up analysis on dust heating over a broader range of spiral galaxies and will discuss additional galaxies where dust may be heated primarily by star forming regions.

Fig. 4 The 100/160, 160/250, 250/350, and 350/500 μm flux ratios compared with the 3.6 and 24 μm emission. Pixels along the major axis shown in blue and off-axis pixels shown in grey. The best-fit linear relation (red line) and the Pearson correlation coefficient are based only on the on-axis pixels.

4.3. SED fitting

Based on the above analysis, it is appropriate to follow previous works (e.g., Smith et al. 2010, 2012b; Verstappen et al. 2013) in fitting the SEDs with a one-component modified blackbody model given by (2)which is a modified version of the equation originally presented by Hildebrand (1983). In this equation, M_dust is the dust mass, B_ν is the Planck function, T_dust is the dust temperature, κ_ν is the dust emissivity and D is the distance to the galaxy (adopting 9.6 Mpc, Strickland et al. 2004). We perform this fit on all data from 100 to 500 μm and, following an approach suggested by several authors (e.g., Bendo et al. 2010a; Smith et al. 2010, 2012b), treat the 70 μm emission, which may originate from small grains with stochastic dust heating, as an upper limit to the fit. Although some of the individual surface brightness ratios for the off-axis regions shown in Fig. 4 exhibited some scatter, the results from the SED fitting across the 100−500 μm should mitigate the effects that noise or PSF mismatch issues in any individual band will have on characterising the overall shape of the SED.

Fig. 5 Left: total FIR/submm SED of NGC 891 obtained from the integrated PACS, SPIRE and SCUBA fluxes (black circles). The red line shows the best-fit one-component model, where the 70 μm point is used as an upper limit. The dashed black line represents the best-fit one-component model found when including the FIR/submm GAUFLUX fluxes reported in the Planck HIFI Catalogue of Compact Sources (blue diamonds; Planck Collaboration XXVIII 2014). For comparison with previous FIR/submm experiments, the flux measurements from IRAS, ISO and IRAM presented in Popescu et al. (2011, see their Table 2) are shown (grey triangles). Right: SEDs measured for different locations in the galaxy, following Dupac et al. (2003a) and Bianchi & Xilouris (2011).

We assume a power-law dust emissivity in the FIR/submm wavelength range, with κ_ν ∝ ν^β. Unfortunately, the absolute value of the dust emissivity and its frequency dependence are still uncertain (e.g., Dupac et al. 2003; Gordon et al. 2010; Planck Collaboration XIX 2011). A variation in the value of κ_ν translates the model SED flux densities, systematically increasing or decreasing the dust property estimates without affecting any correlations that we may find. Following previous works, we adopt κ_ν = 0.192 m² kg^-1 at 350 μm (Draine et al. 2007), but note that this may vary as a function of β (see e.g., Bianchi 2013, and references therein). The choice of value for β is also not so straightforward, because this can potentially create both systematic offsets in the dust property estimates and introduce secondary effects in correlations between the dust mass/temperature estimates, and other quantities we wish to study, such as e.g. the various components of the total gas content. Such effects can lead to an unreliable analysis. There is still much debate about whether the value of β should have a fixed value, or vary as a free parameter, as the dust emissivity spectral index may provide information about the physical properties of the radiating dust grains (see e.g., Smith et al. 2012b). In our analysis, we primarily chose to allow M_dust, T_dust and β to vary as free parameters, and compare our results for the case where β = 1.8 (see Sect. 4.5), the value for dust in the galactic disc found by the Planck Collaboration XIX (2011) and consistent with other studies of nearby galaxies (Boselli et al. 2012; Galametz et al. 2012; Cortese et al. 2014). The fits were done by performing a χ² minimisation using a simple gradient search method.

We estimate errors on the best-fit model parameters via a bootstrap technique. After the best-fit parameters are determined, 200 new sets of data points are created by selecting random flux values from within the errorbars of the observed fluxes. Each new dataset is then fitted again to find the alternative best fitting parameters. We then calculate the 68% interval in both the upper and lower parameter distributions, setting these intervals as the new upper and lower limits. The differences in the parameter values of the original best-fit solution and the extreme values from the bootstrap technique are then taken as our uncertainties in the best-fit model parameters.

4.4. Integrated dust properties

Before proceeding with the resolved dust analysis, we first tested our methodology using integrated quantities. We determined the total flux densities in each FIR (i.e. 70 to 850 μm) waveband via aperture photometry utilizing funcnts, part of the funtools software for DS9. An elliptical source aperture centered on the galaxy was defined based on the extent of the 500 μm emission (shown in Fig. 1). We use identical source apertures for all images except for the SCUBA 850 μm image, which covers a much smaller area than the other observations. To avoid the residual noise seen at the SCUBA field edges, we reduced the size of the source aperture used to measure the total galaxy emission in this band.

Our PACS flux densities are 108 ± 8 Jy, 259 ± 13 Jy and 352 ± 28 Jy at 70, 100 and 160 μm respectively. These fluxes are consistent within the errors with the MIPS photometry presented by Bendo et al. (2012b), reported as 97 ± 10 Jy and 287 ± 35 Jy at 70 and 160 μm. Some of the variation between the MIPS and PACS fluxes may be the result of differences in the filter profiles. We successfully reproduce the SPIRE flux densities found by Bianchi & Xilouris (2011): F₂₅₀ = 155 ± 11 Jy, F₃₅₀ = 66 ± 6 Jy and F₅₀₀ = 23 ± 3 Jy. The discrepancies likely arise from a combination of differences in the choice of apertures in our photometry, updates to the flux calibrations and beam sizes, and sky subtraction. The SCUBA 850 μm image has a total flux density of 5 ± 1 Jy, consistent with values reported by Alton et al. (1998, F₈₅₀ = 4.6 ± 0.6 and Israel et al. (1999, F₈₅₀ = 4.8 ± 0.6, who use an older calibration (see also Haas et al. 2002). Thus, our flux densities are consistent with those previously reported in the literature. These results are summarised in Table 3.

From these results, we construct the integrated SED of NGC 891 and fit the data with our one-component modified black-body model, shown in Fig. 5. For brevity, we focus our discussion on the SED fitting where we allow β to vary as a free parameter, whilst also presenting the results derived from fixing β = 1.8 in Table 4. We find a dust mass and temperature of log M_dust = 7.93 ± 0.05 M_⊙ and T_dust = 23.1 ± 1.2 K, with a β = 1.77 ± 0.17. The total dust mass we uncover (8.5 × 10⁷ M_⊙) is higher than those reported in most previous studies; past results have ranged from M_dust ~ 1.9 to 7 × 10⁷ M_⊙ (Alton et al. 2000; Popescu et al. 2004; Galliano et al. 2008). Our higher dust masses are likely due to some combination of the higher FIR fluxes measured by Herschel compared to the observations from previous FIR experiments and the different assumptions for parameters used in the dust models of these studies.

We briefly note that NGC 891 appears in the Planck HIFI Catalogue of Compact Sources (Planck Collaboration XXVIII 2014). The reported total flux densities measured at FIR/submm wavelengths, i.e. 350, 550 and 850 μm, using the photometry derived from Gaussian fitting, GAUFLUX, are F₃₅₀ = 85.5 ± 0.7 Jy, F₅₅₀ = 19.9 ± 0.3 Jy and F₈₅₀ = 5.6 ± 0.1 Jy. We find these flux values are consistent with those at the respective PACS/SPIRE wavebands. If we include these three data points in our SED fitting, we derive nearly equivalent dust properties as those stated earlier: log M_dust = 7.96 ± 0.07 M_⊙, T_dust = 22.6 ± 1.4 K, with a β = 1.80 ± 0.21 (see Fig. 5). Thus, the difference in the best-fit models is negligible.

As a further check for consistency with the literature, we attempt to reproduce results reported in Dupac et al. (2003a). They traced the FIR SED using observations from ISO, IRAS/HiRes and PRONAOS. The SEDs were constructed from the flux density in three circular apertures with 1.5 arcmin radii (see Fig. 1), where one is centred on the galaxy and other two placed at ±3 arcmin along the galaxy major axis i.e. north-east (NE) and south-west (SW). Their one-component modified blackbody fit to the SEDs yield dust temperatures of 23.5 K, 19.6 K and 18.1 K for the centre, NE and SW regions, respectively, with a range in β ≈ 1.41−1.96. We determine a similar dust temperature for the center region, but the NE and SW regions are warmer by approximately 1.5 to 2 K (see Fig. 5 and Table 4). Although our dust emissivities are slightly higher than Dupac et al. (2003a), ranging from β ≈ 1.76−2.10, they are just within the errors (less than ±0.30). Such agreement in the results is encouraging when one considers that the two datasets came from different instrumentation.

On the SW end of the disc, an asymmetric extension of the disc has been previously detected in ISOPHOT 170 and 200 μm images and found to correspond with the Hi disc (Popescu & Tuffs 2003). The feature was also observed in the SPIRE 250, 350 and 500 μm maps by Bianchi & Xilouris (2011). They calculated the SW extension contained M_dust ~ 6 × 10⁵ M_⊙, albeit with a factor of two error due to the difficulty in determining the dust temperature using only SPIRE flux densities. By adopting a Hi mass of 2.5 × 10⁸ M_⊙ in the SW extension and a dust-to-gas mass ratio of 0.006 (similar to the Milky Way and other galaxies, e.g., Draine et al. 2007), Bianchi & Xilouris (2011) further estimate a total dust mass of 1.5 × 10⁶ M_⊙ for the extension. Combining the new PACS 70, 100 and 160 μm maps with the existing SPIRE observations, we can more accurately trace the peak of the integrated SED of the SW extension. We defined an aperture centred on , (J2000.0) that included all the pixels with detections greater than 3σ in the 500 μm map (see Fig. 1 in Bianchi & Xilouris 2011). We detect the SW extension in the PACS 160 μm map, however there is no detection of the feature in the 70 and 100 μm maps. Thus, we only include the PACS/SPIRE 160 to 500 μm data in the SED for this feature. The integrated SED and best-fit one-component model of the SW extension is shown in Fig. 5. We obtain a dust mass of 8.5 × 10⁵ M_⊙, i.e. 1% of our total dust mass estimate, lying between the two estimates of Bianchi & Xilouris and with an improved error of 15%. However, the dust temperature and emissivity index have larger errors than typically produced from our method, most likely due to the existing difficulty in estimating the dust temperature without the constraint afforded by the 70 and 100 μm data. Yet, our results demonstrate consistency with the previous study.

Having tested our methodology via the reproduction of results from previous studies, we now focus on the central goal of this project: mapping the dust mass and temperature.

4.5. Resolved dust analysis

Using the PACS, SPIRE and SCUBA images, we fit the SEDs of each pixel in the image cube and construct maps of the dust properties. We only consider pixels with a detection greater than 5σ in all maps used for the SED fitting, yielding dust surface density and temperature estimates for 192 pixels. In Fig. 6, we compare the dust surface densities and temperatures derived from fitting the pixel SEDs using a modified blackbody with β varying as a free parameter, against the results we obtain from fixing β = 1.8. For the derived dust masses, allowing β to vary as a free parameter typically yields masses around 15 to 25% higher than if we keep β fixed; the median dust mass per pixel is 3.4 × 10⁵ M_⊙ for a varying β and 2.5 × 10⁵ M_⊙ for β = 1.8 (corresponding to Σ_dust= 1.1 and 0.8 M_⊙ pc^-2, respectively). The dust mass surface density shows good agreement between the two fitting methods, being offset from the 1:1 relationship with a mean difference of 0.32 M_⊙ pc^-2 and a small scatter of σ = 0.20 M_⊙ pc^-2. The corresponding dust temperatures show similar agreement. Allowing β to vary as a free parameter yields temperatures around 10% lower than if we keep β fixed. The median dust temperature is 21.7 K for a varying β and 23.2 K for β = 1.8. The scatter in the temperatures is greater than for dust masses, with a mean difference of 1.1 K and a σ = 0.8 K. We examined the effect of varying the fixed value of the spectral index from β = 1.8 to β = 2.0, finding that the best agreement between the fixed and free β values occurred at β = 1.9. In fact, this result merely reflects that β = 1.9 is the average value obtained when β is allowed to vary as a free parameter (see also Fig. 8). Thus, we find consistent results whether or not β is a free or fixed value.

Fig. 6 Comparison of the dust masses (upper panel) and dust temperatures (lower panel) obtained from allowing β to vary versus fixing the β value. The 1:1 relationship and our linear best fit relationship are indicated by the black and red lines, respectively.

4.5.1. The T_dust–β relation

The fitted dust emissivity index is a parameter that describes the slope and width of the blackbody curve and might, therefore, be related to the physical properties of the dust grains, including composition, grain size and equilibrium temperature. By fitting the observed FIR/submm SEDs of various sources with modified blackbody spectra, Dupac et al. (2001, 2002, 2003b) suggested that β decreases with increasing dust temperatures. They found β ~ 2 in cooler regions (T_dust ≈ 20 K) and β ~ 0.8−1.6 for warmer environments (T_dust ≈ 25−80 K). This anticorrelation has since been observed in numerous experiments probing FIR/submm SEDs. ARCHEOPS data showed a stronger inverse relationship with β ranging from 4 to 1 between T_dust ≈ 7 and 27 K (Désert et al. 2008), and a similar relation between T_dust and β was found in BOOMERanG observations of eight clouds at high Galactic latitude (Veneziani et al. 2010). Recent studies using Herschel observations of the inner regions of the Galactic plane (Paradis et al. 2010) and Galactic cirrus emission (Bracco et al. 2011) also find an inverse relationship.

There is some debate about the nature of such an inverse T_dust–β relationship, which may not have a physical origin, but rather arises from the data and/or methods used to derive the quantities (see e.g., Shetty et al. 2009a; Foyle et al. 2012; Juvela et al. 2013). An inverse correlation between T_dust and β has been shown to arise from least-squares fits due to the uncertainties in data (Shetty et al. 2009a) and line-of-sight temperature variations (Shetty et al. 2009b), since a one-component modified blackbody may not be appropriate for modelling the emission from dust that has a range of temperatures. Such line-of-sight temperature variations may lead to an underestimation of β (e.g., Malinen et al. 2011; Juvela & Ysard 2012). Furthermore, spurious temperature variations may be introduced, such as peaks in the temperature distribution where there are no identifiable heating sources (Galametz et al. 2012).

Yet in a pixel-by-pixel analysis of dust and gas in the Andromeda galaxy, Smith et al. (2012b) observed a radial dependency on the T_dust–β relationship; pixels at radii R > 3.1 kpc produce a steeper relation compared to pixels within the radius R< 3.1 kpc. They argue that neither artefacts in the fitting due to noise nor the need to account for line-of-sight temperature variations can explain this clear separation between the relations from the inner and outer regions of Andromeda. Furthermore, the relation from Smith et al. for the outer region (R > 3.1 kpc) is in good agreement with the T_dust–β relationship observed in the global dust properties derived for galaxies in the Herschel Virgo Cluster Survey (Davies et al. 2014). These results suggest a possible physical origin of the relationship.

Plotting our results (see Fig. 7), we observe a similar relationship between dust emissivity index and temperature. A weak negative correlation exists between the two quantities with a Spearman coefficient of rank correlation of ρ = −0.52, corresponding to a probability P(ρ) > 99.9% that the two variables are correlated. Adopting the same parametrisation for the relation as in previous works (e.g., Smith et al. 2012b), we find a best fit relationship of (3)The distribution of the points are fairly consistent with previous studies, lying offset from the Dupac et al. (2003b) relation and between the two radial relations of Smith et al. (2012b), although their best-fit relation is steeper at lower temperatures compared to our best fit.

As NGC 891 is a nearly fully edge-on spiral (i> 89°, e.g., Kregel & van der Kruit 2005), we cannot investigate the existence in NGC 891 of a radial break in the T_dust–β relation as observed in Andromeda. We attempt to exploit the inclination of NGC 891 by examining any dependence of the T_dust–β relation on the vertical height from the mid-plane. We roughly divide the pixels covering the galaxy into two bins: one containing all pixels with a perpendicular distance from the mid-plane greater than ±1 kpc, and one containing all pixels located within ±1 kpc from the mid-plane. This distance from the disc was chosen based on the extent of the molecular gas in the disc, which has been shown to have a ±1 kpc vertical extension. We remind the reader that at a distance of 9.6 Mpc each 12″ pixel corresponds to approximately 558 pc and, due to the PSF of the SPIRE 500 μm image of 36′′, are not independent. We do not find a significant offset in the mean β of these two regions. The T_dust–β relationships also do not show significant variance between the two bins. This is evident when one inspects the dust temperature and spectral index maps presented in Fig. 8. The temperature distribution is clumpy and asymmetric about the mid-plane, whereas β has a much smoother distribution.

Thus, we do not find a statistical variation of the T_dust–β relation on the vertical height from the disc, suggesting that any change in the T_dust–β relation is predominantly linked to radial variations from the interior to the exterior of a galaxy (Smith et al. 2012b). However, since we are fitting a one-component modified blackbody to dust emission integrated along the line-of-sight through the galaxy disc, covering a range of environments and temperatures, we stress that such a physical interpretation for the variation of the relationship should be regarded with caution.

Fig. 7 Anticorrelation we observe between β and dust temperature from fitting a one-component modified blackbody to the pixel-by-pixel SEDs constructed from our FIR/submm images.

4.5.2. Spatial distributions

Given the plethora of possible problems and caveats that must be considered when attempting to interpret the fitted dust emissivity index, many of which we have discussed above, and, since we find consistent dust properties whether or not β is a free variable or a fixed value, we decide to adopt the dust properties derived from fixing β = 1.8 for the remainder of our analysis.

The resulting maps are presented in Fig. 8. The mean dust mass surface density is 1.0 M_⊙ pc^-2, which peaks at 2.6 M_⊙ pc^-2 in the centre. Either side of the centre and along the semimajor axis are two smaller local maxima that coincide with the surface brightness peaks seen in the WISE, PACS, SPIRE and SCUBA images (see Fig. 1). The distribution shows a fairly smooth decrease away from the galaxy mid-plane. However, the dust temperature has a much more uneven distribution. Ranging from ≈17 to 24 K with a median of 22.9 K, the temperature peaks in the centre and again in a region on the NE end of the disc. The dust in this region appears warmer than the average dust temperature on the opposite end of the disc i.e., towards the SW, an observation which is most prominent in the averaged radial profile of the dust temperature map (see Fig. 3). In fact, the dust temperature profile reflects the same features and asymmetry as seen in the radial profile of the 24 μm emission. This asymmetry in the disc has been previously noted in the literature (e.g., Kamphuis et al. 2007), and we shall discuss this behaviour in detail in Sect. 6.2. We briefly note that the temperature profile exhibits the same break as the FIR/submm profiles, indicating that the break in the radial profiles is most likely due to a change in the ISM properties at radii ±12 kpc. Returning to Fig. 8, we also show the χ² value used to fit the SED of each pixel. The majority of the pixel SEDs have χ² values <1.5. Those pixels with χ²> 1.5 are typically at the outskirts of the disc, where the flux density from the SCUBA 850 μm are less certain with higher errors, as confirmed from a visual inspection of the individual SEDs for these pixels.

Furthermore, we investigated whether the low resolution of the SPIRE 500 μm data (36″) affected our results by omitting this data and performing the same pixel-by-pixel SED fitting at the resolution of the SPIRE 350 μm data (24″). We find that increasing the resolution yields the same distribution of dust properties, albeit with slightly higher mass surface densities and lower temperatures (~6%, i.e. much smaller than the differences seen for treating β as a fixed or variable quantity). However, the SEDs of some individual pixels are less constrained on the Rayleigh-Jeans side, mainly due to the sensitivity of the fitting to the 850 μm flux with respect to the overall FIR SED (see Fig. 5). Since the pixel-by-pixel fits to the SED are better constrained when including rather than excluding the 500 μm data and the enhanced resolution neither improves nor affects our results, we continue the analysis including the SPIRE 500 μm data at the lower resolution.

Fig. 8 Maps of the dust mass surface density (left), dust temperature (middle) and χ2 value (right), obtained from fitting a one-component, modified blackbody model with β = 1.8 to the SED of each pixel. The maps are centred on , (J2000.0), with the beam size indicated by the black circle. The pixel size is 12″ pix-1 corresponding to 0.56 kpc at a distance of 9.6 Mpc.

5. Gas and dust

Having determined the dust properties of NGC 891 in the previous section, we now combine this new information with measurements of the gas content and examine the relationships between gas and dust in the disc.

5.1. Estimating gas content

We first calculate the masses and corresponding surface densities of the various gas components. We base our analysis primarily on the pixels covered by our CO(3−2) map, since it covers a much larger fraction of the galaxy compared to the available Scoville et al. (1993) CO(1−0) emission map (see Fig. 1). We include the latter CO(1−0) observations only as an additional test of the results of our analysis.

5.1.1. CO(3−2) observations

We derive a molecular gas mass estimate for NGC 891 from the new map of ¹²CO(J = 3−2) emission. For all detected pixels, the integrated CO(3−2) line intensity is converted to a H₂ mass, M_H2, according to the equation (4)where I_CO(3−2) is the total integrated line intensity expressed in units of K km s^-1, A represents the surface of the CO(3−2) emitting region and m_H2 is the mass of a hydrogen molecule. The scaling factor to convert an antenna temperature T_A into a main beam temperature T_mb at the JCMT is η_mb = 0.6. We assume a value of 0.3 for the CO(3−2)-to-CO(1−0) line intensity ratio corresponding to the typical ratios found in the diffuse ISM of other nearby galaxies (Wilson et al. 2009). The CO-to-H₂ conversion factor, X_CO, has been difficult to constrain and varies significantly over a range of environments and metallicities (see Bolatto et al. 2013, and references therein). Although we lack the observations to determine the metallicity gradient in NGC 891, we refrain from using an average metallicity gradient derived from observations of other nearby face-on spiral galaxies to estimate the radial variation in the X_CO factor (similar to an approach used in Alton et al. 2000) in order to avoid (i) the large uncertainties in applying such a gradient; and (ii) the uncertainties in the choice of radial relationship between X_CO and metallicity (see e.g. Fig. 9 in Bolatto et al. 2013). We therefore choose to adopt X_CO = 2 × 10²⁰ cm^-2 [K km s^-1]^-1, based on the recommendation of Bolatto et al. (2013) as a conservative choice for the discs of normal, solar metallicity galaxies and the Milky Way (see also Sandstrom et al. 2013). Our adopted values result in a total H₂ mass of 2.9 × 10⁹ M_⊙ contained in the observed area. Lower estimates were obtained by Guelin et al. (1993, 1.5 × 10⁹ M_⊙ and Israel et al. (1999, 1.4 × 10⁹ M_⊙. Our corresponding mean H₂ mass surface density is Σ_H2 = 67 M_⊙ pc^-2, where Σ_H2 ranges from 11 to 144 M_⊙ pc^-2.

5.1.2. Hi observations

We estimate the Hi surface density using the total column density map from Oosterloo et al. (2007), finding a total integrated Hi mass of M_Hi = 4.2 × 10⁹ M_⊙ at our adopted distance of 9.6 Mpc. Previous estimates have ranged from 2.5−8 × 10⁹ M_⊙ (Sancisi & Allen 1979; Rupen 1991; Guelin et al. 1993). It has been shown that up to 30% of the total mass is contained in a large gaseous halo extending ~22 kpc from the disc (Swaters et al. 1997; Oosterloo et al. 2007). We calculate the disc contains 2.5 × 10⁹ M_⊙. The median Hi mass surface density of the 12″ pixels is Σ_Hi = 42 M_⊙ pc^-2 and ranges from 22 to 58 M_⊙ pc^-2 in the disc. The errors are estimated at 10%. Furthermore, we note that we have not accounted for the effects of Hi self-shielding. Whilst we adopt the usual assumption of negligible self-opacity for the transition, Braun et al. (2009) warn that this assumption may overlook the presence of hidden Hi features that are self-opaque in the 21 cm transition and thus underestimate the global atomic gas mass by up to 34% ± 5% (see also Braun 2012). However, any attempt of applying local opacity corrections to our column density map is not straightforward because (i) line-of-sight confusion will tend to overlap discrete features in both position and velocity; and (ii) the physical scales probed in our analysis (~500 pc) are considerably larger than the physical scales of any hidden Hi features (~50−100 pc). We therefore take a conservative approach and refrain from applying this correction, yet caution that our Hi surface density estimates may be underestimated by ~30% and continue our analysis with this caveat in mind.

5.1.3. Total gas content

From these measurements of the atomic and molecular hydrogen masses, we estimate the total gas mass in each pixel as M_gas = M_Hi + M_H2 and include the contribution of helium and heavy elements (i.e., a factor 1.36). Thus, our integrated gas mass for NGC 891 is M_gas = 1.1 × 10¹⁰ M_⊙. For comparison, previous studies found total gas masses ranging from M_gas ≈ 0.3−1.2 × 10¹⁰ M_⊙ (see e.g., Guelin et al. 1993; Israel et al. 1999; Dupac et al. 2003a). We briefly note that our integrated dust, H₂ and Hi masses place NGC 891 very close to the best-fit M_H2–M_dust and M_gas–M_dust relations derived from late-type Virgo cluster galaxies (see Corbelli et al. 2012). The total gas content, Σ_gas, has a mean surface density of 149 M_⊙ pc^-2 and ranges from 56 to 267 M_⊙ pc^-2.

Fig. 9 Pixel-by-pixel relationships between the derived dust properties, Σdust (upper panels) and Tdust (lower panels), and the surface densities of the gaseous ISM components: ΣHi (left), ΣH2 (middle) and Σgas (right). We differentiate between on- and off-axis pixels using solid and open circles, respectively. The best linear fits are indicated by the solid red lines. Error bars are calculated from the standard deviations of the errors on individual pixels and the ρ values are the corresponding Spearman correlation coefficients.

5.2. The dust–gas connection

We now combine all our surface density estimates to study the relationships between dust and gas in the galaxy. We stress that all values related to the gas content reported henceforth in our analysis are derived from H₂ masses from the CO(3−2) map.

We find several strong correlations between the pixel-to-pixel dust and gas surface densities. In Fig. 9, the derived dust surface density and temperature for each pixel are plotted against the corresponding Hi, H₂ and total gas surface densities. As expected, we find only a very weak correlation between Σ_Hi and Σ_dust, as evident both in Fig. 9 and indicated by a low Spearman coefficient of rank correlation (ρ = 0.25). We remind the reader that this is only true for the inner part of the galaxy. However, the clear correlation between Σ_H2 and Σ_dust has a Spearman coefficient of ρ = 0.94, corresponding to a probability P(ρ) > 99.9% that the two variables are correlated. Interestingly, we find a similar correlation when combining the Hi and H₂ gas content; the correlation between Σ_gas and Σ_dust has a Spearman coefficient of ρ = 0.93, with a corresponding P(ρ) > 99.9%. Both the Σ_H2–Σ_dust and Σ_gas–Σ_dust correlations are best fit with the linear relations where all quantities are in units of M_⊙ pc^-2. Such correlations between Σ_dust and Σ_H2 are predicted from the physics of CO formation and dissociation and have already been observed in a number of environments. For example, Roman-Duval et al. (2010) found linear correlations between Σ_gas and Σ_dust on even smaller spatial scales, via Herschel observations of two molecular clouds in the LMC (see also Meixner et al. 2010; Leroy et al. 2011; Galliano et al. 2011, see their Fig. 12). The Σ_gas–Σ_dust correlations have been reproduced in the simulated disc galaxies of Bekki (2013), using star formation histories regulated by the time evolution of interstellar dust, the properties of which control H₂ formation rates. The dust surface density in the simulated discs correlated with Σ_gas than Σ_H2, consistent with our observations.

Fig. 10 Maps of the Hi-to-dust ratio (left), H2-to-dust ratio (middle) and the total gas-to-dust ratio (right). The maps are centred on , (J2000.0), with the beam size indicated by the black circle. The pixel size is 12″ pix-1 corresponding to 0.56 kpc at a distance of 9.6 Mpc.

An immediate observation from these Σ_H2–Σ_dust and Σ_gas–Σ_dust correlations is that they appear to become increasingly bifurcated with increasing surface density. Such bifurcation could indicate problems with offsets between images, misaligned position angles or a mismatching of the PSFs adopted for the image convolution. However, our investigations into the effects of misalignments between the images, by systematically introducing artificial offsets in the position angles and world coordinate systems of the images, failed to confirm that any of these issues were causing the bifurcation. In fact, these artificial misalignments tended to drastically increase the scatter of the observed correlations and thereby hide the bifurcation. The minimum scatter in the correlations actually occurs when adopting the position angle of 22.9°, consistent with our conclusions when extracting the vertical and radial profiles (see Sect. 3). Taking the best fit lines to the Σ_gas–Σ_dust correlations (given by Eqs. (5) and (6)) as references, we note that the highest Σ_H2 and Σ_gas pixels have lower Σ_dust values than these best fit lines would suggest, whereas the pixels with the highest Σ_dust tend to have lower Σ_H2 and Σ_gas values. In other words, there is a clear difference in the gas-to-dust ratios and these results imply that the gas-to-dust ratios in the galaxy center are much lower than the rest of the disc.

To further explore the origin of this birfucation, Fig. 10 presents the mapped ratios of the various gas surface densities to the dust surface density, i.e. the Hi-to-dust ratio (Σ_Hi/Σ_dust), H₂-to-dust ratio (Σ_H2/Σ_dust), and the total gas-to-dust ratio (Σ_gas/ Σ_dust). The Hi-to-dust ratio is lowest in the centre of the disc and increases with increasing radius, reflecting the observations that the dust surface density peaks in the centre whereas the Hi surface density dominates at larger radii. The lack of a clear radial trend lends further evidence to the lack of an overall relationship between Σ_Hi and Σ_dust in the disc (Fig. 9). In contrast, the H₂-to-dust ratio is fairly low (Σ_H2/Σ_dust~ 50) at the disc centre, as predicted from the bifurcation of the Σ_gas–Σ_dust correlations, but peaks in two regions towards the NE and SW of the disc. Most intriguing is the spatial distribution of the total gas-to-dust ratio. The mean Σ_gas/Σ_dust is ~140, slightly lower than the Galactic Σ_gas/Σ_dust value of 158 from Zubko et al. (2004). Most of the disc possesses a Σ_gas/Σ_dust ratio between 110 and 150, but there are two off-centre peaks of Σ_gas/Σ_dust≈ 160 situated on each side of the disc at approximately ±5 kpc from the center. These peaks also roughly coincides with the aforementioned peaks in the dust temperature distribution (see Fig. 8). We briefly note a good agreement between the spatial variation of the gas-to-dust ratio observed in our maps and the variations inferred from the radial gas-to-dust profiles presented in Alton et al. (2000).

Fig. 11 Correlation between the H2-to-dust ratio and the dust temperature. The best linear fits are indicated by the solid red lines, with the mean ΣH2/ Σdust ratio represented by the dashed black line. The cluster of points with ΣH2/ Σdust ratio lower than the mean and possessing a higher dust temperature (Tdust> 23.5 K, shown by the dotted black line) are all central pixels.

The dip in the gas-to-dust ratio in the galaxy center becomes particularly obvious when we consider the relationship between Σ_gas/Σ_dust and the dust temperature (see Fig. 11). We find a tight relationship between these quantities for most of the disc, but the central pixels show a clear separation from this relation due to their higher dust temperatures but on average lower Σ_gas/Σ_dust values. The exact cause of this is uncertain, but may be indicating that we are underestimating the H₂ surface density. Whilst the assumption of a single value for the CO-to-H₂ conversion factor may be a conservative choice for the discs of normal, solar metallicity galaxies and the Milky Way, and so an appropriate choice for the majority of the disc of NGC 891, in this case it may not be appropriate for the nuclear regions. Though the X_CO factor is observed to drop sharply in the central, bright regions of some galaxies, often coincident with bright CO emission and high stellar surface density (see e.g. Bolatto et al. 2013; Sandstrom et al. 2013), our choice of X_CO factor may still be underestimating the Σ_H2 in the center. Another plausible scenario is that the CO(3−2) emission fails to trace the full molecular content in the densest central regions, particularly when integrating along the line-of-sight through the nucleus of an edge-on galaxy, as the CO(3−2) line could become optically thick in the centre of galaxies. The CO(3−2) line becomes optically thick close to a molecular hydrogen column density N_H2 = 2 × 10⁵ cm^-3, which corresponds to a cloud with N_CO > 2 × 10¹⁷ cm^-2 or Tk> 550 K (see e.g. Imai et al. 2012). Furthermore, the geometry of the molecular gas may play a role, as a clumpy distribution of molecular clouds will be locally optically thick and may therefore inefficiently contribute to the CO emission. Undetected CO emission in denser clumps may also explain why Sakamoto et al. (1997) tend to find moderate H₂ densities from the observed CO emission.

Whilst the H₂-to-dust and gas-to-dust ratios show some spatial correlation with the dust temperature, evident from the coincidence of the peaks in their distributions (see also Fig. 11), the gas surface densities typically show weaker correlations with dust temperature compared to the dust surface density (see Fig. 9). Weaker correlations are found between dust temperature and the H₂ and total gas surface densities than compared to the similar relations with dust surface density. The observed correlation between T_dust and Σ_H2 has a Spearman coefficient of rank correlation of ρ = 0.84, corresponding to a probability of correlation of P(ρ) > 99.9%. The T_dust–Σ_gas relation has a Spearman coefficient of ρ = 0.76, also with P(ρ) > 99.9%. The best-fit linear solutions are given by where Σ_H2 and Σ_gas are in units of M_⊙ pc^-2 and T_dust is in K. Although we fit linear relations to the correlations, it is evident that such linear relationships are not really appropriate. For pixels containing Σ_H2< 100 M_⊙ pc^-2, the plotted data appears to follow a linear correlation. This also seems to be reflected in the T_dust–Σ_gas relation for pixels Σ_gas< 200 M_⊙ pc^-2. However, it is not clear at present what causes the pixels at Σ_H2 > 100 M_⊙ pc^-2 (corresponding to Σ_gas > 200 M_⊙ pc^-2) to flatten out towards a peak temperature, rather than continue to follow the linear trend found for the rest of the disc. One potential clue may be gained from the location of these pixels in the galaxy.

Since we have shown that (1) the dust mass surface density typically decreases as a function of radius and | z | from the galaxy centre; (2) the molecular hydrogen surface density distribution is strongly correlated with the dust density distribution; and (3) larger hydrogen surface densities appear to be associated with cooler dust than expected when compared to the dust temperatures for the rest of the galaxy, then it stands to reason that those pixels containing this slightly cooler dust should be located towards the centre of the galaxy. We test this hypothesis, highlighting the location on the temperature map of all pixels possessing H₂ surface densities greater than 100 M_⊙ pc^-2 (see Fig. 12). The cooler-than-expected dust is clearly associated with two regions either side of the nucleus.

Fig. 12 Map of the dust temperature obtained from fitting a one-component, modified blackbody model to the SED of each pixel (colour scale, see Fig. 8). Pixels containing dense regions of molecular hydrogen and cooler-than-average dust are highlighted (purple contours). The WISE F12/F22 flux ratio (black contours) follows the dust temperature.

Finally, we also explored these relations between the dust and gas properties using the Scoville et al. (1993) CO(1−0) emission map⁵, performing the exact same analysis as described above. Although the CO(1−0) map covers a smaller fraction of the galaxy in the vertical direction compared to the CO(3−2) map (see Fig. 1), we were able to recover the same qualitative correlations between the dust properties and gas components as presented above. In fact, the break in the flattening of the T_dust–Σ_gas relation is even more prominent. Our CO(1−0)-based relations only deviated quantitatively from the CO(3−2)-based relations due to a small offset between the two estimates of Σ_H2 and the fewer number of pixels available for the analysis, which were both expected. Also, due to the lower number of pixels in CO(1−0)-based results, it is ambiguous whether we see a similar bifurcation of the Σ_H2–Σ_dust and Σ_gas–Σ_dust correlations. Combining all these observations, an interesting picture emerges.

6. Discussion

In summary, whilst we find only a weak correlation between the Hi gas content and the dust surface density, we find clear correlations between the dust mass surface density and temperature with the H₂ and total gas surface densities. These results are interesting for a number of reasons. Firstly, although the CO lines can become optically thick for an edge-on geometry, the strong spatial correlation between the molecular gas surface density and surface density of dust in the majority of the disc of NGC 891 is also observed in face-on galaxies (e.g., Smith et al. 2010; Foyle et al. 2012; Mentuch Cooper et al. 2012), indicating that CO emission can also be used as a tracer of the molecular gas distribution in edge-on galaxies. However, our observations hint at the possibility that the CO(3−2) line may become optically thick in the central regions where the column density along the line-of-sight is greatest. Secondly, our gas-to-dust ratios seem to show strong variation in the radial direction along the disc. Similar analyses of the gas and dust using Herschel observations of M83 and M51 yielded little or no variation of the gas-to-dust ratio in these face-on systems (Foyle et al. 2012; Mentuch Cooper et al. 2012), whereas many other galaxy studies have found strong gas-to-dust ratio gradients (e.g., Muñoz-Mateos et al. 2009; Bendo et al. 2010b; Magrini et al. 2011; Fritz et al. 2012). Whilst we have yet to fully understand the effects of galaxy structure (e.g. the presence of bars, variations in spirality, etc.) on the gas-to-dust ratio gradients, these features may account for some of the observed radial variations between galaxies. For example, gas tends to be more centrally concentrated in barred galaxies than non-barred galaxies (Sakamoto et al. 1999), possibly due to the presence of a bar enhancing gas flows into the central regions (e.g., Tabatabaei et al. 2013). Finally, the deviations of individual pixel-by-pixel quantities from these overall correlations in the disc are demonstrably indicative of important morphological features of NGC 891, namely a ring-like feature identified from the relationships between the dust properties and Σ_H2, and evidence of disc asymmetry seen in the spatial distribution of the gas-to-dust ratio. For the remainder of this discussion, we focus our attention on the implications of our results in our understanding of these features of NGC 891.

Fig. 13 Pixel-by-pixel relationships between the derived dust properties, Σdust (upper panels) and Tdust (lower panels), and the IRAC 3.6 μm (left panels) and MIPS 24 μm flux densities (right panels). The best linear fits are indicated by the solid red lines. Reported ρ values are the corresponding Pearson correlation coefficients.

6.1. On the molecular ring

One feature observed in the H₂-to-dust ratio and cold dust temperature maps suggest the presence of two regions of dense, relatively warmer dust and molecular hydrogen gas located either side of the disc. The regions occupy radial distances between ~2 and 5 kpc, which appear almost symmetrical in radial size. Such a feature has already been identified by Israel et al. (1999). In their analysis of the CO and dust emission, they describe the disc structure of NGC 891 via three components comprising a very compact central source/circumnuclear disc (see also Garcia-Burillo et al. 1992; Scoville et al. 1993), a “molecular ring” between R = 40′′ and 120′′ (~2−6 kpc) and an extended disc detectable out to R = 200′′ (~9 kpc). This structure in the CO(J = 1−0) emission was also noted by Israel et al. to be tightly correlated with the dust traced via the SCUBA 850 μm emission (see their Fig. 2). The good correlation we observe between the dust surface densities and molecular hydrogen content is thus consistent with these previous results.

We find that the surface densities of dust and H₂ are typically high in the molecular ring compared to the rest of the disc distribution. The gas-to-dust ratio clearly peaks on both the NE and SW side of the disc (~3−6 kpc) as the dust surface density declines faster than the H₂ surface density with increasing radius. The peaks in the H₂ surface density also correlate spatially with these dust temperature peaks, most likely a consequence of the Kennicutt-Schmidt law relating the H₂ surface density to the star formation rate (e.g., Schmidt 1959; Kennicutt 1998), where star formation subsequently heats the dust. This interpretation seems consistent with the fact our derived dust temperature correlates strongly with the star formation traced via the 24 μm emission (see Fig. 13, and also the discussion in Sect. 4.2). Whilst the molecular ring and the outer extended disc appear to follow linear relationships, a flattening of the relation between the dust temperature and the H₂ surface density occurs within R = ± 40′′ (~2−3 kpc). Our T_dust–Σ_H2 diagram suggest that two regions of dense, relatively cooler dust and molecular hydrogen gas reservoirs are located on either side of the galaxy center, corresponding roughly to the circumnuclear disc described in Israel et al. (1999). These regions are on average cooler than the temperature we may predict from the overall correlation between T_dust and Σ_H2. The two regions of cooler dust are also evident in the dust temperature map, lying either side of the central temperature peak and between the NE and SW temperature peaks that are coincident with the Σ_H2 surface density peak. Of course, the T_dust–Σ_H2 relationship is likely more complex than the linear relationship assumed in this work. Nevertheless, this raises the question of why the dust temperature is on average colder within the inner radius of the molecular ring compared with the rest of the disc.

FIR observations of giant molecular clouds and filaments have demonstrated that dust temperatures are typically lower in dense molecular clouds than in diffuse regions (see e.g., Lagache et al. 1998; Stepnik et al. 2003). The dust and H₂ are physically associated with star-forming regions and, in dense molecular clouds, dust shields both the H₂ and CO from dissociation by the ambient interstellar radiation field (ISRF). One explanation for the cooler dust may be that the dust grains in molecular clouds are also opaque to the UV ISRF and so shield the dust embedded deeper within the clouds. An alternative to the possibility of dust shielding in molecular clouds is that there is less dust heating in these regions from star formation. Since both the UV, Hα (see e.g., Kamphuis et al. 2007) and 24 μm maps (Fig. 1) show a similar dip in emission from the central peak to the two NE and SW maxima (see also the radial profiles in Fig. 3), indicating a relative drop in recent SF between these regions, it is possible that the dust is cooler due to a lack of dust heating compared to the rest of the star-forming disc, particularly in the molecular ring. This is perhaps due to lower star formation rates, or even the aforementioned processes of either self-shielding or molecular gas shielding from heating by the UV interstellar radiation field.

Perhaps another possible albeit speculative explanation may lie in the process of dust coagulation, where large aggregates are formed with irregular and “fluffy” shapes (see e.g., Ossenkopf 1993; Köhler et al. 2012, and references therein). The submillimetre emissivity increases with “fluffiness” (e.g., Stognienko et al. 1995), whereas the UV to NIR absorptivity remains constant (e.g., Bazell & Dwek 1990), resulting in fluffy aggregates typically having lower equilibrium temperatures compared to compact dust grains (e.g., Fogel & Leung 1998). Grain-grain coagulation has been invoked to explain the significant dust temperature variations observed in filaments in the Taurus molecular cloud complex (Stepnik et al. 2003; Ysard et al. 2013; see also Paradis et al. 2009). However, given the effects of line-of-sight integration plus the large spatial scales (~0.5 kpc) our observations are probing, we cannot further investigate this scenario. We can only speculate that we are seeing the average effect of dust coagulation or some other process, lowering the observed dust temperature in many unresolved molecular clouds within these regions studied at low resolution. The physical mechanism causing the relationship between the dust temperature and the surface density of molecular hydrogen and total gas remains unknown. Detailed modelling combined with higher resolution data are required to further examine the underlying nature of the correlation.

6.2. On the disc asymmetry

As we previously mentioned in the introduction, NGC 891 displays a NE-SW asymmetry in the star-forming disc. The NE side has more prominent and extended Hα and UV emission than the SW side of the disc (Dettmar 1990; Rand et al. 1990; Kamphuis et al. 2007). Two possible explanations for the nature of this asymmetry have been proposed. Rossa et al. (2004) interpret the NE-SW asymmetry as being due to a higher SFR in the northern part of the disc than in the southern part. Kamphuis et al. (2007) found that star formation tracers affected by dust attenuation (Hα, UV) show greater asymmetry compared to those tracers unaffected by dust attenuation (24 μm, radio continuum). They argue that since the old stellar population and Hi gas distribution are fairly symmetric, the asymmetry in Hα is most likely caused by attenuation by dust in and above the plane. However, the fact that the small asymmetry is also seen in the radio continuum observations (Dahlem et al. 1994), which are not affected by dust attenuation, seems to support the former interpretation.

In order to better understand whether the asymmetry arises due to higher rates of star formation in the north or greater dust obscuration in the south, it is crucial to study the dust surface density and temperature distributions. Firstly, we noted that in the extended disc, beyond the radial extent of the “molecular ring” (i.e. R> 6 kpc), the dust surface density distribution is fairly symmetric between the NE and SW ends of the disc, i.e., both sides contain similar quantities of dust grains available to obscure star-forming regions. Thus, it is unlikely that there is an increase in dust obscuration due to a mere enhancement of the dust distribution in the SW end. However, we cannot completely rule out the scenario due to the presence of a small asymmetry in the distribution. The temperature of the dust in the NE end of the disc is on average warmer (by ~2−3 K) than the dust at corresponding radii on the SW end of the disc (see the dust temperature profile in Fig. 3). This warmer region in the NE is coincident with the peaks in the emission of the Hα and UV images (Kamphuis et al. 2007), and also the peak in the ratio map of 24 μm to 850 μm emission (see Fig. 7 in Whaley et al. 2009; see also our Fig. 13). The latter ratio traces the relative contribution of warm dust associated with star formation and the emission from cold dust. Furthermore, we find that the ISM emission as traced by the WISE F₁₂/F₂₂ flux ratio, i.e. the ratio of the emission from polycyclic aromatic hydrocarbons (PAHs) and dust warmed via the UV radiation field and the old stellar population, is also asymmetric (see the black contours in Fig. 12). The contours tracing the ISM are clearly more radially extended on the NE side than the SW side of the disc, and show some overlap with the peak in the cold dust temperature map.

The combined observational evidence suggests a scenario where the warm and cold dust components are being heated by ongoing or recent (t< 10⁷ yr) star formation in the disc, with the asymmetry at various wavelengths arising from an asymmetry in the SFR. The observed asymmetries in the H₂-to-dust ratio and the total gas-to-dust ratio hint that an enhancement in the SFR may be the result of larger quantities of molecular gas available to fuel star formation in the NE compared to the SW. Furthermore, the substantial Hi plume found by Oosterloo et al. (2007) suggests the possibility of a continual fuelling of star formation on the NE side, by either the accretion of cooling gas in the halo (Fraternali & Binney 2008) or the absorption of a small satellite (e.g. Mouhcine et al. 2010). Recent observations demonstrating solar to supersolar metallicity gas in the halo at 5 kpc from the disc suggests that this halo gas most likely originates from the disc rather than from the accretion of material (Bregman et al. 2013). However, the question of whether gas accretion is fuelling ongoing star formation remains.

Finally, we should also consider the insights afforded by the various attempts to interpret NGC 891 using radiative transfer models (e.g., Kylafis & Bahcall 1987). Smooth, axisymmetric models have so far proved insufficient for accurately fitting the observations of NGC 891. For example, Xilouris et al. (1998) found it necessary to take the average optical emission of the northern and southern halves and use an infinitely-long thin disc of stellar emission to reconcile the two sides of NGC 891, whereas smooth models of NGC 891’s MIR emission have required extra dust and stellar components to fit the data (Popescu et al. 2000; Bianchi 2008). Most recently, Schechtman-Rook et al. (2012) fitted three-dimensional radiative transfer models including dust clumping and a realistic prescription for spirality to Hubble Space Telescope images of NGC 891. Despite some unavoidable degeneracies in their model parameters, the best-fit models determined from a genetic algorithm demonstrated a clear preference for models including spiral structure and clumpy dust (see also the models by Kamphuis et al. 2007; Popescu et al. 2011). Taking these results into account, we suspect that the asymmetry likely arises from dust obscuration due to the geometry of the line-of-sight projection of the spiral arms, but we cannot exclude that there is also an enhancement in the star formation rate in the NE part of the disc.

7. Conclusions

In this work, we investigated the connection between dust and gas in the nearby spiral galaxy NGC 891. High resolution Herschel PACS and SPIRE 70, 100, 160, 250, 350, and 500 μm images are combined with SCUBA 850 μm observations to trace the far-infrared/submillimetre SED. We fit one-component modified blackbody models to the integrated SED, finding a global dust mass of (8.5 ± 2.0) × 10⁷ M_⊙ and temperature of 23 ± 2 K, consistent with results from previous far-infrared experiments. We also fit one-component modified blackbody models to pixel-by-pixel SEDs to produce maps of the dust mass and temperature at physical scales of ~0.5 kpc. The dust mass distribution shows evidence of several peaks along the disc. The derived dust temperature, which ranges from approximately 17 to 24 K, does not correlate with the dust spatial distribution. Allowing the dust emissivity index to vary, we find an average value of β = 1.9 ± 0.3. We confirm an inverse relation between the dust emissivity spectral index and dust temperature, but do not observe any variation of this relationship with vertical height from the mid-plane of the disc. We adopt a fixed value of β = 1.8 for our main analysis comparing the derived dust properties to the gaseous components.

Combining the mapped dust properties with observations of the Hi 21 cm line emission and CO(J = 3 − 2) emission to trace the atomic and molecular hydrogen gas, respectively, we examined the correlations between the dust properties with the gaseous components of the ISM. We found strong spatial correlations between the surface mass densities of dust (Σ_dust) and the molecular hydrogen (Σ_H2) and total gas surface densities (Σ_gas). We observed no correlation between the dust content and the Hi gas content. These observations reveal the presence of regions of dense, cold dust that are coincident with peaks in the gas distribution and appear associated with a molecular ring. Furthermore, the observed asymmetries in the dust temperature, the H₂-to-dust ratio and the total gas-to-dust ratio hint that an enhancement in the star formation rate may be the result of larger quantities of molecular gas available to fuel star formation in the NE compared to the SW. Whilst we suspect that the asymmetry likely arises from dust obscuration due to the geometry of the line-of-sight projection of the spiral arms, we cannot exclude that there is also an enhancement in the star formation rate in the NE part of the disc.

This study raises several intriguing questions. The observed relationship between the dust temperature and the surface density of the molecular hydrogen suggests that dust may be cooler in dense regions of dust and molecular hydrogen, perhaps due to the processes of either self-shielding or molecular gas shielding from heating by the UV interstellar radiation field, or dust coagulation. However, whilst these processes have support from observations at high spatial resolution in nearby molecular clouds, it is unclear whether our lower resolution dataset is capable of uncovering such a correlation between the dust temperature and molecular gas density observed on these large spatial scales. This remains an open issue, which may be resolved by a detailed analysis and modelling of the T_dust–Σ_H2 relation based on larger samples of edge-on galaxies, e.g. HERschel Observations of Edge-on Spirals (HEROES, Verstappen et al. 2013) and the New HErschel Multi-wavelength Extragalactic Survey of Edge-on Spirals (NHEMESES, Holwerda et al. 2012b). Furthermore, it is still an open question whether molecule formation is mainly driven by hydrostatic pressure or UV radiation shielding over different spatial scales and over a large range of metallicities (see e.g., Fumagalli et al. 2010 and references therein). Upcoming facilities such as the Atacama Large Millimeter Array and the Square Kilometer Array should shed some light on this problem. These interferometers will allow the mapping of atomic and molecular gas and dust at high resolution, potentially enabling the study of how the relative quantities of these ISM components vary with metallicity and photodissociating UV radiation intensity in different environments. A formalisation of the processes governing the molecular content in galaxies will undoubtedly aid future cosmological simulations of galaxy formation and evolution.

Acknowledgments

We thank the anonymous referee for comments and suggestions, which improved the quality of the paper. T.M.H. wishes to thank T. Oosterloo for kindly sharing the FITS files of the Hi data, and F. Israel and N. Scoville for their CO observations, all of which were crucial to this work. T.M.H., M.B. and J.F. gratefully acknowledge financial support from the Belgian Science Policy Office (BELSPO) as part of the PRODEX project C90370 (Herschel-PACS Guaranteed Time and Open Time Programs: Science Exploitation). I.D.L. and F.A. are postdoctoral and doctoral fellows of the Flemish Fund for Scientific Research (FWO-Vlaanderen), respectively. M.B. and S.V. also acknowledge financial support from the same Flemish Fund for Scientific Research. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the NASA (USA). This research made use of Montage, funded by the NASA (USA) Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and Caltech, and maintained by the NASA/IPAC Infrared Science Archive. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the JPL/Caltech, funded by NASA (USA).

References

All Tables

All Figures

Fig. 1 Images of NGC 891 at various wavelengths, centred on the galaxy’s coordinates from NED: , (J2000.0). Each image is presented in its native resolution and pixel size, with the open symbols at the bottom left of each panel indicating the beam size at each wavelength. The image characteristics are described in detail in the text (see Sect. 2). The solid white ellipse denotes the aperture for integrated photometry and the dashed white circles represent the apertures to reconstruct the analysis of Dupac et al. (2003b), as described in Sect. 4.4. North is up and east is to the left.

In the text

	Fig. 2 Vertical profiles at 100, 160 and 250 μm, adopting a position angle of 22.9°, are plotted as blue circles connected by a blue line. We fit with an exponential profile (solid red line), which has been convolved with the Herschel beam at the corresponding wavelength (dashed black line).
In the text

Fig. 3 Radial profiles at 100, 160 and 250 μm are compared with the IRAC 3.6 μm and MIPS 24 μm emission. The solid black line represents the optical diameter adopting D25 = 37.7 kpc and the two dashed black lines demark the break radius found at ±12 kpc from Bianchi & Xilouris (2011). The alternative y-axis scale refers only to the dust temperature profile obtained from modelling the SED (see Sect. 4.3), represented by the solid red line.

In the text

	Fig. 4 The 100/160, 160/250, 250/350, and 350/500 μm flux ratios compared with the 3.6 and 24 μm emission. Pixels along the major axis shown in blue and off-axis pixels shown in grey. The best-fit linear relation (red line) and the Pearson correlation coefficient are based only on the on-axis pixels.
In the text

Fig. 5 Left: total FIR/submm SED of NGC 891 obtained from the integrated PACS, SPIRE and SCUBA fluxes (black circles). The red line shows the best-fit one-component model, where the 70 μm point is used as an upper limit. The dashed black line represents the best-fit one-component model found when including the FIR/submm GAUFLUX fluxes reported in the Planck HIFI Catalogue of Compact Sources (blue diamonds; Planck Collaboration XXVIII 2014). For comparison with previous FIR/submm experiments, the flux measurements from IRAS, ISO and IRAM presented in Popescu et al. (2011, see their Table 2) are shown (grey triangles). Right: SEDs measured for different locations in the galaxy, following Dupac et al. (2003a) and Bianchi & Xilouris (2011).

In the text

	Fig. 6 Comparison of the dust masses (upper panel) and dust temperatures (lower panel) obtained from allowing β to vary versus fixing the β value. The 1:1 relationship and our linear best fit relationship are indicated by the black and red lines, respectively.
In the text

	Fig. 7 Anticorrelation we observe between β and dust temperature from fitting a one-component modified blackbody to the pixel-by-pixel SEDs constructed from our FIR/submm images.
In the text

	Fig. 8 Maps of the dust mass surface density (left), dust temperature (middle) and χ2 value (right), obtained from fitting a one-component, modified blackbody model with β = 1.8 to the SED of each pixel. The maps are centred on , (J2000.0), with the beam size indicated by the black circle. The pixel size is 12″ pix-1 corresponding to 0.56 kpc at a distance of 9.6 Mpc.
In the text

Fig. 9 Pixel-by-pixel relationships between the derived dust properties, Σdust (upper panels) and Tdust (lower panels), and the surface densities of the gaseous ISM components: ΣHi (left), ΣH2 (middle) and Σgas (right). We differentiate between on- and off-axis pixels using solid and open circles, respectively. The best linear fits are indicated by the solid red lines. Error bars are calculated from the standard deviations of the errors on individual pixels and the ρ values are the corresponding Spearman correlation coefficients.

In the text

	Fig. 10 Maps of the Hi-to-dust ratio (left), H2-to-dust ratio (middle) and the total gas-to-dust ratio (right). The maps are centred on , (J2000.0), with the beam size indicated by the black circle. The pixel size is 12″ pix-1 corresponding to 0.56 kpc at a distance of 9.6 Mpc.
In the text

	Fig. 11 Correlation between the H2-to-dust ratio and the dust temperature. The best linear fits are indicated by the solid red lines, with the mean ΣH2/ Σdust ratio represented by the dashed black line. The cluster of points with ΣH2/ Σdust ratio lower than the mean and possessing a higher dust temperature (Tdust> 23.5 K, shown by the dotted black line) are all central pixels.
In the text

	Fig. 12 Map of the dust temperature obtained from fitting a one-component, modified blackbody model to the SED of each pixel (colour scale, see Fig. 8). Pixels containing dense regions of molecular hydrogen and cooler-than-average dust are highlighted (purple contours). The WISE F12/F22 flux ratio (black contours) follows the dust temperature.
In the text

	Fig. 13 Pixel-by-pixel relationships between the derived dust properties, Σdust (upper panels) and Tdust (lower panels), and the IRAC 3.6 μm (left panels) and MIPS 24 μm flux densities (right panels). The best linear fits are indicated by the solid red lines. Reported ρ values are the corresponding Pearson correlation coefficients.
In the text