# High-sensitivity high-speed compressive spectrometer for Raman imaging

Benneth Sturm  Département de Physique, Ecole Normale Supérieure / PSL
Research University, CNRS, 24 rue Lhomond, 75005 Paris, France  Fernando
Soldevila  GROC-UJI, Institute of New Imaging Technologies (INIT), Universitat
Jaume I, E12071 Castellò, Spain  Enrique Tajahuerce  GROC-UJI, Institute of
New Imaging Technologies (INIT), Universitat Jaume I, E12071 Castellò, Spain
Sylvain Gigan  Laboratoire Kastler Brossel, CNRS UMR 8552, Ecole Normale
Supérieure, PSL Research University Sorbonne Universités & Université Pierre
et Marie Curie Paris 06, F-75005, Paris, France  Herve Rigneault  Aix
Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, F-13013 Marseille,
France  Hilton B. de Aguiar  Corresponding author<EMAIL_ADDRESS>Département de Physique, Ecole Normale Supérieure / PSL Research University,
CNRS, 24 rue Lhomond, 75005 Paris, France

###### Abstract

Compressive Raman is a recent framework that allows for large data compression
of microspectroscopy during its measurement. Because of its inherent
multiplexing architecture, it has shown imaging speeds considerably higher
than conventional Raman microspectroscopy. Nevertheless, the low signal-to-
noise (SNR) of Raman scattering still poses challenges for high-sensitivity
bio-imaging exploiting compressive Raman: (i) the idle solvent acts as a
background noise upon imaging small biological organelles, (ii) current
compressive spectrometers are lossy precluding high-sensitivity imaging. We
present inexpensive high-throughput spectrometer layouts for high-sensitivity
compressive hyperspectroscopy. We exploit various modalities of compressive
Raman allowing for up to 80X reduction of data storage and 2X
microspectroscopy speed up at a 230-nm spatial resolution. Such achievements
allowed us to chemically image sub-diffraction-limited biological specimens
(lipid bilayers) in few seconds.

##  I Introduction

Raman microspectroscopy allows probing a chemically heterogeneous system with
high spatial resolution and superb molecular selectivity. In particular,
Raman-based imaging holds promising applications for biological specimens
studies Kong  _et al._ (2015), such as for cancer Harmsen  _et al._ (2015);
Jermyn  _et al._ (2015) and bacteria detection Lorenz  _et al._ (2017);
Bodelón  _et al._ (2016), molecular biology imaging (Okada  _et al._ , 2012;
Klein  _et al._ , 2012; Ando  _et al._ , 2015; Donaldson Jr and de Aguiar,
2018), to cite few examples. Nevertheless, hyperspectroscopy (i.e. spatially-
resolved spectroscopy), in general, generates large data sizes which will
eventually become a hurdle when bringing Raman spectroscopic imaging as a
valuable tool for clinical tests, for example, for issues related to data
storage and access. As Raman microspectroscopy acquisition speeds continuously
advance Shipp  _et al._ (2017), it will become increasingly important to
achieve data compression, in particular for dynamic and/or large-scale imaging
at high-resolution.

Single-pixel-based compressive sensing technology has attracted a large
interest recently Duarte  _et al._ (2008). The idea behind compressive sensing
is to obtain the same information as in multi-pixel technologies approaches,
however through a designed optical experiment that will compress the data at
the acquisition stage, by undersampling. In a post-processing step, efficient
computational methods retrieve the desired information accurately.
Interestingly, single-pixel approaches are particularly powerful where multi-
pixel technology is still lacking or challenging, such as in terahertz and
pump-probe spectroscopy methods Berto  _et al._ (2017); Chan  _et al._ (2008).
Single-pixel compressive sensing is a very appealing approach for Raman
imaging as it is inexpensive and may have very fast computational
reconstruction times, allowing for real-time applications. Compressive
microspectroscopy is based on spectral, or spatio-spectral multiplexing August
_et al._ (2013), which is achieved by fast programmable spatial light
modulators, such as digital micromirror devices (DMD), operating at 10’s kHz.

While compressive microspectroscopy Davis  _et al._ (2011); Galvis-Carreño
_et al._ (2014); Thompson  _et al._ (2017); Wilcox  _et al._ (2012, 2013); Wei
_et al._ (2016); Réfrégier  _et al._ (2018); Scotté  _et al._ (2018); Corcoran
(2018); August  _et al._ (2013); Pavillon and Smith (2016); Wadduwage  _et
al._ (2017); Zhang  _et al._ (2018) is a promising technology, it still lacks
sensitivity for bio-imaging when compared to sensitive CCD cameras. High-
sensitivity imaging is of particular relevance in the sub-micron spatial
scale, as it is the size range of many organelle in a cell Klein  _et al._
(2012); Ando  _et al._ (2015); Rieger  _et al._ (2017); Donaldson Jr and de
Aguiar (2018), including lipid droplets, cell membranes, and membrane-less
organelles. The challenges for high-sensitivity compressive microspectroscopy
are two-fold. First, apart from Raman scattering being a weak effect, the
Raman signal ($I_{R}$) scales unfavourably with spatial resolution:
$I_{R}\propto\omega_{0}^{2}$ for objects that fill the focus ($\omega_{0}$
being the focus size of the excitation laser). Therefore, upon low-resolution
imaging, the overall signal increases but the sensitivity decreases for
detecting objects that have a size smaller than the focus, as the surrounding
solvent adds as a strong background. Second, and related to compressive
spectrometer layouts, there is a strong size mismatch (100X) between sensitive
detector sizes and the size of the spectrum spread in the Fourier plane of a
spectrometer (i.e. on the DMD). This mismatch is the main cause of sensitivity
loss in compressive spectrometers.

Here, we exploit novel high-throughput and inexpensive spectrometer designs
for high-resolution and high-sensitivity compressive Raman imaging of
biological specimens. This is achieved by dispersion-controlled layouts to
focus a large spectral bandwidth on a small single-pixel detector. After
discussion of the methods we propose, we showcase two compressive Raman
modalities exploiting different levels of a priori knowledge from the system
being imaged. We then demonstrate chemical imaging of sub-wavelength Raman
sources using the compressive spectrometer.

##  II Methods

### Confocal microscope

The compressive Raman microscope is based on a standard layout of an epi-
detection Raman microscope (Fig. 1). The excitation laser ($\lambda=532$ nm,
150 mW, Oxxius LCX-532) is cleaned up by a bandpass filter (FLH05532-4,
Thorlabs), and expanded to about 7.2 mm, before entering an inverted
microscope (Nikon Eclipse Ti-U). An additional half-waveplate and quarter-
waveplate (not shown) pre-compensate and generate circularly polarized light
Bélanger  _et al._ (2009) after the beam being reflect off from the dichroic
filter (F38-532_ T1, AHF) into a high-numerical aperture (NA), oil-immersion
objective (Nikon 60X, 1.4 NA). Specimens are mounted on a translation-stage
(Physik Instrumente P-545.3R7) for raster scanning. The inelastically
backscattered Raman light (epi configuration) is focussed by the microscope
tube lens on a multimode fiber (M50L02S-A, Thorlabs, 50 $\mu$m core-size),
with a notch filter blocking residual pump light (NF533-17, Thorlabs), thus
guiding the signal to the home-built compressive spectrometer.

### Single-pixel compressive spectrometer

The home-built spectrometer is based on a traditional Czerny-Turner design.
The signal coming out of the fiber is collimated by an achromatic lens L1
(f=19mm, AC127-019-A-ML, Thorlabs), which then impinges on a high-efficiency
transmission grating (FSTG-VIS1379-911, Ibsen Photonics, 1379 l/mm) providing
about 85 % diffraction efficiency around $\lambda$=600 nm. After the grating,
the beam is spectrally dispersed and different colors are focused with an off-
axis parabolic mirror PM (f=101.6 mm, 90∘ off-axis angle, MPD149-P01,
Thorlabs). The DMD module (V-7001, ViALUX, 0.7” XGA resolution, with V4395
controller board) selects the various wavelength to be detected. Typically, we
bin about 8 pixels, thus having a macropixel size that is slightly below the
spectral resolution of the spectrometer.

The signal from the DMD is steered into an avalanche photodiode operating in
photon counting mode (SPCM-AQRH-44, Excelitas Technologies). This module
combines high quantum efficiency from Si avalanche photodiodes (70 %), large
dynamic range (low dark and high maximum count rate), with an active area
diameter of 180 $\mu$m. The small area, while keeping high quantum efficiency,
is of particular importance in the design because focusing the complete
spectrum to the size of the detector-chip efficiently is not straightforward.
We note other single-photon counters exist with mm’s active area size,
however, with lower quantum efficiency.

### Sample preparation

Polystyrene particles (1 $\mu$m and 0.088 $\mu$m-diameter, Polysciences Inc.)
were drop cast on a coverslip, then filled with water. Multilamellar vesicles
(MLV) were prepared by drop-casting a film of lipid solution (10 mg/ml
dipalmitoylphosphatidylcholine in chloroform) on a vial, then adding water and
heat up at $38^{\circ}$ for about 30 min Gasecka  _et al._ (2017). After
spontaneous formation of the vesicles, an aliquot of the dispersion is then
sandwiched between two coverslips.

Figure 1:  Simplified layout of the compressive Raman microscope. A CW laser
source excited Raman scattering process, which is detected in epi
configuration. The sample is raster-scanned with a high-resolution piezo-
scanner. The Raman inelastically scattered light is steered into a home-built
spectrometer designed for compressive single-pixel experiments. The dashed
boxes depicts two strategies (grating or lens telescope) for dispersion
management, introduced by the first grating, to focus large bandwidths in the
high sensitivity single-pixel detector. DMD=digital micromirror device,
PM=off-axis parabolic mirror, L=lens.

## Results

### Dispersion management for single-pixel compressed spectroscopy

The main challenge of a single-pixel compressive spectrometer is as follows.
In a conventional spectrometer, a CCD camera is used to detect the dispersed
input light, which size is typically few mm’s. In contrast, the compressive
spectrometer consists of a DMD, in place of the CCD, which steers wavelength
bins onto a small single-pixel detector (100’s $\mu$m) to achieve high-
sensitivity. Therefore, one has to considerably demagnify the beam reflected
off by the DMD. Hence, we studied two strategies that could keep the
interesting aspects (costs and throughput) of single-pixel detection.

We investigated two approaches that are presented in dashed box insets of Fig.
1. The first approach is based solely on lenses (lens-based), motivated by the
fact that transmission of lenses are highest compared to other optical
elements. To aid the discussion, we recall the Smith-Helmholtz invariant
($y\sigma=y^{\prime}\sigma^{\prime}$) relating planes $y$ (grating 1), ray
angles $\sigma$ (arising from grating dispersion before PM), to $y^{\prime}$
(L3) and $\sigma^{\prime}$ (arising from focusing with L2) of an optical
system. As the grating plane $y$ is imaged onto plane $y^{\prime}$ by a PM and
L2, we have:

$\frac{\sigma}{\sigma^{\prime}}=\frac{y^{\prime}}{y}=\frac{f_{2}}{f_{1}}=m.$
(1)

Eq. 1 shows that dispersions and magnifications can be interchanged. For
instance, $\sigma^{\prime}$ can be reduced by simply increasing the
magnification $m$. However, the maximum magnification is limited by the clear
aperture of L3. The second approach is based on a second grating to cancel the
dispersion of the first one (Fig. 1, grating-based), as typically done in
pulse shapers of ultrafast spectroscopy, if they have the same groove density.
Note that if the second grating is of different dispersion, one can adjust
$f_{2}$ to match the dispersion given by $\sigma^{\prime}$.

