The HASHTAG Project: The First Submillimeter Images of the Andromeda Galaxy from the Ground

Dust itself is interesting for two main reasons. First, it is of great intrinsic interest because it is a vital phase of the interstellar medium (ISM), containing half the heavy metals (James et al. 2002) and being a catalyst in the networks of chemical reactions in the ISM, including the vital one in which atomic hydrogen is transformed into molecular hydrogen. Second, mapping the continuum emission from dust grains is a promising method for both mapping the ISM in galaxies and estimating the total mass of the ISM (Hildebrand 1983; Eales et al. 2012; Magdis et al. 2012; Scoville et al. 2014). The standard tracer of the molecular phase, the CO molecule, has many well-known disadvantages (Bolatto et al. 2013). There is also now the more fundamental problem that one third of the molecular gas in the Galaxy appears to contain no CO (Abdo et al. 2010; Planck Collaboration et al. 2011; Pineda et al. 2013), and there are even galaxies in which the fraction of "CO-dark" gas seems to be close to 100% (Dunne et al. 2021). Some of the advantages of using dust grains rather than CO molecules to trace the ISM are that the dust emission is optically thin, dust grains are robust and not liable to be destroyed by starlight, and the relationship between the gas-to-dust ratio and metallicity seems to be much simpler than that between CO abundance and metallicity (Eales et al. 2012; Sandstrom et al. 2013; Rémy-Ruyer et al. 2014).

The biggest contribution that HASHTAG seems likely to make to our understanding of the dust itself is to show how the properties of dust vary within an individual galaxy. The earlier Herschel observations of Andromeda revealed systematic large-scale spatial variation in the properties of dust. The emission from interstellar dust follows a modified blackbody (S_ν ∝ B_ν (T_d)ν^β ). The Herschel observations revealed that β varies radially within Andromeda's disk (Smith et al. 2012; Draine et al. 2014; Whitworth et al. 2019), and radial variation in β has subsequently been found in the Galaxy (Planck Collaboration et al. 2014a), in M33 (Tabatabaei et al. 2014), and in ≃20 other galaxies (Hunt et al. 2015), although the form of the radial variation varies between galaxies (Hunt et al. 2015). There is also now some evidence that the global value of β varies between galaxies (Lamperti et al. 2019). The variation must be caused by changes in the structure, physics, or chemistry of the dust, although what the key changes are is currently unknown. One clue may be that in Andromeda β does not appear to differ between low-density and high-density gas (Athikkat-Eknath et al. 2021).

HASHTAG will produce measurements of β at ≃70,000 positions in Andromeda's disk, effectively producing a dust atlas for Andromeda. The Hubble observations of one third of the disk have produced estimates, from the optical extinction, of the column density of dust with the same spatial resolution that HASHTAG will provide (Dalcanton et al. 2012). The combination of measurements of the emission properties of dust from submillimeter observations and the absorption properties of dust from optical observations is a powerful one for testing theoretical dust models. The emission and absorption properties of dust, derived from Hubble and Herschel data, are already inconsistent with all existing dust models (Whitworth et al. 2019). The combination of the high-resolution measurements of β from HASHTAG, the Hubble dust measurements, and the maps of the ISM phases, star formation, chemical abundances, and other properties that are available for Andromeda offers at least the prospect of uncovering the physical/chemical causes of the variation in dust.

The use of the dust emission to trace the ISM in Andromeda offers a number of interesting possibilities. By comparing the dust emission to the CO and H i emission it will be possible to search for CO-dark gas in Andromeda (Planck Collaboration et al. 2011; Sandstrom et al. 2013). It will also be possible to produce catalogs of GMCs based on dust emission rather than CO. This technique has already been applied to the Herschel observations of Andromeda, producing a catalog of 326 clouds with masses between 10⁴ M_⊙ and 10⁷ M_⊙ (Kirk et al. 2015). A more recent study suggests that clouds found by the dust method have a much lower CO-to-dust ratio than the clouds found from their CO emission (Athikkat-Eknath et al. 2021), suggesting there is much more variation in the properties of clouds than one would expect from a CO-selected catalog. The clouds in the Herschel catalog are probably associations of GMCs rather than single GMCs, but with the extra resolution of HASHTAG it will be possible for the first time to produce a catalog of dust-selected clouds that are likely to be GMCs rather than GMC associations.

One of the most important properties to measure within a galaxy is the star formation rate, but there is still no gold-standard way of doing this. There are at least 12 different methods, which use different techniques for tracing the obscured and unobscured star formation (Kennicutt & Evans 2012; Speagle et al. 2014; Davies et al. 2016), all of which have limitations, and none is clearly better than the others.

