Study of the Cosmic Rays and Interstellar Medium in Local H i Clouds Using Fermi-LAT Gamma-Ray Observations

Prior to preparing templates of the ISM gas for the γ-ray data analysis, we investigated the properties of their tracers. As neutral gas tracers, we prepared dust maps and a ${W}_{{\rm{H}}{\rm{I}}}$ map stored in a HEALPix (Górski et al. 2005) equal-area sky map of order 9⁶ (with a mean spacing of $6\buildrel{\,\prime}\over{.} 9$ that is commensurate with the $\sim 5^{\prime}$ resolution of the Planck dust maps). We used the Planck dust maps (of R, ${\tau }_{353}$, and T_d) from the public data release 1 (version R1.20)⁷ described by Planck Collaboration XI (2014).⁸ Assuming a uniform dust temperature along the line of sight, they constructed maps by modeling the dust thermal emission with a single modified blackbody (for details of their procedure, see Planck Collaboration XI 2014). As described in Planck Collaboration XI (2014), the dust optical depth is the product of the dust opacity and the total neutral gas column density ${N}_{{\rm{H}}}$. Therefore, if the dust-to-gas ratio and dust cross section at 353 GHz are spatially constant, ${\tau }_{353}$ is proportional to ${N}_{{\rm{H}}}$. The dust radiance R is also expected to trace the total gas column density because it is proportional to ${N}_{{\rm{H}}}$ under the assumption that the dust-to-gas ratio, the optical and UV absorption cross section of the dust, and the interstellar radiation field are spatially uniform. These hypotheses have to be tested because dust evolution models predict that the dust emission properties change across the ISM (e.g., Roy et al. 2013; Ysard et al. 2015); which severely discourages the use of a single conversion factor without verification. For this reason we investigate the change in the ${N}_{{\rm{H}}}$–${D}_{\mathrm{em}}$ relationship with dust properties using γ-ray data (Section 4). We analyzed the northern ROI ($200^\circ \leqslant l\leqslant 260^\circ$ and $22^\circ \leqslant b\leqslant 60^\circ$) and the southern ROI ($210^\circ \leqslant l\leqslant 270^\circ$ and $-60^\circ \leqslant b\leqslant -22^\circ$). The northern ROI is identical to that adopted by Abdo et al. (2009b), but the southern ROI is shifted by $10^\circ$ (see below). To construct the ${W}_{{\rm{H}}{\rm{I}}}$ map, we used the HI4PI survey data (HI4PI Collaboration 2016) integrated over the full velocity range of the survey (from −600 to $600\,\mathrm{km}\,{{\rm{s}}}^{-1}$). In the ${W}_{{\rm{H}}{\rm{I}}}$ map, we identified several bright radio sources and intermediate-velocity clouds (e.g., Wakker 2001). We removed these radio sources from the ${W}_{{\rm{H}}{\rm{I}}}$ map by filling them with the average of the peripheral pixels. We also identified an area where the ${W}_{{\rm{H}}{\rm{I}}}$–${D}_{\mathrm{em}}$ relation is affected by the contamination of an intermediate-velocity cloud and masked the area when studying the ${W}_{{\rm{H}}{\rm{I}}}$–${D}_{\mathrm{em}}$ relation and γ-ray data analysis. We also masked the Orion–Eridanus superbubble (e.g., Ochsendorf et al. 2015) because the CR and ISM properties inside the bubble can be appreciably different from those in other areas (see Appendix A for details). To compensate the loss of photon statistics because of the mask, we shifted the longitude of the southern ROI by $10^\circ$ toward the positive direction (i.e., away from the mask) from that of Abdo et al. (2009b). In the Planck dust maps, we identified several spots with high $R/{\tau }_{353}$ ratios and high R. They are likely infrared sources and were removed from the dust maps by filling them with the average of the peripheral pixels (see Appendix B). Finally, we examined the Planck type 3 CO map that has the highest signal-to-noise ratio among three types of the map (Planck Collaboration XIII 2014) and confirmed that there is no strong CO 2.6 mm emission in our ROI (see Appendix C). We identified a weak spot at (l, b) ∼ ($221\buildrel{\circ}\over{.} 4$, $45\buildrel{\circ}\over{.} 1$). This spot can also be seen in the R and ${\tau }_{353}$ maps and is therefore likely to be an infrared source. Therefore, we removed the source from the dust maps by filling the source area with the average of the peripheral pixels, and we do not consider CO-bright ${{\rm{H}}}_{2}$ in constructing gas models hereafter.

The correlations between ${W}_{{\rm{H}}{\rm{I}}}$ and R, as well as those between ${W}_{{\rm{H}}{\rm{I}}}$ and ${\tau }_{353}$, are shown in Figures 1 and 2, respectively, in which different colors represent different T_d. To match the resolution of the ${W}_{{\rm{H}}{\rm{I}}}$ map, we smoothed the dust maps using a Gaussian kernel with an FWHM of $15\buildrel{\,\prime}\over{.} 4$. In the northern region, we can confirm two trends of the dust–gas relation found by previous studies of the high-latitude sky described in Section 1: (1) we observe in Figure 1(a) a rather good correlation between ${W}_{{\rm{H}}{\rm{I}}}$ and R over a wide range of T_d, as shown by Planck Collaboration XI (2014), which proposed R as a good neutral gas tracer;⁹ and (2) we observe in Figure 1(b) a stronger T_d dependence of the ${W}_{{\rm{H}}{\rm{I}}}$–${\tau }_{353}$ relationship, which Fukui et al. (2015) interpreted as being primarily due to optically thick ${\rm{H}}\,{\rm{I}}$ in low-T_d areas. In the southern region, while both R and ${\tau }_{353}$ show a small T_d dependence in the relation with ${W}_{{\rm{H}}{\rm{I}}}$, a possible nonlinear relationship between ${W}_{{\rm{H}}{\rm{I}}}$ and ${D}_{\mathrm{em}}$ (a break in the ${W}_{{\rm{H}}{\rm{I}}}$/${D}_{\mathrm{em}}$ ratio) can be identified, particularly for the case of ${\tau }_{353}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The correlation between the dust tracers and ${W}_{{\rm{H}}{\rm{I}}}$ alone cannot distinguish which variable (R or ${\tau }_{353}$) is the better tracer of the total dust and total gas column densities, and the correlation with the γ-ray intensity is crucial to reveal the true ${N}_{{\rm{H}}}$ distribution. Therefore, we prepared two types of ${N}_{{\rm{H}}}$ model maps based on R and ${\tau }_{353}$, in addition to the previously described ${N}_{{\rm{H}}}$ map based on ${W}_{{\rm{H}}{\rm{I}}}$, and tested them against the Fermi-LAT γ-ray data, in which the northern and southern regions were analyzed separately.

We converted the ${W}_{{\rm{H}}{\rm{I}}}$ map into an atomic hydrogen column density map by assuming the optically thin approximation (${N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}({\mathrm{cm}}^{-2})=1.82\times {10}^{18}\cdot {W}_{{\rm{H}}{\rm{I}}}({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})$). The obtained ${N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}$ map and the T_d map in our ROI are shown in panels (a) and (b), respectively, in Figures 3 and 4. Under the assumption that all neutral gas is atomic and ${\rm{H}}\,{\rm{I}}$ is optically thin, these ${N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}$ maps will be used as ${N}_{{\rm{H}}}$ template maps in the γ-ray data analysis (Section 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To construct the ${N}_{{\rm{H}}}$ model maps based on the Planck dust models for the northern region, we assumed that ${N}_{{\rm{H}}}$ and ${D}_{\mathrm{em}}$ (R or ${\tau }_{353}$) are proportional and that ${\rm{H}}\,{\rm{I}}$ is optically thin, at least for the areas with high dust temperature (${T}_{{\rm{d}}}\geqslant 21.0\,{\rm{K}}$). This assumption is based on the fact that high-T_d areas show a ${\tau }_{353}/{W}_{{\rm{H}}{\rm{I}}}$ ratio smaller than that in low-T_d areas and are expected to have relatively high gas temperatures, and are therefore likely to be dominated by optically thin ${\rm{H}}\,{\rm{I}}$. The uniformity of the ${N}_{{\rm{H}}}$/${D}_{\mathrm{em}}$ ratio over the whole ROI will be tested in the γ-ray data analysis (Sections 4.1.2 and 4.2.2). First, we made a least-squares fit to the dust–${W}_{{\rm{H}}{\rm{I}}}$ relationship for ${T}_{{\rm{d}}}\geqslant 21.0\,{\rm{K}}$ in Figure 1 using a linear function with an intercept fixed at 0 because the zero level of the dust emission is already subtracted in the Planck dust models (Planck Collaboration XI 2014). We obtained coefficients of $21.1\times {10}^{8}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\left({\rm{W}}{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1}\right)}^{-1}$ and $87.2\times {10}^{6}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for R and ${\tau }_{353}$, respectively. In these fits, ${W}_{{\rm{H}}{\rm{I}}}$ is treated as an independent variable because the signal-to-noise ratio is much greater than that of ${D}_{\mathrm{em}}$. Then, we converted R (or ${\tau }_{353}$) into ${N}_{{\rm{H}}}$ maps using the obtained coefficients and multiplied by $1.82\times {10}^{18}\,{\mathrm{cm}}^{-2}\,{\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}$. This gives a conversion factor for R of $38.4\times {10}^{26}\,{\mathrm{cm}}^{-2}\,{\left({\rm{W}}{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1}\right)}^{-1}$ and for ${\tau }_{353}$ of $159\times {10}^{24}\,{\mathrm{cm}}^{-2}$. The obtained ${N}_{{\rm{H}}}$ template maps are shown in Figure 3 (panels (c) and (d), respectively). By comparing these ${N}_{{\rm{H}}}$ model maps, we can see that the map based on ${W}_{{\rm{H}}{\rm{I}}}$ shows the weakest contrast. The ${\tau }_{353}$-based map predicts the strongest contrast and approximately a factor of 2 higher ${N}_{{\rm{H}}}$ in the dense clouds compared to the ${W}_{{\rm{H}}{\rm{I}}}$-based map, while the R-based map shows a moderate contrast. Thus, they will give different contrasts in the predicted γ-ray map and can be tested by the fit quality to the γ-ray data (Section 4).