We have compared the two methods. For the lens-based approach, we used an
effective 4X demagnification (in the figure it is presented as a single lens
L2). For the grating-based method, we used a second grating with the same
dispersion as the first grating. We have inspected the spectrum of the lens-
based approach (17 nm-bandwidth) and grating-based approach (77 nm-bandwidth).
While both methods coupled large bandwidths (over 100’s cm-1), the lens-based
had too much losses originating from clipping in the edge of lenses.
Therefore, we used the grating-based design. Nevertheless, we believe the
lens-based throughput could overcome the grating-based if a low-dispersion
grating is used.

### Compressive Raman without a priori knowledge

We exploited the spectrometer for two different types of compressive Raman
schemes. Basically, there are two ways of performing compressive Raman: with
(supervised) Wilcox  _et al._ (2012, 2013); Scotté  _et al._ (2018); Wei  _et
al._ (2016); Corcoran (2018) or without (unsupervised) Galvis-Carreño  _et
al._ (2014); Thompson  _et al._ (2017) a priori knowledge about the principal
components (or eigenvectors) of the hyperspectrum. In this section, we present
compressive Raman results without a priori knowledge.

Figure 2:  Compressive Raman imaging without a priori knowledge. Left: Raman
pseudocolored image of polystyrene beads obtained by means of computational
reconstruction. Right: Spectral content of both species present in the sample,
water (green), and polystyrene (red). Compression: 2X.

The measurement process of the hyperspectrum is done by a combined use of the
piezo-scanner and the DMD, described elsewhere Soldevila  _et al._ (2018).
Briefely, the sample is illuminated in a raster scan manner over the image
plane. For each spatial position, the Raman signal generated by the sample is
measured using the compressive spectrometer. Instead of performing a
wavelength scan over the spectral domain, the spectral content is measured
using multiplexed detection with Hadamard functions generated with the DMD.
One key advantage of using Hadamard basis is for increasing the SNR Decker
(1971); DeVerse  _et al._ (2000); August  _et al._ (2013) because half the
pixels are always on. This sampling scheme can be expressed in matrix form as
$\mathbf{y}=\mathbf{S}\mathbf{x}$, where the vector $\mathbf{y}$ contains each
one of the integrated measurements, $\mathbf{x}$ is the hyperspectral object
(in vector form), and the whole sensing process is described by the linear
operator $\mathbf{S}$, which is the Kronecker product of the two sampling
bases used here: $\mathbf{S}=\mathbf{I}\otimes\mathbf{H}$ ($\mathbf{I}$ being
the canonical spatial-scanning basis, and $\mathbf{H}$ the Hadamard spectral-
scanning basis). In order to reduce the number of measurements, we select a
random subset of rows of $\mathbf{S}$ to measure the hyperspectral dataset.
Then, we use compressive sensing algorithms to solve the following
optimization problem:

$\arg\min\mathbf{\hat{x}}=\frac{1}{2}\|\mathbf{y}-\mathbf{S}\mathbf{\hat{x}}\|^{2}_{2}+\tau\|\mathbf{\hat{x}}\|_{TV}$
(2)

based on the Total Variation as the regularizing function, which has been
previously reported to provide good results in multispectral imaging scenarios
Bioucas-Dias and Figueiredo (2007); Rajwade  _et al._ (2013). The
reconstructions can be seen in Fig. 2 of a relatively small region of the
sample (64x64 spatial positions and 128 spectral channels). From the spectral
information, we can distinguish between the water content of the sample
(coloured green), and the polystyrene beads (red).

### Compressive Raman with a priori knowledge

Figure 3:  Compressive Raman imaging with a priori detection. (A) Large scale
imaging with 80X compression (axis are in pixel count). Raman pseudo-color
images of 1 $\mu$m-polystyrene beads (red, upper left inset) and water (green,
lower left inset). Spectra are presented in the inset with the optimized
filters shown on top. Scale bar: 20 $\mu$m. Effective pixel dwell time: 640
$\mu$s. (B) High-sensitivity imaging of 88 nm-polystyrene beads. Only the
polystyrene channel is presented. Scale bar: 1.5 $\mu$m. Effective pixel dwell
time: 3 ms

In supervised compressive Raman, one estimates the proportions of a linear
combination of $M$ eigenspectra ($\mathbf{P}$), from $N$ measurements
($\mathbf{x}$). This is achieved by projecting the input spectrum onto $N$
optimized spectral filters ($\mathbf{F}$) with the DMD. Each spatial pixel can
then be written as $\mathbf{x}=\mathbf{F~{}P^{\top}}~{}\mathbf{c}$, where
proportions $\mathbf{c}$ are estimated from simple inversion of
$\mathbf{F~{}P^{\top}}$. The key step for proportion estimation is to develop
filters that are robust against noise Wilcox  _et al._ (2012); Réfrégier  _et
al._ (2018). We use as a metric the variance of the Cramer-Rao bound (CRB) to
determine the optimized spectral filters. A thorough description of the method
can be found elsewhere Réfrégier  _et al._ (2018); Scotté  _et al._ (2018).

Fig. 3 presents results of large-scale and high-sensitivity chemical imaging.
Fig. 3A shows an example of the compressive detection approach where a large-
scale high-resolution image (0.5 Mpixel) is presented with chemical
selectivity. After generation of the optimized filters (Fig. 3A, optimized
filters are presented above spectra in the top-right inset), we obtain
$\mathbf{x}$ from which we estimate $\mathbf{c}$ in a pixel-wise manner.
Chemically selective images are then presented in pseudo-color images (Fig.
3A) or analysed separately in their respective channels (Fig. 3A insets shows
water and polystyrene proportions). Due to high-resolution imaging, the
pseudo-color images are highly contrasted (lack of yellow color).

To further quantify the resolution and demonstrate high-sensitivity detection,
we present images of 88-nm polystyrene beads, sizes below the diffraction
limit of the system, in Fig. 3B. Despite the low SNR, one observes regions of
aggregated beads with high intensities, and isolated ones with much lower
intensities. To estimate the resolution of the system, we use a simple
deconvolution using the measured line profiles of the low-intensity beads
($\approx$250 nm). We obtain a resolution of $\approx$230 nm, a value close to
the diffraction limit of the system.

### Sensitivity to model membranes

Figure 4:  High-sensitivity detection of lipid membranes in vesicles. (A)
Proportions of lipids (left panel) and water (right panel). The water channel
has been scaled to enhance the appearance of the signal dip. Scale bar: 10
$\mu$m. Total acquisition time: 30 s. Effective pixel dwell time: 2.4 ms. (B)
Normalized (by the image mean) line profiles of the water proportion. The
signal dip is due to a depletion of water from the focal volume, which allows
estimation of the number of bilayers in the focal volume (see text for further
details).

Finally, we use the compressive spectrometer to demonstrate sensitive imaging
of biological specimens with sizes below the diffraction limit. The sample
chosen are vesicles made of lipids as they are model systems in cell membrane
biology Ando  _et al._ (2015); Donaldson Jr and de Aguiar (2018), and also
because the eigenspectrum is not expected to considerably change for different
environments. Fig.4A shows the estimated proportions of the lipid and water
response. In the lipid channel, there are clearly various bilayers wrapped
around each other, a topology common to thin multilamellar vesicles Moreaux
_et al._ (2000). In the water channel, there is a slight decrease in the
overall water response due to the finite thickness of the lipid bilayers. The
drop in intensity in the water channel allow us to estimate the water
depletion thickness, as the Rayleigh length is much smaller than the vesicle
curvature. From the intensity drop of the water channel (Fig. 4B), we estimate
a maximum multilayer thickness of 54 nm (the strongest water proportion
intensity dip), corresponding to 11 bilayers at that location.

## Discussion

We have presented a new spectrometer layout designed for single-pixel
spectroscopy. With this layout, we showcase two types of compressive Raman
modality, either without a priori knowledge or with a priori knowledge of the
eigenspectra of the hyperspectrum. The main purpose of the high-throughput
spectrometer is to enable sensitive measurements of specimens that are smaller
than the microscope resolution; Because the solvent background generates a
strong background, we reduce its intensity by tight focusing, and also
increase the dynamic range of the detection system with low-noise and high-
sensitivity detection.

Remarkably, we could detect with molecular selectivity few bilayers with the
current design. Previously, vesicle systems were mostly imaged by coherent
Raman scattering techniques Müller and Schins (2002); Potma and Xie (2003); Li
_et al._ (2005); Cheng  _et al._ (2003); Chen  _et al._ (2010); Gasecka  _et
al._ (2017); Shen  _et al._ (2017), which are considerably more expensive and
complex than the layout presented here. We estimated that the strongest region
in Fig. 4A contains about 11 bilayers. Given the linearity of the spontaneous
Raman process on the local concentration, we expect single-bilayer Ando  _et
al._ (2015); Rieger  _et al._ (2017); Donaldson Jr and de Aguiar (2018)
imaging with slight improvements in the current design, for instance, by
addition of a second detector Rehrauer  _et al._ (2018) and with larger clear-
aperture gratings to reduce losses due to clipping.

There have been various demonstrations to accelerate Raman imaging for
biological applications. These include line scan methods Okada  _et al._
(2012), spatial focii multiplexing Kong  _et al._ (2016), and widefield
microscopy Müller  _et al._ (2016); Wadduwage  _et al._ (2017). High-
sensitivity in slit-based spectrometers Okada  _et al._ (2012); Kong  _et al._
(2016) is challenging as non-confocal geometries are hard to reach
diffraction-limited performance Schlücker  _et al._ (2003). In that sense,
light-sheet excitation-widefield imaging Müller  _et al._ (2016) could
probably provide highest sensitivity, however the whole image plane has to be
optically clear and static during acquisition because of the spectral/temporal
sampling. The abovementioned methods all use sensitive cameras, contraty to
single-pixel approaches, thus avoiding slow and noisy measurements at high
speeds (mostly read-out noise) Scotté  _et al._ (2018). In addition, the
combination of single-pixel-spectroscopy with tight-focusing (or confocal)
allows for low-power CW lasers to be used. Nevertheless, the main disadvantage
of compressive Raman with a priore information is its need for at least one
long scanning to obtain the eigenspectra ”dictionary”, although a ”blind”
sampling approach to learn the dictionary has been proposed Wilcox  _et al._
(2013). Note that even in biological applications one expects only a few
eigenspectra Klein  _et al._ (2012), which supports that compressive Raman
approaches may have the fastest imaging speeds.

Previous compressive spectrometers used less sensitive detection. Quantum
efficiency and dark counts of the detector are two key parameters for high-
sensitivity detection as they determine the overall throughput and dynamic
range for proportion estimation. The overall throughput scaling of the
grating-based approach is $\phi_{l}^{4}~{}\phi_{\text{\tiny
DMD}}~{}\phi_{g}^{2}~{}\phi_{d}$, where $\phi_{i}$ denotes the throughput of
the lenses ($l$), DMD ($DMD$), gratings ($g$) and detector ($d$) (by detector
throughput we mean quantum efficiency). Given that lenses and gratings may
have negligible losses, the largest penalty is the DMD and detector quantum
efficiency. In our design, we achieve a theoretical throughput of 34%
($\phi_{g}=85\%$, $\phi_{d}=70\%$, $\phi_{\text{\tiny DMD}}=70\%$). A higher
throughput value (50%) is expected if operating in the near-IR as optical
elements are commercially available. Note, however, while the sensitivity
might be higher, the overall signal will be lower as the Raman cross-section
scales as $\lambda^{-4}$. An alternative grating-based approach is to fold the
design (back-reflect the light into the same grating)Davis  _et al._ (2011),
however this would be even more lossy as the DMD has the lowest throughput in
the design. Finally, the lens-based approach may reach overall throughput of
50% after meticulous design of the optical elements.

