# Extragalactic magnetism with SOFIA (SALSA Legacy Program) - III:
First data release and on-the-fly polarization mapping characterization111The
SOFIA Legacy Program for Magnetic Fields in Galaxies (SALSA) provides a
software repository at https://github.com/galmagfields/hawc, and publicly
available data at http://galmagfields.com/

Enrique Lopez-Rodriguez Kavli Institute for Particle Astrophysics & Cosmology
(KIPAC), Stanford University, Stanford, CA 94305, USA  Melanie Clarke  SOFIA
Science Center, NASA Ames Research Center, Moffett Field, CA 94035, USA
Sachin Shenoy  William Vacca  SOFIA Science Center, NASA Ames Research Center,
Moffett Field, CA 94035, USA  Simon Coude  SOFIA Science Center, NASA Ames
Research Center, Moffett Field, CA 94035, USA  Ryan Arneson  SOFIA Science
Center, NASA Ames Research Center, Moffett Field, CA 94035, USA  Peter Ashton
SOFIA Science Center, NASA Ames Research Center, Moffett Field, CA 94035, USA
Sarah Eftekharzadeh  SOFIA Science Center, NASA Ames Research Center, Moffett
Field, CA 94035, USA  Rainer Beck  Max-Planck-Institut für Radioastronomie,
Auf dem Hügel 69, 53121 Bonn, Germany  John E. Beckman  Instituto de
Astrofísica de Canarias, C/ Via Láctea s/n, E-38200 La Laguna, Tenerife, Spain
Departamento de Astrofísica, Universidad de La Laguna, Avda. Astrofísico Fco.
Sánchez s/n, E-38200 La Laguna, Tenerife, Spain  Alejandro S. Borlaff NASA
Ames Research Center, Moffett Field, CA 94035, USA  Susan E. Clark Department
of Physics, Stanford University, Stanford, CA 94305, USA  Kavli Institute for
Particle Astrophysics & Cosmology (KIPAC), Stanford University, Stanford, CA
94305, USA  Daniel A. Dale Department of Physics & Astronomy, University of
Wyoming, Laramie, WY 82071, USA  Sergio Martin-Alvarez Institute of Astronomy
and Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK  Evangelia Ntormousi Scuola Normale Superiore di Pisa,
Piazza dei Cavalieri 7, 56126 Pisa, Italy  Institute of Astrophysics,
Foundation for Research and Technology-Hellas, Vasilika Vouton, GR-70013
Heraklion, Greece  William T. Reach Universities Space Research Association,
NASA Ames Research Center, Moffett Field, CA 94035, USA  Julia Roman-Duval
Space Telescope Science Institute 3700 San Martin Drive Baltimore, MD21218,
USA  Konstantinos Tassis Department of Physics & ITCP, University of Crete,
GR-70013, Heraklion, Greece  Institute of Astrophysics, Foundation for
Research and Technology-Hellas, Vasilika Vouton, GR-70013 Heraklion, Greece
Doyal A. Harper  Department of Astronomy and Astrophysics University of
Chicago, Chicago, IL 60637, USA  Pamela M. Marcum  NASA Ames Research Center,
Moffett Field, CA 94035, USA

###### Abstract

We describe the data processing of the Survey on extragALactic magnetiSm with
SOFIA (SALSA Legacy Program). This first data release presents 33% (51.34h out
of 155.7h, including overheads) of the total awarded time taken from January
2020 to December 2021. Our observations were performed using the newly
implemented on-the-fly mapping (OTFMAP) technique in the polarimetric mode. We
present the pipeline steps to obtain homogeneously reduced high-level data
products of polarimetric maps of galaxies for use in scientific analysis. Our
approach has a general design and can be applied to sources smaller than the
field-of-view of the HAWC+ array in any given band. We estimate that the
OTFMAP polarimetric mode offers a reduction of observing overheads by a factor
2.34, and an improvement in sensitivity by a factor 1.80 when compared to
previously obtained polarimetric observations using the chopping and nodding
mode. The OTFMAP is a significant optimization of the polarimetric mode of
HAWC+ as it ultimately reduces the cost of operations of SOFIA/HAWC+ by
increasing the science collected per hour of observation up to an overall
factor of 2.49. The OTFMAP polarimetric mode is the standard observing
strategy of SALSA. The results and quantitative analysis of this first data
release are presented in Papers IV and V of the series.

††facilities: SOFIA (HAWC+), Herschel††software:  astropy (Astropy
Collaboration et al., 2013), APLpy (Robitaille & Bressert, 2012), matplotlib
(Hunter, 2007)

##  1 Introduction

Infrared (IR) observations from ground-based and suborbital telescopes are
dominated by atmospheric thermal emission and telescope background emission.
Revealing the astrophysical sources requires observing strategies that remove
or minimize the competing contributions of these background emissions. The
general approach is to observe the combined emission from the source and the
background with a detector, and then move the same detector to a part of the
sky with only background emission. The offset between these two positions is
performed on a timescale shorter than the background variation and the
detector response. The subtraction of these two observations reveals the
astrophysical source.

The chopping and nodding (C2N) strategy has been the most commonly used
observing technique in the mid-IR (MIR; $8-30$ $\mu$m) wavelength range (e.g.
Burtscher et al., 2020, and references therein) and sometimes used in the far-
IR (FIR; $30-300$ $\mu$m) wavelength range (e.g. Hildebrand et al., 2000). The
secondary mirror alternates between an on-axis and an adjacent off-axis
position on a fast cadence. The angular separations (i.e. chop-throw) between
on- and off-axis positions are usually within two to three fields-of-view
(FOV) of the array in a given band, yielding $\sim 5-30$″ for MIR in 8-m class
telescopes and up to $500$″ for FIR wavelengths. The cadence of $<1$s
facilitates subtraction of fast sky variations. This telescope movement is
known as chopping. Each pair of on-axis and off-axis images are subtracted to
eliminate the background emission from the science image. However, as the
secondary mirror is not aligned with the optical axis of the telescope while
observing the nearby sky location, a residual background emission (also named
radiative offset) is still present in the subtracted image. To minimize the
radiative offset, the telescope is moved on its optical axis every $\leq 60$s
with the same, or different, amplitude and direction as the chopping
configuration. A pair of subtracted chop-nod observations produces a positive
image of the science object.

The C2N observing strategy has several disadvantages. First, although the
residual offset is subtracted, final reduced images still have some level of
background that generates emission patterns and/or residual levels of
radiation that are not entirely removed. These residuals may arise from
radiation of structures in the dome, sky variations, detector response
variations, and/or dependencies of radiation offset as the telescope rotates.
Second, chop-throws are usually $\leq 10-500$″, which can complicate the
background subtraction for objects with extended diffuse thermal emission
(e.g. the Galactic center, star-forming regions, molecular clouds). Finally,
this observing mode spends approximately half of the total observing time
pointing to a nearby sky location away from the science object, which makes
the observations inefficient with large overheads.

The on-the-fly mapping (OTFMAP) strategy is an alternative observing mode in
which the telescope on its optical axis scans the science object and
surrounding areas. This technique is commonly used in the FIR and sub-mm
wavelength ranges (e.g. Tegmark, 1997; Hildebrand et al., 2000; Reichertz et
al., 2001; Weferling et al., 2002; Waskett et al., 2007; Chapin et al., 2013).
Single-dish radio telescopes also use this observing technique (Haslam, 1974;
Müller et al., 2017). This technique has recently gained attention in the MIR
wavelength range (e.g. Ohsawa et al., 2018) due to its prospects for the next
generation of 30-m class telescopes. The telescope on its optical axis
performs a scanning strategy (e.g. Kovács, 2008a) described by a parametric
curve and a mapping speed. The mapping speed is chosen such that the beam size
(full width at half maximum, FWHM) of the instrument is oversampled at the
sampling rate of the detector readout. For this configuration: a) the sampling
rate has to be at least twice the maximum frequency of the astrophysical
and/or sky signal (i.e. spatial and temporal Nyquist sampling), and b) the
system noise must be comparable to or lower than the power spectrum of the sky
at the desired spatial frequency to be recovered (e.g. Emerson & Graeve,
1988).

The OTFMAP observing strategy presents advantages in terms of observing
overheads as at any given time the scan is integrating over a deeper image.
Also, the telescope can be moved to farther adjacent sky regions to reach
‘true’ zero-level background emission. For bolometers, the signal is dominated
by the $1/f$ noise, which is characterized by having higher amplitudes at
longer timescales. This behavior results in signal drifts that can affect all
spatial structures in the final image. Most of the efforts of the OTFMAP
strategy focus on the map-making algorithms and signal filtering to recover
diffuse and extended emission from the astrophysical object. A map-making
algorithm is required to reconstruct the sky image while subtracting the sky
and detector variation responses (i.e. Roussel, 2013; Kovács, 2006; Waskett et
al., 2007; Patanchon et al., 2008; Kovács, 2008b; Cantalupo et al., 2010;
Chapin et al., 2013).

Chop-nod has been the standard observing strategy for the polarimetric mode of
the High-resolution Airborne Wideband Camera-plus (HAWC+; Vaillancourt et al.,
2007; Dowell et al., 2010; Harper et al., 2018) onboard the 2.7-m
Stratospheric Observatory For Infrared Astronomy (SOFIA). The OTFMAP observing
mode had previously only been used for total intensity observations (Harper et
al., 2018). Recent observations of galaxies, Centaurus A (Lopez-Rodriguez,
2021) and NGC 1097 (Lopez-Rodriguez et al., 2021), and the molecular cloud
Taurus/L1495 (Li et al., 2021) have successfully used the OTFMAP polarimetric
mode with HAWC+. However, a full characterization of this new observing mode
for HAWC+ has not been presented until now.

This manuscript describes the data reduction scheme of the OTFMAP polarimetric
mode of HAWC+ for sources smaller than the FOV in any given band. In Section
2, we present the individual steps from raw data to high-level data products
for use in scientific analysis. We present approaches for background
subtraction, precipitable water vapour correction, and flux and polarization
calibration. In Section 3, we show a quantitative comparison between C2N and
OTFMAP polarization modes and estimate the sensitivities for the OTFMAP
observing mode. First results of this data release will be presented in Paper
IV and V of the series.