For the southern region, we also assumed that ${W}_{{\rm{H}}{\rm{I}}}$ and ${D}_{\mathrm{em}}$ are proportional and followed the same procedure as that for the northern region. The obtained coefficients for the dust–${W}_{{\rm{H}}{\rm{I}}}$ relationship are $17.6\times {10}^{8}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\left({\rm{W}}{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1}\right)}^{-1}$ and $66.9\times {10}^{6}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and the conversion factors to construct the ${N}_{{\rm{H}}}$ model are $32.0\times {10}^{26}\,{\mathrm{cm}}^{-2}\,{\left({\rm{W}}{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1}\right)}^{-1}$ and $122\times {10}^{24}\,{\mathrm{cm}}^{-2}$ for R and ${\tau }_{353}$, respectively. The obtained ${N}_{{\rm{H}}}$ template maps are shown in Figure 4; the contrast of the gas density is weakest for the map based on ${W}_{{\rm{H}}{\rm{I}}}$, while that based on ${\tau }_{353}$ is the strongest.

Another possible source of diffuse γ-ray emissions is CR interactions with ionized gas. To estimate this contribution, we referred to Casandjian (2015) and used the free–free intensity map at a frequency of 22.7 GHz extracted from 9 yr of WMAP observations, namely, wmap_K_mem_freefree_9yr_v5.fits (Bennett et al. 2013), as a template for the γ-ray emission correlated with the ionized hydrogen (${\rm{H}}\,{\rm{II}}$). We used the scaling factor adopted by Casandjian (2015), $1.3\times {10}^{20}\,{\mathrm{cm}}^{-2}\,{\mathrm{mK}}^{-1}$, and constructed an ${\rm{H}}\,{\rm{II}}$ column density (${N}_{{\rm{H}}{\rm{II}}}$) model map. We confirmed that the average column density is only a few percent of the neutral gas column density estimated by ${W}_{{\rm{H}}{\rm{I}}}$ shown in Figures 3(a) and 4(a). Therefore, the contribution of the ionized gas to γ-ray emission is small, and the ${N}_{{\rm{H}}{\rm{II}}}$-related term will be fixed in the γ-ray data analysis (Section 3).

The LAT on board the Fermi Gamma-ray Space Telescope, launched in 2008 June, is a pair-tracking γ-ray telescope, detecting photons in the energy range from ∼20 MeV to more than 300 GeV. Details of the LAT instrument and the pre-launch performance expectations can be found in Atwood et al. (2009), and the on-orbit calibration is described in Abdo et al. (2009a). Past studies of the Galactic diffuse emissions by Fermi-LAT can be found in, e.g., Ackermann et al. (2012c) and Casandjian (2015).

Routine science operations with the LAT started on 2008 August 4. We used a data-taking interval from 2008 August 4 to 2016 August 1 (i.e., 8 yr) to study diffuse γ-rays in our ROI ($200^\circ \leqslant l\leqslant 260^\circ$ and $210^\circ \leqslant l\leqslant 270^\circ$ for the northern and southern regions, respectively, with $22^\circ \leqslant | b| \leqslant 60^\circ$ for both). During most of this time, the LAT was operated in sky survey mode, obtaining complete sky coverage every two orbits (which corresponds to $\sim 3$ hr) and a relatively uniform exposure. Because the diffuse γ-ray emission from the ISM is greatly extended and the intensity is rather weak at high latitudes, having a clean (low background) sample of photons is crucial. Therefore, we used the latest release of the Pass 8 data (P8R3) recently made public by the LAT Collaboration. The previous data set (P8R2) is known to suffer from a residual background with a peak along the ecliptic plane that is inside our ROI in the northern region. We used the standard LAT analysis software, Fermi Science Tools¹⁰ version v10r00p05, to select events above 100 MeV¹¹ because good angular resolution is essential for examining the correlation between γ-rays and gas distribution, and we selected events of the SOURCE class. The instrument response function (IRF) that matches our data set and event selection, P8R3_SOURCE_V2, was used in the following analysis. We also analyzed data using the cleanest ULTRACLEANVETO class and corresponding IRF and found that a decrease of the isotropic component (see Section 3.2), which includes the residual background, was marginal (within $\leqslant 10 \%$). Therefore, we keep using the SOURCE class to maximize the photon statistics. In addition, we required that the reconstructed zenith angles of the arrival direction of the photons be less than $100^\circ$ and $90^\circ$ for energies above and below 200 MeV, respectively, to reduce contamination by photons from Earth's atmosphere. To accommodate the rather poor angular resolution at low energy, below 200 MeV, we used events and the responses of point-spread function (PSF) event types 2 and 3; meanwhile, above 200 MeV, we did not apply selections based on PSF event types and used events and the responses of Front + Back. As described in Section 3.3, we carried out a bin-by-bin likelihood fitting, and data below and above 200 MeV are analyzed separately. We stopped at 25.6 GeV because of poor photon statistics above that energy. We used the gtselect command to apply the selections described above.

In addition, we excluded periods of time during which the LAT detected bright γ-ray bursts or solar flares (by using the gtmktime command); the integrated period of time excluded in this procedure was negligible (less than 1%) compared to the total observation time. The data count maps in the northern and southern regions of our ROI are given in Figure 5.

Zoom In Zoom Out Reset image size

We modeled the γ-ray emission above 100 MeV observed by Fermi-LAT as a linear combination of the ${N}_{{\rm{H}}}$ model map constructed from the HI4PI survey data (${W}_{{\rm{H}}{\rm{I}}}$) or the Planck dust map (R or ${\tau }_{353}$; see Section 2.2), the template map for ${N}_{{\rm{H}}{\rm{II}}}$, the IC emission, the isotropic component, the emission from the Sun and Moon (e.g., Abo et al. 2011; Ackermann et al. 2016), and the γ-ray point sources. The use of the gas column density map as a template is based on the assumption that γ-rays are generated via interactions between CRs and the ISM gas and that the CR intensity does not vary significantly over the scale of the interstellar complexes in this study. This assumption is simple but very plausible, particularly in high-latitude regions (where most of the gas is expected to be local), such as the one studied here. We started with a single ${N}_{{\rm{H}}}$ map (in Sections 4.1.1 and 4.2.1) and then tested multiple ${N}_{{\rm{H}}}$ maps based on the Planck dust model maps (in Sections 4.1.2 and 4.2.2). To model the γ-rays produced via IC scattering, we used GALPROP¹² (e.g., Strong & Moskalenko 1998; Strong et al. 2007), a numerical code that solves the CR transport equation within the Galaxy and predicts the γ-ray emission produced via the interactions of CRs with interstellar gas and low-energy photons. The IC emission is calculated from the distribution of the propagated electrons and the interstellar radiation field developed by Porter et al. (2008). Here we adopted the IC model map produced in the GALPROP run 54_Yusifov_z4kpc_R30kpc_Ts150K_EBV2mag,¹³ which was developed by Ackermann et al. (2012c) and has been used in other LAT Collaboration publications¹⁴ such as Planck Collaboration Int. XXVIII (2015) and Remy et al. (2017). To model individual γ-ray sources, we referred to the preliminary Fermi-LAT 8 yr point-source list (FL8Y),¹⁵ which is based on the first 8 yr of the science phase of the mission and includes more than 5000 sources detected at a significance of ≥4σ. We note that the list is built using the same lowest-energy threshold (0.1 GeV) as the present study. For our analysis, we considered 151 and 142 FL8Y sources (detected at a significance of ≥4σ) in the northern and southern regions, respectively, and 27 bright sources just outside the ROIs (≥50σ within $10^\circ$ and ≥20σ within $1^\circ$) to take into account possible contaminations from their emission caused by overlap because of the breadth of the PSF. We modeled a (quasi-)isotropic component to represent the contributions to our ROI due to extragalactic diffuse emission and the residual instrumental background from misclassified CR interactions in the LAT detector by the uniform emission in our ROI. A template of the ionized gas was constructed based on the WMAP free–free intensity map as described in Section 2.2. We included the templates for the γ-ray emission from the Sun and Moon that were used in the initial sky models for the fourth LAT source catalog (4FGL; Abdollahi et al. 2019), which is based on the first 8 yr of LAT observations. Specifically, we used template_SunDiskMoonv3r2_8years_zmax105.fits and template_SUNICv2_8years_zmax105.fits.¹⁶

Then, the γ-ray intensity ${I}_{\gamma }(l,b,E)\,(\mathrm{ph}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}\,{\mathrm{MeV}}^{-1})$ can be modeled as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\gamma }(l,b,E) & = & {q}_{\gamma }(E)\left[{c}_{1}(E)\cdot {N}_{{\rm{H}}}(l,b)+{N}_{{\rm{H}}{\rm{II}}}(l,b)\right]\\ & & +{c}_{2}(E)\cdot {I}_{\mathrm{IC}}(l,b,E)+{I}_{\mathrm{iso}}(E)\\ & & +{c}_{3}(E)\cdot {I}_{\mathrm{SM}}(l,b,E)+\sum _{j}{P}_{j}(l,b,E),\end{array}\end{eqnarray} \tag{ 1 }$$

where ${N}_{{\rm{H}}}(l,b)$ is the total neutral gas column density model (based on ${W}_{{\rm{H}}{\rm{I}}}$, R, or ${\tau }_{353}$) in ${\mathrm{cm}}^{-2}$; ${N}_{{\rm{H}}{\rm{II}}}(l,b)$ is the gas column density model from ionized gas in ${\mathrm{cm}}^{-2}$; ${q}_{\gamma }(E)$ ($\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,{\mathrm{MeV}}^{-1}$) is the model of the differential γ-ray yield or γ-ray emissivity per H atom; ${I}_{\mathrm{IC}}(l,b,E)$, ${I}_{\mathrm{iso}}(E)$, and ${I}_{\mathrm{SM}}(l,b,E)$ are the IC model, (quasi-)isotropic background intensities, and model of the emission from the Sun and Moon (each in units of $\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}\,{\mathrm{MeV}}^{-1}$), respectively; and ${P}_{j}(l,b,E)$ represents the point-source contributions. We used the γ-ray emissivity model ${q}_{\gamma }(E)$ of the local interstellar spectrum (LIS) of the CRs (protons and electrons) adopted in Abdo et al. (2009b), specifically the model curve for the nuclear enhancement factor ${\epsilon }_{{\rm{M}}}$ (a scale factor to take into account the effect of heavy nuclei in both CRs and the target matter) of 1.84 (Mori 2009). In other words, ${q}_{\gamma }(E)$ is decomposed as $1.84\cdot {q}_{\gamma }(\mathrm{pp})+{q}_{\gamma }(\mathrm{brems})$, where ${q}_{\gamma }(\mathrm{pp})$ and ${q}_{\gamma }(\mathrm{brems})$ are the γ-ray emissivity models due to the proton–proton interaction and electron bremsstrahlung, respectively. To accommodate the uncertainties in the LIS and ${\epsilon }_{{\rm{M}}}$, we included a scale factor (${c}_{1}(E)$ in Equation (1)) as a free parameter. This parameter is 1.0 if the measured γ-ray emissivity agrees with the LIS and ${\epsilon }_{{\rm{M}}}$ we adopted. Because the estimated ionized gas column density is small (at the maximum $\sim 1\times {10}^{20}\,{\mathrm{cm}}^{-2}$), we fixed the scale factor at 1.0 for the ${N}_{{\rm{H}}{\rm{II}}}(l,b)$ template to obtain a stable fitting. The IC emission model is also uncertain, and we included another scale factor (${c}_{2}(E)$ in Equation (1)) as a free parameter. The extragalactic diffuse emission and the residual background are modeled by the isotropic term, I_iso, as a uniform template. The intensity of this component could exhibit large-scale fluctuations due to the uncertainty of the LAT response and/or possible nonuniformity of the background. Therefore, we include the isotropic term in each energy band rather than adopt the template provided by the LAT Collaboration¹⁷ or determine it by the fit outside of our ROI. We note that the isotropic term is a dominant component (see Section 4); therefore, we should be able to constrain its contribution using the data in our ROI. The scale factor for ${I}_{\mathrm{SM}}(l,b,E)$, ${c}_{3}(E)$, was taken to be a free parameter in the northern region where the ecliptic circle goes through our ROI, even though it was fixed to 1.0 in the southern region. The point-source contributions were also taken to be free parameters as a function of energy. The positions of the sources were fixed to the values in FL8Y. To model the contamination from outside the ROI, we used the ${N}_{{\rm{H}}}$, I_IC, and I_SM maps including peripheral regions.