## Conclusion

High-sensitivity Raman imaging means detecting faint signals over strong idle
backgrounds. To address that, we have presented two inexpensive layouts of a
high-throughput spectrometer for single-pixel compressive Raman
microspectroscopy, and demonstrated its application with two different levels
of a priori information. Such layout allowed for high-sensitivity imaging of
few-bilayers vesicle in few seconds, therefore paving the way of compressive
Raman for cell biology applications.

## Acknowledgements

We thank Tom Sperber for fruitful discussions. H.B.A. was supported by LabEX
ENS-ICFP: ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. This works was supported
by European Research Council (ERC) (724473), Universitat Jaume I
(PREDOC/2013/32), Generalitat Valenciana (PROMETEO/2016/079), Ministerio de
Economia y Competitividad (MINECO, FIS2016-75618-R), Agence Nationale de la
Recherche (ANR) France Bio Imaging (ANR-10-INSB-04-01) and France Life Imaging
(ANR-11-INSB-0006) infrastructure networks and Plan cancer INSERM (PC201508).
S.G. is a member of the Institut Universitaire de France.

## References

  * Kong  _et al._ (2015) K. Kong, C. Kendall, N. Stone, and I. Notingher, Advanced drug delivery reviews 89, 121 (2015).
  * Harmsen  _et al._ (2015) S. Harmsen, R. Huang, M. A. Wall, H. Karabeber, J. M. Samii, M. Spaliviero, J. R. White, S. Monette, R. O’Connor, K. L. Pitter, S. A. Sastra, M. Saborowski, E. C. Holland, S. Singer, K. P. Olive, S. W. Lowe, R. G. Blasberg, and M. F. Kircher, Science Translational Medicine 7, 271ra7 (2015), http://stm.sciencemag.org/content/7/271/271ra7.full.pdf .
  * Jermyn  _et al._ (2015) M. Jermyn, K. Mok, J. Mercier, J. Desroches, J. Pichette, K. Saint-Arnaud, L. Bernstein, M.-C. Guiot, K. Petrecca, and F. Leblond, Science Translational Medicine 7, 274ra19 (2015), http://stm.sciencemag.org/content/7/274/274ra19.full.pdf .
  * Lorenz  _et al._ (2017) B. Lorenz, C. Wichmann, S. Stöckel, P. Rösch, and J. Popp, Trends in microbiology 25, 413 (2017).
  * Bodelón  _et al._ (2016) G. Bodelón, V. Montes-García, V. López-Puente, E. H. Hill, C. Hamon, M. N. Sanz-Ortiz, S. Rodal-Cedeira, C. Costas, S. Celiksoy, I. Pérez-Juste, _et al._ , Nature materials 15, 1203 (2016).
  * Okada  _et al._ (2012) M. Okada, N. I. Smith, A. F. Palonpon, H. Endo, S. Kawata, M. Sodeoka, and K. Fujita, Proceedings of the National Academy of Sciences 109, 28 (2012).
  * Klein  _et al._ (2012) K. Klein, A. M. Gigler, T. Aschenbrenner, R. Monetti, W. Bunk, F. Jamitzky, G. Morfill, R. W. Stark, and J. Schlegel, Biophysical journal 102, 360 (2012).
  * Ando  _et al._ (2015) J. Ando, M. Kinoshita, J. Cui, H. Yamakoshi, K. Dodo, K. Fujita, M. Murata, and M. Sodeoka, Proc. Natl. Acad. Sci. U. S. A. 112, 4558 (2015), http://www.pnas.org/content/112/15/4558.full.pdf .
  * Donaldson Jr and de Aguiar (2018) S. H. Donaldson Jr and H. B. de Aguiar, The journal of physical chemistry letters 9, 1528 (2018).
  * Shipp  _et al._ (2017) D. W. Shipp, F. Sinjab, and I. Notingher, Advances in Optics and Photonics 9, 315 (2017).
  * Duarte  _et al._ (2008) M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. E. Kelly, R. G. Baraniuk, _et al._ , IEEE Signal Processing Magazine 25, 83 (2008).
  * Berto  _et al._ (2017) P. Berto, C. Scotte, F. Galland, H. Rigneault, and H. B. de Aguiar, Opt. Lett. 42, 1696 (2017).
  * Chan  _et al._ (2008) W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, Applied Physics Letters 93, 121105 (2008).
  * August  _et al._ (2013) Y. August, C. Vachman, Y. Rivenson, and A. Stern, Applied optics 52, D46 (2013).
  * Davis  _et al._ (2011) B. M. Davis, A. J. Hemphill, D. Cebeci Maltas, M. A. Zipper, P. Wang, and D. Ben-Amotz, Analytical chemistry 83, 5086 (2011).
  * Galvis-Carreño  _et al._ (2014) D. F. Galvis-Carreño, Y. H. Mejía-Melgarejo, and H. Arguello-Fuentes, Dyna 81, 116 (2014).
  * Thompson  _et al._ (2017) J. V. Thompson, J. N. Bixler, B. H. Hokr, G. D. Noojin, M. O. Scully, and V. V. Yakovlev, Optics Letters 42, 2169 (2017).
  * Wilcox  _et al._ (2012) D. S. Wilcox, G. T. Buzzard, B. J. Lucier, P. Wang, and D. Ben-Amotz, Analytica chimica acta 755, 17 (2012).
  * Wilcox  _et al._ (2013) D. S. Wilcox, G. T. Buzzard, B. J. Lucier, O. G. Rehrauer, P. Wang, and D. Ben-Amotz, Analyst 138, 4982 (2013).
  * Wei  _et al._ (2016) D. Wei, S. Chen, Y. H. Ong, C. Perlaki, and Q. Liu, Optics Letters 41, 2783 (2016).
  * Réfrégier  _et al._ (2018) P. Réfrégier, C. Scotté, H. B. de Aguiar, H. Rigneault, and F. Galland, J. Opt. Soc. Am. A 35, 125 (2018).
  * Scotté  _et al._ (2018) C. Scotté, H. B. de Aguiar, D. Marguet, E. Green, P. Bouzy, S. Vergnole, P. Winlove, N. Stone, and H. Rigneault, Anal. Chem. 90, 7197 (2018).
  * Corcoran (2018) T. C. Corcoran, Applied spectroscopy 72, 392 (2018).
  * Pavillon and Smith (2016) N. Pavillon and N. Smith, Optics Express 24, 30038 (2016).
  * Wadduwage  _et al._ (2017) D. N. Wadduwage, V. R. Singh, H. Choi, Z. Yaqoob, H. Heemskerk, P. Matsudaira, and P. T. So, Optica 4, 546 (2017).
  * Zhang  _et al._ (2018) S. Zhang, Z. Song, G. D. P. Godaliyadda, D. H. Ye, A. U. Chowdhury, A. Sengupta, G. T. Buzzard, C. A. Bouman, and G. J. Simpson, Analytical chemistry 90, 4461 (2018).
  * Rieger  _et al._ (2017) S. Rieger, D. Grill, V. Gerke, and C. Fallnich, Journal of Raman Spectroscopy 48, 1264 (2017), jrs.5205.
  * Bélanger  _et al._ (2009) E. Bélanger, S. Bégin, S. Laffray, Y. D. Koninck, R. Vallée, and D. Côté, Opt. Express 17, 18419 (2009).
  * Gasecka  _et al._ (2017) P. Gasecka, A. Jaouen, F.-Z. Bioud, H. B. De Aguiar, J. Duboisset, P. Ferrand, H. Rigneault, N. K. Balla, F. Debarbieux, and S. Brasselet, Biophys. J. 113, 1520 (2017), http://biorxiv.org/content/early/2017/02/03/105965.full.pdf .
  * Soldevila  _et al._ (2018) F. Soldevila, J. Dong, E. Tajahuerce, S. Gigan, and H. B. de Aguiar, Submitted (2018).
  * Decker (1971) J. A. Decker, Applied optics 10, 510 (1971).
  * DeVerse  _et al._ (2000) R. DeVerse, R. Hammaker, and W. Fateley, Applied Spectroscopy 54, 1751 (2000).
  * Bioucas-Dias and Figueiredo (2007) J. M. Bioucas-Dias and M. A. Figueiredo, IEEE Transactions on Image processing 16, 2992 (2007).
  * Rajwade  _et al._ (2013) A. Rajwade, D. Kittle, T.-H. Tsai, D. Brady, and L. Carin, SIAM Journal on Imaging Sciences 6, 782 (2013).
  * Moreaux  _et al._ (2000) L. Moreaux, O. Sandre, M. Blanchard-Desce, and J. Mertz, Opt. Lett. 25, 320 (2000).
  * Müller and Schins (2002) M. Müller and J. M. Schins, The Journal of Physical Chemistry B 106, 3715 (2002), http://pubs.acs.org/doi/pdf/10.1021/jp014012y .
  * Potma and Xie (2003) E. O. Potma and X. S. Xie, J. Raman Spectrosc. 34, 642 (2003).
  * Li  _et al._ (2005) L. Li, H. Wang, and J.-X. Cheng, Biophysical journal 89, 3480 (2005).
  * Cheng  _et al._ (2003) J. Cheng, S. Pautot, D. Weitz, and X. Xie, Proc. Natl. Acad. Sci. U. S. A. 100, 9826 (2003).
  * Chen  _et al._ (2010) B.-C. Chen, J. Sung, and S.-H. Lim, The Journal of Physical Chemistry B 114, 16871 (2010), http://pubs.acs.org/doi/pdf/10.1021/jp104553s .
  * Shen  _et al._ (2017) Y. Shen, Z. Zhao, L. Zhang, L. Shi, S. Shahriar, R. B. Chan, G. Di Paolo, and W. Min, Proceedings of the National Academy of Sciences 114, 13394 (2017).
  * Rehrauer  _et al._ (2018) O. G. Rehrauer, V. C. Dinh, B. R. Mankani, G. T. Buzzard, B. J. Lucier, and D. Ben-Amotz, Applied spectroscopy 72, 69 (2018).
  * Kong  _et al._ (2016) L. Kong, M. Navas-Moreno, and J. W. Chan, Analytical Chemistry 88, 1281 (2016), pMID: 26654100, http://dx.doi.org/10.1021/acs.analchem.5b03707 .
  * Müller  _et al._ (2016) W. Müller, M. Kielhorn, M. Schmitt, J. Popp, and R. Heintzmann, Optica 3, 452 (2016).
  * Schlücker  _et al._ (2003) S. Schlücker, M. D. Schaeberle, S. W. Huffman, and I. W. Levin, Analytical Chemistry 75, 4312 (2003).


# The data reduction pipeline for the Hi-GAL survey

A. Traficante1,9, L. Calzoletti2,8, M. Veneziani3,9, B. Ali5, G. de Gasperis1,
A.M. Di Giorgio6, F. Faustini2, D. Ikhenaode4, S. Molinari6, P. Natoli1,2,7,
M. Pestalozzi6, S. Pezzuto6, F. Piacentini3, L. Piazzo4, G. Polenta2,8 and E.
Schisano6
1Dipartimento di Fisica, Università di Roma “Tor Vergata”, Italy
2ASI Science Data Center, I-00044 Frascati (Rome)
3Dipartimento di Fisica, Università di Roma “La Sapienza”, Italy
4DIET - Dipertimento di Ingegneria dell’ Informazione, Elettronica e
Telecomunicazioni, Università di Roma “La Sapienza”, Italy
5Nasa Herschel Science Center, Caltech, Pasadena, CA
6INAF-Istituto Fisica Spazio Interplanetario I-00133 Rome
7INFN Sezione di Tor Vergata, Rome, Italy
8INAF, Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monte Porzio
Catone, Italy
9Spitzer Science Center, Caltech, Pasadena, CA E-mail:
<EMAIL_ADDRESS>
(Submitted)