##  2 Observations and Data Reduction

###  2.1 Observations

New multi-wavelength polarimetric observations of 10 galaxies were performed
as part of the SOFIA Legacy Program222SOFIA Legacy Program on extragalactic
magnetism can be found at http://galmagfields.com/ (ID: 08_0012, PI: Lopez-
Rodriguez, E. & Mao, S. A.) using HAWC+ on the 2.7-m SOFIA telescope. This
first data release presents 33% (51.34h out of 155.7h, including overheads) of
the total awarded time taken from January 2020 to December 2021. Table 1 shows
the details of the new observations presented in this manuscript. In addition,
we also include observations of Centaurus A (Lopez-Rodriguez, 2021) at 89
$\mu$m and Circinus at 53, 89, and 214 $\mu$m from program ID 07_0032 (PI:
Lopez-Rodriguez, E.), M51 (Borlaff et al., 2021) at 154 $\mu$m under programs
70_0509 (Guaranteed Time Observations by the HAWC+ Team), 76_0003
(Discretionary Director Time), and 08_0260 (PI: Dowell, D.), and NGC 1097
(Lopez-Rodriguez et al., 2021) at 89 and 154 $\mu$m from program ID 07_0034
(PI: Lopez-Rodriguez, E.). These additional observations comprised a total
exposure time of 11.84h (including overheads) taken from 2017 to 2021. A total
of 14 galaxies are presented in this data release.

All galaxies presented in this data release were homogeneously reduced using
the steps presented below. For those galaxies already published, i.e.
Centaurus A and NGC 1097, we reprocessed these datasets to produce data
products within a homogeneous reduction scheme.

Table 1: Summary of new OTFMAP polarimetric observations. _Columns, from left
to right:_ a) Object name. b) Central wavelength of the band used for
observations. c) Observation date (YYYYMMDD). d) Flight ID. e) Sea-level
altitude during the observations (kft.). f) Speed of the scan (″/sec.). g)
Amplitudes in elevation (EL) and cross-elevation (XEL) of the scan (″). h)
Time per scan (s). i) Number of observation sets used (and rejected). j) On-
source observation time (s).

Object | Band | Date | Flight | Altitude | Scan Rate | Scan Amplitude | Scan Duration | #Sets (bad) | Obs. Time
---|---|---|---|---|---|---|---|---|---
|  | (YYYYMMDD) |  | (kft) | (″/sec.) | (EL $\times$ XEL; ″) | (s) |  | (s)
(a) | (b) | (c) | (d) | (e) | (f) | (g) | (h) | (i) | (j)
Circinus⋆ | A (53 $\mu$m) | 20190716 | F596 | 42 | 100 | 50$\times$50 | 100 | 1(2) | 400
|  |  |  |  |  |  |  | Total time | 400 (0.11h)
| C (89 $\mu$m) | 20190716 | F596 | 41-42 | 100 | 70$\times$70 | 100 | 2 | 800
|  |  |  |  |  | 90$\times$90 | 100 | 6 | 2400
|  |  |  |  |  |  |  | Total time | 3200 (0.89h)
| E (214 $\mu$m) | 20190716 | F596 | 43 | 100 | 120$\times$120 | 100, 60 | 2, 1 | 1040
|  |  |  |  |  |  |  | Total time | 1040 (0.29h)
M 82 | A (53 $\mu$m) | 20210505 | F726 | 41-42 | 100 | 80$\times$80 | 78, 60 | 1(2), 8(1) | 2232
|  | 20210513 | F731 | 43 | 100 | 80$\times$80 | 78 | 11 | 3432
|  |  |  |  |  |  |  | Total time | 5664 (1.6h)
| C (89 $\mu$m) | 20210518 | F733 | 43 | 100 | 90$\times$90 | 78 | 11 | 3432
|  | 20210521 | F736 | 42-43 | 100 | 90$\times$90 | 78 | 8(3) | 2496
|  |  |  |  |  |  |  | Total time | 5928 (1.65h)
| D (153 $\mu$m) | 20210520 | F735 | 41-43 | 200 | 100$\times$100 | 78 | 18 | 5616
|  | 20210521 | F736 | 42 | 200 | 100$\times$100 | 78 | 5 | 1560
|  |  |  |  |  |  |  | Total time | 7176 (1.99h)
| E (214 $\mu$m) | 20210506 | F727 | 39-40 | 200 | 180$\times$180 | 78 | 3(3) | 936
|  |  |  |  |  | 210$\times$210 | 78 | 7(1) | 2184
|  | 20210507 | F728 | 39 | 200 | 210$\times$210 | 78 | 0(1) | 0
|  | 20210511 | F729 | 39 | 200 | 210$\times$210 | 78 | 2 | 624
|  | 20210513 | F731 | 43 | 200 | 180$\times$180 | 78 | 1(2) | 312
|  |  |  |  |  |  |  | Total time | 4056 (1.13h)
M 83 | D (154 $\mu$m) | 20210512 | F730 | 38-41 | 200 | 100$\times$100 | 120 | 15 (1) | 6720
|  | 20210513 | F731 | 41 | 200 | 100$\times$100 | 120 | 8 | 3840
|  | 20210520 | F735 | 39-41 | 200 | 100$\times$100 | 120, 60 | 16, 1 | 7920
|  | 20210521 | F736 | 39-42 | 200 | 100$\times$100 | 120 | 10 | 4800
|  |  |  |  |  |  |  | Total time | 23760 (6.60h)
NGC 253 | C (89 $\mu$m) | 20210901 | F776 | 43 | 100 | 70$\times$70 | 100 | 9 | 3600
|  |  |  |  |  |  |  | Total time | 3600 (1.00h)
| D (154 $\mu$m) | 20210915 | F784 | 43 | 100 | 90$\times$90 | 100 | 7(1) | 3200
|  |  |  |  |  |  |  | Total time | 3200 (0.89h)
NGC 1068 | A (53 $\mu$m) | 20210828 | F774 | 43 | 100 | 60$\times$60 | 120 | 3 | 1440
|  | 20211103 | F786 | 43 | 100 | 60$\times$60 | 120, 45 | 6(1), 1 | 3060
|  |  |  |  |  |  |  | Total time | 4500 (1.06h)
| C (89 $\mu$m) | 20211207 | F800 | 40-42 | 100 | 70$\times$70 | 120, 60 | 11, 1 | 5520
|  | 20211209 | F802 | 42-43 | 100 | 70$\times$70 | 120, 60 | 3(1), 1 | 1680
|  |  |  |  |  |  |  | Total time | 7200 (2.00h)
NGC 1097⋆ | D (154 $\mu$m) | 20200128 | F654 | 38 | 200 | 90$\times$90 | 120 | 2 | 960
|  |  |  |  |  |  |  | Total time | 960 (0.27h)
NGC 2146 | A (53 $\mu$m) | 20200125 | F653 | 41-42 | 100, 60 | 30$\times$30, 60$\times$60 | 100 | 2, 1 | 1200
|  |  |  |  | 100 | 60$\times$60 | 60 | 2 | 480
|  | 20200131 | F657 | 40-42 | 100 | 30$\times$30 | 100 | 3(2) | 1200
|  |  |  |  |  |  | 60 | 7 | 1680
|  | 20210828 | F774 | 43 | 100 | 30$\times$30 | 100 | 0(12) | 0
|  | 20210831 | F775 | 41 | 100 | 30$\times$30 | 100 | 3 | 1200
|  |  |  |  |  |  |  | Total time | 5760 (1.60h)
| C (89 $\mu$m) | 20200130 | F656 | 41 | 100 | 50$\times$50 | 100 | 5(1) | 2000
|  | 20200131 | F657 | 42 | 100 | 50$\times$50 | 80 | 4(1) | 1280
|  | 20210914 | F783 | 40-43 | 100 | 50$\times$50 | 100 | 10(1) | 4000
|  |  |  |  |  |  |  | Total time | 7280 (2.02h)
| D (154 $\mu$m) | 20200131 | F657 | 42-43 | 100 | 70$\times$70 | 80 | 2 | 640
|  |  |  |  |  | 90$\times$90 | 80, 60 | 4(1), 1 | 1200
|  | 20210512 | F730 | 42-43 | 200 | 100$\times$100 | 100 | 12(3) | 3600
|  | 20210911 | F782 | 41-42 | 200 | 100$\times$100 | 100 | 4(3) | 1600
|  | 20210914 | F783 | 40 | 200 | 100$\times$100 | 100 | 4(1) | 1600
|  |  |  |  |  |  |  | Total time | 8640 (2.40h)
| E (214 $\mu$m) | 20210512 | F730 | 43 | 200 | 120$\times$120 | 100 | 4(2) | 1600
|  | 20210519 | F734 | 39 | 200 | 120$\times$120 | 100 | 2(3) | 800
|  | 20210831 | F775 | 41-43 | 200 | 120$\times$120 | 100 | 15 | 6000
|  |  |  |  |  |  |  | Total time | 8400 (2.33h)
NGC 3627 | D (154 $\mu$m) | 20210511 | F729 | 40-41 | 200 | 110$\times$110 | 120 | 12(1) | 5760
|  | 20210514 | F732 | 39-41 | 200 | 110$\times$110 | 120 | 11(3) | 5280
|  | 20210519 | F734 | 38 | 200 | 110$\times$110 | 120 | 11(1) | 5280
|  |  |  |  |  |  |  | Total time | 16320 (4.53h)
NGC 4736 | D (154 $\mu$m) | 20211207 | F800 | 43 | 200 | 110$\times$100 | 102, 90 | 2, 12 | 5136
|  | 20211208 | F8001 | 43 | 200 | 170$\times$170 | 102 | 7 | 2856
|  |  |  |  |  |  |  | Total time | 7992 (2.22h)
NGC 4826 | C (89 $\mu$m) | 20210520 | F735 | 43 | 200 | 110$\times$110 | 100 | 4 | 1600
|  | 20210521 | F736 | 43 | 200 | 110$\times$110 | 100 | 7 | 2800
|  |  |  |  |  |  |  | Total time | 4400 (1.22h)
NGC 6946 | D (154 $\mu$m) | 20200909 | F683 | 41-42 | 200 | 120$\times$120 | 100 | 1 | 400
|  |  |  | 42-43 | 200 | 180$\times$180 | 80 | 17 | 5440
|  |  |  | 43 | 200 | 220$\times$220 | 80 | 1 | 360
|  | 20200910 | F684 | 43 | 200 | 120$\times$120 | 80 | 1 | 320
|  |  |  |  |  | 200$\times$200 | 80, 60 | 2(4), 5(2) | 1840
|  | 20200911 | F685 | 43 | 200 | 120$\times$120 | 80 | 2 | 640
|  |  |  |  |  | 180$\times$180 | 100 | 14 | 5600
|  | 20200922 | F688 | 41-43 | 200 | 200$\times$200 | 60, 90 | 12, 3 | 3960
|  | 20200923 | F689 | 42-43 | 200 | 180$\times$180 | 80 | 10(1) | 3200
|  | 20200924 | F690 | 41 | 200 | 200$\times$200 | 80 | 13(1) | 4160
|  |  |  |  |  |  |  | Total time | 25920 (7.20h)
NGC 7331 | D (154 $\mu$m) | 20210511 | F729 | 43 | 200 | 120$\times$120 | 120 | 9 | 4320
|  | 20210513 | F731 | 43 | 200 | 120$\times$120 | 120 | 7 | 3360
|  | 20210514 | F732 | 42-43 | 200 | 120$\times$120 | 120, 90 | 12(1), 1 | 6120
|  | 20210518 | F733 | 43 | 200 | 120$\times$120 | 120 | 2 | 960
|  | 20210519 | F734 | 40-43 | 200 | 120$\times$120 | 120, 60 | 9, 1 | 4560
|  | 20210903 | F778 | 43 | 200 | 120$\times$120 | 120, 60 | 3, 1 | 1680
|  | 20210910 | F781 | 43 | 200 | 120$\times$120 | 120 | 2(4) | 960
|  | 20211103 | F786 | 42 | 200 | 120$\times$120 | 120 | (4) | 0
|  |  |  |  |  |  |  | Total time | 21960 (6.10h)