Among the model components described above, the γ-ray emission model from Sun and Moon was adopted from the work for the 4FGL while it was in the development phase. The model will be outdated when the work to construct the 4FGL is complete. However, we confirmed that the choice of the Sun and Moon template does not affect the fit for the neutral gas component as described in Sections 4.1.1 and 4.1.2. Because the IC model has uncertainty, we also examined the effect as described in Sections 4.1.2 and 4.2.2.

We divided the γ-ray data from 0.1 to 25.6 GeV into eight energy ranges using logarithmically equally spaced energy bands, and each band was divided into four sub-bins. Then, we fit Equation (1) to data in each energy band in $0\buildrel{\circ}\over{.} 5\times 0\buildrel{\circ}\over{.} 5$ bins below 400 MeV (where the PSF (68% contaminant radius) is $\geqslant 2^\circ$) and $0\buildrel{\circ}\over{.} 25\times 0\buildrel{\circ}\over{.} 25$ bins above 400 MeV using the binned maximum likelihood method with Poisson statistics implemented in the Science Tools. In each narrow energy band, c₁, c₂, and c₃ were modeled as energy-independent normalization factors. Because ${I}_{\mathrm{iso}}(E)$ is the most dominant component over the entire energy range, it was modeled via a power-law function with both photon index and normalization as free parameters in each energy band. ${P}_{j}(l,b,E)$ were modeled via separate power-law functions in each energy band with only the normalization allowed to vary; the photon index was fixed at 2.2 as a representative value of the high-latitude LAT sources (FSRQ and BL Lac; Ackermann et al. 2015).

When modeling the point sources, we iteratively included them in several groups at a time in order of decreasing significance. First, we included and fitted the brightest sources detected in FL8Y at more than 35σ; then, we added and fit a second group detected at 20σ–35σ, freezing the already-included source parameters. In this way, we worked down to the sources detected at more than 4σ in FL8Y. For each step, the parameters of the diffuse emission model (c₁, c₂, c₃, and I_iso) were always left free to vary. After we reached the least significant sources (more than 4σ), we went back to the brightest ones (more than 35σ significance) and let them and diffuse emission models free to vary, while parameters of other sources are kept fixed to those already determined. We repeated the process until the increment of the log-likelihoods, $\mathrm{ln}L$,¹⁸ was less than 0.1 over one loop in each energy band. To model the contamination from outside the ROI, we took into account sources (with the model parameters fixed to those of FL8Y) detected above 50σ located at a distance $\leqslant 10^\circ$ from the region boundaries, or detected above 20σ and located within $\leqslant 1^\circ$ from the boundaries. Other sources outside of our ROIs are not considered.

4.1.1. Initial Modeling with a Single Gas Map

First, we analyzed the northern region and started our data analysis using individual ${N}_{{\rm{H}}}$ model maps based on the HI4PI survey data or the Planck dust model maps described in Section 2.2: we used ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})=1.82\times {10}^{18}\cdot {W}_{{\rm{H}}{\rm{I}}}({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})$, or ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})=38.4\,\times {10}^{26}\cdot R({\rm{W}}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1})$, or ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})\,=159\times {10}^{24}\cdot {\tau }_{353}$.

We fit the γ-ray data as described in Section 3.3; c₁, c₂, c₃, I_iso, and P_j are free parameters in each energy band. As described in Section 3.3, point sources are included iteratively. We remind the reader that the dust-based templates are equivalent to ${N}_{{\rm{H}}}=\left(\overline{{N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}/{D}_{\mathrm{em}}}\right)\cdot {D}_{\mathrm{em}}$, where $\overline{{N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}/{D}_{\mathrm{em}}}$ is the average of the ${N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}$/${D}_{\mathrm{em}}$ ratio calibrated in the high-T_d regions. Therefore, in this initial modeling, we implicitly assume that the ${N}_{{\rm{H}}}$/${D}_{\mathrm{em}}$ ratio is constant over the whole T_d range.

The obtained values of $\mathrm{ln}L$ summed over the entire energy range with the R-based and ${\tau }_{353}$-based ${N}_{{\rm{H}}}$ maps are 62.2 and −455.5, respectively, when compared to that of the ${W}_{{\rm{H}}{\rm{I}}}$-based ${N}_{{\rm{H}}}$ map. Therefore, the R-based ${N}_{{\rm{H}}}$ map is preferred by the γ-ray data and the ${\tau }_{353}$-based map is disfavored. This conclusion is unchanged even if we do not include weak sources (detected at less than 5σ), or if we increase the energy threshold (200 or 400 MeV) to reduce the coupling with sources. By looking at the residual maps summarized in Figure 6, we can identify extended positive residuals covering the range of l in 245°–260° and b in 30°–55° in the ${\tau }_{353}$-based map, which makes the $\mathrm{ln}L$ significantly worse. Although the difference of the $\mathrm{ln}L$ values is smaller between ${W}_{{\rm{H}}{\rm{I}}}$-based analysis and R-based ones, we can also identify positive residuals in the ${W}_{{\rm{H}}{\rm{I}}}$-based map at around $(l,b)=(259^\circ ,24^\circ )$, ($256^\circ$, $33\buildrel{\circ}\over{.} 5$), and ($236^\circ$, $37\buildrel{\circ}\over{.} 5$), which corresponds to a coherent low-T_d and high-${W}_{{\rm{H}}{\rm{I}}}$ (and ${D}_{\mathrm{em}}$) area. These are the areas where the ${W}_{{\rm{H}}{\rm{I}}}$-based ${N}_{{\rm{H}}}$ map predicts smaller gas column density. We note that a region of apparently flat residuals at $(l,b)\sim (236\buildrel{\circ}\over{.} 5,38\buildrel{\circ}\over{.} 5)$ (where a difference of predicted ${N}_{{\rm{H}}}$ is the largest) is visible in all three analyses, even though the predicted ${N}_{{\rm{H}}}$ is rather different there. This is likely due to the interplay with a weak (≤8σ) γ-ray source FL8Y J0946.2+0104 located at (l, b) ∼ ($235\buildrel{\circ}\over{.} 37$, $38\buildrel{\circ}\over{.} 56$). The averages of the normalizations (weighted inversely by the square of the error in each band) for the neutral gas component, c₁ in Equation (1), are 1.028 ± 0.011, 0.964 ± 0.012, and 0.615 ± 0.008 for the ${W}_{{\rm{H}}{\rm{I}}}$-based, R-based, and ${\tau }_{353}$-based maps, respectively.

Zoom In Zoom Out Reset image size

The emission from the Sun and the Moon is subdominant in the region, but it could still have a small effect on the results. To test this effect, we redid the analysis with the model of the Sun and Moon emission used for the LAT 3FGL catalog (Acero et al. 2015). The change in the scale factors of the ${N}_{{\rm{H}}}$ map is negligible ($\leqslant 1 \%$), and the results are thus not affected by this component.

We also note that the normalization of the IC term (${c}_{2}(E)$ in Equation (1)) is nearly 0 for the R-based analysis, likely (at least partially) because of the interplay with the isotropic term (the dominant component¹⁹ in our ROI). Because the data prefer the fit with ${c}_{2}(E)\sim 0$, we maintain the results but will examine the systematic uncertainty by employing alternative IC models and fixing ${c}_{2}(E)$ to 1.0 (Sections 4.1.2 and 5.2).

4.1.2. Dust-temperature-sorted Modeling

As we saw in Section 2.1 (Figure 1), the correlation between ${W}_{{\rm{H}}{\rm{I}}}$ and ${D}_{\mathrm{em}}$ depends on T_d, and the temperature dependence is significantly different in the cases of R or ${\tau }_{353}$. Although the model with ${N}_{{\rm{H}}}$ proportional to R is preferred to the one proportional to ${\tau }_{353}$ (and the one proportional to ${W}_{{\rm{H}}{\rm{I}}}$) in terms of $\mathrm{ln}L$ by the γ-ray data analysis as described in Section 4.1, the true ${N}_{{\rm{H}}}$ distribution could be appreciably different from either of them; see Section 2.1 for the prerequisites for ${D}_{\mathrm{em}}$ to be proportional to ${N}_{{\rm{H}}}$ and also Mizuno et al. (2016), who reported the apparent T_d dependence of ${N}_{{\rm{H}}}$/R and ${N}_{{\rm{H}}}/{\tau }_{353}$ in the MBM 53, MBM 54, and MBM 55 clouds and the Pegasus loop.