###### Abstract

We present the data reduction pipeline for the Hi-GAL survey. Hi-GAL is a key
project of the Herschel satellite which is mapping the inner part of the
Galactic plane ($|l|\leqslant$ $70^{\circ}$ and $|b|\leqslant 1^{\circ}$),
using 2 PACS and 3 SPIRE frequency bands, from $70$$\mu$m to $500$$\mu$m. Our
pipeline relies only partially on the Herschel Interactive Standard
Environment (HIPE) and features several newly developed routines to perform
data reduction, including accurate data culling, noise estimation and minimum
variance map-making, the latter performed with the ROMAGAL algorithm, a deep
modification of the ROMA code already tested on cosmological surveys. We
discuss in depth the properties of the Hi-GAL Science Demonstration Phase
(SDP) data.

###### keywords:

instrumentation – data reduction

††pagerange: The data reduction pipeline for the Hi-GAL
survey–References††pubyear: 2002

##  1 Introduction

The Herschel Space Observatory was launched from Kourou in May 2009 aboard an
Ariane 5 rocket. Two of three scientific instruments on the focal plane (PACS
and SPIRE) are capable to observe the infrared sky with unprecedented angular
resolution and sensitivity, providing photometric observations in 6 different
bands ($70$$\mu$m, $100$$\mu$m, $160$$\mu$m, $250$$\mu$m, $350$$\mu$m and
$500$$\mu$m: Pilbratt et al., 2010, and reference therein).

The PACS photometer is composed of two bolometer arrays: a $64\times 32$ pixel
matrix arranged from 8 monolithic subarrays of $16\times 16$ pixels each
centered on the $70\mu$m and $100\mu$m wavelength (blue and green bands), and
a $32\times 16$ pixel matrix organized in two subarrays for the band centered
on $160\mu$m (red band), see Poglitsch et al. (2010).

SPIRE comprises a three band photometer, operating in spectral bands centered
on $250\mu$m, $350\mu$m and $500\mu$m. Every band uses a matrix of germanium
bolometers (139, 88 and 43 respectively) coupled to hexagonally packed conical
feed horns (Griffin et al., 2010).

In order to handle science data provided by the Herschel instruments,
including the data retrieval from the Herschel Science Archive, the data
reduction through the standard pipelines and the scientific analysis, an
official software environment called HIPE (Herschel Interactive Processing
Environment, Ott et al., 2010) is available from ESA.

The raw data provided by the satellite are reduced in HIPE to generate
scientific data (so-called Level 2) and intermediate products of the data
reduction process (Level 0 and Level 1 data).

In this paper, we describe the dedicated pipeline created to obtain maps for
Hi-GAL (Herschel Infrared Galactic Plane Survey, Molinari et al., 2010a). Hi-
GAL aims to homogeneously cover with observations in 5 contiguous IR bands
between 70$\mu$m and 500$\mu$m a 2 degrees wide stripe of galactic plane
between $l=-70^{\circ}$ and $l=70^{\circ}$.

The Galactic plane shows emission varying from point-like sources to large-
scale structures and with intensity varying over a wide dynamic range. In this
work we show that the existing standard reduction strategy (which is based on
HIPE version 4.4.0, released in November, 11th 2010) is not optimized to
reduce Hi-GAL data and that a dedicated pipeline can enhance the quality of
the Level 2 products.

After Herschel successfully passed the Performance Verification Phase (PV
Phase), two fields of the Hi-GAL survey were acquired during the Science
Demonstration Phase (SDP): 2x2 square degree areas of the Galactic plane
centered on 30∘ of longitude (hereafter, $\textit{l}=30^{\circ}$) and on 59∘
($\textit{l}=59^{\circ}$).

We describe the data reduction tools used to obtain high quality maps from SDP
data, with the aim to provide a reliable environment for the Routine Phase
(RP) data. The maps provided by our pipeline are successfully used for several
works like, e.g. , Molinari et al. (2010b), Martin et al. (2010), Peretto et
al. (2010).

The paper is organized as follows: in Section 2 we describe the acquisition
strategy for Hi-GAL data; in Section 3 we describe the pre-processing steps of
the data reduction pipeline, necessary to prepare data for the map making and
the tools that we have developed to that purpose. In Section 4 we describe the
ROMAGAL map making algorithm used in order to estimate the maps. ROMAGAL is
used in place of the MadMap code which is the map making algorithm offered in
HIPE. The quality of the ROMAGAL maps for both PACS and SPIRE instruments,
related to the SDP observations, will be analyzed in Section 5; in Section 6
we draw our conclusions.

##  2 The Hi-GAL acquisition strategy

Hi-GAL data are acquired in PACS/SPIRE Parallel
mode111http://herschel.esac.esa.int/Docs/PMODE/html/parallel_om.html, in which
the same sky region is observed by moving the satellite at a constant speed of
60$\arcsec$/sec and acquiring images simultaneously in five photometric bands:
$70$$\mu$m and $160$$\mu$m for PACS and $250$$\mu$m, $350$$\mu$m and $500\mu$m
for SPIRE.

The whole data acquisition is subdivided in $2^{\circ}\times 2^{\circ}$ fields
of sky centered on the Galactic plane. Every Hi-GAL field is composed of the
superposition of two orthogonal AOR (Astronomical Observation Requests). Each
of them is based on a series of consecutive, parallel and partly overlapped
scan legs covering $2^{\circ}\times 2^{\circ}$ square degrees. The scanning
strategy adopted for Hi-GAL is fully described in (Molinari et al., 2010a).
The superposition is performed in order to obtain a suitable data redundancy
and for better sampling the instrumental effect like the high-frequency
detector response (Molinari et al., 2010a).

The acquisition rate for the parallel mode is 40 Hz for PACS and 10 Hz for
SPIRE, although the PACS data are averaged on-board for an effective rate of 5
Hz and 10 Hz for the $70$$\mu$m and $160$$\mu$m array respectively. The
implications of the PACS data compression are detailed in Section 5.3.

An example of the scanning strategy of the Hi-GAL survey is shown in Figure 1.
The coverage map of the PACS blue array is shown on the left panel and on the
right panel we highlight the superposition of one scan leg to the following
ones, by enlarging the bottom right corner of the left image. Two calibration
blocks for each AOR, during which the instrument observes two internal
calibration sources located on the focal plane, were provided during the SDP
observations. They appear as higher than mean coverage areas and are marked by
black and green circles for the 2 AORs in Figure 1. Higher coverage zones are
also clearly visible in the slewing region at the end of each scan leg, where
the telescope decelerates and then accelerates before initiating the next scan
leg.

Figure 1: Left: the coverage map of PACS blue array. Right: a zoom of the
bottom right corner where is clear the effect of the superposition from one
scan-leg to the next. Black and green circles highlight the calibration blocks
that, during PV and SDP phase, were observed twice during which the internal
calibration sources are observed.

##  3 Map making pre-processing

The aim of the pre-processing is to prepare Hi-GAL data for mapmaking.

While map making is performed using a Fortran parallel code borrowed from
cosmological observations (see Section 4), the preprocessing is done through a
series of IDL and jython tools to be run on the data one after the other.
After having tried to map Hi-GAL data using the standard tools provided within
HIPE, we decided to develop our own routines that we tailored specifically for
the reduction of data affected by bright and irregular background, as in the
Galactic plane. In fact, high-pass filtering used in HIPE to cure the long-
time drift also removes a large part of the diffuse Galactic emission.
Furthermore, the standard deglitching embedded in HIPE , the MMT
(Multiresolution Median Transform, Starck et al., 1995), generates false
detections in correspondence to the very bright sources when we apply this
task to the PACS data. On the other hand, the deglitching procedure based on
wavelet analysis used by SPIRE does not affect the data, given also the lower
spatial resolution compared to the PACS one. We therefore use the HIPE task
for SPIRE data only.

Herschel data are stored in subsequent snapshots of the sky acquired by the
entire detector array, called frames. In a first processing step, the standard
modules within HIPE are used to generate Level 1 data for both PACS (except
the deglitching task) and SPIRE. Thus, the data are rearranged into one time
series per detector pixel, called Time Ordered Data (TOD), in which the
calibrated flux ( Jy beam-1 for SPIRE and Jy sr-1 for PACS) and the celestial
coordinates are included. At the end of this step, TODs are exported outside
HIPE in fits format. In the subsequent processing steps, TODs are managed by a
series of IDL tools, in order to produce final TODs free of systematic effects
due to electronics and of glitches and corrupted chunks of data due to cosmic
rays. To each TOD a flag file is attached to keep track of any flagging done
during data reduction steps.

Preprocessing includes identification of corrupted TODs (or part of them),
drift removal and deglitching. The following steps will be the Noise
Constraint Realization (NCR) and ROMAGAL mapmaking; they will be both
described in detail in the next Sections. The summary of the entire pipeline
is shown in the diagram in Figure 2.

Figure 2: Schematic representation of the Hi-GAL data reduction pipeline

###  3.1 Corrupted TODs

A TOD can be partially corrupted by the random hiting of charged particles
(cosmic rays) which produce strong and spiky signal variations called
glitches. Two different effects on the TODs can be identified: the glitch
corrupts a single or few consecutive samples, generating spiky changes along
the TOD. This is the most common effect and in Section 3.3 we describe how to
detect and mask the data for PACS, as well as to mask any possible residual
glitches for SPIRE. Very powerful glitches, on the other hand, make the
detector signal unstable for a considerable interval of time (see, e.g.,
Billot et al., 2010). This effect depends on the glitch amplitude and on the
bolometers time response. These events affect a much larger part of the TOD
that cannot be used for mapmaking.

Their impact results in a bias on the bolometer timeline with a low-frequency
drift which can involve a considerable number of samples, as shown in Figure
3.

In that Figure, blue crosses represent the observed timeline of one bolometer
of the blue array.

Automatic identification of the (partially) corrupted TODs exploits the first
derivative of the TOD to detect extraordinary “jumps” in the signal. In order
to determine what portion of the TOD is to be flagged, the flagging tool
exploits the fact that the detector pixels response that have been hit by a
cosmic ray is mostly exponential. Data samples ranging from the jump to the
time sample at which an exponential fit reaches 98% of the signal level before
the event are identified as bad data and stored as such in the corresponding
flag file. In case the exponential fit does not reach 98% of the starting
value before the end of the TOD, then all data starting from the hit will be
flagged as is the case in Figure 3. This procedure is applied both in the
cases of a changing in responsivity or a detector offset alteration. In the
latter, we estimate the fit with an exponent equal to 0. The described
procedure was adopted to process both PACS and SPIRE data.

Figure 3: Timeline of a PACS blue bolometer. The exponential decay illustrates
the change in responsivity after frame 40000 due to the impact of a powerful
glitch.

###  3.2 Drift removal

After having identified corrupted data we proceed to the elimination of
changes in responsivity over time. The procedures are in principle identical
for PACS and SPIRE data, the only differences account for different
composition of the detector arrays and the data acquisition of the two
instruments.

The signal in PACS TODs that exit HIPE does not represent the true sky but is
dominated by the telescope background and the (unknown) zero level of the
electronics. The electronics further introduce significant pixel-to-pixel
offsets. For each TOD, we mitigate the effect of pixel-to-pixel offset by
calculating and then subtracting the median level for each pixel from each
readout. This ensures that all pixels have median level equal to 0. The median
value is preferred over mean because the median is a much better
representation of the sky+telescope background flux, and is much less
sensitive to signal from astrophysical sources.