⋆⋆footnotetext: Circinus observations from program ID 07_0032 (PI: Lopez-
Rodriguez, E.) and NGC 1097 observations from program ID 07_0034 (PI: Lopez-
Rodriguez, E.).

###  2.2 Data reduction of the on-the-fly-map (OTFMAP) polarimetric mode

Figure 1: On-the-fly-map (OTFMAP) polarimetry data reduction flow-chart. The
raw data is combined in sets of four half-wave plate position angles (HWP PA)
within a sky rotation tolerance of $3^{\circ}$. The sets are reduced using
crush, which generates a scan per $\mathcal{R}$ and $\mathcal{T}$ arrays and
per HWP PA. Then, these scans are zero-level background corrected and
combined, from which Stokes parameters are derived. After the Stokes
parameters are corrected for instrumental polarization, the scans are sky
rotated and flux calibrated. Finally, all Stokes parameters from different
legs, flights, and/or cycles are merged to compute the final images, which
contain Stokes $IQU$, polarization fraction, $P$, position angle, $PA$, of
polarization, and polarized fluxes, $PI$. The pipeline computes uncertainties
and covariances, which are taken into account during all stages of the data
reduction.

We performed observations using the OTFMAP polarimetric mode with HAWC+. This
technique was an experimental observing mode used during the SOFIA Cycle 7
observations as part of engineering time (IDs: 71_0019 and 71_0023, PI: Lopez-
Rodriguez, E.) to optimize the polarimetric observations of HAWC+. This
observing mode has been routinely used during SOFIA Cycles 8 and 9 for further
characterization as part of the programmatic responsibilities of this SOFIA
Legacy Program and offered in ‘shared-risk’ mode for those proposals that
desired to test this new observing mode. This technique was successfully
applied to both galaxies (i.e. Centaurus A (Lopez-Rodriguez, 2021) and NGC
1097 (Lopez-Rodriguez et al., 2021)) and Galactic star-forming regions (the
filamentary cloud Taurus/L1495 (Li et al., 2021), and OMC-3/4; Li et al.,
submitted to ApJ). These early results were part of the SOFIA Cycle 7 and 8
observations. Here, we describe the details of the OTFMAP polarimetric mode as
the new observing strategy for HAWC+ on the study of magnetic fields in
galaxies.

####  2.2.1 Processing raw data using crush

We reduced the raw data using the Comprehensive Reduction Utility for SHARP II
v.2.50-1 (crush; Kovács, 2006, 2008b) and the hawc_drp_v2.7.0 pipeline
developed by the data reduction pipeline group at the SOFIA Science Center333A
full description of the most recent pipeline by the SOFIA Science center can
be found at https://www.sofia.usra.edu/data/data-analysis. SOFIA provided a
version, hawc_drp_v3.0.0, that replaces crush, written in java, with a new
version written in python. This new version444SOFIA data pipeline in python
can be found at: https://github.com/SOFIA-USRA/sofia_redux and
https://www.sofia.usra.edu/data/data-pipelines is a one-for-one drop-in
replacement with all crush’s options included. The analysis presented here
still holds for this new version.

HAWC+ polarimetric observations simultaneously measure two orthogonal
components of linear polarization arranged in two arrays of $32\times 40$
pixels each, labeled as reflective ($\mathcal{R}$) and transmissive
($\mathcal{T}$). Table 2 shows the characteristics of the HAWC+ bands, pixel
scales, beam sizes (FWHM), and FOVs of the array for the polarimetric
observations. For more information about HAWC+, we refer the reader to Harper
et al. (2018).

Table 2: HAWC+ configuration. _Columns, from left to right:_ a) Band name. b) Central wavelength of the band in $\mu$m. c) FWHM bandwidth in $\mu$m. d) Pixel scale of the band in ″. e) FWHM in ″. and f) FOV for polarimetric observations in ′.  Band | $\lambda_{c}$ | $\Delta\lambda$ | Pixel scale | Beam size | Polarimetry FOV
---|---|---|---|---|---
| ($\mu$m) | ($\mu$m) | (″) | (″) | (′)
(a) | (b) | (c) | (d) | (e) | (f)
A | $53$ | $8.7$ | $2.55$ | $4.85$ | $1.4\times 1.7$
C | $89$ | $17$ | $4.02$ | $7.80$ | $2.1\times 2.7$
D | $154$ | $34$ | $6.90$ | $13.6$ | $3.7\times 4.6$
E | $214$ | $44$ | $9.37$ | $18.2$ | $4.2\times 6.2$

The OTFMAP polarimetric mode performs observations in a sequence of four
Lissajous scans, where each scan has a different halfwave plate (HWP) position
angle (PA) in the following sequence: $5^{\circ}$, $50^{\circ}$,
$27.5^{\circ}$, and $72.5^{\circ}$. A starting angle of $5^{\circ}$ is chosen
to avoid the reference position, $0^{\circ}$, of the HWP (Harper et al.,
2018). This sequence is called a ‘set’ hereafter (Figure 1-first row). In this
HAWC+ observing mode, the telescope is driven to follow a parametric curve
with a nonrepeating period whose shape is characterized by the relative phases
and frequency of the motion. Each scan is characterized by its amplitude,
rate, angles, and duration. Table 1 shows the Lissajous parameters for each
galaxy (examples of Lissajous patterns are shown in Figures 10 and 15).

Each set was reduced using crush, which estimates and removes the correlated
atmospheric and instrumental signals, solves for the relative detector gains,
and determines the noise weighting of the time streams in an iterative
pipeline scheme. Each reduced scan produces two images associated with each
array, $\mathcal{R}$ and $\mathcal{T}$. Both images are orthogonal components
of linear polarization at a given HWP PA. At this stage, eight different scans
are computed (Figure 1-second row).

One of the limitations of the OTFMAP observing mode lies in the recovery of
large-scale diffuse and faint emission from the astrophysical objects. This
challenge is a result of the finite size of the array, variable atmosphere
conditions, variable detector temperatures, scan rates, and the applied
filters in the reduction steps used to recover extended emission. We applied
several filters using crush to recover large-scale emission structures of the
galaxies while paying close attention to any change that may compromise the
intrinsic polarization pattern of the astrophysical source. A total of 16
different pipeline schemes (e.g. faint, extended, bright, number of
iterations, several filtering options) were performed per galaxy per band. For
each of these pipeline schemes, the full data reduction steps were performed
and statistical comparisons between the Stokes parameters $IQU$ were computed.
Specifically, flux conservation was compared between pipeline schemes. The
morphologies between Stokes $IQU$ per pipeline scheme and between Stokes $I$
and PACS/Herschel observations at the closest wavelength were cross-
correlated, and Anderson-Darling tests were conducted using the polarization
fraction and angles between pipeline schemes. In addition, final products
using these different schemes were compared to the galaxies with available C2N
observations (M82, NGC 1068, and Circinus). We find that any additional
options to the pipeline configuration described above introduce artificial
polarization patterns to the galaxies. Some examples of the artificial
polarization patterns include: a) constant PA of polarization across the FOV,
b) changes of the PA of polarization larger than $30^{\circ}$ in regions of
high signal-to-noise ratio in Stokes $I$ (SNR${}_{I}\geq 500$), c) additional
extended structures in Stokes $I$ that were not present in Herschel images, d)
flux calibration incompatible with Herschel images, and e) differences in the
PA of polarization of up to $50^{\circ}$ when compared with the C2N observing
mode. Thus, the standard pipeline using the nominal configuration for crush
provided by the SOFIA Science Center was used.