Therefore, we performed an analysis with the T_d-sorted ${N}_{{\rm{H}}}$ template maps. We split the ${N}_{{\rm{H}}}$ template map (constructed from R or ${\tau }_{353}$) into four templates based on T_d, for ${T}_{{\rm{d}}}\leqslant 19\,{\rm{K}}$, ${T}_{{\rm{d}}}=19\mbox{--}20\,{\rm{K}}$, ${T}_{{\rm{d}}}=20\mbox{--}21\,{\rm{K}}$, and ${T}_{{\rm{d}}}\geqslant 21\,{\rm{K}}$,²⁰ and fit the γ-ray data with Equation (1) but using ${\sum }_{i}{c}_{1,i}(E)\cdot {N}_{{\rm{H}},i}(l,b)$ instead of ${c}_{1}(E)\cdot {N}_{{\rm{H}}}(l,b)$, where ${c}_{1,i}(E)$ and ${N}_{{\rm{H}},i}(l,b)$ represent the scale factor and the template gas map for each of the four templates. Under the assumption of a uniform CR intensity, ${c}_{1,i}(E)$ should not show T_d dependence if the ${N}_{{\rm{H}}}$/${D}_{\mathrm{em}}$ ratio is constant. To keep the fit stable while accommodating more free parameters, in this analysis (and that in Section 4.2.2) we omitted the lowest-energy band (where the angular resolution is poor) and the highest-energy band (where the photon statistics is low) so that the energy bands are restricted to the range 0.2–12.8 GeV. The improvement in the fits, the likelihood test statistics $\mathrm{TS}\equiv 2{\rm{\Delta }}\mathrm{ln}L$, were 41.4 and 496.6 with 18 more degrees of freedom (giving statistical significances of 3.2σ and 20σ) for the R-based and ${\tau }_{353}$-based fit, respectively. Therefore, the fit improvement was significant in both cases;²¹ however, the R-based analysis was still preferred.

Tables 1 and 2 show how the scaling factor ${c}_{1,i}(E)$ depends on T_d, and the averages over 0.2–12.8 GeV are summarized in Figure 7. The ${\tau }_{353}$-based analysis reveals a positive correlation (a lower scaling factor in ${T}_{{\rm{d}}}\leqslant 20\,{\rm{K}}$), implying an overestimation of ${N}_{{\rm{H}}}/{\tau }_{353}$ in the low-T_d area. Even though a negative correlation might be seen in the R-based analysis, the dependence on T_d is less clear and the fit improvement over the analysis with the single ${N}_{{\rm{H}}}$ template map is moderate. Therefore, careful examination of the systematic uncertainties is required.

Zoom In Zoom Out Reset image size

Table 1.  Results with the R-based ${N}_{{\rm{H}}}$ Maps Sorted by T_d in the Northern Region

Energy	${c}_{\mathrm{1,1}}$	${c}_{\mathrm{1,2}}$	${c}_{\mathrm{1,3}}$	${c}_{\mathrm{1,4}}$	${c}_{2}$	${c}_{3}$	${I}_{\mathrm{iso}}$	${I}_{\mathrm{iso}}$
(GeV)	(${T}_{{\rm{d}}}\leqslant 19\,{\rm{K}}$)	(19–20 K)	(20–21 K)	(${T}_{{\rm{d}}}\geqslant 21\,{\rm{K}}$)			(Intensitya)	(Index)
0.2–0.4	0.89 ± 0.04	0.77 ± 0.04	0.96 ± 0.03	0.78 ± 0.06	≤0.20	1.14 ± 0.06	40.7 ± 0.4	2.28 ± 0.01
0.4–0.8	0.99 ± 0.04	0.89 ± 0.03	0.99 ± 0.02	0.87 ± 0.05	≤0.10	1.31 ± 0.10	16.1 ± 0.2	2.27 ± 0.02
0.8–1.6	1.00 ± 0.05	0.96 ± 0.05	1.05 ± 0.04	0.92 ± 0.07	≤0.31	1.38 ± 0.12	6.02 ± 0.17	2.44 ± 0.02
1.6–3.2	1.02 ± 0.06	0.97 ± 0.05	1.05 ± 0.04	0.91 ± 0.07	≤0.17	1.20 ± 0.17	2.23 ± 0.06	2.31 ± 0.04
3.2–6.4	0.99 ± 0.10	0.89 ± 0.08	1.03 ± 0.07	0.86 ± 0.13	≤0.38	1.12 ± 0.26	0.99 ± 0.04	2.25 ± 0.06
6.4–12.8	1.44 ± 0.21	1.17 ± 0.18	1.23 ± 0.16	1.04 ± 0.27	≤1.16	1.73 ± 0.43	0.34 ± 0.04	2.44 ± 0.10

Notes. The errors are the 1σ statistical uncertainties. Each of the four scale factors (${c}_{\mathrm{1,1}}$, ${c}_{\mathrm{1,2}}$, ${c}_{\mathrm{1,3}}$, and ${c}_{\mathrm{1,4}}$) gives the normalization for a specified range of T_d of the neutral gas template in each energy bin. c₂ and c₃ give the normalizations for the IC and the γ-ray emission from the Sun and Moon, respectively, in each energy bin. I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters.

^aThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

Table 2.  Results with the ${\tau }_{353}$-based ${N}_{{\rm{H}}}$ Maps Sorted by T_d in the Northern Region

Energy	${c}_{\mathrm{1,1}}$	${c}_{\mathrm{1,2}}$	${c}_{\mathrm{1,3}}$	${c}_{\mathrm{1,4}}$	${c}_{2}$	${c}_{3}$	${I}_{\mathrm{iso}}$	${I}_{\mathrm{iso}}$
(GeV)	(${T}_{{\rm{d}}}\leqslant 19\,{\rm{K}}$)	(19–20 K)	(20–21 K)	(${T}_{{\rm{d}}}\geqslant 21\,{\rm{K}}$)			(Intensitya)	(Index)
0.2–0.4	0.54 ± 0.03	0.49 ± 0.03	0.70 ± 0.02	0.76 ± 0.06	0.39 ± 0.14	0.88 ± 0.07	41.0 ± 0.7	2.25 ± 0.01
0.4–0.8	0.60 ± 0.02	0.58 ± 0.02	0.77 ± 0.02	0.83 ± 0.05	$\leqslant 0.24$	0.86 ± 0.10	17.4 ± 0.3	2.26 ± 0.01
0.8–1.6	0.61 ± 0.03	0.65 ± 0.02	0.82 ± 0.02	0.91 ± 0.06	0.29 ± 0.24	1.00 ± 0.13	6.25 ± 0.22	2.43 ± 0.02
1.6–3.2	0.60 ± 0.04	0.63 ± 0.03	0.80 ± 0.03	0.85 ± 0.08	$\leqslant 0.34$	0.87 ± 0.17	2.46 ± 0.07	2.33 ± 0.03
3.2–6.4	0.61 ± 0.06	0.62 ± 0.06	0.83 ± 0.06	0.85 ± 0.14	$\leqslant 0.48$	0.89 ± 0.25	1.04 ± 0.04	2.27 ± 0.05
6.4–12.8	0.91 ± 0.14	0.80 ± 0.12	0.96 ± 0.13	1.01 ± 0.28	$\leqslant 1.34$	1.53 ± 0.43	0.35 ± 0.04	2.46 ± 0.10

Notes. The errors are the 1σ statistical uncertainties. Each of the four scale factors (${c}_{\mathrm{1,1}}$, ${c}_{\mathrm{1,2}}$, ${c}_{\mathrm{1,3}}$, and ${c}_{\mathrm{1,4}}$) gives the normalization for a specified range of T_d of the neutral gas template in each energy bin. c₂ and c₃ give the normalizations for the IC and the γ-ray emission from the Sun and Moon, respectively, in each energy bin. I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters.

^aThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

A possible systematic effect that might affect the T_d dependence is the uncertainty of the IC model. Even if we adjust the IC spectrum by scaling it in each energy range, the spatial distribution of our IC model might not be accurate. This could affect the results shown in Figure 7 in two ways: by changing the T_d dependence of the scaling factor (i.e., the measured T_d dependence of ${N}_{{\rm{H}}}$/${D}_{\mathrm{em}}$), and by changing the values of the scaling factor (i.e., the measured γ-ray emissivity or CR intensity). If there is a small spatial variation of the (quasi-)isotropic component that is dominant in our ROI, and this uncertainty is not absorbed by the IC model, the scale factors of the ${N}_{{\rm{H}}}$ templates will be affected as well. In order to investigate such possibilities, we tested several alternative IC models. As described in Section 3.2, we used the IC model produced in the GALPROP run 54_Yusifov_z4kpc_R30kpc_Ts150K_EBV2mag as our baseline model. This configuration assumes a CR source distribution proportional to the pulsar distribution of Yusifov & Kücük (2004) and a CR halo size z_h of 4 kpc. As described in Ackermann et al. (2011, 2012c) and de Palma et al. (2012), the CR source distribution and the halo height typically most strongly affect the propagated CR spatial distribution (in the Galactocentric distance and the height from the Galactic plane) and therefore the IC spatial distribution (in l and b). Therefore, we tested two additional CR source distributions, the pulsar-based distribution by Lorimer et al. (2006) and the supernova remnant (SNR) based distribution by Case & Bhattacharya (1998), and one additional CR halo height, 10 kpc; these are also available in Ackermann et al. (2012c), where the diffusion coefficient was adjusted when changing the CR source distribution and halo height to match the direct CR measurements at Earth's surface. The obtained IC maps show the smallest gradient in the Galactic longitude direction with the SNR-based CR source distribution (the flattest distribution of our three choices) and the smallest gradient in the Galactic latitude direction with ${z}_{h}=10\,\mathrm{kpc}$ (the larger halo height of our two choices). These tests result in a very small variation in the scaling factors (comparable to the statistical errors) in the case of the R-based ${N}_{{\rm{H}}}$ template maps. This is due to the small normalization of the IC terms, which is nearly 0, as exemplified by Table 1, and likely underestimates the systematic uncertainty. Therefore, we also employed a fit using our baseline IC model in which the normalization of the IC term was fixed to 1.0. While the scale factors were found to be robust against fixing the IC normalization to 1.0 for the case of the ${\tau }_{353}$-based analysis, the scale factors for the R-based analysis decreased by ∼10%. The systematic uncertainty due to the choice of IC model, evaluated as described above, is shown by the shaded bands in Figure 7. We found that the T_d dependence seen for the ${\tau }_{353}$-based analysis is robust against systematic uncertainties examined here, while that for the R-based analysis is less significant (scale factors of each T_d range are roughly the same within errors) if not zero. We also tested the alternative Sun/Moon emission model (see Section 4.1.1) and confirmed that the change in the T_d dependence of the scale factors for the R-based analysis is $\leqslant 1 \%$ and negligible compared to the systematic uncertainty due to the IC modeling.