Also, the subtraction of this offset from each TOD does not alter the signal
in the final map, but it introduces only a global offset constant over the
entire area covered in the observation. However it should be kept in mind that
bolometers are inherently differential detectors which bear no knowledge of an
absolute signal value; besides, any optimized map making methods like the one
we employ (see Section 4) produce maps with an unknown offset value which
needs to be independently calibrated. So it is important to reduce all the
bolometers to the same median value, regardless of its amount. The pixel-to-
pixel median subtraction has the effect seen in Figure 4. Diffuse emission and
compact sources are clearly visible in the frame.

Figure 4: Blue PACS frame after the median subtraction on each pixel. Diffuse
emission and compact source are visible in the frame.

Still, when plotting a detector pixel timeline we see that the signal
decreases systematically from the start to the end of the observation. This
trend is likely due to a combination of small changes in the thermal bath of
the array and to small drifts in the electronics. The former affects the
entire array, while the latter affects subunits of the detector array (in
fact, PACS blue is divided into 8 units, PACS red into 2, electronically
independent units (Poglitsch et al., 2008)). The drift is then a combination
of these two effects: drifts of the entire array and drift of a single
subunit. These effects are dominant with respect to the $1/f$ noise pattern,
which will be described in the next Sections.

Figure 5: Median behavior computed on the whole array for each frame. The
slow-drift behavior is due to the electronics and the thermal bath.

It is in principle not a trivial task to decide which drift has to be
subtracted first: the drift from the thermal bath (affecting the entire
detector array) or the drift from the readout electronics (affecting sub-
arrays differently)? Ideally both should be subtracted, if only it were
possible to separate each component accurately, as the net effect in the data
is the sum of both.

Our methodology for removing the correlated signal drifts (on both the
bolometer module/unit level and the array level) is based on tracing the low
signal envelope of the unit or array median levels. In Figure 5, this envelope
is the curve defined by the lowest signal values. and estimated as follows:

  * i

We compute the median value of the entire bolometer array/unit for each
bolometer readout. Figure 5 shows one example for the entire array.

  * ii

The median values thus obtained are segmented and grouped by scan legs. Each
scan leg is composed of $\sim 1000$ frames and we observed 54 scan leg for
each $2^{\circ}\times 2^{\circ}$ Hi-GAL field.

  * iii

For each scan leg we compute the minimum value of the array/unit medians.

  * iv

The resulting set of minimum median values for all scan legs are fit with a
polynomial.

The median value for each array/unit readout is chosen because it is closest
to the actual sky+telescope background. However, as clearly seen in Figure 5,
in the presence of strong astrophysical sources the median value is incorrect
for our purposes. The strong sources appear as signal spikes in Figure 5.
Hence, we take the next step of finding the minimum value from the set of
medians belonging to a single scan leg. The idea is that at some point during
the scan leg the median was indeed a true representation of the local
sky+telescope and relatively free of source emission. This step allows us to
reject the sources at the expense of degrading our time-resolution to scan-leg
duration ($\sim 240$ sec). The polynomial fit allows us to estimate the drift
behavior at the same time resolution as the PACS signal (5Hz and 10Hz for
70$\mu$m and160$\mu$m band respectively). We further note that the correlated
signal drift is relatively flat over a single scan leg; hence, the minimum
value is not significantly affected by the presence of the monotonic signal
drift itself in the first place.

The minimum median method discussed above removes background levels from
spatial emission structures that are of the order or larger than the scan legs
yet preserves the spatial structures smaller than scan legs. In essence,
information about the absolute calibration zero-point (Miville-Deschênes &
Lagache, 2005) is lost but all spatial structures within our map boundaries
are preserved.

In Figure 6 we reported the minimum median values of each subarray. The common
downward trend is due to the (common) thermal bath drift, while the different
readout electronics are responsible for the differences in the subarray
slopes.

We therefore decide to subtract the subarray drifts in order to consider both
the thermal bath and the readout electronics behaviors, but separately for
each subunit.

Figure 6: Interpolation of the minima of the median evaluated on every scan
leg. Each curve refers to a PACS subarray. The curves mimic the same behavior
with slopes due to the different subarray electronics.

Once the individual subarray drift is removed, the remaining dispersion on the
whole array is only a residual scatter due to the intrinsic efficiency of the
removal tool, as shown in Figure 7.

Figure 7: Minima of the median evaluated on the whole PACS blue array after
the subarray drift subtraction. The dispersion is due to the intrinsic
efficiency of the drifts subtraction tool. There is no residual behavior and
the scatter is one order of magnitude under the value of the initial drift.

SPIRE array detectors are not divided into subarrays, so every procedure that
has to be run 8 or 2 times for PACS data is only performed once per SPIRE
band. SPIRE uses blind bolometers (5 in total) as thermistors to evaluate the
most relevant correlated noise component: the bath temperature fluctuations. A
standard pipeline module uses this information to perform an effective removal
of the common drift present along the scan observation. HIPE also corrects for
the delay between the detector data and the telescope position along the scan,
using an electrical low-pass filter response correction. But despite these
initial (and very effective) corrections, we apply the drift removal tool to
SPIRE data in the same way as for PACS data: we fit a polynomial to the
minimum of the median of each scan leg (calculated over the entire detector
array), that we then subtract to all TODs. Experience on the data shows that a
residual long time drift is often present in SPIRE data.

Finally, when removing drifts it is important to know how the observational
scans are oriented. In fact, as the Galactic plane is very bright, scans
across the plane will give rise to an increase of signal, on top of the
general drift. On the other hand, when the scan is almost parallel to the
plane of the Galaxy, the signal can dominated by its bright emission, also on
the evaluation of the minima of the median.

In this case, the curve to fit is estimated only in the scan legs where the
median is not affected by the signal.

Since the procedure is not automatic, care has to used when choosing what
polynomial to fit and subtract from the data, in order not to eliminate
genuine signal from the sky. From our experience the best choice is a first or
a second degree polynomial, depending on the signal behavior observed.

Higher polynomial degree can be necessary when part of the drift has to be
extrapolated in order to avoid signal contamination.

An example for one subarray is shown in Figure 8.

Figure 8: Blue curves: interpolation of the minima of the median for one
subarray of PACS blue channel when the scan direction is almost parallel to
the Galactic plane. Red line: the fit that we choose to evaluate the drift
without considering the minima affected by the Galactic emission (central
bump).

Each scan leg is bigger then the overlapping region of the two scan directions
for the scanning strategy adopted by Hi-GAL (see Figure 1). Since the Hi-GAL
fields are squared regions, the slowly-traversed direction of the AOR within
the overlapping region have a length comparable with the scan leg. Thus, we
assume that even if there is a signal gradient along the slowly-traversed
direction of the AOR, it is not filtered out by the array medians subtraction.

###  3.3 Deglitching

To remove outliers in the time series of the bolometers we exploit the spatial
redundancy provided by the telescope movement which ensures that each sky
pixel of the final map is observed with different bolometers. Outliers
detection is done with the standard sigma-clipping algorithm: given a sample
of $N$ values, first estimates for the mean and for the standard deviations
are derived; then all the values that differ from the mean by more than $n$
standard deviations are considered outliers and removed.

For this algorithm the choice of $n$, the parameter that defines the threshold
above which a value is considered an outlier, is usually arbitrary: a certain
$n$ is chosen, very often equal to 3, without any statistical justification.
Recently Pezzuto (2010) has derived a formula that starting from the
properties of the error function for a Gaussian distribution and exploiting
the discreteness of the number of available measures, relates $n$ to the size
of the sample. The formula is

$n=-0.569+\sqrt{-0.072+4.99\log(N)}$ (1)

As a consequence, in the central region of the map, where the coverage (and so
$n$) is high, the number of standard deviation is larger than in the outskirt
of the map where the coverage is low. For instance, if a sky pixel has been
observed with 40 bolometers, the above formula gives $n=2.25$; so, once we
have estimated the mean $m$ and the standard deviations $\sigma$, all the
values $x_{i}$ such that ABS$(x_{i}-m)>2.25\sigma$ are flagged as outliers. If
a pixel has been observed with 20 bolometers the threshold lowers to
1.96$\sigma$.

This procedure is automatically iterated until outliers are no longer found.
However, the procedure converges within 1 iteration in $\sim$98% of the cases
in which we have applied the analysis.

The outliers detection is done in this way for both instruments, however for
SPIRE, as explained before, we also make use of the standard deglitching
algorithm (wavelet based) implemented in the official pipeline. But we found
some weak glitches left in the SPIRE TODs so that we decided to run our
deglitching algorithm also on SPIRE data.

The number of glitches found is on the average about 15%, a value which is
likely larger than the real percentage. For PACS we are now working on a
different way to associate each bolometer to the sky pixels, taking into
account the finite size of the sky pixels. For the first test cases we run,
the percentage of detected glitches is now around 5-6%.

##  4 The ROMAGAL Map Making algorithm

The ROMAGAL algorithm is based on a Generalized Least Square (GLS) approach
(Lupton, 1993). Since the TOD is a linear combination of signal and noise, we
can model our dataset $\textbf{d}_{k}$ for each detector $k$ as (Wright,
1996):

$\textbf{d}_{k}=P\textbf{m}+\textbf{n}_{k}$ (2)

where $P$ is the pointing matrix, which associates to every sample of the
timeline a direction in the sky, m is our map estimator of the “true” sky and
$\textbf{n}_{k}$ is the noise vector.

The observed sky, $P\textbf{m}$, is the map estimator of the “true sky”
convolved with the instrumental transfer function and the optical beam.
However, in case of circularly symmetric beam profile, $\mathbf{m}$ is a beam
smeared, pixelised image of the sky.

In this case the pointing matrix has only a non-zero entry per row
corresponding to the sky pixel observed at a given time. Since the beam
profiles for PACS (Poglitsch et al., 2008) and SPIRE (Griffin et al., 2008)
are only weakly asymmetric we can put ourselves in this simple case. Note that
the transpose of the $P$ operator performs a data binning (without averaging)
into the sky pixels.

Equation 2 holds only if the noise vector of each detector $\textbf{n}_{k}$ is
composed of statistical random noise, with Gaussian distribution and null
average. All the relevant systematic effects (offset, glitches) have then to
be removed with an accurate data preprocessing before map production, as
explained in Section 3.

The formalism can be easily extended in the case of multidetector analysis. In
this case the vector d contains the data relative to each detector. Rather,
one has to take care to upgrade also the noise vector n, accordingly with the
correct noise value for each detector.

The GLS algorithm produces minimum noise variance sky maps. Noise properties
for each detector have to be previously estimated and provided in input to the
algorithm as described in Section 4.1.

The GLS estimate for the sky, $\tilde{\textbf{m}}$, is (Natoli et al., 2001)

$\tilde{\textbf{m}}=(P^{T}\textbf{N}^{-1}P)^{-1}P^{T}\textbf{N}^{-1}\textbf{d}$
(3)

where $\textbf{N}=\langle\textbf{nn}^{T}\rangle$ is the noise covariance
matrix, which takes into account noise time correlation between different
samples. Such correlation is particularly high at low frequencies because of
the $1/f$ (or long memory) noise. In case of uncorrelated noise (or white
noise) the N matrix becomes diagonal and the problem is greatly simplified. If
we further assume stationary uncorrelated noise, Equation 3 reduces to:

$\tilde{\textbf{m}}=(P^{T}P)^{-1}P^{T}\textbf{d}.$ (4)

$P^{T}P$ is the number of observations of a pixel of the map, so we are
averaging the different TOD values into that pixels assigning the same weight
to each sample. We will refer to this map estimate as “naive” or “binned” in
the following.