The angular extension of all galaxies are within the polarimetry FOVs of the
array in any given band (Table 2). For all observations, the final FOV of the
images have zero-level background which enables a measure of the true sky
background without contribution from the astrophysical object. Thus, the
standard pipeline recovers the large-scale extended emission of the galaxies
in our sample. In contrast, observations of molecular clouds or the Galactic
center may have extended emission much larger than the polarimetry FOV in any
given band. We point the reader to Taurus/L1495 observed with HAWC+ at 214
$\mu$m, where further pipeline steps were required to recover those large-
scales (Li et al., 2021). For all the galaxies in our sample, observations
with a sky rotation smaller than $3^{\circ}$ were reduced within the same set,
which increases the signal-to-noise ratio (SNR) per pixel to optimize the
correlation between the astrophysical, sky, and instrumental signals. This
approach optimizes the recovery of large-scale and low-surface brightness of
the science object while keeping the rotation of the polarization angle vector
within the intrinsic angular uncertainty of $3^{\circ}$ of HAWC+ (Harper et
al., 2018).

Table 1 shows several rejected sets for some galaxies, where $13$% of the data
was removed from the full sample. In general, these sets were removed due to
tracking issues during the observations, producing WCS offsets in the
timestreams of the detector pixels that affect the map-making process. The
worst case scenario was found during the observations of NGC 2146 at $89$
$\mu$m, flight F774. The object was found to have an offset from the position
of the boresight of the array for each HWP PA during the full length of
observations. Post-correction of the WCS as a function of the HWP PAs: a) did
not fully correct the angular offsets, and b) produced a final product with
polarization maps incompatible with the other flights. Given these issues, NGC
2146 observations in Band A during flight F774 were removed from the final
data products.

###  2.3 Zero-level background correction

HAWC+ measures the power of the emissive and variable atmosphere and the
astrophysical object. The data reduction scheme described above may produce
regions of negative flux in areas of extended and low surface brightness due
to the similar levels of noise and astrophysical signal. Thus, characterizing
and estimating the zero-level background across the FOV of the observations of
galaxies is imperative in mitigating the need to approximating and adding the
lost flux back to the full image later.

We have determined and corrected the zero-level background of our observations
as follows. Since all galaxies are smaller than the HAWC+ FOV for a given band
and they are isolated objects without large-scale extended thermal emission,
the scan amplitudes (Table 1) were selected to have at least 1/3 of the final
FOV with ‘true’ zero-level background counts. In addition, the polarization
skydips (Section 2.5) show that the sky is unpolarized. Thus, the sky
background around the galaxies represents the ‘true’ zero-level, where flux
and polarimetric calibration can be performed. Specifically, pixels with a SNR
in total intensity $\geq 3$ were masked for each of the scans per HWP PA and
$\mathcal{R}$ and $\mathcal{T}$ arrays produced by crush (Figure 1-second
row). The masked background was fitted with a second-order surface and then
added to the unmasked image, which ensures a positive and flat background
across the full FOV. The second-order surface has six free parameters,
$f(x,y)=Ax^{2}+By^{2}+Cxy+Dx+Ey+F$, that were fit to the Nyquist-sampled
images using several hundreds pixels. Using the Nyquist sampling or an
individual pixel per beam did not affect the final fitting. A second-order
surface provides better results than a first-order (i.e. flat) surface. The
background has similar curvature as the exposure maps, which have an exposure
time that varies with the distance from the center of the image. This radial
variation and the weighting computed by crush per iteration produce the
observed curvature in the background. As a test, this approach was also
applied to unpolarized objects (i.e. planets), and the measured fluxes and
instrumental polarization were found to be compatible with the C2N and skydips
observing modes (Section 2.5). We estimated that the total flux from the
second-order surface is $\leq 1\sigma$ of the pixel-to-pixel variation within
each of individual scans, and the polarization associated with it is lower
($\leq 0.3\%$) than the instrumental polarization ($\sim 1.6-2.1$%, Section
2.5) of HAWC+ in any given band.

A similar approach was used in Taurus/L1495 by Li et al. (2021), where a
simpler approach to estimating the background in a C2N observation was
performed. For these observations, an adjacent region equal to the FOV of the
HAWC+ array close to the molecular cloud was identified using Herschel images.
Using the HAWC+ observations, the mean of the background was estimated and
considered as the zero-level background. Finally, the mean was added to the
full FOV of the observations. These authors found that this technique
contributed $\sim 14$% to the polarized flux in their science observations.
Here, we use a second-order surface to fit the masked image after the removal
of the astrophysical object, yields lower artificial polarized flux and an
optimal correction of the background structure. At this stage, eight different
scans with positive and flat background are computed (Figure 1-third row).

###  2.4 Precipitable water vapour correction

The standard pipeline corrects for the precipitable water vapour (PWV) based
on the altitude of the aircraft at the time of the observations. As real-time
PWV data are not available, this correction is performed using a model
developed by the SOFIA Science Center that provides flux calibrations within a
maximum of $20$% uncertainty. If PWV variability and/or bad weather conditions
were present during observations, further correction may be required to
minimize flux variations that can affect the final flux calibration and/or
polarization measurements.

We examined the flux variations of the scans per HWP PA and array after
background correction for all galaxies as a function of the altitude and
weather conditions. All galaxies show observed flux variations $\leq 15\%$ at
all bands independently of the altitude and weather conditions. Figure 2 shows
an example of flux variations as a function of time and aircraft altitude for
all HAWC+ bands for M82. This figure demonstrates that flux variations are
negligible during climbing in the flight profile. Also, the relative flux
contributions from HWP PAs and arrays are mostly constant, although the
relative mean flux of the full scan may vary within the maximum 20% flux
uncertainty. This result implies that the intrinsic polarization of the source
is conserved relative to the altitude and weather conditions within a set of
four HWP PAs.

Figure 2: Example of flux variation as a function of time and aircraft
altitude for M82 at all HAWC+ bands. Peak fluxes (first and third row) of the
scans per HWP PA and $\mathcal{R}$ and $\mathcal{T}$ arrays at 53 (top left),
89 (top right), 154 (bottom left) and 214 (bottom right) $\mu$m of M82. The
aircraft altitude (second and fourth rows) at the time of data acquisition is
shown. the vertical lines show the boundaries between different flights.

The only exceptions are M83 and NGC 6946. We found that these galaxies have
strong flux variations during several days of observations that cannot be
explained by altitude changes or instrumental performance. Figure 3 shows the
measured flux variations for M83 at 154 $\mu$m. Flights F735 and F736 have
flux variations that exceed the nominal $20$% flux uncertainties provided by
HAWC+. These flux variations are not directly related to changes in the
aircraft altitude (Figure 3-A vs. 4-Bottom). Instead, changes in the PWV
during the data acquisition explain the flux variations. Specifically, Flight
F731 has a constant mean flux across the 1.5h of observations with also a
constant PWV of $\sim 11$ $\mu$m. However, the flux increases by a factor of
1.8 in flight F735, while the altitude is constant across two thirds of the
observations with a climb of 2000ft. and 1000ft. at the beginning and end of
the 2.5h of observations. The expected PWV for flight F735 also decreases by a
factor 1.8 during the observations (Figure 4-Top) and shows a bump at
approximately halftime of the observations. This behavior causes an increase
in flux as well as a bump in the flux measured by HAWC+, which is in agreement
with the measured fluxes and the factor of the increased fluxes (Figure 3-A
vs. 4-Top, as well as Figure 5).

We correct the PWV variations for each flight as follows. Fluxes for flights
F735 and F736 were fitted with a polynomial function of order 3 and rescaled
to the mean flux from flight F731. Figure 5-left shows the measured
(uncorrected) fluxes with the fitted polynomial function. To conserve the
relative contribution of the fluxes between $\mathcal{R}$ and $\mathcal{T}$
scans for each HWP PA per set, the polynomial function was fit using all
measured fluxes as shown in Figure 5. The fit was normalized to the mean
measured flux from F731, which shows a constant flux within the flight at the
highest altitude from the full set of observations. Finally, each measured
flux from each set was multiplied by the normalization factor to correct for
the PWV variation (Figure 5-right). Final corrected fluxes are shown in Figure
3-B. In average, fluxes vary by $\sim 15$% across the full observations with
some outlier sets at $\sim 20$%. The same methodology was applied to NGC 6946
(Appendix A). At this stage, eight different scans with positive background
and corrected by PWV varations are computed (Figure 1-third row).

Figure 3: Flux variation and aircraft altitude as a function of time for the
observations of M83 at 154 $\mu$m. A) Measured peak fluxes of the scans per
HWP PA and $\mathcal{R}$ and $\mathcal{T}$ arrays; B) PWV corrected fluxes
using the fits shown in Figure 5. The mean (black solid line) and $2\sigma$
uncertainty (black dotted line) are shown. C) Aircraft altitude at the time of
data acquisition. The PWV during the observations of M83 (green for F731 and
cyan for F735 shadowed regions) is shown in Figure 4.  Figure 4: PWV variation
as a function of time for the flights associated with M83 at 154 $\mu$m (see
Figure 3). Top panel: expected PWV for flight F735. The PWV during the
observations of M83 (cyan shadowed region) is shown. Bottom panel: expected
PWV for flight F731. The PWV during the observations of M83 (green shadowed
regions) is shown.  Figure 5: Fluxes uncorrected and corrected due to PWV
variations during observations of M83. Measured fluxes (left column) are fit
with a polynomial function of order 3 (black solid line) for flights F0735
(top), F0736 (bottom). Corrected fluxes (left column) by the normalized fit to
the mean flux of flight F731 are shown.

###  2.5 Stokes parameters and instrumental polarization