In Figure 7, we observed a clear positive T_d dependence for ${N}_{{\rm{H}}}$/${\tau }_{353}$. This trend implies an underestimation of the ${N}_{{\rm{H}}}/{\tau }_{353}$ ratio in low-temperature areas and cannot be interpreted as being due to the properties of CRs, because the physical parameters that determine ${\tau }_{353}$ (e.g., dust-to-gas ratio or the dust cross section at 353 GHz) do not affect the CR intensity. The only possible explanation in terms of the CR properties is the exclusion of charged particles in dense clouds with large magnetic fields expected in areas of low T_d. However, CRs have been confirmed to penetrate into dense cloud cores with ${W}_{\mathrm{CO}}\geqslant 10\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (e.g., Abdo et al. 2010; Ackermann et al. 2011, 2012a), which corresponds to densities much larger than those of the clouds studied here with conventional values of X_CO. With ${X}_{\mathrm{CO}}=\left(1\mbox{--}2\right)\times {10}^{20}\,{\mathrm{cm}}^{-2}\,{{\rm{K}}}^{-1}\,{\left(\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}$ inferred in nearby clouds (e.g., Grenier et al. 2015), we would have ${N}_{{\rm{H}}}=\left(2\mbox{--}4\right)\times {10}^{21}\,{\mathrm{cm}}^{2}$ for ${W}_{\mathrm{CO}}=10\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1};$ this is larger than for the densest clouds in the region studied (see Figures 1 and 2). Indeed, ${c}_{1,i}(E)$ in Table 2 shows that the trend in Figure 7 is seen in both the low- and high-energy bands. Therefore, the main cause of the T_d dependence found here is not the properties of the CRs. Instead, it suggests that the dust opacity increases as T_d decreases.

4.1.3. Final Modeling

As shown in Section 4.1.2, while the γ-ray data analysis reveals that the ${N}_{{\rm{H}}}/{\tau }_{353}$ ratio has a positive T_d dependence, the ${N}_{{\rm{H}}}$/R ratio is rather constant. Because an R-based analysis was still preferred in terms of $\mathrm{ln}L$ in T_d-sorted modeling, we used R to construct our best estimate of the ${N}_{{\rm{H}}}$ distribution. Within the systematic uncertainties, no dependence on T_d is indicated (see Figure 7); therefore, applying a simple correction such as a linear function of T_d to the conversion factor for R ($38.4\times {10}^{26}\,{\mathrm{cm}}^{-2}\,{\left({\rm{W}}{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1}\right)}^{-1}$) is not necessary. The most noticeable feature in Figure 7 is an apparent decrease in the scaling factor for ${T}_{{\rm{d}}}\geqslant 21\,{\rm{K}}$ where the conversion from R to ${N}_{{\rm{H}}}$ was calibrated (Section 4.1). To examine whether this marked decrease is robust, we performed a further analysis using two T_d-sorted ${N}_{{\rm{H}}}$ maps, one with ${T}_{{\rm{d}}}\leqslant 21\,{\rm{K}}$ and the other with ${T}_{{\rm{d}}}\geqslant 21\,{\rm{K}}$. The obtained scale factors for the neutral gas template agree within 1% for the low- and high-T_d maps when averaged over the entire energy band. We also tested other thresholds of T_d (19.0 and 20.0 K) and confirmed that the scaling factors agree within 1% for the low- and high-T_d maps. These results indicate that there is no strong statistical evidence for a T_d dependence, and therefore the single R-based ${N}_{{\rm{H}}}$ map is preferred. Therefore, we adopt it as our best estimate of the ${N}_{{\rm{H}}}$ distribution. Finally, we fit the γ-ray data using this map (in the same way as in Section 4.1.1) but with finer energy bins, in order to study the spectral shape in more detail; each band was divided into two sub-bands. We summarize the best-fit model parameters and the obtained spectral components in Table 3 and Figure 8, respectively.

Zoom In Zoom Out Reset image size

Table 3.  Results of the Fit with the R-based Single ${N}_{{\rm{H}}}$ Map in the Northern Region

Energy	c1	(${E}^{2}\cdot {c}_{1}\cdot {q}_{\gamma }$)a	c2	c3	Iiso	Iiso
(GeV)					(Normb)	(Index)
0.10–0.14	0.93 ± 0.10	1.11	$\leqslant 0.29$	0.88 ± 0.09	48.1 ± 1.4	1.93 ± 0.05
0.14–0.20	0.87 ± 0.06	1.41	0.27 ± 0.23	1.17 ± 0.09	32.7 ± 1.0	2.06 ± 0.05
0.20–0.28	0.84 ± 0.04	1.71	0.18 ± 0.24	1.15 ± 0.08	24.1 ± 0.5	2.26 ± 0.03
0.28–0.40	0.91 ± 0.03	2.12	0.28 ± 0.22	1.26 ± 0.10	14.8 ± 0.4	2.30 ± 0.04
0.40–0.57	0.93 ± 0.03	2.22	$\leqslant 0.18$	1.15 ± 0.12	9.93 ± 0.17	2.25 ± 0.04
0.57–0.80	1.00 ± 0.03	2.37	$\leqslant 0.26$	1.47 ± 0.15	5.97 ± 0.14	2.28 ± 0.05
0.80–1.13	1.04 ± 0.05	2.30	0.18 ± 0.49	1.34 ± 0.18	3.56 ± 0.21	2.28 ± 0.06
1.13–1.60	0.99 ± 0.05	2.00	$\leqslant 0.52$	1.45 ± 0.20	2.29 ± 0.17	2.52 ± 0.08
1.60–2.26	1.05 ± 0.04	1.84	$\leqslant 0.26$	1.06 ± 0.22	1.33 ± 0.05	2.30 ± 0.10
2.26–3.20	0.98 ± 0.05	1.42	$\leqslant 0.39$	1.43 ± 0.26	0.87 ± 0.04	2.23 ± 0.11
3.20–4.53	0.91 ± 0.08	1.05	$\leqslant 0.80$	0.99 ± 0.32	0.62 ± 0.04	2.48 ± 0.13
4.53–6.40	1.12 ± 0.11	0.99	$\leqslant 0.54$	1.34 ± 0.42	0.36 ± 0.02	2.20 ± 0.18
6.40–9.05	1.36 ± 0.16	0.92	$\leqslant 1.10$	1.34 ± 0.53	0.22 ± 0.02	2.52 ± 0.23
9.05–12.8	0.97 ± 0.24	0.50	1.54 ± 1.05	2.57 ± 0.73	0.11 ± 0.03	1.87 ± 0.38
12.8–18.1	1.17 ± 0.36	0.46	1.08 ± 1.38	1.22 ± 0.91	0.08 ± 0.02	2.01 ± 0.39
18.1–25.6	0.87 ± 0.60	0.27	$\leqslant 2.64$	1.40 ± 1.24	0.06 ± 0.03	2.36 ± 0.39

Notes. The errors are the 1σ statistical uncertainties. Parameters c₁, c₂, and c₃ give the normalizations for the neutral gas template, IC, and the γ-ray emission from the Sun and Moon, respectively, in each energy bin. I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters. For convenience, the best-fit value of the emissivity multiplied by E² is also tabulated.

^aThe emissivity (${c}_{1}\times {q}_{\gamma }$) multiplied by E², where $E=\sqrt{{E}_{\min }{E}_{\max }}$ in each energy bin in units of ${10}^{-24}\,{\mathrm{MeV}}^{2}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}\,{\mathrm{MeV}}^{-1}$. ^bThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

4.2.1. Initial Modeling with a Single Gas Map

We then proceeded to model the southern region, in which the same analysis procedures were used as for the northern region. Since the ecliptic plane does not run through the ROI, we fixed the scale factors for the Sun and Moon emission template to 1.0. The ${N}_{{\rm{H}}}$ model maps were prepared as described in Section 2.2: we used ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})\,=1.82\times {10}^{18}\cdot {W}_{{\rm{H}}{\rm{I}}}({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})$, or ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})=32.0\,\times \,{10}^{26}\cdot R({\rm{W}}\,{{\rm{m}}}^{-2}{\mathrm{sr}}^{-1})$, or ${N}_{{\rm{H}}}({\mathrm{cm}}^{-2})=122\,\times \,{10}^{24}\cdot {\tau }_{353}$. As a side effect of masking the Orion–Eridanus superbubble, the ${N}_{{\rm{H}}}$ template map has a similar spatial distribution (larger intensity toward the Galactic center) to that of the IC model, and they are degenerate with each other. As we progressed in the iterative method, we observed that the ${N}_{{\rm{H}}}$ template is given a progressively lower scale factor (c₁) while the IC component is given a higher one (c₂) in low-energy bands, although the contribution of point sources to the total γ-ray flux is almost unchanged. To mitigate this, we employed a (semi)global fitting as a preparatory stage. We first adopted wider energy bins allowing overlaps (70.7–282.8 MeV, 141.4–565.7 MeV, etc.), and we included point sources iteratively until the fit improvement is saturated as we did for the analysis of the northern region. Allowing the energy bins to overlap makes the fit results stabler by encouraging spectral interbin consistency. We then fixed the IC normalization to the best-fit value and repeated the analysis with the original energy bins (100–200 MeV, 200–400 MeV, etc.). This "two-stage" analysis is employed hereafter for the southern region.

The obtained log-likelihoods summed over individual energy ranges in 0.1–25.6 GeV with the R-based and ${\tau }_{353}$-based ${N}_{{\rm{H}}}$ maps are 19.4 and −74.6, respectively, when compared to that of the ${W}_{{\rm{H}}{\rm{I}}}$-based ${N}_{{\rm{H}}}$ map. Therefore, the R-based ${N}_{{\rm{H}}}$ map is preferred compared to the ${\tau }_{353}$-based map. Like the northern region, this conclusion is unchanged against the significance threshold of the point-source model or the lowest-energy threshold. The averages of the normalization for the neutral gas component, c₁ in Equation (1), are 0.946 ± 0.008, 0.946 ± 0.007, and 0.690 ± 0.005 for the ${W}_{{\rm{H}}{\rm{I}}}$-based, R-based, and ${\tau }_{353}$-based maps, respectively. The residual maps are summarized in Figure 9.

Zoom In Zoom Out Reset image size

4.2.2. Dust-emission-sorted Modeling