When non negligible noise correlation is present, as in the case of PACS
(Poglitsch et al., 2008) and SPIRE (Schulz et al., 2008), Equation 3 must be
solved. This is a challenging computational task since it requires, in
principle, the inversion of the large (of the order of the number of pixels in
the map) matrix $P^{T}\textbf{N}^{-1}P$, which is the covariance matrix of the
GLS estimator (Lupton, 1993). One key simplifying assumption is to invoke that
the noise is stationary. In this case, the N matrix has a Toeplitz form which
can be approximately treated as circulant, ignoring boundary effects (Natoli
et al., 2001). A circulant matrix is diagonal in Fourier space and its inverse
is also circulant, so the product between $\textbf{N}^{-1}$ and a vector is a
convolution between the same vector and a filter provided by any of the rows
of the matrix. In the following we will refer to any of these rows as a noise
filter. Its Fourier transform is the inverse of the noise frequency power
spectrum.

Considering the conditions listed above, the GLS map making algorithm performs
the following operations, starting with rewriting the Equation 3 in the form

$(P^{T}\textbf{N}^{-1}P)\textbf{m}_{0}-P^{T}\textbf{N}^{-1}\textbf{d}=\textbf{r}$
(5)

where $\textbf{m}_{0}$ is the starting map used at the first iteration,
generally the naive map.

The $P\textbf{m}_{0}$ product projects the map onto a timeline. Application of
$\textbf{N}^{-1}$: this is a convolution which can be performed in Fourier
space. Application of $P^{T}$: this step projects the convolved timeline back
into a map.

The second term performs the convolution with the filter (applying
$\textbf{N}^{-1}$ to the data vector d in Fourier space) and then the
projection of the convolved timeline into a map (applying $P^{T}$ to the
product $\textbf{N}^{-1}\textbf{d}$).

Then, we need to evaluate the residual r. If the residual is higher than a
fixed threshold, it is used to produce a new map, $\textbf{m}_{1}$, as
described in Hestenes & Stiefel (1952). This map will be considered instead of
$\textbf{m}_{0}$ for evaluating again the Equation 5, until convergence. This
is achieved by running a Conjugate Gradient algorithm, an iterative method
useful to obtain a numerical solution of a system (Hestenes & Stiefel, 1952),
until convergence is reached with the residual lower then the threshold.

The algorithm outlined is the same described in (“unroll, convolve and bin”:
Natoli et al., 2001) and is implemented in the ROMAGAL code. The next section
explains the strategy employed to estimate the noise filters used by ROMAGAL
directly from in-flight data.

###  4.1 Noise estimation

In order to estimate the noise filters for ROMAGAL we need to investigate the
noise statistical properties in the timelines. Data are mostly affected by two
kind of statistical noise: $1/f$ noise due both to the electronics and thermal
background radiation from the telescope or the instruments, and photon noise
(see Poglitsch et al. 2008; Schulz et al. 2008).

The detector $1/f$ noise arises in the electronic chain and it impacts
particularly regions with low signal-to-noise ratio (SNR), where only diffuse
emission is present. In those regions it can be of the same order of magnitude
of the signal or even higher. In these cases GLS treatment is particularly
effective.

Photon noise is due to statistical fluctuation in the photon counts. This
process follows poissonian statistic, so the SNR is proportional to the square
root of counts. Since poisson distribution tends to be Gaussian for large
numbers, we can approximate photon noise as Gaussian on the map if the number
of counts is large enough.

Since bolometers are organized in matrix and sub-matrix, the signal of a
bolometer can be correlated with the signal of another, generally adjacent,
bolometer. These effects could be both statistical and deterministic. We
already described how to remove the deterministic common mode from TOD (like
the thermal bath variations, see Section 3).

One possible source of statistical cross-correlated noise is the crosstalk
between bolometers: the signal of one pixel may contaminate the signal of its
neighbors through the capacitive inductive coupling, generating a common mode
called “electrical crosstalk”. On the contrary, “optical crosstalk” is due to
diffraction or aberrations in the optical system which could drive an
astronomical source to fall on inappropriate detectors (Griffin et al., 2008).

We then analyze the residual contribution of the statistical component of the
correlated noise. We found that the residual correlated noise level in each
pixel is negligible with respect to the intrinsic detector noise level for
both PACS and SPIRE instruments, as described in the following.

In principle, the noise properties vary significantly across the array and we
had to estimate the noise power spectrum for each bolometer. To do that we
have processed “blank” sky mode (i.e. filled with negligible contribution from
sky signal) data acquired during the PV Phase.

In Figure 9 we show a typical noise spectrum estimated for a pixel of the
$160\mu$m PACS band (black) and the cross-spectrum between two adjacent
bolometers (red). The cross-spectrum evaluates the impact of the cross-
correlated noise in the frequency domain between two different bolometers. The
level of the cross-correlated noise is at least 4 order of magnitude below the
auto-correlated noise power spectrum of each pixel. Note that this means we do
not see any relevant cross-correlated noise, despite the fact that crosstalk
can be present into the timeline.

In Figure 10 we show noise power spectra of the 250$\mu$m SPIRE band
bolometers. Also in this case the cross-spectrum is negligible.

Noise spectra of both PACS and SPIRE display low-frequency noise excess
($1/f$). In case of SPIRE spectra (Figure 10) a high frequency rise is also
evident, which is due to the deconvolution of the bolometric response
function. PACS spectra do not show this behavior because the bolometer
transfer function is not deconvoluted by the standard pipeline.

Figure 9: Black line: typical noise spectrum of a PACS $160$$\mu$m detector,
estimated on blank sky data. Red line: cross spectrum between two detectors of
the same subarray. The level of the cross-correlated noise is significantly
under the noise level of each single bolometer, so we can reasonably neglect
it. Figure 10: Same as Figure 9 for a SPIRE $250$$\mu$m bolometer. For SPIRE
also the noise level of cross-spectrum is reasonably negligible with respect
to the auto spectrum level.

###  4.2 From ROMA to ROMAGAL

The ROMAGAL GLS code has been optimized to recover the features in the Hi-GAL
fields with high accuracy.

Hi-GAL observes the Galactic plane where the dynamic range of sources spans
over several orders of magnitudes. This poses strong constraints on the map
making algorithm: both the weak diffuse emission and the bright compact
sources in e.g., star forming regions have to be recovered with high accuracy.
The signal often exhibits steep gradients that are hard to follow for the GLS
solver, which relies on the assumption that the sky signal does not vary
significantly within a resolution element (see below 5.3). At the same time,
several systematics affect the dataset. As explained above, many of them are
cured at the preprocessing level. However, their removal generates a
conspicuous amount of transient flagging, that must be correctly handled by
the GLS code.

The core of ROMAGAL is the same of the ROMA code (de Gasperis et al., 2005)
where the input-output routines have been deeply modified to adapt to the HIPE
generated dataset. ROMAGAL inputs are the TOD generated by HIPE, pointing and
transient flagging. These have the same format for both PACS and SPIRE.
ROMAGAL outputs are fits file containing the optimal map, arranged in standard
gnomonic projection routines. The code is written in FORTRAN 95 and relies on
the MPI library for parallel calls. It runs on the Hi-GAL dedicated machine,
called BLADE, a cluster of 104 processors at 2.5 GHz each and 208 Gb RAM
total. Its nodes are interconnected with MPI-infiniBAND. The machine is
located at IFSI-Rome.

As explained in the previous Section, the computation of Equation 3 is
unfeasible due to the size of the map’s covariance matrix. However, we assume
the noise of each Hi-GAL field to be stationary to set up an FFT based solver
built upon a conjugate gradient iterative algorithm (see Section 4). Such a
scheme can estimate the final maps with a precision of order of
$\epsilon=10^{-8}$ in $\sim$ 150 iterations for Hi-GAL . ROMAGAL computational
time scales linearly with the size of the dataset and only weakly with the
number of pixels in the output maps. The scaling with the number of processors
is highly optimal in the range of cores required for the Hi-GAL analysis
($<50$). For the largest channels (PACS blue band), a final GLS map of a
$2^{\circ}\times 2^{\circ}$ field requires about 16 Gbytes of RAM and 1400 sec
on 8 cores. Due to the high number of array pixels (2048), this channel is the
largest dataset to analyze as well as the most demanding in terms of
computational resources. Further information on resource consumptions can be
found in Table 1.

Band | Total Time (sec) | RAM (Gb)
---|---|---
$70\mu$m |  $\sim$ 1400 | 16
$160\mu$m |  $\sim$ 1000 | 8
$250\mu$m |  $\sim$ 180 | 4
$350\mu$m |  $\sim$ 130 | 1
$500\mu$m |  $\sim$ 100 | 1
Table 1: Time and minimum RAM amount required from ROMAGAL for each PACS and
SPIRE band using 8 BLADE processors.

###  4.3 Optimal treatment of transient flagging

As mentioned above, the timelines are inspected for bad data samples that must
be excluded from map making as part of the preprocessing pipeline. Bad data
can arise due to a variety of reasons. They are generally caused by transient
events, either unexpected (e.g., glitches, anomalous hardware performance) or
expected (e.g., detectors saturating because of a bright source, observation
of a calibrator). Once identified, a flag record is generated and stored for
these anomalous events, so that their contribution can be safely excluded from
map making. Flags in the TOD pose a potential problem for ROMAGAL because its
solver is based on the FFT as discussed in the previous section. The FFT
requires the timeline to be continuous and uniformly sampled. Since noise in
the PACS and SPIRE data is correlated, just excising the flagged sampled to
fictitiously create a continuous timeline would interfere with noise
deconvolution, and is thus not a safe option. Instead, we advocate using a
suitable gap-filling technique. The rest of this section is mostly devoted to
defining which, among the various options for gap filling, is best suited for
the Hi-GAL maps.

It is important to realize that the output map will depend on the content of
the flagged sections of the TOD even if these values are not associated to any
(real) map pixel. This is due to the convolution performed within the solver
that “drags” out data from the flagged section of the timeline, even if, when
the data are summed back into a map, the $P$ operator is applied only to the
unflagged samples. Since one is not in control of the content of the flagged
section of the timeline, a kind of gap-filling must be performed in any case.
We have tested different recipes for gap-filling making extensive use of
simulations. We have treated separately the signal and noise components of the
timelines, running separately noise dominated and signal dominated cases,
because the behavior towards flags of the two components is different, as it
will be shown below.

The simplest form of gap filling is to remove the content of the flags
altogether, replacing the flagged sections with a nominal constant (usually
null) value. This works very well on a signal-only simulation of the Hi-GAL
field. However, it fails dramatically when noise is present, as evident from
Figure 11 (left panel), where a markedly striped pattern in the reconstructed
map is seen (in this simulation, the Hi-GAL noise has been amplified by a
factor 100 to make its pattern more evident). The reason for this failure is
readily understood: The GLS map making employed requires noise stationarity
(see Section 4.2 above), which is obviously not preserved by zeroing the gaps.
A less obvious effect is that even if gaps are filled with a noise realization
with the correct statistical properties, but unconstrained, the GLS map making
is bound to fail as well, as shown in the middle panel of Figure 11. A noise
realization is said to be constrained when it respects certain boundary
conditions (Hoffman & Ribak, 1991), which in our case are represented by the
noise behavior outside the gap. Unconstrained noise inside the gap, despite
having the correct statistical properties, creates a border discontinuity that
causes the GLS map maker to behave sub-optimally (Stompor et al., 2002). We
have employed a code to create Gaussian noise constrained (NCR) realizations,
based on the algorithm set forth by Hoffman & Ribak (1991). The code uses in
input the noise pattern and statistical properties, as measured from the
timelines. The results on noisy simulations are excellent, as set forth by the
third (rightmost) panel in Figure 11. Note, however, that Figure 11 refers to
a noise dominated simulation. We now turn to discuss the effect of a non
negligible signal contribution (a far more realistic case for Hi-GAL ).