The Stokes parameters $IQU$ were estimated using the double difference method
(equations shown in Figure 1) in the same manner as the standard C2N
observations carried by HAWC+ described in Section 3.2 by Harper et al.
(2018). In Fig.1, $N$ is the number of sets, and $\mathcal{R}_{\phi}$ and
$\mathcal{T}_{\phi}$ are the total intensity images of the $\mathcal{R}$ and
$\mathcal{T}$ arrays per HWP PA. The instrumental polarization (IP) was
corrected using the polarization skydips. Specifically, the polarization
skydips were performed by moving the telescope from $57^{\circ}$ to
$23^{\circ}$ in elevation while the HWP rotates at a constant speed (Harper et
al., 2018). These observations provide an estimate of the Stokes $qu$ per
detector pixel across the full array (Appendix B, Figure 14). The final IP was
estimated as the weighted median of all pixels and the associated uncertainty
was estimated as the weighted standard deviation of the median. Figure 6 shows
the estimated IP as a percentage of the normalized Stokes $qu$ in the
instrument reference frame. The numerical values and their uncertainties are
shown in Table 3. We estimate an IP variation of $\leq 0.2$% across the FOV
with a median IP of $1.6-2.1$% within the wavelength range of $53-214$ $\mu$m.

For all galaxies, the IP was corrected by subtraction a constant Stokes $qu$
from the measured Stokes $qu$ of each galaxy. We use the normalized Stokes
$qu$ from the polarization skydips. The IP correction was performed in the
instrument reference frame before the rotation to the sky coordinates. This
correction has a polarization uncertainty of $\leq 0.2$% across the FOV in any
band (Harper et al., 2018, and Section 3.3, Figure 14). We include an
uncertainty of $0.2$% in the polarization fraction in quadrature to the
uncertainties estimated in the scientific analysis of this project. Finally,
for each set, we rotate the Stokes $QU$ from the instrument to the sky
coordinates. At this stage, the Stokes $QU$ corrected by IP are in sky
coordinates for each set (Figure 1-fourth row).

To quantify the feasibility of our data reduction scheme, we performed OTFMAP
polarimetric observations of planets (i.e. Neptune and Uranus) in all bands
during 2021. We reduced these observations using the same approach as
described above. Specifically, for each set of four HWP PAs, we reduced the
raw data before estimating the Stokes $IQU$ images. Then, Stokes $QU$ were
normalized such that $q=Q/I$ and $u=U/I$. We computed the median and standard
deviation of $q$ and $u$ within the beam size centered at the peak intensity
of Stokes $I$ for each band. The final IP was estimated as the weighted median
of all measurements and the associated uncertainty was estimated as the
weighted standard deviation of the median. As an example, Appendix B (Figure
15) shows the Stokes $I$ and normalized Stokes $qu$ for Neptune at $53$
$\mu$m. The central beam (white circle) used to estimate the Stokes parameters
is shown. After we apply the same IP correction to the planets, we estimate a
residual polarization of $\leq 0.3$%.

Table 3: Instrumental polarization from OTFMAP, C2N, and skydips. _Columns,
from left to right:_ a) Band name, b) percentage Stokes $qu$ for the OTFMAP
observing mode, c) percentage Stokes $qu$ for the C2N observing mode, and d)
percentage Stokes $qu$ for polarization skydips.

Band | OTFMAP | C2N | Skydips
---|---|---|---
|  $q$ (%) |  $u$ (%) |  $q$ (%) |  $u$ (%) | q (%) | u (%)
(a) | (b) | (c) | (d) | (e) | (f) | g
| 2021† | 2017 | Oct-Nov 2017
A | $-1.57\pm 0.43$ | $-0.36\pm 0.40$ | $-1.54\pm 0.30$ | $-0.30\pm 0.30$ | $-1.60\pm 0.07$ | $-0.38\pm 0.02$
C | $-1.59\pm 0.35$ |  +$1.14\pm 0.24$ | $-1.51\pm 0.30$ |  +$0.90\pm 0.30$ | $-1.64\pm 0.08$ |  +$0.82\pm 0.07$
D |  +$0.13\pm 0.41$ |  +$2.12\pm 0.52$ |  +$0.28\pm 0.30$ |  +$1.91\pm 0.30$ |  +$0.14\pm 0.14$ |  +$1.89\pm 0.11$
E | $-1.30\pm 0.34$ | $-1.24\pm 0.34$ | $-1.29\pm 0.30$ | $-1.11\pm 0.30$ | $-1.09\pm 0.16$ | $-1.41\pm 0.11$

††footnotetext: The year in which the data were obtained, with polarization
measured using one of the three listed techniques.

Figure 6: Instrumental polarization for the C2N and OTFMAP observing modes at
all bands. IP for the C2N (circles), OTFMAP (stars), and polarization skydips
(squares) for bands A (blue), C (red), D (orange), and E (grey) are shown.
Definitions and details of the estimation of the IP are described in Section
2.5 and 3.3.

###  2.6 Flux calibration and data merging

The flux is calibrated using observations of standard sources (i.e. planets,
asteroids). The SOFIA Science Center routinely observes standard sources
through each observing run and derives the flux calibration factors to be
applied to each band. Harper et al. (2018) describes the flux calibration
processes in more detail. We used the standard flux calibration factors
(associated uncertainties are $\leq 15$%) provided by the SOFIA Science Center
for each observation. At this stage, Stokes $IQU$ in sky coordinates for each
set are flux calibrated (Figure 1-fifth row).

Multi-flight and multi-cycle observations were merged. Stokes $IQU$ were
merged using an adaptive weighted average555A detailed description of the
adaptive weighting algorithm used by SOFIA can be found at https://sofia-
usra.github.io/sofia_redux/sofia_redux/toolkit/resampling/parameters.html#adaptive-
weighting within a common grid with a specified pixel size, where North is up
and East is left. The pixelization of this grid (i.e. final data product) is
set to match the detector pixel size in a given band (Table 2), a value
equivalent to Nyquist sampling. For each pixel, the weighted average within
the FWHM of the associated band is estimated. Thus, the final data product has
a Nyquist sampling pixelization with correlated fluxes within the FWHM of the
band. The error maps are estimated from the input variances for the pixels
involved in each weighted average.

The polarization fraction, $P$, position angle of polarization, $PA$, and
polarized flux, $PI$, are derived from the merged Stokes $IQU$. The
polarization fraction is then debiased,
$P^{\prime}=\sqrt{P^{2}-\sigma_{p}^{2}}$, and corrected for polarization
efficiency, $P^{\prime\prime}=P^{\prime}/P_{eff}$, where $\sigma_{p}$ is the
polarization uncertainty, and $P_{eff}$ is the fractional polarization
efficiency of $0.842,0.939,0.975,0.978$ at $53$, $89$, $154$, and $214$
$\mu$m, respectively. HAWC+ has an absolute error of $3^{\circ}$ in
polarization angle and $0.4$% in polarization fraction (Harper et al., 2018).

Final Stokes $I$ values were spatially cross-correlated with PACS/Herschel
images at the closest wavelength to correct the WCS of the HAWC+ images. Small
corrections of $1-2$ pixels were required for all HAWC+ observations to match
with the central core of the galaxies observed with Herschel. These small
offsets were applied to the maps of the Stokes $IQU$, $P$, $PA$, $PI$, and
their associated uncertainties. At this stage, Stokes $IQU$, $P$, $PA$, $PI$
values and their associated uncertainties are computed and ready for
scientific analysis (Figure 1-sixth row).

###  2.7 Data products

The data products generated within this data release are publicly available on
the SOFIA Legacy Program website via http://galmagfields.com/. The database
contains 14 galaxies observed at multiple wavelengths as shown in Table 4. A
total on-source time of 47.54h (51.34h including overheads) is presented,
which corresponds to 33% of the total awarded time of 155.70h (including
overheads) for this SOFIA Legacy Program. An additional 11.84h (including
overheads) of observations for Centaurus A (Lopez-Rodriguez, 2021), Circinus,
M51 (Borlaff et al., 2021), and NGC 1097 (Lopez-Rodriguez et al., 2021) are
included in this data release. Each file contains Stokes $IQU$, $P$, $PA$,
$PI$ values and their associated uncertainties with a pixel scale equal to the
half-beam size for a given band. The released files have the same file format
as table 2 of Gordon et al. (2018).

In this manuscript, we focus on the characterization of the new observing
technique and show the overall status of the observations and the polarization
maps of each object. A detailed quantitative analysis of the polarization
fraction and magnetic field orientation will be presented in Papers IV and V
of the series. Figures 7, 8, and 9 show the polarization maps of the released
data. Polarization measurements are shown with a constant length to display
the magnetic field orientation. The polarization measurements (E-vector) were
rotated by $90^{\circ}$ to show the B-field orientation. We selected
polarization measurements with $P/\sigma_{P}\geq 3$, $PI/\sigma_{PI}\geq 3$,
$P\leq 30$%, and $I/\sigma_{I}\geq 50$ for all galaxies, except for M82, NGC
253, and NGC 2146 with $I/\sigma_{I}\geq 80$, where $\sigma_{I}$ and
$\sigma_{PI}$ are the uncertainties in Stokes $I$ and polarization intensity,
respectively.

Table 4: Status of released galaxy sample. _Columns, from left to right:_ a)
Galaxy name. b) Central wavelength of the band. C) On-source requested time.
D) Observed on-source time. E) completed fraction per object per band.

Galaxy | Band | Requested | Observed | Completed⋆
---|---|---|---|---
|  | Time† | Time† |
| ($\mu$m) | (h) | (h) |
(a) | (b) | (c) | (d) | (e)
Centaurus A | 89 | - | 0.89 | -
Circinus | 53 | - | 0.11 | -
| 89 | - | 0.89 | -
| 214 | - | 0.29 | -
M51 | 154 | - | 2.78 | -
M82 | 53 | 2.00 | 1.60 | ✓
| 89 | 2.00 | 1.65 | ✓
| 154 | 2.00 | 1.99 | ✓
| 214 | 2.00 | 1.13 | ✓
M83 | 154 | 6.80 | 6.60 | ✓
NGC 253 | 89 | 3.00 | 1.00 | 33%
| 154 | 5.00 | 0.89 | 18%
NGC 1068 | 53 | 3.07 | 1.06 | 35%
| 89 | 7.07 | 2.00 | 28%
NGC 1097 | 53 | - | 1.87 | -
| 154 | - | 0.27 | -
NGC 2146 | 53 | 3.00 | 1.60 | 53%
| 89 | 3.00 | 2.02 | 67%
| 154 | 3.00 | 2.40 | 80%
| 214 | 3.00 | 2.33 | ✓
NGC 3627 | 154 | 6.80 | 4.53 | 67%
NGC 4736 | 154 | 6.80 | 2.22 | 33%
NGC 4826 | 89 | 6.80 | 1.22 | 18%
NGC 6946 | 154 | 6.80 | 7.20 | ✓
NGC 7331 | 154 | 6.80 | 6.10 | ✓