As we saw in Section 2.1 (Figure 2), the correlation between ${W}_{{\rm{H}}{\rm{I}}}$ and ${D}_{\mathrm{em}}$ shows a change of the slope in particular for ${\tau }_{353}$, while the T_d dependence is small; unlike the case of the northern region, the ${W}_{{\rm{H}}{\rm{I}}}$–${D}_{\mathrm{em}}$ relationship depends on ${D}_{\mathrm{em}}$ rather than T_d. Even though the R-based ${N}_{{\rm{H}}}$ map is preferred for the three models in terms of $\mathrm{ln}L$, the true ${N}_{{\rm{H}}}$ distribution could be appreciably different, and a possible nonlinear relationship between ${N}_{{\rm{H}}}$ and ${D}_{\mathrm{em}}$ should be examined. Therefore, we performed an analysis with R-sorted and ${\tau }_{353}$-sorted ${N}_{{\rm{H}}}$ maps (instead of T_d-sorted ${N}_{{\rm{H}}}$ maps applied for the northern region). We split the ${N}_{{\rm{H}}}$ template map (constructed from R) into three templates based on R, for $R\leqslant 12$, $R=12\mbox{--}20$, and $R\geqslant 20$ in units of ${10}^{-8}\,{\rm{W}}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1}$,²² and fit the γ-ray data with Equation (1), using ${\sum }_{i}{c}_{1,i}(E)\cdot {N}_{{\rm{H}},i}(l,b)$ instead of ${c}_{1}(E)\cdot {N}_{{\rm{H}}}(l,b)$, where ${c}_{1,i}(E)$ and ${N}_{{\rm{H}},i}(l,b)$ represent the scale factor and template gas map for each of the three templates. For the ${\tau }_{353}$-based ${N}_{{\rm{H}}}$ template, we split it into three as ${\tau }_{353}\leqslant 4$, ${\tau }_{353}=4\mbox{--}6$, and ${\tau }_{353}\geqslant 6$ in units of 10⁻⁶.²³ With a plausible assumption of a uniform CR intensity, ${c}_{1,i}(E)$ is expected to trace the ${N}_{{\rm{H}}}$/${D}_{\mathrm{em}}$ ratio. We used data in the range of 0.2–12.8 GeV to avoid a possible unstable fit in the lowest (0.1–0.2 GeV) and highest (12.8–25.6 GeV) energy bands; then, we obtained a value of TS of 22.0 and 101.8 for R and ${\tau }_{353}$, respectively, with 18 more degrees of freedom. This indicates a 7.4σ improvement when the ${\tau }_{353}$ model for ${N}_{{\rm{H}}}$ is used, while the improvement of the fit is not significant for R. We also note that ${\tau }_{353}$-sorted modeling is still not favored compared to the analysis using the single ${N}_{{\rm{H}}}$ map based on R in terms of $\mathrm{ln}L$. The R and ${\tau }_{353}$ dependence of the scaling factors is summarized in Tables 4 and 5, and the averages over 0.2–12.8 GeV are summarized in Figure 10.

Zoom In Zoom Out Reset image size

Table 4.  Results with the R-based ${N}_{{\rm{H}}}$ Maps Sorted by R in the Southern Region

Energy	${c}_{\mathrm{1,1}}$	${c}_{\mathrm{1,2}}$	${c}_{\mathrm{1,3}}$	${c}_{2}$	${I}_{\mathrm{iso}}$	${I}_{\mathrm{iso}}$
(GeV)	(Ra ≤ 12)	(12–20)	($R\geqslant 20$)		(Normb)	(Index)
0.2–0.4	1.00 ± 0.04	0.88 ± 0.02	0.88 ± 0.03	2.72 ± 0.10	23.6 ± 0.3	2.35 ± 0.02
0.4–0.8	1.04 ± 0.03	1.01 ± 0.02	0.91 ± 0.03	2.32 ± 0.14	9.97 ± 0.16	2.33 ± 0.03
0.8–1.6	1.11 ± 0.04	1.09 ± 0.02	1.02 ± 0.04	2.53 ± 0.19	3.41 ± 0.08	2.46 ± 0.05
1.6–3.2	1.07 ± 0.06	1.01 ± 0.04	1.02 ± 0.05	2.56 ± 0.26	1.18 ± 0.05	2.25 ± 0.08
3.2–6.4	1.01 ± 0.11	1.09 ± 0.07	0.87 ± 0.09	2.37 ± 0.38	0.57 ± 0.03	2.14 ± 0.11
6.4–12.8	1.36 ± 0.22	1.15 ± 0.13	0.96 ± 0.17	3.07 ± 0.58	0.19 ± 0.02	2.83 ± 0.20

Notes. The errors are the 1σ statistical uncertainties. Each of the three scale factors (${c}_{\mathrm{1,1}}$, ${c}_{\mathrm{1,2}}$, and ${c}_{\mathrm{1,3}}$) gives the normalization for a specified range of R of the neutral gas template in each energy bin. Parameter c₂ is the IC template normalization for each energy bin obtained at the first stage of the fitting (see text). I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters.

^aR is given in units of ${10}^{-8}\,{\rm{W}}\,{{\rm{m}}}^{-1}\,{\mathrm{sr}}^{-1}$. ^bThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

Table 5.  Results with the ${\tau }_{353}$-based ${N}_{{\rm{H}}}$ Maps Sorted by ${\tau }_{353}$ in the Southern Region

Energy	${c}_{\mathrm{1,1}}$	${c}_{\mathrm{1,2}}$	${c}_{\mathrm{1,3}}$	${c}_{2}$	${I}_{\mathrm{iso}}$	${I}_{\mathrm{iso}}$
(GeV)	(τ353a ≤ 4)	(4–6)	(${\tau }_{353}\geqslant 6$)		(Normb)	(Index)
0.2–0.4	0.84 ± 0.04	0.67 ± 0.03	0.64 ± 0.02	3.01 ± 0.10	22.7 ± 0.3	2.33 ± 0.02
0.4–0.8	0.83 ± 0.03	0.82 ± 0.03	0.63 ± 0.02	2.85 ± 0.13	9.35 ± 0.15	2.31 ± 0.03
0.8–1.6	0.87 ± 0.03	0.85 ± 0.03	0.71 ± 0.02	3.15 ± 0.18	3.15 ± 0.08	2.47 ± 0.05
1.6–3.2	0.89 ± 0.05	0.78 ± 0.04	0.70 ± 0.03	3.06 ± 0.25	1.07 ± 0.04	2.27 ± 0.09
3.2–6.4	0.80 ± 0.09	0.91 ± 0.07	0.60 ± 0.05	2.83 ± 0.37	0.54 ± 0.03	2.16 ± 0.12
6.4–12.8	1.00 ± 0.18	0.87 ± 0.15	0.63 ± 0.10	3.54 ± 0.54	0.18 ± 0.02	2.33 ± 0.21

Notes. The errors are the 1σ statistical uncertainties. Each of the three scale factors (${c}_{\mathrm{1,1}}$, ${c}_{\mathrm{1,2}}$, and ${c}_{\mathrm{1,3}}$) gives the normalization for a specified range of ${\tau }_{353}$ of the neutral gas template in each energy bin. Parameter c₂ is the IC template normalization for each energy bin obtained at the first stage of the fitting (see text). I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters.

^a ${\tau }_{353}$ is given in units of 10⁻⁶. ^bThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

As for the northern region, we examined the systematic uncertainty due to the choice of the IC model and plotted the results in Figure 10; the outer polygonal area shows the full uncertainty evaluated by trying all six IC models, and the inner shaded area shows the variation where the model with the worst $\mathrm{ln}L$ was excluded. One may also argue that the obtained scale factors of the IC term and the isotropic emission intensity are greatly different between the northern and southern ROIs (Tables 1, 2, 4, and 5), and that the IC normalization is nearly 0 (which is not physical) with the R-based ${N}_{{\rm{H}}}$ template fit for the northern region. This indicates that our model does not completely describe the γ-ray data. For example, the IC spatial template may not agree with the true distribution, or there may be a small variation of the (quasi-)isotropic component. If the IC spatial template is not representing the underlying distribution of the IC emission observed by the LAT, it will most likely be ingested in the isotropic component contributions. Therefore, it is the total IC and isotropic that matters for these particular ROIs, and the sum of the IC term and isotropic emission is similar between two ROIs (Figures 8 and 11). Therefore, most of uncertainties of the IC term and the isotropic emission are mutually absorbed, and we believe that the effect on the neutral gas component is properly examined by employing several IC models as described above and in Section 4.1.2. We also note that while the differences of log-likelihoods are very small (19.4) between the R-based model and the ${W}_{{\rm{H}}{\rm{I}}}$-based one, the specific choice of the template does not affect the gas (and CR) properties very much; the normalizations of the neutral gas component are almost identical between the two gas models as described in Section 4.2.1.

Zoom In Zoom Out Reset image size

In Figure 10, we observe a clear negative ${\tau }_{353}$ dependence of the scale factor that is robust against the systematics owing to the choice of IC model. This trend implies an overestimation of the ${N}_{{\rm{H}}}/{\tau }_{353}$ ratio in the high-density area. As discussed in Section 4.1.2, this cannot be interpreted as being due to the properties of CRs. Instead, it is likely due to the dust properties, such as the change of the dust cross section by dust grain evolution (e.g., Roy et al. 2013). For the case of R-sorted analysis, there seems to be a negative dependence. However, the change of ${N}_{{\rm{H}}}$/R ratio (inversely proportional to the scale factor) is much smaller than that for ${\tau }_{353}$ whatever IC model is employed, and the fit improvement is not significant. Therefore, an R-based single ${N}_{{\rm{H}}}$ map is suggested to reproduce the ${N}_{{\rm{H}}}$ distribution inferred by γ-ray data.

4.2.3. Final Modeling

We employed the γ-ray data as a robust tracer of the total neutral gas distribution. We therefore can apply a correction to the ${N}_{{\rm{H}}}$ model template based on γ-ray data in principle. To prove this concept, and to see if this affects the choice of dust tracer (R or ${\tau }_{353}$), we started with the uncorrected ${N}_{{\rm{H}}}$ map that is proportional to ${\tau }_{353}$ (denoted as ${N}_{{\rm{H}},{\tau }_{353}}$) and modified the gas column density to take into account the observed ${\tau }_{353}$ dependence (Figure 10). We assumed that the ${N}_{{\rm{H}}}$ is proportional to ${\tau }_{353}$ up to a particular value (${\tau }_{\mathrm{bk}}$) and deviates from the proportionality linearly above that. Then, we can apply the correction to the ${N}_{{\rm{H}}}$ model using the empirical function below:

