Extragalactic magnetism with SOFIA (SALSA Legacy Program) -- III: First data release and on-the-fly polarization mapping characterization

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.


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 µm) wavelength range (e.g. Burtscher et al. 2020, and references therein) and sometimes used in the far-IR (FIR; 30 − 300 µ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 ∼ 5 − 30 for MIR in 8-m class telescopes and up to 500 for FIR wavelengths. The cadence of < 1s 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 ≤ 60s 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 ≤ 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 , 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.

Observations
New multi-wavelength polarimetric observations of 10 galaxies were performed as part of the SOFIA Legacy Program 1 (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 µm and Circinus at 53, 89, and 214 µm from program ID 07 0032 (PI: Lopez-Rodriguez, E.), M51 (Borlaff et al. 2021) at 154 µm under programs 70 0509 (Guaranteed Time Observations by the HAWC+ Team), 76 0003 (Discretionary Director Time), and 08 0260 (PI: NGC 1097 (Lopez-Rodriguez et al. 2021) at 89 and 154 µ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

Data reduction of the on-the-fly-map (OTFMAP) polarimetric mode
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 of- 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 • . The sets are reduced using CRUSH, which generates a scan per R and 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, P A, of polarization, and polarized fluxes, P I. The pipeline computes uncertainties and covariances, which are taken into account during all stages of the data reduction. fered 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.

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 2006Kovács , 2008b and the HAWC DRP V2.7.0 pipeline developed by the data reduction pipeline group at the SOFIA Science Center 2 . SOFIA provided a version, HAWC DRP V3.0.0, that replaces CRUSH, written in JAVA, with a new version written in PYTHON. This new version 3 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 × 40 pixels each, labeled as reflective (R) and transmissive (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). 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 • , 50 • , 27.5 • , and 72.5 • . A starting angle of 5 • is chosen to avoid the reference position, 0 • , 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, R and 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 • in regions of high signal-to-noise ratio in Stokes I (SNR I ≥ 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 • 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 µ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 • 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 • 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 µ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.

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 ≥ 3 were masked for each of the scans per HWP PA and R and 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 ≤ 1σ of the pixel-to-pixel variation within each of individual scans, and the polarization associated with it is lower (≤ 0.3%) than the instrumental polarization (∼ 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 ∼ 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).

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 ≤ 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. 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 µ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 ∼ 11 µ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 R and 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 ∼ 15% across the full observations with some outlier sets at ∼ 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). 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 R φ and T φ are the total intensity images of the R and 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 • to 23 • 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 ≤ 0.2% across the FOV with a median IP of 1.6 − 2.1% within the wavelength range of 53 − 214 µm.

Stokes parameters and instrumental polarization
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 ≤ 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 µ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 ≤ 0.3%.

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 ≤ 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 average 4 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, P A, and polarized flux, P I, are derived from the merged Stokes IQU . The polarization fraction is then debiased, P = P 2 − σ 2 p , and corrected for polarization efficiency, P = P /P ef f , where σ p is the polarization uncertainty, and P ef f is the fractional polarization efficiency of 0.842, 0.939, 0.975, 0.978 at 53, 89, 154, and 214 µm, respectively. HAWC+ has an absolute error of 3 • 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 , P A, P I, and their associated uncertainties. At this stage, Stokes IQU , P , P A, P I values and their associated uncertainties are computed and ready for scientific analysis (Figure 1-sixth row).

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 , P A, P I 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 • to show the B-field orientation. We selected polarization measurements with P/σ P ≥ 3, P I/σ P I ≥ 3, P ≤ 30%, and I/σ I ≥ 50 for all galaxies, except for M82, NGC 253, and NGC 2146 with I/σ I ≥ 80, where σ I and σ P I are the uncertainties in Stokes I and polarization intensity, respectively.

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

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 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 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 '-'. 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 − 50s (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 • , 50 • , 27.5 • , and 72.5 • . As the yield of operational detectors is    boresight in a given HAWC+ band. Before and after the sequence of four HWP PAs, internal calibration images (i.e. INTCAL) of 15s 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 on-source,C2N = nodtime 2 × n dithers × n HWPPAs (1) where nodtime is the nodding time where half of the time is spent on the science object, n dithers is the number of dither positions, and n HWPPAs is the number of HWP PAs. The total on-source time for the OTFMAP polarimetric mode is estimated as t on-source,OTFMAP = t ScanDuration × n HWPPAs (2) where t ScanDuration is the scan duration per scan as shown in Table 1 Table 1 were reduced to have on-source times similar to the C2N observations. OTFMAP polarimetric observations with an on-source time of 2880s and 400s 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 3041s and 428s 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 µ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 t C2N = 2.53 × t on-source,C2N (3) t OTFMAP = 1.08 × t on-source,OTFMAP The OTFMAP overhead arises from: a) the time spent to rotate the four HWP PA, and b) the waiting time of ∼ 2 − 4s 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 15s 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.

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 µ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+ 5 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 3600s for an extended source; and the Minimum Detectable Continuum Polarized Flux (MDCPF) is given as the total flux needed to reach a 4σ detection in 900s for point sources.
Using the OTFMAP and C2N polarimetric observations of NGC 1068 and Circinus at 89 µm, Figures 10 and 15 clearly show an increase in the number of polarization measurements with the same statistical significance (i.e. P/σ P ≥ 3), and a larger spatial Table 5. OTFMAP vs. C2N: Overheads and sensitivities using observations of Circinus and NGC 1068 at 89 µ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 900s (including overheads), d) minimum detectable continuum polarized flux (MDCPF) for a 4σ detection in 900s (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. extension of the polarized flux with the same P I/σ P I ≥ 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 µ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.

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.

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 ∼ 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.

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 µ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.     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.  C2N and OTFMAP within the same FOV as above. Contours start at 32σI and increase in steps of 2 n σ, where n = 5, 5.5, 6, . . . and σ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 × 120 sqarcsec. Polarization measurements (white lines) are shown for P/σP ≥ 3, P I/σP I ≥ 3, P < 30%, and I/σI ≥ 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σP I and increase in steps of 1σP I , where σP I = 0.45 mJy/sqarcsec.