††footnotetext: On-source times. Overhead of 1.08 for OTFMAP observations
(Table 5), except for M51 with a specific overhead of 2.59 for C2N
observations. Observations from other SOFIA programs are labeled as ‘-’.

⋆⋆footnotetext: the final observations are considered as completed (i.e. ✓)
when a) completion $>85$% or b) a desired SNR was reached. Observations from
other SOFIA programs are labeled as ‘-’.

Figure 7: Active galactic nuclei. Centaurus A at 89 $\mu$m, Ciricinus at 53,
89, and 214 $\mu$m, NGC 1068 at 53 and 89 $\mu$m, and NGC 1097 at 89 and 154
$\mu$m are shown. Total intensity (colorscale) in logarithmic scale with
contours increasing in steps of $2^{n}\sigma$, with $n=5,5.5,6,\dots$ The
measurements of polarization fraction (white lines) are shown with constant
length to display the B-field orientations. Only measurements of polarization
fraction with $P/\sigma_{P}\geq 3$ and $I/\sigma_{I}\geq 80$ were selected. A
physical scale (black line) of $0.5$, $1.0$ or $2.0$ kpc and the beam size
(red circle) are shown.  Figure 8: Starburst galaxies. M82 at 53, 89, 154, and
214 $\mu$m, NGC 2146 at 53, 89, 154, and 214 $\mu$m, and NGC 253 at 89 and 154
$\mu$m are shown. Total intensity (colorscale) is shown in logarithmic scale
with contours increasing in steps of $2^{n}\sigma$, with $n=5,5.5,6,\dots$ The
measurements of polarization fraction (white lines) are shown with constant
length to display the B-field orientations. Only measurements of polarization
fraction with $P/\sigma_{P}\geq 3$ and $I/\sigma_{I}\geq 80$ were selected. A
physical scale (black line) of $1.0$ kpc and the beam size (red circle) are
shown.  Figure 9: Spiral galaxies. M51, M83, NGC 3627, NGC 4736, NGC 6946, and
NGC 7331 at 154 $\mu$m, and NGC 4826 at 89 $\mu$m are shown. Total intensity
(colorscale) is shown in logarithmic scale with contours increasing in steps
of $2^{n}\sigma$, with $n=5,5.5,6,\dots$ The measurements of polarization
fraction (white lines) are shown with constant length to display the B-field
orientations. Only measurements of polarization fraction with
$P/\sigma_{P}\geq 3$ and $I/\sigma_{I}\geq 80$ were selected. A physical scale
(black line) of $1.0$ kpc and the beam size (red circle) are shown.

##  3 OTFMAP and C2N comparison

This section shows a quantitative comparison of the overheads and
sensitivities between the C2N and OTFMAP observation modes of HAWC+.

###  3.1 Observing overheads

The OTFMAP observational strategy for polarimetry has been described in
Section 2.2.1. Here, we briefly describe the specifics of the C2N polarimetric
observations using HAWC+ (for further details see Harper et al., 2018). The
telescope with the secondary mirror on the optical axis points to the science
object. Then, the secondary mirror is moved at a frequency of $10.2$ Hz
between the science object and an adjacent position on the sky. This movement
is characterized by its amplitude (i.e. chop-throw) and direction (i.e. chop-
angle). Each pair of science and sky images are subtracted to eliminate the
background emission from the science image. To minimize the radiative offset,
the telescope is moved every $30-50$s (i.e. nod time) with the same amplitude
and direction as the chop. Asymmetric chop-nod is not performed with HAWC+ due
to the associated large overheads. The final image is computed after a pair of
subtracted chop-nod observations, which produces a positive image of the
science object at the location of the boresight within the FOV of HAWC+.
Chopping within the array is not performed due to the small FOVs (Table 2).
This procedure is repeated per HWP PA in a sequence of $5^{\circ}$,
$50^{\circ}$, $27.5^{\circ}$, and $72.5^{\circ}$. As the yield of operational
detectors is $\sim 70$% (Harper et al., 2018), these sets of observations are
repeated in a 4-point dither sequence with a 3 pixel offset around the
boresight in a given HAWC+ band. Before and after the sequence of four HWP
PAs, internal calibration images (i.e. INTCAL) of $15$s each are taken to
calibrate the C2N observations and compute observing flatfields.

The total observing time is estimated using the start and end time of
observations for each file, and the total on-source times are computed as
follows. The total on-source time for the C2N polarimetric mode is estimated
as

$t_{\mbox{on-source,C2N}}=\frac{\mbox{nodtime}}{2}\times
n_{\mbox{dithers}}\times n_{\mbox{HWPPAs}}$ (1)

where nodtime is the nodding time where half of the time is spent on the
science object, $n_{\mbox{dithers}}$ is the number of dither positions, and
$n_{\mbox{HWPPAs}}$ is the number of HWP PAs.

The total on-source time for the OTFMAP polarimetric mode is estimated as

$t_{\mbox{on-source,OTFMAP}}=t_{\mbox{ScanDuration}}\times n_{\mbox{HWPPAs}}$
(2)

where $t_{\mbox{ScanDuration}}$ is the scan duration per scan as shown in
Table 1.

We compare the observing overheads of the C2N and OTFMAP polarimetric
observations of NGC 1068 and Circinus at $89$ $\mu$m. The C2N observations of
NGC 1068 were presented by Lopez-Rodriguez et al. (2020) with a total
observing time of $7615$s and a total on-source time of $2910$s. The chop-
throw is $180$″, and the chop-angle is $90^{\circ}$ to always be along the
short axis of the array. The C2N observations of Circinus were performed
during the same flight, F596, as those of the OTFMAP polarimetric observations
with a total observing time of $615$s and a total on-source time of $252$s.
The chop-throw is $180$″, and the chop-angle is $322^{\circ}$ counter-
clockwise in the East of North direction

A subgroup of the OTFMAP observations presented in Table 1 were reduced to
have on-source times similar to the C2N observations. OTFMAP polarimetric
observations with an on-source time of $2880$s and $400$s were reduced
following the steps described in Section 2 for NGC 1068 and Circinus,
respectively. The associated total observing times of these subsets of
observations are $3041$s and $428$s for NGC 1068 and Circinus, respectively.
Figures 10 and 15 show the C2N configuration and Lissajous curves from one of
the scans over the Herschel image at $70$ $\mu$m for NGC 1068 and Circinus,
respectively.

The observing overhead factor for the C2N observing mode is estimated to be
$2.53$, while an overhead of $1.08$ is estimated for the OTFMAP observing mode
(Table 5). Then, the total observing time is computed as

$\displaystyle t_{\mbox{C2N}}$ $\displaystyle=$ $\displaystyle 2.53\times
t_{\mbox{on-source,C2N}}$ (3) $\displaystyle t_{\mbox{OTFMAP}}$
$\displaystyle=$ $\displaystyle 1.08\times t_{\mbox{on-source,OTFMAP}}$ (4)

The OTFMAP overhead arises from: a) the time spent to rotate the four HWP PA,
and b) the waiting time of $\sim 2-4$s to ensure tracking is correctly
performed between scans. The C2N overhead arises from a) chopping on the sky
positions and nodding the telescope, b) INTCALs of $15$s each before and after
the sequence of four dither positions, and c) the rotation of the HWP PAs. We
estimate that the OTFMAP polarimetric mode provides a total improvement in
observing time overhead of a factor of $2.53/1.08=2.34$ with respect to the
C2N polarimetric observations. The main reason for this improvement is that
the galaxies are always within the FOV of the HAWC+ array during the whole
integration time.

Figure 10: Comparison between C2N (left) and OTFMAP (right) observations of
NGC 1068 at $89$ $\mu$m. First row: Herschel observations at $70$ $\mu$m
(colorscale) of NGC 1068 within a FOV of $500\times 500$ sqarsec with overlaid
C2N (left) and OTFMAP (right) configurations. For C2N, on-axis (white) and
off-axis (yellow) positions with the FOV of HAWC+ at $89$ $\mu$m are shown.
For OTFMAP, the Lissajous curve (white solid line) of one scan is shown.
Second row: HAWC+ total flux observations for C2N and OTFMAP within the same
FOV as above. Contours start at $32\sigma_{I}$ and increase in steps of
$2^{n}\sigma$, where $n=5,5.5,6,\dots$ and $\sigma_{I}=0.3$ mJy/sqarcsec. The
beam size is shown as a red circle in the bottom left. Third row: Same as
above within a FOV of $120\times 120$ sqarcsec. Polarization measurements
(white lines) are shown for $P/\sigma_{P}\geq 3$, $PI/\sigma_{PI}\geq 3$,
$P<30$%, and $I/\sigma_{I}\geq 50$. A $5\%$ polarization measurement is shown
at the bottom right. Fourth row: HAWC+ polarized flux observations
(colorscale) within the same FOV and polarization as above. Contours start at
$3\sigma_{PI}$ and increase in steps of $1\sigma_{PI}$, where
$\sigma_{PI}=0.46$ mJy/sqarcsec.  Table 5: OTFMAP vs. C2N: Overheads and
sensitivities using observations of Circinus and NGC 1068 at $89$ $\mu$m.
_Columns, from left to right:_ a) Observing mode, b) observing time overhead
(i.e. total exposure time / o-source time), c) minimum intensity flux for
polarimetry (MIfP) to achieve a polarization uncertainty of $0.3$% in $900$s
(including overheads), d) minimum detectable continuum polarized flux (MDCPF)
for a $4\sigma$ detection in $900$s (including overheads), e) same as c) but
only taken into account on-source time, f) same as d) but only taken into
account on-source time.