$$\begin{eqnarray}&&{N}_{{\rm{H}},\mathrm{mod}}=\left\{\begin{array}{l}{N}_{{\rm{H}},{\tau }_{353}}\,({\tau }_{353}\lt {\tau }_{\mathrm{bk}}),\\ {N}_{{\rm{H}},\mathrm{bk}}+(1-0.1\cdot C)\cdot ({N}_{{\rm{H}},{\tau }_{353}}-{N}_{{\rm{H}},\mathrm{bk}})\,({\tau }_{353}\geqslant {\tau }_{\mathrm{bk}}),\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where ${N}_{{\rm{H}},\mathrm{bk}}$ is the (uncorrected) gas column density (proportional to ${\tau }_{353}$) at ${\tau }_{\mathrm{bk}}$. C = 1 corresponds to a 10% decrease in ${N}_{{\rm{H}}}$ above ${N}_{{\rm{H}},\mathrm{bk}}$. We carried out a grid scan (${\tau }_{\mathrm{bk}}$ = 2, 3, 4, 5, 6 in units of 10⁻⁶ and C = 2, 3, 4, 5, 6) and found that ${\tau }_{\mathrm{bk}}=4$ and C = 4 give the best representation of the Fermi-LAT data. This configuration increases the scale factor of the neutral gas component by 20% and makes it agree with that from the R-based one within 15%.

The corrected ${N}_{{\rm{H}}}$ model based on ${\tau }_{353}$, however, still gives smaller $\mathrm{ln}L$ compared to the single ${N}_{{\rm{H}}}$ model based on R. Because we found that the R dependence of the scaling factors of the neutral gas component is small, if not zero, and the fit improvement is not significant over the single, R-based ${N}_{{\rm{H}}}$ template, we adopted this R-based ${N}_{{\rm{H}}}$ model as our best estimate of the ${N}_{{\rm{H}}}$ distribution. We therefore fit the γ-ray data using this map with finer energy bands to study the spectral shape in more detail as we did for the northern region. The best-fit model parameters and the obtained spectral components are summarized in Table 6 and Figure 11, respectively.

Table 6.  Results of the Fit with the R-based Single ${N}_{{\rm{H}}}$ Map in the Southern Region

Energy	c1	(${E}^{2}\cdot {c}_{1}\cdot {q}_{\gamma }$)a	c2	Iiso	Iiso
(GeV)				(Normb)	(Index)
0.10–0.14	0.91 ± 0.04	1.08	2.38 ± 0.10	26.7 ± 0.4	1.81 ± 0.08
0.14–0.20	0.86 ± 0.03	1.40	2.38 ± 0.10	20.1 ± 0.3	2.00 ± 0.07
0.20–0.28	0.86 ± 0.02	1.76	2.65 ± 0.10	15.6 ± 0.1	2.38 ± 0.05
0.28–0.40	0.91 ± 0.02	2.11	2.65 ± 0.10	9.63 ± 0.10	2.28 ± 0.06
0.40–0.57	0.97 ± 0.02	2.34	2.34 ± 0.14	6.39 ± 0.08	2.28 ± 0.07
0.57–0.80	1.00 ± 0.02	2.36	2.34 ± 0.14	3.96 ± 0.06	2.33 ± 0.09
0.80–1.13	1.08 ± 0.02	2.39	2.63 ± 0.19	2.16 ± 0.04	2.39 ± 0.12
0.13–1.60	1.05 ± 0.03	2.13	2.63 ± 0.19	1.35 ± 0.03	2.64 ± 0.15
1.60–2.26	1.00 ± 0.03	1.76	2.56 ± 0.26	0.78 ± 0.02	2.60 ± 0.20
2.26–3.20	1.08 ± 0.04	1.57	2.56 ± 0.26	0.48 ± 0.02	2.11 ± 0.25
3.20–4.53	0.97 ± 0.06	1.12	2.44 ± 0.37	0.36 ± 0.01	2.33 ± 0.27
4.53–6.40	1.16 ± 0.09	1.03	2.44 ± 0.37	0.22 ± 0.01	2.52 ± 0.34
6.40–9.05	1.08 ± 0.13	0.72	3.06 ± 0.57	0.14 ± 0.01	2.05 ± 0.44
9.05–12.8	1.18 ± 0.19	0.61	3.06 ± 0.57	0.08 ± 0.01	1.82 ± 0.57
12.8–18.1	0.81 ± 0.27	0.32	3.44 ± 0.90	0.06 ± 0.01	1.72 ± 0.65
18.1–25.6	0.63 ± 0.37	0.19	3.44 ± 0.90	0.04 ± 0.01	3.92 ± 0.78

Notes. The errors are the 1σ statistical uncertainties. Parameters c₁ and c₂ give the normalization for the neutral gas template and IC, respectively, in each energy bin (the best-fit values obtained at the first stage of the fitting are given for the latter). I_iso is modeled using a power law with the integrated intensity and the photon index as free parameters. For convenience, the best-fit value of the emissivity (${c}_{1}\times {q}_{\gamma }$) multiplied by E² is also tabulated.

^aThe emissivity (${c}_{1}\times {q}_{\gamma }$) multiplied by E², where $E=\sqrt{{E}_{\min }{E}_{\max }}$ in each energy bin in units of ${10}^{-24}\,{\mathrm{MeV}}^{2}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}\,{\mathrm{MeV}}^{-1}$. ^bThe integrated intensity (${10}^{-7}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}$) in each band.

Download table as:  ASCIITypeset image

In Section 4, we used the GeV γ-rays observed by Fermi-LAT as robust tracers of the ISM gas under the assumption of a uniform CR intensity and obtained the ${N}_{{\rm{H}}}$ distributions inferred by the γ-ray data. Trends of the scale factor for ${N}_{{\rm{H}}}$ templates (T_d dependence and ${\tau }_{353}$ dependence in the northern and southern regions, respectively) are commonly seen between low- and high-energy bands (see Tables 1, 2, 4, and 5); this supports the uniformity of the CR intensity in each ROI. We found that the ${N}_{{\rm{H}}}$ template based on the R data best matches the γ-ray observations in both northern and southern regions, and in the following discussion we assume that R is a good tracer of ${N}_{{\rm{H}}}$. The obtained relationships between ${W}_{{\rm{H}}{\rm{I}}}$ and ${N}_{{\rm{H}}}$ are shown in Figure 12, together with maps of the excess gas column density above ${N}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{thin}}$. We point out that the ${W}_{{\rm{H}}{\rm{I}}}$/${\tau }_{353}$ ratio (and ${N}_{{\rm{H}}}$/${\tau }_{353}$ ratio) strongly depends on T_d in the northern region, while in the southern region this dependence is weaker. The dust optical depth ${\tau }_{353}$ depends on T_d and the dust spectral index, β; the two properties are tightly connected (Planck Collaboration XI 2014). The anticorrelation between T_d and β is apparent in the northern region, while the southern region presents more dispersion. Differences in these dust properties suggest different grain evolution (e.g., Jones et al. 2013; Köhler et al. 2015). How this grain evolution is related to the observed ${\tau }_{353}$ dependence of the ${N}_{{\rm{H}}}$/${\tau }_{353}$ ratio in the southern region is not clear. In the following, we will focus on discussing implications of Figure 12.

Zoom In Zoom Out Reset image size

In the northern region (panel (a)), the ${N}_{{\rm{H}}}$/${W}_{{\rm{H}}{\rm{I}}}$ ratio in the area of low T_d (18–19 K) and high ${W}_{{\rm{H}}{\rm{I}}}$ ($\geqslant 300\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$) is greater than those in other areas. This area corresponds to the excess gas column density at around $(l,b)=(236^\circ ,38\buildrel{\circ}\over{.} 5)$ seen in panel (c). It also corresponds to the positive residual at around $(l,b)=(236^\circ ,37\buildrel{\circ}\over{.} 5)$ seen in the ${W}_{{\rm{H}}{\rm{I}}}$-based analysis (see also Figure 6(a) and Section 4.1.1). Because our ROI does not include strong CO emission (see Appendix C) and the data for the area of interest span a wide range of ${W}_{{\rm{H}}{\rm{I}}}$, these gas-related emissions are likely due to optically thick ${\rm{H}}\,{\rm{I}}$, and we will consider this scenario hereafter. The brightness temperature for the ${\rm{H}}\,{\rm{I}}$ emission at the velocity v is given by

$$\begin{eqnarray}&&{T}_{{\rm{b}}}(v)=[{T}_{{\rm{S}}}(v)-{T}_{\mathrm{bg}}]\cdot [1-\exp (-{\tau }_{{\rm{H}}{\rm{I}}}(v))],\end{eqnarray} \tag{ 3 }$$

where ${T}_{\mathrm{bg}}$ is the background continuum radiation temperature and ${T}_{{\rm{S}}}(v)$ and ${\tau }_{{\rm{H}}{\rm{I}}}(v)$ are, respectively, a harmonic mean of the spin temperature at velocity v on the line of sight and an integration of the optical depth at this velocity. Then, if we approximate the ${\rm{H}}\,{\rm{I}}$ emission spectrum by a single boxcar spectrum on the line of sight with a spectral width of ${\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}$, ${T}_{{\rm{S}}}(v)$ and ${\tau }_{{\rm{H}}{\rm{I}}}(v)$ can be expressed by single values independent of v, and thus the ${W}_{{\rm{H}}{\rm{I}}}$ can be correlated to ${N}_{{\rm{H}}}$ as a function of T_S such that (e.g., Fukui et al. 2015)

$$\begin{eqnarray}\begin{array}{rcl}{W}_{{\rm{H}}{\rm{I}}}({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}) & = & [{T}_{{\rm{S}}}({\rm{K}})-{T}_{\mathrm{bg}}({\rm{K}})]\cdot {\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}(\mathrm{km}\,{{\rm{s}}}^{-1})\\ & & \cdot \,[1-\exp (-{\tau }_{{\rm{H}}{\rm{I}}})]\end{array}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{\tau }_{{\rm{H}}{\rm{I}}}=\displaystyle \frac{{N}_{{{\rm{H}}}_{\mathrm{tot}}}({\mathrm{cm}}^{-2})}{1.82\times {10}^{18}}\cdot \displaystyle \frac{1}{{T}_{{\rm{S}}}({\rm{K}})}\cdot \displaystyle \frac{1}{{\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}(\mathrm{km}\,{{\rm{s}}}^{-1})},\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}$ is defined as ${W}_{{\rm{H}}{\rm{I}}}$/(peak ${\rm{H}}\,{\rm{I}}$ brightness temperature). In Figure 12(a), assuming that all of the gas is atomic, we overlaid the model curves for several choices of T_S with ${\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}=18\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (the median line width in the northern region). As inferred from Figure 12(a), while most of the region is compatible with being optically thin, the area with ${T}_{{\rm{d}}}=18\mbox{--}19\,{\rm{K}}$ gives, on average, ${T}_{{\rm{S}}}\sim 40\,{\rm{K}}$.