HASHTAG will produce high-resolution maps of the bolometric dust emission, which is a direct measurement of the emission from the obscured OB stars, although it is really an upper limit since much of the bolometric dust emission is reradiated emission from the older stellar population (Bendo et al. 2012, 2015; Kennicutt & Evans 2012; Ford et al. 2013; Viaene et al. 2017). However, since Andromeda is close enough for us to detect individual stars, there is at least the possibility of combining Hubble observations of the unobscured OB stars and the HASHTAG observations of the emission from the obscured stars to provide a direct measurement of the star formation rate, rather than the measurements produced by the current methods, which mostly rely on indirect tracers of the star formation. The PHAT team made a first attempt to do this (Lewis et al. 2017), using optical extinction measurements to correct for obscuration, but their method was unable to account for the OB stars that are still deep in GMCs and are completely hidden by dust. If it is possible to correct for the part of the bolometric dust emission that is reradiated emission from the older stellar population, HASHTAG will provide estimates of numbers of these missing OB stars.

As soon as one of our observations was made, we processed it using a quick-look script that produced images at the two wavelengths (450 and 850 μm) and saved the cleaned data from the individual bolometers. We used these images to check for any severe problems (e.g., array failures or low sensitivity), which for SCUBA-2 data is extremely rare (only two out of 235 observations had issues). We used the saved bolometer data from this initial processing as the input for the next stage, which resulted in a considerable rationalization of the data, since the initial processing converted the hundreds of raw data files into just eight (one per array).

To reduce SCUBA-2 data, the observatory provides an iterative data-reduction procedure called makemap, which at the end of every iteration provides a better estimate of the astronomical signal and the noise until no further improvement is made. Here we provide a brief summary of the makemap algorithm but refer the reader to Chapin et al. (2013) for full details. makemap is provided as part of the starlink software package (Currie et al. 2014), and throughout this paper we use a recent development build of starlink (version 9072c4434 from 2020 April), because we required some updates since the last 2018A stable release.

Initially makemap splits the raw data into individual observations (subdivided into chunks if enough RAM is not available, although for HASHTAG "chunking" was not required). There is an initial cleaning step for each observation in which bad bolometers are masked and artifacts (e.g., glitches) are removed from the timelines.

makemap then starts the iterative process. At the beginning of each iteration, the common-mode signal (common to all bolometers) is removed. The data are then corrected for atmospheric extinction, and a high-pass filter is applied to the timeline data to remove any residual, slowly varying signal. Since in the SCUBA-2 arrays are moving across the sky, the removal of a slowly varying signal is equivalent to removing emission on a large angular scale. An image is then made from the data. Any real astronomical signal is then identified in the image, and this astronomical signal is then removed from the timelines (an optional signal-to-noise cut or mask can be applied). The process is then repeated on the new timeline data, with the astronomical signal being updated in each iteration. The process stops when the pixel variations in the map at the end of each iteration fall below a set threshold (i.e., when the map has "converged"). If we had used the standard implementation of makemap, our final image of Andromeda would have been a mosaic of the images made by makemap from each individual observation.

However, a weakness of makemap is that a single ∼43 minute Pong observation does not have the sensitivity to detect the low-surface-brightness emission in Andromeda's disk. Recognising this limitation, the starlink team created a script called skyloop, ⁴² which runs makemap on all the observations, one iteration at a time, combining the individual images at the end of each iteration to produce the best estimate of the astronomical signal, thus maximizing the signal-to-noise ratio in the faint extended structure. However, the volume of our data is so large that skyloop would take too long to run. We built on our previous work (Smith et al. 2019) to develop a new version of skyloop that we nicknamed "scubaDuperSkyloop". The main difference from the old script is that while skyloop is given all the observations and internally processes them individually, mosaicking the images at the end of every iteration, our modified script calls makemap separately for each iteration and each observation. We found that the new script was more stable on our system (possibly due to more regular memory clearance); however, the main advantage of our approach is that each observation can be processed in parallel, or even on separate machines. Although makemap can process individual observations using multiple threads, there are diminishing improvements in processing time.

To make the final images shown in this paper, we reduced all the data using a new processing machine with 768 GB of RAM and 32 cores/64 threads, which allowed us to run five invocations of makemap in parallel. The main factor limiting the speed of "scubaDuperSkyloop" was the speed of our hard disk drives. We increased the speed of the simulations described in Section 5 by running the script in parallel on two separate smaller machines (with common storage access).

Ground-based submillimeter surveys of sources with extended emission (greater than a few arcminutes) face the challenge of slow variations both in the atmospheric emission and within the camera, which need to be removed from the data. This is overcome using harsh high-pass filtering, which also removes real astronomical signal on large angular scales. However, there are submillimeter observations of many objects with the space telescopes Planck and Herschel. The images made with these telescopes do preserve the large-scale structure; however, because of the small sizes of the telescopes' mirrors they do not have as good resolution as images made with telescopes on the ground. In principle, we can now produce high-fidelity images of submillimeter sources with large angular sizes by combining observations with telescopes on the ground, which provide the small-scale structure, with observations made with telescopes in space, which provide the missing large-scale structure. This technique has often been used in radio astronomy to combine single-dish and interferometer measurements; as the latter sparsely samples the UV plane, the case here is potentially simpler.