Obs. Mode | Overhead | with Overheads | without Overheads
---|---|---|---
|  | MIfP† | MDCPF† | MIfP | MDCPF
|  | with $\sigma_{p}=0.3$% in 3600s | with $4\sigma$ in 900s | with $\sigma_{p}=0.3$% in 3600s | with $4\sigma$ in 900s
|  | [MJy/sr] | [% Jy] | [MJy/sr] | [% Jy]
(a) | (b) | (c) | (d) | (e) | (f)
C2N | 2.53 | 7010 | 59 | 4333 | 37
OTFMAP | 1.08 | 2819 | 21 | 2410 | 21
Improvement${}_{OTFMAP}^{\star}$ | 2.34 | 2.49 | 2.81 | 1.80 | 1.76

⋆⋆footnotetext: Improvement of the OTFMAP shown as a ratio between OTFMAP and
C2N of the quantities shown for each column.

††footnotetext: Harper et al. (2018) quote a MIfP and MDCPF, including
overheads, of $6000$ MJy/sr and $50$ % Jy, respectively.

###  3.2 Sensitivities

The sensitivities of the total intensity mode using OTFMAP and polarimetry
mode using C2N have been already characterized (Harper et al., 2018). We
present the sensitivity estimations of the OTFMAP polarimetric mode at 89
$\mu$m (Band C) as an example of the improvement provided by this new
observing mode. Further characterization in different bands will be performed
and released by the SOFIA Science Center. The SOFIA observer’s handbook for
HAWC+666Observer’s Handbook for HAWC+ can be found at
https://www.sofia.usra.edu/instruments/hawc. defines the sensitivities as
follows: the Minimum Intensity flux for Polarimetry (MIfP) is given as the
flux needed to reach $0.3$% polarization uncertainty in $3600$s for an
extended source; and the Minimum Detectable Continuum Polarized Flux (MDCPF)
is given as the total flux needed to reach a $4\sigma$ detection in $900$s for
point sources.

Using the OTFMAP and C2N polarimetric observations of NGC 1068 and Circinus at
$89$ $\mu$m, Figures 10 and 15 clearly show an increase in the number of
polarization measurements with the same statistical significance (i.e.
$P/\sigma_{P}\geq 3$), and a larger spatial extension of the polarized flux
with the same $PI/\sigma_{PI}\geq 3$. Figure 11 and 16 show the histograms of
the total and polarized fluxes, as well as a 1:1 plot between the polarization
fraction and polarization angle for NGC 1068 and Circinus, respectively. For
the 1:1 plot, measurements associated with the same WCS were plotted. We
conclude that the total flux, polarized flux, polarization fraction, and
position angles are in agreement between both OTFMAP and C2N polarimetric
observations. In addition, after the standard deviation of the observations
are normalized using their associated on-source time, we find that the
standard deviation of the polarized flux decreases for the OTFMAP polarimetric
observations (middle right panels of Figures 11 and 16), demonstrating a
significant improvement in sensitivity by the OTFMAP polarimetric mode.

Using the observing overheads presented in Section 3.1, we estimate a MIfP of
$7010$ MJy/sr and MDCPF of $59$ % Jy in C2N mode at $89$ $\mu$m. These values
are in agreement, within the nominal $20$% flux uncertainty, with the quoted
MIfP of $6000$ MJy/sr and $50$ % Jy by Harper et al. (2018) and the HAWC+
observer’s handbook. The estimated MIfP and MDCPF for the OTFMAP polarimetric
mode are shown in Table 5. Finally, we estimate a sensitivity improvement of
$1.80$ in MIfP, without accounting for the observing overheads. After
accounting for observing time overheads, the OTFMAP polarimetric mode offers
total improvement in MIfP of $2.49$ times over the C2N observing mode.

Figure 11: Sensitivity comparison between C2N (red) and OTFMAP (blue)
observations of NGC 1068 at $89$ $\mu$m. Histograms of the total intensity
(top left), polarized intensity (top right), standard deviations of the total
flux (middle left) and polarized flux (middle right) are shown. The median and
$1\sigma$ uncertainty of the distribution of the standard deviations are shown
in the middle panels. The polarization fraction (bottom left) and polarization
angle (bottom right) for C2N and OTFMAP observations are shown. A 1:1 line
(black dashed line) is shown.

###  3.3 Instrumental polarization

The IP was estimated with unpolarized standard objects using C2N, OTFMAP, and
polarization skydips. The IP from C2N and polarization skydips were collected
in 2017 and computed by Harper et al. (2018). The polarization from OTFMAP of
planets were presented and computed in Section 2.5. Figure 6 shows the
estimated IP as a percentage of the normalized Stokes $qu$ in the instrument
reference frame for the C2N, OTFMAP, and polarization skydips. The numerical
values and their uncertainties are shown in Table 3. We find that all three
methods produce consistent and reproducible results. The uncertainties are
larger in C2N and OTFMAP than in the polarization skydips due to the number of
measurements used for the statistics. We take the measurements from the
polarization skydips for the IP correction in our observations (Section 2.5).

We find a rotation in the IP with wavelength. The HAWC+ filters are metal
grids (Harper et al., 2018), which probably introduce some level of intrinsic
polarization. However, the polarization properties of these filters were not
measured before they were installed in the instrument. The important result
here is that the IP is reproducible and consistent throughout the different
methods, and time and conditions of observations. The IP is known to arise
from the tertiary mirror (Harper et al., 2018). Thus, our results show that
the IP remains constant, within the uncertainties, from 2017 to 2021. This
result suggests that the physical conditions of the tertiary mirror produce a
negligible difference in IP during five years of observations.

###  3.4 Further implementations of the OTFMAP polarization for HAWC+

Although we have shown that the OTFMAP polarimetric mode has substantially
improved the overheads and sensitivities when compared to the C2N polarimetric
mode, further developments on the OTFMAP polarimetric mode are still required.
As mentioned above, the data processing presented in Section 2 works for
objects with thermal emission within the FOV of HAWC+ in a given band.
However, this procedure may not be as efficient for large-scale and diffuse
thermal emission that cover areas much larger than the FOV of the HAWC+ array.
This well-known deficiency of the OTFMAP observing mode can be mitigated by
application of alternative OTFMAP strategies in the recovery of large-scale
emission, as discussed below.

For FIR polarimetric observations, Hildebrand et al. (2000) suggested an
observing strategy involving the telescope chopping to an adjacent sky
location at the same time as scanning. As in C2N, the chop is then subtracted
from the closest scan in time. While this technique may remove the sky
fluctuations more efficiently than the procedure we presented earlier, it will
increase the observing overheads and reduce the sensitivities as there may be
a residual radiative offset after the chop subtraction. An alternative
strategy consists of a continuous rotation of the HWP to modulate signal at a
faster rate, usually $\sim 2$ Hz, than the atmospheric fluctuations Note that
HAWC+ can continuously rotate the HWP to $0.5$ Hz, providing a polarization
modulation frequency of $2$ Hz (Harper et al., 2018). The time streams for
each pixel can then be linearly decomposed to IP, astrophysical polarization,
detector gains, detector temperature fluctuations, and atmospheric
transparency fluctuations. This technique is commonly used in sub-mm
polarimetric observations, for example by POL-2 on the James Clerk Maxwell
Telescope (JCMT; Ward-Thompson et al., 2017), and cosmic microwave background
(CMB) experiments (e.g., Johnson et al., 2007). While this observing approach
may be the most efficient, its implementation requires new software and
hardware development. These common issues are also well-known problems in
radio polarimetric observations using scanning techniques. Müller et al.
(2017) applied a technique called ’basket-weaving’ to maps with orthogonal
scanning on the sky, and Emerson & Graeve (1988) applied a technique called
’PLAIT’ for scans oriented at any angle. Although evaluating these observing
modes is outside of the scope of this manuscript, we have mentioned them as
options for future implementations that could reduce the observational
overheads and, as a result, SOFIA’s operational cost.

##  4 Conclusions

We present the first data release (DR1) of the SOFIA Legacy program (PI:
Lopez-Rodriguez, E. & Mao, S. A.) on extragalactic magnetism. The DR1 consists
of 14 galaxies, including active galactic nuclei, starbursts, and spirals,
from 53 to 214 $\mu$m, which comprises 33% (51.43h out of 155.7h) of the total
awarded time of this program. We release homogeneously reduced high-level data
products, ready for scientific analysis.

We have presented the data processing of the OTFMAP polarization mode for
HAWC+. The pipeline steps were applied to the newly acquired observations of
galaxies from January 2020 to December 2021 within the SOFIA Legacy program on
extragalactic magnetism. This new observing mode had been successfully applied
to objects smaller than the FOV of the HAWC+ array in any given band. The data
processing includes: Zero-level background subtraction, PWV correction,
instrumental polarization subtraction, and merging. Comparison with C2N
observations show that the OTFMAP has greatly improved the polarimetric mode
of HAWC+. Specifically, we estimated that OTFMAP polarimetric mode provides an
improvement in observing time overhead of a factor of $2.34$, and an
improvement of $1.80$ in MIfP without accounting for the observing overheads
from the C2N polarimetric observations. Including these two factors, we
estimate that the OTFMAP polarimetric mode offers a total improvement in MIfP
of $2.49$ times over the C2N observing mode. The OFTMAP is a significant
optimization of the polarimetric mode of HAWC+ as it ultimately reduces the
cost of operations of SOFIA/HAWC+–more data can be collected per hour of
flight, more than doubling the number of programs/papers per hour of
observation.

Although the polarimetric mode of HAWC+ has been improved, we emphasize that
further development of this observing mode for objects with large-scale
diffuse polarized emission is still required. We briefly described several
observing strategies to improve this technique for HAWC+ based on the
approaches from sub-mm observations and CMB experiments.