In the southern region (Figure 12(b)), while most of the data lie along a mildly curved line, those in regions with ${T}_{{\rm{d}}}\leqslant 18\,{\rm{K}}$ show a flat profile with ${W}_{{\rm{H}}{\rm{I}}}\sim 100\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in the plot. This corresponds to a spot seen in the dust data at (l, b) ∼ ($230^\circ$, $-28\buildrel{\circ}\over{.} 5$) (see also Figure 12(d)); it is also seen in the residual map for the ${W}_{{\rm{H}}{\rm{I}}}$-based analysis (Figure 9(a)). A plausible interpretation is that the spot consists of CO-dark ${{\rm{H}}}_{2}$ (e.g., Smith et al. 2014) because the flat profile means that the column density of ${\rm{H}}\,{\rm{I}}$ is nearly constant. Finally, we overlaid the model curves for several choices of T_S with ${\rm{\Delta }}{V}_{{\rm{H}}{\rm{I}}}=21\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (the median velocity of the southern region with several areas masked); the majority of data agree, on average, with the model curve for ${T}_{{\rm{S}}}=125\,{\rm{K}}$, supporting the choice of T_S by Abdo et al. (2009b). This value also agrees well with the average value found in the local ISM, ${T}_{{\rm{S}}}=140\,{\rm{K}}$ (Casandjian 2015). We also note that there is a scatter around the model curve for ${T}_{{\rm{S}}}=125\,{\rm{K}}$. Because ${W}_{{\rm{H}}{\rm{I}}}$ and the optically thin approximation give the lower limit of ${N}_{{\rm{H}}}$, the spread is likely because of the uncertainty of the ${D}_{\mathrm{em}}$/${N}_{{\rm{H}}}$ ratio and thus could introduce over- and underestimation of the ${N}_{{\rm{H}}}$. For example, negative values of the excess gas column density are seen around $(l,b)=(240^\circ ,-35^\circ )$ in Figure 12(d), and they correspond to the residual (i.e., underestimate of ${N}_{{\rm{H}}}$) in Figures 9(b) and (c). This is a drawback of employing ${D}_{\mathrm{em}}$ as a tracer of the total gas column density, and the small scatter around model curves in Figure 12 should not be taken at face value.

We point out that the ${N}_{{\rm{H}}}$ in Figure 12 is proportional to R. Although this is our best estimate of ${N}_{{\rm{H}}}$ based on the correlation between the γ-ray data and gas templates, we did not measure the ${N}_{{\rm{H}}}$ distribution on pixel scales. Accordingly, overinterpretation of Figure 12 (e.g., estimating the T_S and the excess gas column density on very small scales) should be avoided.

The average column density of the neutral gas ($\overline{{N}_{{\rm{H}}}}$) in the northern region is obtained as $\sim 3.4\times {10}^{20}\,{\mathrm{cm}}^{-2}$ based on either ${W}_{{\rm{H}}{\rm{I}}}$ (with optically thin approximation) or R (with the conversion factor determined in Section 2.2). On the other hand, that based on ${\tau }_{353}$ is $\sim 4.3\times {10}^{20}\,{\mathrm{cm}}^{-2}$, indicating that an ${N}_{{\rm{H}}}\propto {\tau }_{353}$ model (not favored by γ-ray analysis) overestimates the gas column density by $\sim 30 \%$. In the southern region, again, while the values of $\overline{{N}_{{\rm{H}}}}$ inferred by ${W}_{{\rm{H}}{\rm{I}}}$ and R are similar ($\sim 2.2\times {10}^{20}\,{\mathrm{cm}}^{-2}$), that based on ${\tau }_{353}$ is $\sim 15 \%$ higher. This supports our earlier statement that the use of the γ-ray data is crucial to accurately determine the ${N}_{{\rm{H}}}$ distribution.

Next, we discuss the ${\rm{H}}\,{\rm{I}}$ emissivity spectra obtained in this study, which are summarized in Figures 13(a) and (b). To investigate possible systematic uncertainties due to the choice of the ${N}_{{\rm{H}}}$ template and the IC model, in the northern region (where the normalization of the IC model is $\sim 0$ for all six models) we bracketed the spectrum with that obtained using the ${N}_{{\rm{H}}}$ model based on ${W}_{{\rm{H}}{\rm{I}}}$ (i.e., the pure optically thin ${\rm{H}}\,{\rm{I}}$ scenario) and with that obtained using the ${N}_{{\rm{H}}}$ model based on R but with the normalization of the IC model fixed to 1.0. For the southern region, we bracketed the spectrum with those obtained using all six IC models. We also took into account the LAT effective area uncertainty;²⁴ we assumed the uncertainty to be 10% below 200 MeV (where we used only events of PSF classes 2 and 3) and 5% above 200 MeV. The fractional uncertainties of the spectrum due to the modeling (${N}_{{\rm{H}}}$ and IC for the northern and southern regions, respectively) and that due to the effective area uncertainty are summed in quadrature and shown as shaded bands; they are 11%–13% below 200 MeV and $\lesssim 10 \%$ above 200 MeV. For comparison, we plotted the model curve for the LIS that we adopted with ${\epsilon }_{{\rm{M}}}$ of 1.84 in the same figure. In order to approximately indicate the uncertainty in the emissivity model (primarily due to the uncertainty in the elemental composition of the CRs and the cross sections other than for proton–proton (p–p) collisions), we also plotted the model curve for ${\epsilon }_{{\rm{M}}}=1.45$ (the lowest value referred to in Mori 2009), which gives 15%–20% lower emissivity. We also plotted the emissivity spectrum of the same ROI measured by Abdo et al. (2009b) using 6 months of Fermi-LAT data. Our results favor the model curve with ${\epsilon }_{{\rm{M}}}=1.84$ and agree well with those by Abdo et al. (2009b), but here we cover a wider energy range and investigate northern and southern regions separately. In other words, the analysis presented here shows for the first time that the ${\rm{H}}\,{\rm{I}}$ emissivities are consistent between the northern and southern regions at the 10% level, supporting the hypothesis that the CR intensity is uniform in the vicinity of the solar system. The integral emissivities above 100 MeV are $(1.58\pm 0.04)\,\times {10}^{-26}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,{\rm{H}} \mbox{-} {\mathrm{atom}}^{-1}$ and $(1.59\,\pm 0.02)\times {10}^{-26}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,{\rm{H}} \mbox{-} {\mathrm{atom}}^{-1}$ in the northern and southern regions, respectively, and those above 300 MeV are $(0.68\pm 0.01)\times {10}^{-26}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,{\rm{H}} \mbox{-} {\mathrm{atom}}^{-1}$ and $(0.69\pm 0.01)\times {10}^{-26}\,\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,{\rm{H}} \mbox{-} {\mathrm{atom}}^{-1}$ in the northern and southern regions, respectively, with an additional systematic error of $\sim 10 \%$ due to the modeling and the effective area uncertainties (see above). We also note that our emissivities agree (within $\leqslant 10 \%$) with the results of Shen et al. (2019), where the authors analyzed the same northern region employing a template-fitting method with the assumption of a uniform T_S for the atomic gas phase and discussed the local CR spectrum based on recent p–p interaction models and AMS–02 data. In other words, our analysis supports their findings by showing that ${T}_{{\rm{S}}}=125\,{\rm{K}}$ is compatible with most of the ${N}_{{\rm{H}}}$ distribution and that the CR spectrum is uniform.

Zoom In Zoom Out Reset image size

As in the case of discussing the ISM gas densities (Section 5.1), evaluating the gas model using the γ-ray data is crucial to accurately constrain the ${\rm{H}}\,{\rm{I}}$ emissivity and CR intensity. If we use the ${N}_{{\rm{H}}}\propto {\tau }_{353}$ models, the scale factors of the ${\rm{H}}\,{\rm{I}}$ emissivity (∝CR intensity) are 30%–40% lower (see Sections 4.1.1 and 4.2.1). Other sources of uncertainty on the CR intensity are the hadronic interaction cross section and the elemental composition of CRs as indicated by two curves (${\epsilon }_{{\rm{M}}}=1.84$ and 1.45) in Figure 13(a). If we adopt ${\epsilon }_{{\rm{M}}}=1.45$, we would need $\sim 25$% higher proton LIS flux, which might be incompatible with the proton flux directly measured at Earth. Given the uncertainty and the fact that the directly measured CRs do not necessarily represent the LIS, we do not deny such a possibility. See discussions in, e.g., Strong (2015), Orlando (2018), and Shen et al. (2019). In Figure 13(a), one may also recognize that the model overestimates the data below a few × 10² MeV while it predicts lower flux above 1 GeV. This might indicate a possible spectral break of the proton LIS. For example, Strong (2015) reported a break at a few GeV based on Casandjian (2015) that gives a similar spectrum to ours (see also Figure 13(b)). To reach a robust conclusion, constraining the electron LIS using radio synchrotron emission (e.g., Orlando 2018) and an accurate determination of the emissivity spectrum below a few ×10² MeV is crucial. Because the analysis suffers from coupling with the point sources through the IC model in low-energy bands (see Section 4.1.1), we defer such a study to future projects using gas-rich areas.

Finally, we compare our ${\rm{H}}\,{\rm{I}}$ emissivity spectrum (the average of the northern and southern regions) with several other Fermi-LAT studies of nearby clouds: the average spectrum found in the local ISM in $10^\circ \leqslant | b| \leqslant 70^\circ$ by Casandjian (2015), that toward the Orion molecular clouds by Ackermann et al. (2012b), that toward the Chamaeleon molecular clouds by Planck Collaboration Int. XXVIII (2015), that toward the South Taurus cloud by Remy et al. (2017), and that toward the MBM 53, MBM 54, and MBM 55 clouds and the Pegasus loop by Mizuno et al. (2016), as summarized in Figure 13(b). Although the spectral shape does not change significantly over the samples examined here, the peak-to-peak variation of the normalization is by a factor of ∼2 even in nearby clouds. Given the diffusive nature of the Galactic CRs and the lack of significant change of the spectral shape, the variation is mostly attributable to uncertainties of the gas column density, particularly due to assumptions of the value of T_S. For example, as discussed by Planck Collaboration Int. XXVIII (2015), the emissivity toward the Chamaeleon clouds agrees with that found by Casandjian (2015) within $\sim 20 \%$ if we assume ${T}_{{\rm{S}}}=140\,{\rm{K}}$ for ${\rm{H}}\,{\rm{I}}$ clouds, although the γ-ray fit favors higher T_S. A small gradient at the 20% level could be possible and of interest to understand the CR generation and propagation in the vicinity of the solar system, and a systematic study of nearby clouds is therefore important. In such future studies, one should overcome the uncertainty on the ${N}_{{\rm{H}}}$ distribution by using the ${\rm{H}}\,{\rm{I}}$, CO, and dust data together with GeV γ-rays as a robust tracer of the ISM gas.