We have produced the high-fidelity image of Andromeda at 850 μm by combining the SCUBA-2 image at 850 μm and the Planck one at 353 GHz, which is virtually the same wavelength. At 450 μm, we have combined the SCUBA-2 image at this wavelength with the Herschel image at 500 μm (Smith et al. 2012). To combine the low- and high-resolution images we have written a Python module that applies a "feathering" technique (Bajaja & van Albada 1979). The module performs the following steps:

• 1.  
  The low-resolution FITS image is reprojected so that the pixel scale and the celestial coordinates of the pixels are the same as in the high-resolution image. The units of both images are converted into Jy beam⁻¹ using information contained in the header or supplied by the user.
• 2.  
  We apply color corrections to the high- and low-resolution images to correct for the effect of different instrumental filters, calibration schemes (e.g., different reference spectra), and differing central frequencies. To apply this correction the user can specify a fixed dust temperature and β, provide maps of the dust parameters, or a cube created by the PPMAP algorithm (Marsh et al. 2015) with the surface density of dust for a grid of dust temperatures and β's. A precomputed grid specific to each far-infrared/submillimeter instrument is used to perform the corrections, and is provided with the task.
• 3.  
  The median value in each image is subtracted from the image and "NaN" pixels are replaced with zeros to avoid artifacts.
• 4.  
  Both images are Fourier transformed and shifted so a spatial frequency of zero is assigned to the center. The values in the Fourier transform (FT) of the low-resolution image are scaled by the ratio of the beam areas of the high-resolution and low-resolution images.
• 5.  
  A filter is then created to weight the FT images by the selected amount when they are combined. The standard filter in the module, which we used for HASHTAG, is a Gaussian filter in Fourier space. We chose the value of the filter's standard deviation using the simulations described in the next section. The FT of the high-resolution image is multiplied by one minus the Gaussian filter. The FT of the low-resolution image by default is not weighted, because the image is effectively already weighted by the point-spread function of the image; this is the same method employed by casa (McMullin et al. 2007). In our simulations we also try the alternative where the low-resolution FT image is weighted by the Gaussian filter, which is applicable when the filter is on significantly larger spatial scales than the resolution of the images. For a full discussion of combining images in the Fourier plane see Stanimirovic (2002). The weighted FT images are then added together. There is also an option in the module to use either a Butterworth filter (Csengeri et al. 2016) or sigmoid filter, in which case both FTs are multiplied by the filter. We tried these filters for HASHTAG but found they did not produce appreciable benefit over the Gaussian filter.
• 6.  
  The combined FT image is then inverse Fourier transformed (with appropriate inverse shifts applied). The NaN pixels are restored, the median background from the original low-resolution image is added back to the new image, and the keywords in the image header are updated. ⁴³
• 7.  
  As many spectral energy distribution (SED) fitting procedures apply a color-correction step in their processing, we remove the color corrections performed in step 2, so the flux densities in the final images are based on the same assumptions as regular SCUBA-2 images.

One key parameter in the "feathering" algorithm is the scale of the Gaussian filter. As the feathering step is a relatively quick process (compared to the SCUBA-2 pipeline) we keep this as a free parameter, which we optimize in Sections 5 and 6. For the color corrections we assume the SED in each pixel estimated by Whitworth et al. (2019), who applied the PPMAP algorithm to the Herschel data set to generate SEDs with the angular resolution of the highest-resolution Herschel image. For the Herschel and Planck images, we calculated the color corrections using these SEDs and the filter curves available on the observatories' websites, after removing the standard SED used to estimate the Planck and Herschel flux densities (F_ν ∝ ν⁻¹). SCUBA-2 flux densities, on the other hand, are calibrated relative to Mars and Uranus, which means that they are based on a very different assumption about SEDs (roughly F_ν ∝ ν^1.7, Lellouch & Amri 2008; Orton et al. 2014). We calculated the color corrections for the SCUBA-2 images after removing this assumption and then using the PPMAP SED in each pixel. For SCUBA-2, an additional complication is that the effective filter function is the product of the actual filter function and the atmospheric transmission. We calculated the effective filter function from the filter function available on the observatory's website and a model for the transmission of the atmosphere ⁴⁴ with τ_{225 GHz} = 0.065 to match the weather for our survey. The color corrections are small (≤3%) for the two SCUBA-2 filters and for the Herschel 500 μm filter, apart from the large correction needed to change the flux at 500 μm to one at 450 μm, which ranges from a factor of ∼1.5 in the center to ∼1.34 in the ring. The color correction for the Planck filter is ≃10%.