Based on observations made with the NASA/DLR Stratospheric Observatory for
Infrared Astronomy (SOFIA) under the 07_0034, 08_0012 Program. SOFIA is
jointly operated by the Universities Space Research Association, Inc. (USRA),
under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under
DLR contract 50 OK 0901 to the University of Stuttgart. KT has received
funding from the European Research Council (ERC) under the European Unions
Horizon 2020 research and innovation programme under grant agreement No.
771282. EN is supported by the ERC Grant ”Interstellar” (Grant agreement
740120) and has received funding from the Hellenic Foundation for Research and
Innovation (H.F.R.I., project No 224). SMA is supported by the ERC Starting
Grant 638707 “Black holes and their host galaxies: co-evolution across cosmic
time” and by STFC

##  Appendix A PWV correction of NGC 6946

NGC 6946 shows flux variations in several flights due to PWV fluctuations
during the observations. Following the approach described in Section 2.4, the
measured fluxes (Figure 12-top) were corrected by fitting a linear function
(Figure 13) normalized to the mean fluxes of flights F684 and F685. The
corrected fluxes (Figure 12-middle) show a flux variation of $15$% across the
full observations.

Figure 12: Flux variation as a function of time and aircraft altitude for
NGC6946 at 154 $\mu$m. Top panel: Measured peak fluxes of the scans per HWP PA
and $\mathcal{R}$ and $\mathcal{T}$ arrays. Middle panel: PWV corrected fluxes
using the fits shown in Figure 13. The mean (black solid line) and $2\sigma$
uncertainty (black dotted line) are shown. Bottom panel: aircraft altitude at
the time of data acquisition.  Figure 13: Fluxes uncorrected and corrected by
a linear fit due to PWV variations during observations of NGC 6946. Measured
fluxes (left column) are fit with a linear function (black solid line) for
flights F0683 (top), F0689 (middle), and F0690 (bottom). Corrected fluxes
(right column) by the normalized fit to the mean flux of flights F684 and F685
are shown.

##  Appendix B Instrumental polarization

This section shows the supporting figure described in Section 2.5 and 3.3.

Figure 14: OTFMAP polarimetric observations of a planet (Neptune) and
polarization skydip (poldip) at $53$ $\mu$m. Top panels: Polarization
observations of a planet. Stokes I (left) in $\log_{10}$ scale, normalized
Stokes $q=Q/U$ (middle) and $u=U/I$ (right) are shown. The FOV is $40\times
40$ pixels ($48\times 48$ sqarcsec, with pixelscale of $1\farcs 21$). Bottom
panels: Polarization observations of a skydip. Normalized Stokes $q=Q/U$
(middle) and $u/I$ (right) are shown. The FOV is $32\times 40$ pixels
($82\times 102$ sqarcsec, with a pixelsxale of $2\farcs 55$). Definitions and
details of the estimation of the IP are described in Section 2.5 and 3.3.

##  Appendix C OTFMAP vs. C2N for Circinus

OTFMAP and C2N comparisons for the observations of Circinus at 89 $\mu$m.
Figure 15 shows the instrumental configuration and polarization maps for both
observing modes. Figure 16 presents the histograms of the total and polarized
flux, and their standard deviations. The 1:1 plots of the polarization
fraction and polarization angle are shown.

Figure 15: Comparison between C2N (left) and OTFMAP (right) observations of
Circinus at $89$ $\mu$m. First row: Herschel observations at $70$ $\mu$m
(colorscale) of Circinus within a FOV of $350\times 350$ sqarsec with overlaid
C2N (left) and OTFMAP (right) configurations. For C2N, on-source (white) and
off-source (yellow) positions with the FOV of HAWC+ at $89$ $\mu$m are shown.
For OTFMAP, the Lissajous curve (white solid line) of the scan is shown.
Second row: HAWC+ total flux observations for C2N and OTFMAP within the same
FOV as above. Contours start at $32\sigma_{I}$ and increase in steps of
$2^{n}\sigma$, where $n=5,5.5,6,\dots$ and $\sigma_{I}=0.3$ mJy/sqarcsec. The
beam size is shown as a red circle in the bottom left. Third row: Same as
above within a FOV of $120\times 120$ sqarcsec. Polarization measurements
(white lines) are shown for $P/\sigma_{P}\geq 3$, $PI/\sigma_{PI}\geq 3$,
$P<30$%, and $I/\sigma_{I}\geq 50$. A $5\%$ polarization measurement is shown
at the bottom right. Fourth row: HAWC+ polarized flux observations
(colorscale) within the same FOV and polarization as above. Contours start at
$3\sigma_{PI}$ and increase in steps of $1\sigma_{PI}$, where
$\sigma_{PI}=0.45$ mJy/sqarcsec.  Figure 16: Sensitivity comparison between
C2N (red) and OTFMAP (blue) observations of Circinus at $89$ $\mu$m.
Histograms of the total intensity (top left), polarized intensity (top right),
standard deviations of the total flux (middle left) and polarized flux (middle
right) are shown. The median and $1\sigma$ uncertainty of the distribution of
the standard deviations are shown in the middle panels. The polarization
fraction (bottom left) and polarization angle (bottom right) for C2N and
OTFMAP observations are shown. A 1:1 line (black dashed line) is shown.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Borlaff et al. (2021) Borlaff, A. S., Lopez-Rodriguez, E., Beck, R., et al. 2021, arXiv e-prints, arXiv:2105.09315.  https://arxiv.org/abs/2105.09315
  * Burtscher et al. (2020) Burtscher, L., Politopoulos, I., Fernández-Acosta, S., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114477L, doi: 10.1117/12.2576271
  * Cantalupo et al. (2010) Cantalupo, C. M., Borrill, J. D., Jaffe, A. H., Kisner, T. S., & Stompor, R. 2010, ApJS, 187, 212, doi: 10.1088/0067-0049/187/1/212
  * Chapin et al. (2013) Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545, doi: 10.1093/mnras/stt052
  * Dowell et al. (2010) Dowell, C. D., Cook, B. T., Harper, D. A., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Proc. SPIE, 77356H, doi: 10.1117/12.857842
  * Emerson & Graeve (1988) Emerson, D. T., & Graeve, R. 1988, A&A, 190, 353
  * Gordon et al. (2018) Gordon, M. S., Lopez-Rodriguez, E., Andersson, B. G., et al. 2018, arXiv e-prints, arXiv:1811.03100.  https://arxiv.org/abs/1811.03100
  * Harper et al. (2018) Harper, D. A., Runyan, M. C., Dowell, C. D., et al. 2018, Journal of Astronomical Instrumentation, 7, 1840008, doi: 10.1142/S2251171718400081
  * Haslam (1974) Haslam, C. G. T. 1974, A&AS, 15, 333
  * Hildebrand et al. (2000) Hildebrand, R. H., Davidson, J. A., Dotson, J. L., et al. 2000, PASP, 112, 1215, doi: 10.1086/316613
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 3
  * Johnson et al. (2007) Johnson, B. R., Collins, J., Abroe, M. E., et al. 2007, ApJ, 665, 42, doi: 10.1086/518105
  * Kovács (2006) Kovács, A. 2006, PhD thesis, Caltech
  * Kovács (2008a) Kovács, A. 2008a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7020, Millimeter and Submillimeter Detectors and Instrumentation for Astronomy IV, ed. W. D. Duncan, W. S. Holland, S. Withington, & J. Zmuidzinas, 702007, doi: 10.1117/12.790272
  * Kovács (2008b) Kovács, A. 2008b, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7020, Proc. SPIE, 70201S, doi: 10.1117/12.790276
  * Li et al. (2021) Li, P. S., Lopez-Rodriguez, E., Ajeddig, H., et al. 2021, MNRAS, doi: 10.1093/mnras/stab3448
  * Lopez-Rodriguez (2021) Lopez-Rodriguez, E. 2021, Nature Astronomy, doi: 10.1038/s41550-021-01329-9
  * Lopez-Rodriguez et al. (2020) Lopez-Rodriguez, E., Dowell, C. D., Jones, T. J., et al. 2020, ApJ, 888, 66, doi: 10.3847/1538-4357/ab5849
  * Lopez-Rodriguez et al. (2021) Lopez-Rodriguez, E., Beck, R., Clark, S. E., et al. 2021, ApJ, 923, 150, doi: 10.3847/1538-4357/ac2e01
  * Müller et al. (2017) Müller, P., Krause, M., Beck, R., & Schmidt, P. 2017, A&A, 606, A41, doi: 10.1051/0004-6361/201731257
  * Ohsawa et al. (2018) Ohsawa, R., Sako, S., Miyata, T., et al. 2018, ApJ, 857, 37, doi: 10.3847/1538-4357/aab6ae
  * Patanchon et al. (2008) Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2008, ApJ, 681, 708, doi: 10.1086/588543
  * Reichertz et al. (2001) Reichertz, L. A., Weferling, B., Esch, W., & Kreysa, E. 2001, A&A, 379, 735, doi: 10.1051/0004-6361:20011227
  * Robitaille & Bressert (2012) Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python.  http://ascl.net/1208.017
  * Roussel (2013) Roussel, H. 2013, PASP, 125, 1126, doi: 10.1086/673310
  * Tegmark (1997) Tegmark, M. 1997, Phys. Rev. D, 56, 4514, doi: 10.1103/PhysRevD.56.4514
  * Vaillancourt et al. (2007) Vaillancourt, J. E., Chuss, D. T., Crutcher, R. M., et al. 2007, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6678, Proc. SPIE, 66780D, doi: 10.1117/12.730922
  * Ward-Thompson et al. (2017) Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66, doi: 10.3847/1538-4357/aa70a0
  * Waskett et al. (2007) Waskett, T. J., Sibthorpe, B., Griffin, M. J., & Chanial, P. F. 2007, MNRAS, 381, 1583, doi: 10.1111/j.1365-2966.2007.12327.x
  * Weferling et al. (2002) Weferling, B., Reichertz, L. A., Schmid-Burgk, J., & Kreysa, E. 2002, A&A, 383, 1088, doi: 10.1051/0004-6361:20011617