Figure 11: Shown are the results obtained in a noise dominated regime (normal
Hi-GAL noise is amplified by a factor 100). Left panel is the map obtained
replacing flagged data samples with null values (clearly it does not work),
middle panel is the map obtained replacing data samples with unconstrained
noise realization (does not work either), right panel shows the map obtained
using our NCR code (does work).

We have verified that the presence of non negligible signal in the timelines
does not affect the results provided that the NCR is performed using the
underlying noise field as a baseline. Measuring the latter in presence of
signal is however impractical. It would be significantly simpler if the NCR
could be run directly on the timelines themselves, constraining thus the
(fake) noise realization within the gap to the outside noise plus signal
(true) values. Unfortunately, this poses a problem for Hi-GAL data: large
signal jumps are present in the field and the resulting gap filling
realizations are affected in a non negligible manner by the boundary
conditions, at least with the present version of ROMAGAL. This behavior is
different from what happens for experiments aimed at the Cosmic Microwave
Background (see e.g. Masi et al. (2006)) where NCR codes are routinely run on
the timelines as they are. In order to find a workaround that would spare us
the inconvenience of estimating the underlying noise field to serve as a NCR
input, we have modified the flag treatment of the ROMAGAL code as explained in
the following.

The original version of ROMAGAL makes use of a single extra pixel (dubbed
“virtual pixel”) to serve as junk bin where the contents of the gaps are sent
when applying the $P^{T}$ operator within the ROMAGAL solver. This approach,
as stated above, works excellently in presence of both signal and noise,
irrespective of their relative amplitude, provided the NCR code assumes the
underlying noise field as a baseline to perform the realization. In order to
relax this assumption, we have modified ROMAGAL to take into account not a
single virtual pixel but an entire virtual map. In other words, we introduce a
virtual companion for each pixel of the map, and use it as a junk bin to
collect the output from the gaps they correspond to. The hope is to
redistribute the content of the flagged sections more evenly, preventing
artifacts. This approach obtains satisfactory results when the NCR code is run
on the original (signal plus noise) timelines, as shown in Figure 12.

Figure 12: Top row: For a signal dominated case, shown are the relative
differences between the input simulated map and the output obtained with
ROMAGAL (in the “virtual map” mode) with NCR performed on the underlying noise
(left), on the original signal plus noise TOD (middle), and without NCR after
replacing the gaps with null values. The latter case is clearly striped, but
no signal related artifacts are present. Bottom row: For a noise dominated
case, shown again are the relative difference versus the input map obtained by
ROMAGAL with “virtual map”, assuming NCR on underlying noise only (left),
without NCR (middle) and with NCR on the signal plus noise TOD (right). As
expected, the left panel shows the best residuals, but the right one appears
as a good compromise (see also text).

To summarize our findings:

  * •

Using a NCR realization code is always necessary in order to avoid artifact
striping effect into the map

  * •

If the NCR code is run on the underlying noise-only timelines (which are
cumbersome to estimate) we obtain the best quality output map, with no signal
related artifacts and a noise standard deviation which is lower (by a factor
$\sim 2$) with respect to the case in which no NCR is performed.

  * •

Running the NCR on the original timelines is possible with the “virtual map”
approach. No signal artifacts are detected in the difference maps and the
advantage in terms of noise standard deviation with respect to the no NCR case
is still present, and of order 10% to 20% on average.

We have therefore chosen this latter approach as the baseline for the Hi-GAL
pipeline.

##  5 ROMAGAL maps

In this section we analyze the final map obtained running our dedicated
pipeline. We analyze the Point Spread Function (PSF) for the five bands in
order to fix the resolution of the ROMAGAL maps. We compare the final GLS map
with the naive map and we point out the differences and the capability to
recover the diffuse emission from the GLS map. Finally we discuss about the
noise residuals on the maps.

###  5.1 Point Spread Function and pixel size

The angular resolution (pixel size) of the final map is a free parameter. Its
choice is a compromise between competing requirements. A small pixel size
assures a correct sampling of the PSF;

indeed, assuming a Gaussian profile for the PSF (which is reasonable as
discussed in the following), the Nyquist theorem imposes that to better sample
a 2-d image we need to set a pixel size which is at most one third of its FWHM
value.

On the other hand, a too small pixel size can cause the loss of redundancy,
useful to reduce the white noise level, and (even) some non-observed pixels in
the final map.

The diffraction limited beam of PACS at $70\mu$m is $5.2\arcsec$. Thus, we
should build the map with a pixel size of at least $1.8\arcsec$. However, due
to the limited bandwidth available for transmission, especially in PACS/SPIRE
Parallel mode, PACS frames are coadded on-board the satellite before
broadcasting. For the $70$$\mu$m channel, a coaddition of 8 consecutive frames
is applied by on-board software. Since the acquisition rate is 40Hz and the
scanning speed for Hi-GAL is set to $60\arcsec/s$, two close frames are
snapshots of the sky acquired $1.5\arcsec$ apart. Due to coaddition, the
satellite provides one averaged frame every $12\arcsec$; in spatial
coordinates, this is twice the beam width of the PACS blue channel. The
measured PSF then is not the diffraction limited beam but it results in an
elongation along the scan direction due to the convolution of the coaddition
with the beam. As shown in Lutz (2009), the observations of Vesta and
$\alpha$Tau with the blue channel evidenced a FWHM equal to $5.86\arcsec\times
12.16\arcsec$ as a result of a 2-d Gaussian fitting, elongated in the scan
direction.

PACS $160\mu$m is also affected by the averaging on-board, but only 4 frames
are coadded together. The nominal instrumental beam is $12.0\arcsec$, while
the measured is $11.64\arcsec\times 15.65\arcsec$ (Lutz, 2009), elongated
along the scan direction. However, in this case we can sample the beam without
issues, and the effect of coaddition on the final map is negligible.

For the Hi-GAL fields the scanning strategy consists of two orthogonal AORs,
therefore the redundancy regularizes the PSF, resulting in approximately 2-d
symmetric Gaussian profile, as shown in Table 2.

According to the values reported in Table 2, we observe quasi-symmetric beams
with an averaged ellipticity of less than $15\%$ for both blue and red channel
and the axis oriented randomly with respect to the scan direction.

We choose a pixel size of $3.2\arcsec$ for the PACS $70\mu$m band which
samples the observed beam at Nyquist frequency. Below this threshold, the
diffuse emission areas become too noisy due to the low SNR. Similarly, we can
choose a pixel size of $4.5\arcsec$ for red band without loosing redundancy.

SPIRE does not suffer from on-board coadding and the detectors were built to
reach the band diffraction limit. In-flight data show that the SPIRE beam is
well approximated by a 2-d asymmetric Gaussian curve with the major axis
orientation independent of the scan direction, with an ellipcticity not bigger
then $10\%$ (see Sibthorpe et al. (2010)). We set the pixel size for each
SPIRE band equal to one third of the nominal beam. In Table 2 we also report
the beam measured in the SPIRE maps.

The average ellipticity we observe agrees with found by Sibthorpe et al.
(2010). On the contrary, while the FWHM for the two axis found by Sibthorpe et
al. (2010) are in agreement with the nominal (within the error), our measured
beam results in a FWHM larger then the nominal of $\sim 25\%$.

Band | Nominal Beam (arcsec) | Measured Beam (arcsec) | Ellipticity | Pixel Size (arcsec)
---|---|---|---|---
$70\mu$m | $5.2\times 5.2$ | $\sim 9.7\times\sim 10.7$ | 14.6% | 3.2
$160\mu$m | $12.0\times 12.0$ | $\sim 13.2\times\sim 13.9$ | 14.7% | 4.5
$250\mu$m | $18\times 18$ | $\sim 22.8\times\sim 23.9$ | 8.3% | 6.0
$350\mu$m | $24\times 24$ | $\sim 29.3\times 31.3$ | 8.8% | 8.0
$500\mu$m | $34.5\times 34.5$ | $\sim 41.1\times\sim 43.8$ | 9.7% | 11.5
Table 2: Nominal (2nd column), map-measured beam (two AOR, 3rd column) and
ellipticity (4th column) of each band.

###  5.2 Hi-GAL SDP results

The quality of the outcome can be qualified by the comparison between the igls
maps with the naive maps. In fact, the naive map is the simple averaging of
signal recorded by every spatial pixel and it represents “the least
compromised view of the sky”.

Since the TOD are created at the end of the preprocessing steps, when the data
are a combination of only signal and $1/f$ noise, we expect the $1/f$
residuals in the naive map as well as a “pure” sky map produced by ROMAGAL.

In Figure 13 a comparison between the naive map and the ROMAGAL map of
$\textit{l}=30^{\circ}$ field at $70\mu$m is shown. The GLS code is capable to
remove the $1/f$ residuals without loosing any signal both on bright sources
and on diffuse emission.

Figure 13: Left: particular of the ROMAGAL map of the red PACS array,
$\textit{l}=30^{\circ}$ field. Right: particular of the naive map, same band
and field. The $1/f$ noise is evident on the naive map, as well as its
minimization is evident on the GLS map without loosing the signal of the
diffuse component.

In particular, we choose three main proxies:

  * •

the difference between naive and igls should show only a pattern due to the
$1/f$ noise residuals in the binned map. The pattern of this low-frequency
noise is recognizable as stripes superimposed on the sky signal in the naive
map. The stripes are the consequence of the $1/f$ noise due to the adopted
scanning strategy.
In Figure 14 we show the difference between igls map and naive map. The $1/f$
noise is removed in the igls map but not in the naive. The residual stripes
due to the low-frequency noise are clearly visible.

Figure 14: Particular of map difference between PACS $160\mu$m igls and naive,
same region of the previous Figure. The stripes due to $1/f$ removal made in
the igls are evident.

  * •

the source fluxes should be the same between the igls and naive maps.
This item is quantified by the map difference, where the pattern is due only
to the noise without any residual signal, except across the very brilliant
source. This last effect is discussed on the next Section.

  * •

a statistical analysis of the background noise level should show a decrease of
the rms value in the igls with respect to the naive map.
In Tables 3 and 4 we report the rms residuals of the PACS and SPIRE maps
respectively, calculated in a diffuse emission area of each map. Since the
flux between the maps is conserved, a decreasing of the rms noise level
assures an increasing of S/N ratio in the ROMAGAL maps.

PACS $\textit{l}=30^{\circ}$ field
---
Band | rms igls (MJy/pixel) | rms naive(MJy/pixel) | ratio
$70\mu$m | 0.0085 | 0.026 |  $\sim$ 3.1
$160\mu$m | 0.047 | 0.102 |  $\sim$ 2.2
PACS $\textit{l}=59^{\circ}$ field
Band | rms igls (MJy/pixel) | rms naive(MJy/pixel) | ratio
$70\mu$m | 0.004545 | 0.02208 |  $\sim$ 4.9
$160\mu$m | 0.01899 | 0.03586 |  $\sim$ 1.9
Table 3: rms compared from GLS and naive map on both the SDP observations, for
PACS bands, measured on a background region of $50\times$50 pixels. In the
last column is reported the ratio between the naive and GLS rms.

The ratio between the naive and igls rms shows an improvement of a factor
$\sim 2-5$ in the PACS ROMAGAL maps, and a factor $\sim 1-2$ in the SPIRE
case. The difference is mostly due to an intrinsically different $1/f$ noise
level.