Figure 3 shows the results of the feathering technique when applied to SCUBA-2 simulation (see Section 5), illustrating how it is effective at restoring the structure on all spatial scales. The red line shows the power spectrum of a simulated "true" image of the sky at 850 μm. The green and blue lines show models of what Planck and SCUBA-2 would see, respectively, in one case missing the high-k and in the other case the low-k Fourier components. The orange line shows the power spectrum of our reconstruction of the sky by applying our feathering technique to the artificial Planck and SCUBA-2 images.

Zoom In Zoom Out Reset image size

We performed a sanity-check of the method by applying it to a Herschel 250 μm image of the Milky Way (Molinari et al. 2016), in which the large-scale emission is detected with high signal-to-noise ratio. We created an artificial image at the Planck resolution by smoothing the Herschel image, and we created a rough simulation of what a camera such as SCUBA-2 would see by using the Nebuliser algorithm developed by the Cambridge Astronomical Survey Unit ⁴⁵ to remove the structure on large scales. When we produced our combined image by applying our feathering technique to the two artificial images, we found a good agreement with the original image (to within a few percent), although in the brightest regions there were differences up to the 10% level.

There are also some parameters to tune in makemap, in particular the scale of the high-pass filter used to remove the signal from the atmosphere and the camera itself. To produce a reliable map of Andromeda it is crucial to get these values right. This is a particular challenge at the longer of the two SCUBA-2 wavelengths because of the need to ensure that emission on the SCUBA-2 images is preserved on all angular scales up to the angular resolution of Planck (≃5', FWHM). In the next section, we describe the simulations of the sky that we used to determine the best values of the parameters.

We chose to carry out a simulation of a single 30' Pong observation with a similar sensitivity to our pilot field. While ideally we would have simulated observations of the whole galaxy, the processing time would have been too long and there was no suitable SCUBA-2 data we could use for a simulation of the entire galaxy at our depth. For our simulation we used data from the SCUBA-2 Cosmology Legacy Survey (Geach et al. 2017). We chose to use the data from the survey of the Lockman Hole because it was carried out in similar weather conditions and consisted of 35 30' Pongs, the same as we used in our observations of our pilot field (Section 3), reaching a similar sensitivity.

We made our "true" image out of the Herschel 250 μm image, which has a resolution (18'', FWHM) that is not very different from that of SCUBA-2 at 850 μm (135, FWHM). We first reprojected the Herschel image onto a 4'' pixel grid and then multiplied the intensity values in each pixel by the ratio of the global 250 and 850 μm fluxes (Planck Collaboration et al. 2015). The Lockman Hole data are actually slightly less sensitive than the data for our pilot field (3.4 versus 3.0 mJy beam⁻¹), so to make the signal-to-noise ratio of our artificial image the same we multiplied the intensity values by the ratio of the noises. We also applied color corrections. We then converted the intensity values into instrumental units of picowatts (pW) and ran makemap with the options "fakemap" and "exportclean" set so that our artificial 850 μm image was added to the SCUBA-2 timelines for the Lockman Hole survey. We used the data files exported by this process, which also did the initial cleaning of the timelines, to carry out our simulations. We produced our artificial Planck image by convolving our input image to the Planck resolution.

We performed five different sets of simulations, to optimize different aspects of the method. While in an ideal world every possible combination of parameters in the SCUBA-2 data-reduction pipeline would be tested, each individual simulation required ∼8–10 hr to run. Therefore, we varied one parameter at a time. However, our reconstruction process in which we combined the processed SCUBA-2 image with our artificial Planck image was quite quick to run, so at the end of every run of the pipeline, we applied our reconstruction technique many times, each time with a different value for the scale of the Gaussian filter, with scales ranging from 160'' to 840''.

We measured the success of a simulation by the differences between our final image and the original true image. We assessed the significance of the differences using the noise image produced by the pipeline, but since changes to the pipeline also change the noise in the final image we also used a reference noise image (from a run using typical values of all the pipeline parameters). For each final SCUBA-2 image, we created several "difference maps" of the difference between the image and the true image: (1) a basic residual image, i.e., the final image minus the "true" image; (2) the residual image divided by the noise image; (3) the residual image divided by the reference noise image; (4) the residual image divided by the true image.

To assess the agreement, we measured statistics in two different regions of these difference maps: (1) the full-depth region of the Pong; (2) the full-depth region of the Pong but only for pixels where the flux in the pixel in the Herschel 500 μm image is above a critical value (500 μm is used as it is the closest in wavelength to both SCUBA-2 bands). The point of the second region was to stop the more numerous pixels outside the galaxy with little emission biasing the results, because a method producing a flat map (e.g., harsh Fourier filtering) would be preferred. We inspected all of these methods of assessing the agreement at one time or another. The statistic that we found most useful was the mean of the absolute difference between the final image and the true image divided by the reference noise image, for pixels in the deep region above the 500 μm threshold. This statistic using a reference noise map has the advantage that changes in the noise map do not bias the estimate of how well we recover the galaxy, while still accounting for variations in the sensitivity across the map.

The most important parameter in the SCUBA-2 data-reduction pipeline is the scale of the high-pass filter used to remove residual emission from the atmosphere or the instrument. We started our simulations with the expectation that we would need to set the scale to roughly the angular resolution of Planck. In our early results we found that with a filter scale of 340'' we were able to reproduce the true image well.

In 2019 April, however, the observatory released a new mode for the SCUBA-2 pipeline in which principal component analysis (PCA) is used to remove residual atmospheric and instrumental noise. The advantage of PCA is that it makes it possible to increase the angular scale of the high-pass filter, reducing the attenuation of the emission on the SCUBA-2 images on large angular scales. Of course, if one allows PCA to remove too many components, it is also possible to remove real astronomical emission from the image. In the SCUBA-2 pipeline the default number of PCA components to remove is 20 components per subarray, but after inspecting the final images made with this setting we decided it led to the removal of some real emission. We therefore decided to run a suite of simulations in which we varied both the scale of the high-pass filter and the number of PCA components.

In our initial simulations without PCA we had found an optimum filter scale of 340''. We realized that with PCA we should be able to increase this scale. We therefore ran simulations with filter scales between 340'' and 560'' and the number of PCA components between 0 and 20, running 91 simulations to cover this 2D parameter space (20 PCA components is the default in the new pipeline mode).

Figure 4 shows the effect of changing the number of PCA components while keeping a constant filter scale of 520''. The top left panel shows what happens if no PCA components are removed. The filter may not remove much large-scale emission from the galaxy but much of the emission visible in the picture is clearly spurious. Increasing the number of PCA components leads to better removal of these artifacts, but the removal of too many PCA components leads to the removal of real astronomical signal, producing the negative regions around the image in the bottom right panel, which has had 20 PCA components removed.

Zoom In Zoom Out Reset image size

An important point to note is that, although the eye tells us that the patches in the top left panel in Figure 4 are artifacts, the SCUBA-2 data-reduction pipeline treats these as real emission, which means that the values in the noise image produced by the pipeline are too low. Figure 5 shows how the average noise on an image depends on the number of PCA components and on the scale of the high-pass filter. The solid lines show the results if the noise is measured from empty areas of the final image; the dashed lines show the result if the noise is measured from the noise image produced by the SCUBA-2 data-reduction pipeline. The pipeline estimates are much lower, showing the effect of the pipeline treating the artifacts as real astronomical signal (users of the SCUBA-2 pipeline beware!). The more reliable noise values, measured from the images themselves, show that the noise in an image can be reduced either by decreasing the scale of the high-pass filter or by increasing the number of PCA components. In either case, of course, one also runs the risk of removing real astronomical signal.

Zoom In Zoom Out Reset image size

We chose the best combination of filter scale and number of PCA components based on a combination of visual inspection of the residual maps and the statistical estimates of the difference between the recovered image and the true image. Figure 6 shows the mean absolute residual for the difference map made by dividing the residual image by the reference noise image (number three in the list in Section 5.1). The statistic has been calculated for the pixels that are in the full-depth region and are above the 500 μm threshold value (see above).

Zoom In Zoom Out Reset image size

The figure shows that there is an advantage in increasing the scale of the high-pass filter from the 340'' that we had originally considered (see above) if removal of PCA components is included in the analysis. The best agreement between the recovered and true image is obtained for a filter scale of ≃520''. The figure also shows that there is an optimum number of PCA components of ≃8—removing more increases the difference between the recovered and true images. Based on these results, we visually inspected the recovered images for a filter scale of 480'' and five or six PCA components and for a filter scale of 520'' with seven or eight PCA components, concluding that the best results came with a filter scale of 520'' with seven PCA components. However, in the next stage of the simulations we also tested the method with a filter scale of 480'' and five PCA components, in order to check whether the optimum values shifted if other parameters in the method were varied.

In the simulations, we also checked the number of iterations in the pipeline required for each set of parameters, since the computer processing time is directly proportional to the number of iterations. We found, as expected, that there is a processing cost to setting a larger filter scale (from 8 to 19 iterations) but including PCA analysis can reduce this.

An important element in the SCUBA-2 data-reduction pipeline is a mask provided by the user as their best guess of the area in which real astronomical signal will be found. This greatly helps the convergence of the iterative procedure because it makes it easier for the algorithm to distinguish between real astronomical emission and extended structures that are actually the result of atmospheric emission or noise in the camera (see Figures 4 and 5). This does not mean that a pixel in the final image that is not within the mask will necessarily contain no astronomical signal. The software only subtracts astronomical signal from samples in the timelines that will contribute to image pixels within the mask, but once all data-reduction stages are completed (subtraction of PCA components etc.) even pixels outside the mask, which are the averages of many samples in the timeline, may still contain astronomical signal.

The mask we used was created from the Herschel 500 μm image of Andromeda (Smith et al. 2012), the Herschel image closest in wavelength to our 850 μm image. We defined the mask as all pixels above a 500 μm flux-density threshold, with this threshold providing another knob we could twiddle in our analysis. Too high a threshold makes it possible for the algorithm to treat real astronomical signal as noise and remove it from the timelines; too low a threshold leads to slower convergence of the algorithm and higher noise. In the simulations described in the previous section, we set the 500 μm flux-density threshold at 200 mJy beam⁻¹, which seemed a reasonable compromise because the mask then included most of the disk and some inner regions of the galaxy while excluding some regions between the rings where the emission is faint.

Once we had identified the best combinations of filter scale and number of PCA components (520'' and seven PCA components or 480'' and five PCA components), we used these to test the effect of changing the flux threshold used to construct the mask. We ran simulations with values of the 500 μm flux threshold between 120 and 520 mJy beam⁻¹. Figure 7 shows masks created with a range of flux thresholds that are representative of the ones we used in the simulations. Figure 8 shows the results of the simulations. The agreement between the input and output images is clearly best for a 500 μm flux threshold of 280 mJy beam⁻¹ (this includes ∼25% of the total flux of M31), which is the threshold we adopted to create the mask for the real 850 μm observations of Andromeda. We found that changing the flux threshold and thus the mask had a negligible effect on the noise of the final image, a very slight increase (<1%) for masks generated with a 500 μm flux threshold below 200 mJy beam⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The next parameter we investigated was the map-tolerance parameter, which the SCUBA-2 data-reduction pipeline uses to decide whether the algorithm has converged or whether more iterations are required. The default tolerance value is 0.05, which means that the iterations stop when the average change in the flux in a pixel from the last iteration is less than 0.05σ, σ being the noise in that pixel (calculated from the distribution of instrument samples contributing to that pixel). There have, however, been some studies (Mairs et al. 2015; Smith et al. 2019) that suggest a lower tolerance value might improve results—something we wanted to explore in our simulations.

Given our huge volume of data, another important consideration was processing time, which increases if the tolerance value is reduced because a greater number of iterations are then needed to reach a lower tolerance value (processing time scales linearly with number of iterations). We therefore carried out simulations over a fairly small range of tolerance values: 0.0075–0.05. In Section 5.3 we found that we achieved the best results with a mask generated with a 500 μm threshold of 280 mJy beam⁻¹. In the simulations described in this section, we also experimented with masks created with 500 μm thresholds of 200 and 240 mJy beam⁻¹ to see whether the choice of best mask changed if we also changed the tolerance parameter. We also tried both winner and runner-up from the competition between filter-scale/PCA combinations of Section 5.2 to see whether the order might be reversed with a different value of the tolerance parameter.

Figure 9 shows that decreasing the value of the tolerance parameter does improve the agreement between the true and recovered images. It also shows that the best choices for mask and filter-scale/PCA combination generally remain the best choices at all values of the tolerance parameter. The improvement with decreasing tolerance parameter is fairly slow, and there is a high price in increased processing time. Therefore, as a trade-off, we adopted a tolerance parameter of 0.03 for the real observations. Even with this tolerance parameter, reducing the current HASHTAG 850 μm data set, which is only 70% of the final data set, required 7.5 days of computer processing time.

Zoom In Zoom Out Reset image size

The final parameter we investigated was the scale of the Gaussian filter, the "feather scale," used to combine the SCUBA-2 and Planck images. In this final round of simulations, we changed the pixel scale from 40 to 45, which makes the pixels in the final 850 μm image close to one third of the FWHM of the point-spread function, the value that was eventually adopted after some experimentation for making Herschel images. Figure 10 shows the results from the simulations. The agreement between the "true" input image and the recovered output image is best when the scale of the filter is 320'', which is roughly what we expected; the resolution of Planck is 48 (290'') (Planck Collaboration et al. 2014b) and so the Planck image should supply all Fourier components on angular scales larger than this. In all the 850 μm simulations the feathering method where the low-resolution image is not weighted is preferred (see Section 4.3).

Zoom In Zoom Out Reset image size

As the result of these simulations, we adopted the following values for all of the parameters when reducing the real SCUBA-2 data.

• 1.  
  We set the scale of the high-pass filter in the SCUBA-2 data-reduction pipeline to 520'' (flt.filt_edge_largescale = 520).
• 2.  
  We set the number of PCA components per array in the SCUBA-2 data-reduction pipeline to seven (pca.pcathresh = −7) with the values of all the other parameters in the PCA analysis having their default values.
• 3.  
  We use an input mask to define the region where there is likely to be astronomical signal created from the Herschel 500 μm image (Smith et al. 2012). This is used by the SCUBA-2 data-reduction pipeline (the AST model in the pipeline). We defined the mask as all pixels with a 500 μm flux greater than 280 mJy beam⁻¹.
• 4.  
  We set the tolerance parameter in the SCUBA-2 data-reduction pipeline to 0.03.
• 5.  
  We use a pixel scale for the final 850 μm image of 45.
• 6.  
  We used a Gaussian filter with a scale (the "feathering scale") of 320'' to combine the final SCUBA-2 image with the Planck image.

These pipeline parameters are optimized for the observing strategy, source properties, and weather for HASHTAG and M31. For other SCUBA-2 data sets with extended structure, we would recommend performing a similar simulation to optimize the processing; however, these results should provide a useful initial guess.

The final stage in the simulations was to assess the fidelity of the final image produced with the values of the parameters listed in the previous section. How close is the structure in the final image to the structure in the original true image? This analysis gives us a useful estimate of the fidelity of our real image. Note, however, that analysis will yield an upper limit to the errors on the real image because the simulations have been carried out for only one Pong (Section 3); the spatial overlaps of the many Pong fields for the real observations (Figure 2) should improve the fidelity of the final real image.

Figure 11 shows the input "true" image, the recovered output image, and the difference between the two divided by the noise image produced by the SCUBA-2 data-reduction pipeline. If our recovery method were perfect, the final panel should simply show random noise, but in fact there is some faint structure in the noise that is clearly correlated with bright structures. We therefore need to assess the importance of these systematic errors.

Zoom In Zoom Out Reset image size

We estimated the random and systematic errors in the flux densities from a plot of D = (F_o − F_i)/σ_pipe versus F_i, in which F_o is the flux in a pixel in the output image, F_i is the flux in that pixel in the input image, and σ_pipe is the estimate from the SCUBA-2 data-reduction pipeline for the noise in that pixel. The top panel of Figure 12 is a surface-density plot showing how the number of pixels depends on D and F_i. We have only included pixels in the full-sensitivity central circular region of the image (diameter of 30'). If the fidelity of the final image were perfect, and if our noise estimates were correct, D should have a Gaussian distribution around zero with a standard deviation of one. In reality, the plot shows that there is a clear bias in D, which is systematically higher than zero at high flux densities in the input true image. The standard deviation of D is also slightly higher than one, showing that the estimate of the noise produced by the pipeline is too low.

Zoom In Zoom Out Reset image size

Given the discrepancy between the actual distribution of D and its predicted distribution, we have assumed that the true uncertainty of the flux in each pixel is given by three uncertainties added in quadrature:

$$\begin{eqnarray}&&{\sigma }_{{\rm{tot}}}^{2}={a}^{2}+{\left(b{\sigma }_{{\rm{pipe}}}\right)}^{2}+{\left(f{S}_{{\rm{o}}}\right)}^{2}\end{eqnarray} \tag{ 1 }$$

in which a is a constant, b is a multiplicative factor for the error given by the SCUBA-2 data-reduction pipeline (σ_pipe), and f is a multiplicative factor for the flux in the output image (S_o). The third error in Equation (1) is effectively a photometric calibration error, which exists for all astronomical observations, although in this case it is an error on top of the standard SCUBA-2 photometric calibration error, which we have not included in the equation. We estimated the values of a, b, and f by applying the minimization package lmfit (Newville et al. 2016) so that the standard deviation of D approached as closely as possible a value of 1. We carried out the minimization on a rolling group of 800 pixels ranked in flux. We found a = 0, b = 1.04, showing the SCUBA-2 data-reduction pipeline had slightly underestimated the random errors in the fluxes, and that f = 0.12, showing that there is a systematic error that depends on the brightness of the emission, confirming the qualitative impression produced by Figure 11.

The bottom panel of Figure 12 shows the same surface-density plot of pixels as the top panel but with the noise value predicted by the pipeline (σ_pipe) replaced by the noise value calculated from Equation (1) (σ_tot). The distribution shows that while the width of the distribution of D now is roughly correct, the distribution is still distorted, in the sense that as the output flux density (F_o) increases, it becomes progressively higher than the input flux density (F_i), although there is a suggestion above 30 mJy beam⁻¹ that this reduces. We could have corrected the flux densities in the real HASHTAG image using the red curve in the bottom panel of the figure. We decided not to do this for two reasons. First, the real HASHTAG image is made of a large number of spatially overlapping data sets, so it is possible this effect is less for the real image. Second, this systematic effect is fairly small compared with the statistical error: only 0.6σ at F_o = 21 mJy beam⁻¹.

The SCUBA-2 Makemap routine produces maps in instrumental units of picowatts, and so a flux conversion factor (FCF) is used to convert these units into either Jy beam⁻¹ or Jy arcsec⁻². There have recently been two changes at the JCMT, which have adjusted the standard FCF values used at the observatory. First, in 2016 November a filter set in SCUBA-2 was changed, which predominantly affected the 850 μm FCF. Second, in 2018 May the maintenance of the secondary mirror resulted in a change in the value of the FCF for 450 μm. We have observations both before and after these changes. Based on an interim analysis by observatory staff (private communication), we assume FCF values of 3.62 and 2.14 $\mathrm{Jy}\ {\mathrm{pW}}^{-1}\ {\mathrm{arcsec}}^{-2}$ at 450 and 850 μm, respectively, for observations post 2018 May. ⁴⁶ For our older observations, we multiply the "cleaned" data (see Section 4.1) by a correction factor, so that the final data can be calibrated with the same FCF. We adopted values for this correction factor of 1.21053 and 1.06481, for 450 and 850 μm, respectively.

The final factor we have to consider is that the standard JCMT calibration scheme is not designed for such extended objects as Andromeda. The JCMT calibration scheme (Dempsey et al. 2013) uses the flux of a calibrator source within a 30'' radius aperture, after subtracting a background measured in an annulus around the calibrator between radii of 45'' and 60''. This scheme is fine for calibrating images that contain point sources. But we are trying to calibrate very extended emission, and the beam of the telescope extends to much larger radii than the radii used in the standard calibration scheme.

We have adopted the calibration scheme we devised for the JINGLE Large Program (Smith et al. 2019), where we multiplied the FCFs so that the flux densities in the images matched the convention of other telescopes, that an aperture centered on a galaxy would only include the entire flux of the galaxy if the radius were increased to infinity. By integrating the beam model in Dempsey et al. (2013), we calculated that the standard 850 μm FCF should be multiplied by 0.91, which agreed with the curve of growth found by Dempsey et al. (2013). Doing the same calculation at 450 μm, we obtain a correction factor of 0.99. The curve of growth given in Dempsey et al. (2013), however, suggests a slightly smaller factor of 0.97. We suspect the difference is caused by extra non-Gaussian features in the beam at large radii. We therefore adopt the smaller value of 0.97. The final values of the FCF used to create these HASHTAG images were therefore the FCF values given above multiplied by these correction factors. These are 3.51 and 1.95 Jy pW⁻¹ arcsec⁻² at 450 and 850 μm, respectively.

The real HASHTAG data were reduced using the method outlined in the previous subsections, using the best values of the parameters found from the simulations. Our final maps (Data Release 1, DR1) are composed of two multi-extension fits files that contain the flux density, uncertainty, and sensitivity maps for the 450 and 850 μm images, respectively.

The uncertainty map contains the true uncertainty values we derived in Section 5.7, which include both the statistical uncertainty in each pixel and the systematic uncertainty as the result of the flux density in that pixel (the third term in Equation (1)). The sensitivity map is the same except that this map does not include the systematic term (f in Equation (1) is set to zero). At both wavelengths, the sensitivity of our final images exceeds our targets. In the 10 kpc ring, the typical sensitivity is ∼2.0 and ∼30 mJy beam⁻¹ at 850 and 450 μm, respectively, with peak sensitivities in the center of 1.5 and 20.6 mJy beam⁻¹ at 850 and 450 μm, respectively. As a rough comparison the Herschel 500 μm observations have an instrumental sensitivity of ∼11 mJy beam⁻¹ (Smith et al. 2017), and the point source sensitivity of Planck at 850 μm is ∼69 mJy beam⁻¹ (Planck Collaboration et al. 2014c).

The flux-density and sensitivity maps are shown in Figure 16. The sensitivity is not quite as uniform at 450 μm, which we attribute to variations in submillimeter opacity during the periods we took the data, since opacity variations have a bigger effect at the shorter wavelength.

Zoom In Zoom Out Reset image size

As well as the images shown in Figure 16, we also provide versions that have three different levels of Gaussian smoothing (for example, Figure 1), so that users can balance resolution versus signal-to-noise ratio. In these smoother images, the raw images have been smoothed with Gaussians with FWHMs of 4'', 5'', and 79 at 450 μm and 7'', 10'', and 13'' at 850 μm (this equates to effective resolutions of 89, 93, and 112 at 450 μm, and 152, 168, and 191 at 850 μm).