SPIRE $\textit{l}=30^{\circ}$ field
---
Band | rms igls (MJy/beam) | rms naive(MJy/beam) | ratio
$250\mu$m | 0.1749 | 0.2868 |  $\sim$ 1.6
$350\mu$m | 0.1569 | 0.2302 |  $\sim$ 1.5
$500\mu$m | 0.2659 | 0.4065 |  $\sim$ 1.5
SPIRE $\textit{l}=59^{\circ}$ field
Band | rms igls (MJy/beam) | rms naive(MJy/beam) | ratio
$250\mu$m | 0.09857 | 0.1123 |  $\sim$ 1.1
$350\mu$m | 0.0734 | 0.08164 |  $\sim$ 1.1
$500\mu$m | 0.1073 | 0.2101 |  $\sim$ 1.9
Table 4: rms compared from GLS and naive map on both the SDP observations, for
SPIRE bands, measured on a background region of $50\times$50 pixels. In the
last column is reported the ratio between the naive and GLS rms. The $1/f$
noise is less evident in the SPIRE bolometers with respect to the PACS one,
but its effect is still remarkable.

###  5.3 Consistency of data analysis assumptions

One of the assumptions of ROMAGAL, as well as of all Fourier based GLS map
making codes, is that the underlying sky signal is constant within a chosen
resolution element. If this is not the case, artifacts (stripes) will be
generated in the final map, contributing to the so called pixelization noise
(Poutanen et al., 2006). In the case of Hi-GAL the situation is complicated by
several effects:

  * •

On-board coaddition of samples: each PACS 70$\mu$m (160 $\mu$m) frame is the
result of an on-board average of eight (four) consecutive frames reducing the
effective sampling frequency of the instrument (see Section 5.1). Thus, sky
signal is low-pass filtered by only partially effective data-windowing, rather
then telescope response, leaving room for signal aliasing.

The map making code is quite sensitive to aliasing since it works in Fourier
space. The situation is worsened by the large dynamic range of the Hi-GAL
fields, especially when scanning across bright sources.

  * •

Time bolometer response induces signal distortions along the scan. While
within HIPE the SPIRE detector response is deconvoluted from the data (Griffin
et al., 2008), the same is not true for PACS. Redundancy in each pixel is
obtained by scans coming from different directions, thus the effect
contributes further signal mismatch at the pixel level.

  * •

Pointing error: as analyzed in detail in Garcia Lario et al. (2007) the
pointing performance of Herschel, which mean the capability of assign
coordinates to each sample in a given reference frame, can be affected by
several pointing error effects; the main contributor is due to the angular
separation between the desired direction and the actual instantaneous
direction, based on the position-dependent bias within the star-trackers.

The Herschel AOCS goal is that the mismatch between real coordinates and
assigned coordinates along the scan-leg is smaller than $1.2\arcsec$ at 1
sigma (Garcia Lario et al., 2007). So that a 2$\sigma$ event becomes
significant compared to the PSF of the PACS blue band.

Figure 15: Zoom around a compact sources for the PACS blue optimal map. The
striping dragged along the scan directions are remarkable.

All of the above effects challenge the basic assumptions (no signal aliasing
and no sub-pixel effect) under which ROMAGAL works. Our simulations suggest
that signal aliasing contribute significantly more than the other two. The net
result on maps is the striping of bright sources in the GLS maps. An example
is shown in Figure 15 for the PACS blue band.

It is important to notice that the stripes are present only around the point-
like sources (where - of course - the signal aliasing is more evident),
regardless of their magnitude. However, magnitude influences the amplitude of
the stripes. For a source peak which is only 10 times higher that the rms
background, the intensity of the stripes are within $1\sigma$ from the
background dispersion. For the more intense sources, the stripes magnitude in
the pixels surrounding the PSF shape can be $100\sigma$ times away from the
background value.

Since these stripes are produced within the GLS solver, which performs
repetitive convolutions along the scan directions, but do not affect the naive
map the obvious workaround is to implement dedicated map making which
considers a different analysis around the bright sources.

However, the detailed accounting of the above effects and the enhanced map
making strategies to address them will be the subject of a forthcoming paper
Piazzo et al. (2011).

##  6 Summary and conclusions

This paper describes in detail all steps of the processing of Herschel data,
from originally downloaded frames to high-quality maps, used in the framework
of Hi-GAL survey. Hi-GAL data are taken in fast scan mode ($60\arcsec$/sec )
and simultaneously by PACS and SPIRE (Parallel mode). We test our pipeline
reducing data from the Science Demonstration Phase and present results, taking
as proxy for the quality of the final images their comparison with naive maps.

We divided data processing into two distinct phases: preprocessing and
mapmaking. Pre-processing aims to accurately remove systematics and random
effects in the data, in order to prepare them for the ROMAGAL map making
algorithm, which implements the minimum variance GLS approach in order to
minimize the noise into the data. It turns out that NCR is a fundamental step
in the pre-processing because ROMAGAL, as an FFT mapmaking code needs
continuous and regular data time series as input.

Noise residuals in the diffuse emission of the two test fields (SDP Hi-GAL
data, two 2${}^{\circ}\times 2^{\circ}$ tiles centered on the Galactic plane
at l = 30∘ and l = 59∘) show that we obtain optimal maps, getting rid of
glitches and systematic drifts, as well as minimizing the 1/ f and white noise
component. The remaining effects, which do not affect the overall quality of
the maps except across the bright source on the PACS 70$\mu$m maps, are under
investigation that will appear in a dedicated publication to be available
shortly.

## Acknowledgements

Special thanks to Göran Pilbratt for allowing the using of PV Phase blank
data.

## References

  * Billot et al. (2010) Billot N., Sauvage M., Rodriguez L., Horeau B., Kiss C., Aussel H., Okumura K., Boulade O., Altieri B., Poglitsch A., Agnèse P., 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7741 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, CEA bolometer arrays: the first year in space
  * de Gasperis et al. (2005) de Gasperis G., Balbi A., Cabella P., Natoli P., Vittorio N., 2005, A&A, 436, 1159
  * Garcia Lario et al. (2007) Garcia Lario P., Heras A. M., Sanchez-Portal M., 2007, Herschel Observer Manual
  * Griffin et al. (2008) Griffin M., Dowell C. D., Lim T., Bendo G., Bock J., Cara C., Castro-Rodriguez N., Chanial P., Clements D., Gastaud R., Guest S., Glenn J., Hristov V., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7010 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, The Herschel-SPIRE photometer data processing pipeline
  * Griffin et al. (2008) Griffin M., Swinyard B., Vigroux L., Abergel A., Ade P., André P., Baluteau J., Bock J., Franceschini A., Gear W., Glenn J., Huang M., Griffin D., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7010 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Herschel-SPIRE: design, ground test results, and predicted performance
  * Griffin et al. (2010) Griffin M. J., Abergel A., Abreu A., Ade P. A. R., André P., Augueres J., Babbedge T., Bae Y., Baillie T., Baluteau J., Barlow M. J., Bendo G., Benielli D., et al. 2010, A&A, 518, L3+
  * Hestenes & Stiefel (1952) Hestenes M. R., Stiefel E., 1952, Journal of Research of the National Bureau of Standards, 49, 409
  * Hoffman & Ribak (1991) Hoffman Y., Ribak E., 1991, ApJL, 380, L5
  * Lupton (1993) Lupton R., 1993, Statistics in theory and practice. Princeton University Press, Princeton
  * Lutz (2009) Lutz D., 2009, PACS photometer point spread function
  * Martin et al. (2010) Martin P. G., Miville-Deschênes M., Roy A., Bernard J., Molinari S., Billot N., Brunt C., Calzoletti L., et al. 2010, A&A, 518, L105+
  * Masi et al. (2006) Masi S., Ade P. A. R., Bock J. J., Bond J. R., Borrill J., Boscaleri A., Cabella P., Contaldi C. R., Crill B. P., de Bernardis P., de Gasperis G., de Oliveira-Costa A., de Troia G., et al. 2006, A&A, 458, 687
  * Miville-Deschênes & Lagache (2005) Miville-Deschênes M., Lagache G., 2005, ApJS, 157, 302
  * Molinari et al. (2010b) Molinari S., Swinyard B., Bally J., Barlow M., Bernard J., Martin P., Moore T., Noriega-Crespo A., Plume R., Testi L., Zavagno A., et al. 2010b, A&A, 518, L100+
  * Molinari et al. (2010a) Molinari S., Swinyard B., Bally J., Barlow M., Bernard J., Martin P., Moore T., Noriega-Crespo A., Plume R., Testi L., Zavagno A., et al. 2010a, PASP, 122, 314
  * Natoli et al. (2001) Natoli P., de Gasperis G., Gheller C., Vittorio N., 2001, A&A, 372, 346
  * Ott et al. (2010) Ott S., Science Centre H., Space Agency E., 2010, ArXiv e-prints
  * Peretto et al. (2010) Peretto N., Fuller G. A., Plume R., Anderson L. D., Bally J., Battersby C., Beltran M. T., Bernard J., Calzoletti L., Digiorgio A. M., et al. 2010, A&A, 518, L98+
  * Pezzuto (2010) Pezzuto S., 2010, ”A new approach to the detection of outliers with the sigma-clipping method applied to PACS photometry observation”. Poster presented at ”The impact of Herschel surveys on ALMA Early Science” Garching, 16-19 November 2010
  * Piazzo et al. (2011) Piazzo L., Ikhenaode D., Natoli P., Pestalozzi P., Piacentini F., Traficante A., 2011, ”Artifacts removal for GLS map makers”, in preparation
  * Pilbratt et al. (2010) Pilbratt G. L., Riedinger J. R., Passvogel T., Crone G., Doyle D., Gageur U., Heras A. M., Jewell C., Metcalfe L., Ott S., Schmidt M., 2010, A&A, 518, L1+
  * Poglitsch et al. (2008) Poglitsch A., Waelkens C., Bauer O. H., Cepa J., Feuchtgruber H., Henning T., van Hoof C., Kerschbaum F., Krause O., Renotte E., Rodriguez L., Saraceno P., Vandenbussche B., 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7010 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, The Photodetector Array Camera and Spectrometer (PACS) for the Herschel Space Observatory
  * Poglitsch et al. (2010) Poglitsch A., Waelkens C., Geis N., Feuchtgruber H., Vandenbussche B., Rodriguez L., Krause O., Renotte E., van Hoof C., Saraceno P., Cepa J., Kerschbaum F., et al. 2010, A&A, 518, L2+
  * Poutanen et al. (2006) Poutanen T., de Gasperis G., Hivon E., Kurki-Suonio H., Balbi A., Borrill J., Cantalupo C., Doré O., Keihänen E., et al. 2006, A&A, 449, 1311
  * Schulz et al. (2008) Schulz B., Bock J. J., Lu N., Nguyen H. T., Xu C. K., Zhang L., Dowell C. D., Griffin M. J., Laurent G. T., Lim T. L., Swinyard B. M., 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7020 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Noise performance of the Herschel-SPIRE bolometers during instrument ground tests
  * Sibthorpe et al. (2010) Sibthorpe B., Ferlet M., Schulz B., 2010, SPIRE Beam Model Release Note. ICC, London
  * Starck et al. (1995) Starck J., Murtagh F., Louys M., 1995, in R. A. Shaw, H. E. Payne, & J. J. E. Hayes ed., Astronomical Data Analysis Software and Systems IV Vol. 77 of Astronomical Society of the Pacific Conference Series, Astronomical Image Compression Using the Pyramidal Median Transform. pp 268–+
  * Stompor et al. (2002) Stompor R., Balbi A., Borrill J. D., Ferreira P. G., Hanany S., Jaffe A. H., Lee A. T., Oh S., Rabii B., Richards P. L., Smoot G. F., Winant C. D., Wu J., 2002, PRD, 65, 022003
  * Wright (1996) Wright E. L., 1996, paper presented at the IAS CMB Data Analysis Workshop (astro-ph/9612006)