FERMI-LAT OBSERVATIONS OF HIGH- AND INTERMEDIATE-VELOCITY CLOUDS: TRACING COSMIC RAYS IN THE HALO OF THE MILKY WAY

The LAT is a pair-conversion γ-ray telescope for the energy range 20 MeV to $\gt 300$ GeV. The instrument is described in Atwood et al. (2009), and its on-orbit performance discussed in Ackermann et al. (2012c). We use LAT observations accumulated from the beginning of the Science phase of the mission, MET⁹ 239500801 (2008 August 4), up to MET 431481603 (2014 September 4), for a total of 73 months. We use the P7REP data set (Bregeon et al. 2013) and select "Clean" class events in order to limit the contamination from backgrounds due to residual CRs misclassified as γ rays in the study of large-scale diffuse features.

We apply further restrictions to the data set used:

• 1.  
  We select only time intervals when the LAT was operated in normal Science mode (e.g., excluding calibration data) and when the data quality was considered good (e.g., excluding solar flares which caused pile-up phenomena in the detector).
• 2.  
  We retain only candidate γ rays detected at inclinations from the Earth zenith $\lt 95^\circ$, in order to limit the contamination from the bright γ-ray emission produced by CR interactions in the Earth's atmosphere that forms extended features in the sky.
• 3.  
  We consider only candidate γ rays with measured energies between 300 MeV and 10 GeV. Below 300 MeV the LAT angular resolution rapidly deteriorates (Ackermann et al. 2012c), which makes separating different components of diffuse γ-ray emission based on their morphologies more difficult. Additionally, below ∼300 MeV one would need to model the residual emission from CR interactions in the upper atmosphere. This can introduce a faint diffuse glow at moderately high declinations, which would need to be modeled, and would further complicate the study of faint extended features especially near the North Celestial pole (NCP) where our targets are located. Above 10 GeV the photon counts are sparse and the data are not very useful to constrain CR interactions in clouds owing to their spectrum being softer than inverse-Compton (IC) emission from Galactic CR electrons and isotropic backgrounds (e.g., Ackermann et al. 2012a). Limiting the energy range also helps to attenuate the systematic uncertainties due to the spectral models assumed for the various components of the γ-ray sky.

For the γ-ray analysis described in this section and elsewhere in the article we used the official LAT Science Tools, version v09r34p01.¹⁰

If the CR densities are uniform within a cloud, the γ-ray intensities are simply proportional to the target gas column densities. Thus, in order to model γ-ray emission from IVCs and HVCs and foreground emission from local gas we rely on a set of ISM tracers described here.

We employ as tracer of atomic gas the Leiden-Argentine-Bonn (LAB) all-sky survey of the 21 cm ${\rm{H}}\;{\scriptsize{\rm I}}$ emission line by Kalberla et al. (2005). The survey provides the ${\rm{H}}\;{\scriptsize{\rm I}}$ line brightness temperatures over a grid of 05 (with an effective angular resolution of ∼06), and covers velocities with respect to the local standard of rest from −450 to +400 km s⁻¹, at a resolution of 1.3 km s⁻¹. Data were corrected for stray radiation reaching a residual noise contamination $\lt 0.09$ K for most directions in the sky. We use data from Kalberla et al. (2005) to derive $N({\rm{H}}\;{\scriptsize{\rm I}})$ column-density maps as described in Section 4.1.

Molecular hydrogen cannot be traced directly in its most common cold state. The most widely used surrogate tracer is the 2.6 mm line of ¹²CO. For this work we use a map of CO brightness temperature integrated over Doppler velocities, ${W}_{\mathrm{CO}}$, derived from the Center for Astrophysics composite CO survey (Dame et al. 2001) supplemented with higher-latitude observations of the Ursa Major region (de Vries et al. 1987). The CO data, denoised using the moment-masking technique (Dame 2011), were used to construct the ${W}_{\mathrm{CO}}$ map on a $0\buildrel{\circ}\over{.} 25$ grid. The ${W}_{\mathrm{CO}}$ intensities are assumed to be proportional to the molecular hydrogen column densities $N({{\rm{H}}}_{2})$.

Only ROI A contains significant CO emission (Figure 1). The CO emission detected is associated with the Ursa Major molecular clouds, located at a few 100 pc from the solar system (e.g., de Vries et al. 1987). Note that, according to the results from the component separation of the all-sky survey performed using the Planck-HFI instrument (Planck Collaboration XIII 2014), there is no significant CO emission in our ROIs besides what is included in this CO map.

Zoom In Zoom Out Reset image size

Various observations indicate that a linear combination of $N({\rm{H}}\;{\scriptsize{\rm I}})$ and ${W}_{\mathrm{CO}}$ maps does not properly account for the totality of the ISM neutral gas, but large fractions of the gas can be in a DNM phase, also known as dark gas (e.g., Magnani et al. 2003; Grenier et al. 2005; Abdo et al. 2010; Langer et al. 2010; Planck and Fermi Collaborations 2014). The DNM may be composed of molecular gas not associated with CO, revealed, e.g., by alternative tracers like CH (e.g., Magnani et al. 2003), that is likely to exist because the CO molecule is more prone to photodissociation than ${{\rm{H}}}_{2}$ in the external layers of molecular clouds (e.g., Wolfire et al. 2010; Smith et al. 2014). There can also be a contribution from atomic gas overlooked because of the approximations applied in deriving column densities from the ${\rm{H}}\;{\scriptsize{\rm I}}$ line intensities (Fukui et al. 2015). On large scales the DNM in the Milky Way is best traced by dust thermal emission or extinction and γ rays, which reveal correlated excesses on top of the components traced by ${\rm{H}}\;{\scriptsize{\rm I}}$ and CO that can be interpreted as a contribution from undetected, but otherwise normal, neutral interstellar matter (e.g., Grenier et al. 2005; Abdo et al. 2010; Planck and Fermi Collaborations 2014).

In order to trace the DNM in this work we rely on the all-sky model of thermal dust emission by the Planck Collaboration XI (2014) based on Planck and IRAS data (Planck public data release R1.20). We test both the map of optical depth at 353 GHz, ${\tau }_{353}$, and the map of total dust radiance, R. The angular resolution of these maps is 5'. They are used to construct maps of the DNM as described in 4.2.

The column densities of atomic gas, $N({\rm{H}}\;{\scriptsize{\rm I}})$, in IVCs and HVCs and those in the local foreground gas can be separated based on Doppler velocities in order to separate contributions from CR interactions with gas in different locations in the Galaxy. To do so we adapt the procedure described in Abdo et al. (2010). The procedure used in this work has three steps.

• 1.  
  For each ROI we define some preliminary velocity boundaries that separate different components along the line of sight, namely low-velocity (local) gas and gas in HVCs/IVCs (Table 1).
• 2.  
  For each line of sight on the grid of the LAB survey we adjust the velocity boundaries so that they fall at the closest minimum (or inflection point) in the brightness temperature profile. The column density $N({\rm{H}}\;{\scriptsize{\rm I}})$ associated with each component is proportional to the integral of the brightness temperature over velocity within the adjusted velocity boundaries.
• 3.  
  For each line of sight on the ${\rm{H}}\;{\scriptsize{\rm I}}$ grid we fit to the brightness temperature profile a combination of Gaussian functions centered at the profile maxima (or inflection points if a maximum is not found for a velocity component) and use the results of the fit to correct the column densities of each component for the spillover from/into the adjacent components.

The procedure is illustrated for an example line of sight in Figure 2. Technical details of the procedure are illustrated by (Tibaldo 2011, 2.1.2, Appendix B).

Zoom In Zoom Out Reset image size

In region A we set the preliminary velocity boundaries at −40 km s⁻¹ to separate local gas from the low-latitude IV Arch, and at −90 km s⁻¹ to separate the low-latitude IV Arch from Complex A. At latitudes $b\lt 24^\circ$ there is a significant contribution at higher velocities from gas in the disk of the Galaxy belonging to the Perseus and outer arms. Because Abdo et al. (2010) found that the γ-ray emissivity of the gas in these arms is compatible with the local emissivity, we assign all of this gas to the local component. We note that in any case this low-latitude region is not included in the ROI, but only used as part of the border to ensure a proper convolution of the model with the LAT PSF.

In region B we set the preliminary velocity boundaries at −15 km s⁻¹ to separate local gas from the lower IV Arch, and at −60 km s⁻¹ to separate the lower from upper IV Arch. In the case of region C we set the preliminary velocity boundaries at −20 km s⁻¹ to separate local gas from the IV Spur, and at −90 km s⁻¹ to separate the IV Spur from high-velocity gas; no significant emission is present at high velocities and the corresponding map is not used in the construction of the γ-ray model.

The Gaussian fits reproduce the measured brightness temperature profiles within ∼5% on average (for the ratio of the integral of the absolute values of the differences between measured and fitted temperature along the line of sight to the integral of measured temperatures). The fit is significantly worse along a few lines¹¹ of sight with differences up to ∼25%. The magnitude of the correction for the spillover is on average 3% of the uncorrected column densities. (The average is the largest for ROI B, where it is ∼6%.) Even though the correction is only approximate, owing to the imperfect modeling of the temperature profiles by the Gaussian functions, the average magnitudes of the errors are $\lt 0.5\%$ of the total column densities. In a small subset of the lines of sight (∼0.5%) the fits do not converge properly, mostly due to the faintness of ${\rm{H}}\;{\scriptsize{\rm I}}$ emission in these high-latitude regions. In such cases we do not apply any correction for spillover.

We have calculated the $N({\rm{H}}\;{\scriptsize{\rm I}})$ column densities under the approximation of small optical depth. An approximate correction for ${\rm{H}}\;{\scriptsize{\rm I}}$ optical depth based on the assumption of a uniform spin temperature has often been used in the literature (e.g., Abdo et al. 2010; Planck and Fermi Collaborations 2014), but the low column densities in the regions studied make the effect of the correction small. Assuming a uniform spin temperature of 80 K (as in Planck Collaboration XXIV 2011) would increase the column densities on average at most by 5% for the local component in region A, and less for other components and regions. We also note that, by using dust maps to trace gas missed by the $N({\rm{H}}\;{\scriptsize{\rm I}})$ and ${W}_{\mathrm{CO}}$ maps as described in 4.2, we are effectively accounting also for ${\rm{H}}\;{\scriptsize{\rm I}}$ densities underestimated due to the approximation of small optical depth.

We replace the $N({\rm{H}}\;{\scriptsize{\rm I}})$ column density seen in ROI A at l = 1385, b = 37° with the average of neighboring lines of sight owing to an anomalous temperature profile indicative of an artifact in the LAB data.

The resulting $N({\rm{H}}\;{\scriptsize{\rm I}})$ maps are shown in Figures 3–5. Figures 4(a) and (b) show some level of anticorrelation between higher/lower $N({\rm{H}}\;{\scriptsize{\rm I}})$ intensity regions. This is due to some genuine anticorrelation between the intensities of ${\rm{H}}\;{\scriptsize{\rm I}}$ emission in velocity space. At the same time, the lower IV Arch shows a bright feature at $\sim -30$ km s⁻¹ that blends into the local ${\rm{H}}\;{\scriptsize{\rm I}}$ components, which makes the separation particularly challenging. This is addressed by the spillover correction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the work of Grenier et al. (2005) on the construction of models of interstellar γ-ray emission, dust residuals, i.e., total dust opacity or extinction minus the best-fit linear combination of ${\rm{H}}\;{\scriptsize{\rm I}}$ and ${W}_{\mathrm{CO}}$ maps, have been used to trace the DNM, neutral interstellar matter not properly traced by the linear combination of $N({\rm{H}}\;{\scriptsize{\rm I}})$ and ${W}_{\mathrm{CO}}$. However, the original implementation of this method has two limitations: (1) the fit of the dust map is biased by the missing DNM component; (2) the effect of the assumptions on the stated errors in the fit of the dust maps is overlooked.

Recently, the Planck and Fermi Collaborations (2014) performed an analysis of Planck and LAT data of the Chameleon interstellar complex where they refined the original method by: (1) developing a "circular" fitting method, where dust residuals are used to model the DNM in the γ-ray analysis, and γ-ray residuals to model the DNM in the dust analysis iteratively until a stable solution is found and (2) exploring several assumptions on the errors used in the fits and folding the differences into the errors on the final results.

The method developed by the Planck and Fermi Collaborations (2014) derives a map of the DNM column densities from γ rays to use in fitting the dust maps. Owing to the low gas column densities (on average one order of magnitude lower than in the Chameleon complex) and the limited statistics of the γ-ray data set, it was not applicable to our ROIs. We therefore developed an alternative method to acknowledge and address the issue of the missing DNM component in the dust fit that does not rely on γ-ray data and is described in the rest of the section.

Owing to the limited angular resolution of γ-ray data and the other ISM tracers, and to the input projections accepted by the LAT Science Tools, we reproject the maps from the Planck Collaboration XI (2014) onto the grid chosen for the γ-ray analysis in each ROI.

Under the assumptions that dust grains are well mixed with gas in the cold and warm phases of the ISM (e.g., Bohlin et al. 1978), and that the properties of dust grains and photon fields are uniform within each phase and complex, we build a model for the dust map D (i.e., the map of optical depth at 353 GHz, ${\tau }_{353}$, or radiance, R) by assuming that either quantity is proportional to the atomic and molecular column densities as traced by $N({\rm{H}}\;{\scriptsize{\rm I}})$ and ${W}_{\mathrm{CO}}$, respectively, plus an isotropic term representing any zero-level shift, for example, due to residual contamination from the cosmic infrared background.

$$\begin{eqnarray}{M}^{0}(l,b) & = & \underset{{\imath }=1}{\overset{3}{\displaystyle \sum }}{y}_{{\rm{H}}\;{\rm{I}},{\imath }}^{0}\cdot N({\rm{H}}\;{\scriptsize{\rm I}}){(l,b)}_{i}\\ & & +{y}_{\mathrm{CO}}^{0}\cdot {W}_{\mathrm{CO}}(l,b)+{D}_{\mathrm{iso}}^{0}\end{eqnarray} \tag{ 1 }$$

Here and in the following equations we indicate with ${y}_{{\rm{H}}\;{\rm{I}},{\imath }}$ the dust specific power (opacity) per hydrogen atom, and with y_CO the specific power (opacity) per ${W}_{\mathrm{CO}}$ unit in the molecular phase. We use indices ${\imath }=1,2,3$ to label the complexes separated kinematically along the line of sight as described in Section 4.1. The model M⁰ is fit to the dust map D using the minimum ${\chi }^{2}$ method. The ${\chi }^{2}$ is defined, for a generic model M, as

$$\begin{eqnarray}&&{\chi }^{2}=\underset{l,b}{\displaystyle \sum }\displaystyle \frac{{[D(l,b)-M(l,b)]}^{2}}{{\sigma }^{2}(l,b)}.\end{eqnarray} \tag{ 2 }$$

Following the work of the Planck and Fermi Collaborations (2014), we consider different possible assumptions on the uncertainties $\sigma (l,b)$. Our baseline assumption is that the uncertainties are mainly due to imperfections of the model because of the assumption of proportionality between radiance/opacity and gas column densities. In the absence of a reliable estimate of the uncertainties on M, we assume the uncertainties to be $\sigma (l,b)\propto M(l,b)$. We empirically consider two alternative hypotheses to assess the impact of this choice:

• 1.  
  $\sigma (l,b)\propto D(l,b)$, which tests how the agreement between M and D affects the results of the fits and
• 2.  
  $\sigma (l,b)\propto {\sigma }_{\tau }(l,b)$, where ${\sigma }_{\tau }(l,b)$ is the map of statistical uncertainties on ${\tau }_{353}$ provided by the Planck Collaboration XI (2014) (for fits of ${\tau }_{353}$ only).

We initially assume that the proportionality coefficients are 30% for the model or the dust map itself, and 1 when using ${\sigma }_{\tau }(l,b)$. These are adjusted later as a part of the iterative procedure to obtain a reduced ${\chi }^{2}=1$. The ${\chi }^{2}$ is minimized to determine the best-fit coefficients, ${\tilde{y}}^{0}$ and ${\tilde{D}}_{\mathrm{iso}}^{0}$.

From the results we derive a dust residual map

$$\begin{eqnarray}{D}_{\mathrm{res}}^{0}(l,b) & = & D(l,b)-\left[\underset{{\imath }=1}{\overset{3}{\displaystyle \sum }}{\tilde{y}}_{{\rm{H}}\;{\rm{I}},{\imath }}^{0}\cdot N({\rm{H}}\;{\scriptsize{\rm I}}){(l,b)}_{i}\right.\\ & & +\left.{\tilde{y}}_{\mathrm{CO}}^{0}\cdot {W}_{\mathrm{CO}}(l,b)+{\tilde{D}}_{\mathrm{iso}}^{0}\phantom{\rule[-0000]{0ex}{3.40ex}}\right].\end{eqnarray} \tag{ 3 }$$

The map is then filtered to extract the significant positive residuals associated with the DNM. The residual map is convolved with a Gaussian kernel with standard deviation of 1° in order to suppress statistical fluctuations while retaining large-scale real features. Then we build a histogram of the residuals and we fit a Gaussian curve to the core of the distribution, within ±2 standard deviations of the distribution itself. We set a threshold at 2 standard deviations of the fitted Gaussian curve, and produce the filtered residual map ${\bar{D}}_{\mathrm{res}}^{0}$ by retaining pixels of the original residual map when the values in the convolved residual map are above the threshold and setting other pixels to 0.

Then we proceed iteratively to refine the model. At each iteration $\alpha =1,2,\ldots$:

• 1.  
  The model is defined:

  $$\begin{eqnarray}{M}^{\alpha }(l,b) & = & \underset{{\imath }=1}{\overset{3}{\displaystyle \sum }}\;{y}_{{\rm{H}}\;{\rm{I}},{\imath }}^{\alpha }\cdot N({\rm{H}}\;{\scriptsize{\rm I}}){(l,b)}_{i}+{y}_{\mathrm{CO}}^{\alpha }\\ & & \times {W}_{\mathrm{CO}}(l,b)+{y}_{\mathrm{res}}^{\alpha }\cdot {\bar{D}}_{\mathrm{res}}^{\alpha -1}(l,b)+{D}_{\mathrm{iso}}^{\alpha },\end{eqnarray} \tag{ 4 }$$
• 2.  
  We fit the model ${M}^{\alpha }$ to the dust map D minimizing the ${\chi }^{2}$ as defined in Equation (2);
• 3.  
  The dust residual map is built:

  $$\begin{eqnarray}{D}_{\mathrm{res}}^{\alpha }(l,b) & = & D(l,b)-[\underset{{\imath }=1}{\overset{3}{\displaystyle \sum }}{\tilde{y}}_{{\rm{H}}\;{\rm{I}},{\imath }}^{\alpha }\cdot N({\rm{H}}\;{\scriptsize{\rm I}}){(l,b)}_{i}\\ & & +{\tilde{y}}_{\mathrm{CO}}^{\alpha }\cdot {W}_{\mathrm{CO}}(l,b)+{\tilde{D}}_{\mathrm{iso}}^{\alpha }]\phantom{\rule[-0000]{0ex}{2.95ex}}\end{eqnarray} \tag{ 5 }$$

  where ${\tilde{y}}^{\alpha }$ and ${\tilde{D}}_{\mathrm{iso}}^{\alpha }$ are the best-fit coefficients.
• 4.  
  We filter ${D}_{\mathrm{res}}^{\alpha }$ to extract the significant positive residuals ${\bar{D}}_{\mathrm{res}}^{\alpha }$ following the procedure described for iteration 0.

The iterative procedure is stopped when the two following conditions are simultaneously satisfied: (1) the differences between all the best-fit coefficients in iteration α and iteration $\alpha -1$ are compatible with 0 within $1\sigma$ statistical uncertainties from the fit and (2) the difference between minimum ${\chi }^{2}$ in iteration α and iteration $\alpha -1$ is $\lt 1$.

After the first pass just described, the scaling coefficients for the uncertainties $\sigma (l,b)$ are re-defined, so that the reduced ${\chi }^{2}=1$ in the latest iteration. Requiring the reduced ${\chi }^{2}$ to equal 1 ensures that the statistical uncertainties from the fit are truly representative of the differences between dust maps and models. This, in turn, also makes the criteria to stop the iteration of the fitting procedure more meaningful, even if only in a data-driven fashion. Using the new scaling coefficients we perform a second pass, after which the final reduced ${\chi }^{2}$ values are found to be 1 within 1% in all ROIs and variants of the maps. The distribution of the residuals in the final iteration does not show any significant negative values, which demonstrates how the procedure is appropriately addressing the initial bias due to the missing template for the DNM (see Appendix A for further details). The scaling coefficients for the errors on ${\tau }_{353}$ provided by the Planck Collaboration XI (2014) are found to be systematically greater than 1. The error map represents only statistical uncertainties from a fit where the dust thermal emission spectral energy distribution was fitted using a modified blackbody spectrum. Hence, a scaling coefficient greater than 1 may capture the presence of systematic errors, or discrepancies between the distribution in the sky of dust optical depth and our simple model based on the assumption that the optical depth scales with gas column densities.

Because the Planck Collaboration XI (2014) recommends the use of R as the best dust column-density tracer at high latitudes, we take as our baseline DNM templates those obtained from R assuming that the error is proportional to the model. These maps are shown in Figure 6. We note that the DNM component is particularly abundant in ROI A, where it follows a known interstellar structure dubbed the NCP Loop (Meyerdierks et al. 1991), possibly related to the Ursa Major molecular complex (e.g., Pound & Goodman 1997). Further results from the analysis described in this section are discussed in Appendix A.

Zoom In Zoom Out Reset image size

The γ-ray emission measured using the LAT can be modeled as the sum of diffuse components and discrete sources. Diffuse emission is produced in the interstellar space of the Milky Way by energetic interactions of CRs with gas, via nuclear collisions leading to γ-ray production mainly through neutral pion decay and by electron Bremsstrahlung, and with low-energy radiation fields, due to IC scattering by electrons.

Multiple studies of interstellar γ-ray emission have shown that CR fluxes can be approximated as uniform on the scales of interstellar complexes located away from potential CR sources (e.g., Abdo et al. 2010; Planck and Fermi Collaborations 2014). Under this hypothesis, the γ-ray intensities from interactions with gas are proportional to the target gas column densities. We can therefore model the diffuse intensity as a combination of the gas tracer maps described in Section 4. We assume the spectral shape to be that measured using LAT data for the local interstellar space (Casandjian 2015) in order to perform the fit over a broad energy range. This assumption will be critically evaluated against the data in the assessment of systematic uncertainties (Section 5.4). For each tracer map we allow a free scaling coefficient in the model, to define the overall CR flux in that component.¹²

The IC component is modeled based on the work by de Palma et al. (2013). They considered eight models of IC emission calculated using the GALPROP¹³ code (Moskalenko & Strong 2000; Strong et al. 2000; Porter et al. 2008) for different values of some model parameters, and tuned them to the LAT data over the whole sky. The work by de Palma et al. (2013), originally based on the LAT P7 data, was recently updated for a P7REP data selection consistent with our data set. For the IC component we allow in the fit a free scaling coefficient, as well as a log-parabolic spectral correction to mitigate the effect of uncertainties in the local electron spectra. Among the models in de Palma et al. (2013), we arbitrarily select that with the CR source distribution based on the work by Lorimer (2004), a maximum height of the CR confinement halo of 10 kpc and a ${\rm{H}}\;{\scriptsize{\rm I}}$ spin temperature of 150 K. The other models presented in de Palma et al. (2013) will be used later in Section 5.4 to assess the impact of this choice on the results.

Finally, there is a diffuse component with almost isotropic distribution over the sky interpreted as the superposition of extragalactic diffuse γ-ray emission and residual backgrounds from CR interactions in the LAT misclassified as γ rays. This is modeled using a spectral shape derived from the LAT data (see de Palma et al. 2013 for the details of the procedure) consistently with our IC models, separately for γ rays converting in the front and back sections of the LAT tracker, which are characterized by different residual background rates relative to their acceptances. Owing to limitations in the assumption that the residual CR background can be treated as an isotropic γ-ray term, we allow a free scaling coefficient in the fit also for this term. Note that, due to the limited size of the ROIs and the quasi-isotropic distribution of IC γ-ray emission in the models, both IC and isotropic components have to be considered as effective background templates in our fits and the respective scaling factors are not the object of any interpretation in this work.

We incorporate point sources from the 3FGL list (Fermi-LAT Collaboration 2015). For each ROI we include in the model all the sources within 10° from the ROI border. The source positions are fixed at their catalog values. We assume for each source the spectral shape given in the catalog as well. For sources outside the ROI but within 10° from the ROI border all the spectral parameters are fixed at their catalog values. For sources inside the ROI we leave all the spectral parameters free if the average significance S of the source in 3FGL is $\gt {\theta }_{1}$, just the global normalization free if ${\theta }_{2}\lt S\lt {\theta }_{1}$, and all the parameters fixed otherwise. The thresholds ${\theta }_{1}$ and ${\theta }_{2}$ are chosen for each ROI so that the number of free parameters associated with sources is $\lt 35$ to prevent convergence problems in the fitting procedure. The pairs of thresholds $({\theta }_{1},{\theta }_{2})$ chosen are $(15\sigma ,9\sigma )$, $(30\sigma ,13\sigma )$, and $(10\sigma ,5\sigma )$, for regions A, B, and C, respectively. We verified that none of the 3FGL sources are positionally coincident with structures in the $N({\rm{H}}\;{\scriptsize{\rm I}})$ maps in Figures 3–5. Following a visual inspection of the residual maps after a preliminary analysis we left free in addition to those sources also 3FGL J0809.5+5342 in ROI A and 3FGL J1159.2+2914 in ROI C. We will show later in Figures 7–9 that additional two years of data in our analysis with respect to 3FGL do not introduce any additional sources that significantly affect the results of the analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To summarize, for each ROI the model for the γ-ray intensities is described by Equation (6).

$$\begin{eqnarray}{I}_{\gamma }(l,b,E) & = & \underset{{\imath }=1}{\overset{3}{\displaystyle \sum }}\;\left[{x}_{{\rm{H}}\;{\rm{I}},{\imath }}\cdot N({\rm{H}}\;{\scriptsize{\rm I}}){(l,b)}_{{\imath }}\;{q}_{\mathrm{LIS}}(E)\right]\\ & & +{x}_{\mathrm{CO}}\cdot {W}_{\mathrm{CO}}(l,b)\;{q}_{\mathrm{LIS}}(E)\\ & & +{x}_{\mathrm{DNM}}\cdot \bar{D}(l,b)\;{q}_{\mathrm{LIS}}(E)\\ & & +{x}_{\mathrm{IC}}\cdot \mathrm{IC}(l,b,E){\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{{\alpha }_{\mathrm{IC}}+{\beta }_{\mathrm{IC}}\mathrm{ln}\frac{E}{{E}_{0}}}\\ & & +{x}_{\mathrm{Iso}}\cdot {I}_{\mathrm{iso}}(E)+\mathrm{SRC}\end{eqnarray} \tag{ 6 }$$

There, ${q}_{\mathrm{LIS}}(E)$ stands for the tabulated emissivity spectrum of ${\rm{H}}\;{\scriptsize{\rm I}}$ in the local interstellar space from Casandjian (2015). $N{({\rm{H}}\;{\rm{I}})}_{{\imath }}$ are the ${\rm{H}}\;{\scriptsize{\rm I}}$ column density maps that we constructed in Section 4.1, so that, for each complex ${x}_{{\rm{H}}\;{\rm{I}}}$ is the scaling factor with respect to the local emissivity. Kinematically separated complexes along the line of sight are numbered as (1) low-velocity, i.e., local gas, (2) and (3) IVCs and HVCs. ${W}_{\mathrm{CO}}$ is the CO intensity map described in Section 3.2 (used only for ROI A); hence x_CO can be interpreted as the scaling factor that transforms the local ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity into γ-ray emissivity per ${W}_{\mathrm{CO}}$ intensity unit in the region studied. Similarly, $\bar{D}$ is one of the templates of the DNM constructed in Section 4.2, and x_DNM is the scaling factor that transforms the local ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity into γ-ray emissivity per dust radiance (or optical depth) unit. IC represents the IC emission model: x_IC is an overall global scaling factor, ${\alpha }_{\mathrm{IC}}$ and ${\beta }_{\mathrm{IC}}$ are the parameters of the log-parabolic spectral correction, while E₀ is a reference energy which was set to 700 MeV (the choice is arbitrary because the parameter is degenerate with the others). Finally, I_iso is the tabulated isotropic background spectrum (different for front- and back-converting events), with x_Iso as the common scaling factor, and SRC represents the model for individual sources described above.

Using the data set described in Section 3.1, for each ROI in Table 1 we build binned count cubes with $0\buildrel{\circ}\over{.} 25$ spacing in angle and 15 logarithmically spaced energy bins from 300 MeV to 10 GeV for front- and back-converting events. The model in Equation (6) is fit to the data by using the binned likelihood analysis procedure implemented in the Science Tools that takes into account the LAT exposure and PSF. We employ a joint likelihood technique to independently treat the front- and back-converting events in the fit in order to exploit the better angular resolution of the former. We use the P7REP_V15_CLEAN LAT instrument response functions.

The maximum likelihood values of the coefficients of the interstellar emission models for the baseline analysis are reported in Table 2. In Figures 7–9 we show the count maps and residual maps summed over the entire energy range. The residual maps do not show any large-scale structures correlated with the templates used to model the γ-ray emission. Some localized residuals may hint at missing point sources not included in 3FGL, which was developed using the first 48 months of LAT observations, as opposed to the 73 months used in our analysis.

Table 2.  Interstellar Best-fit Coefficients from the Baseline Analysis

Parameter	ROI A	ROI B	ROI C
${x}_{{\rm{H}}\;{\rm{I}},1}$	0.999 ± 0.017	0.94 ± 0.07	1.08 ± 0.04
${x}_{{\rm{H}}\;{\rm{I}},2}$	0.94 ± 0.16	1.08 ± 0.04	0.68 ± 0.06
${x}_{{\rm{H}}\;{\rm{I}},3}$	${0.00}_{-0.0}^{+0.06}$	${0.00}_{-0.0}^{+0.10}$	⋯
xCOa	2.3 ± 0.2	⋯	⋯
xDNMb	0.24 ± 0.03	${0.00}_{-0.0}^{+0.11}$	0.34 ± 0.09
xIC	1.19 ± 0.05	1.31 ± 0.05	1.82 ± 0.09
${\alpha }_{\mathrm{IC}}$	−0.158 ± 0.014	−0.06 ± 0.02	0.02 ± 0.05
${\beta }_{\mathrm{IC}}$	0.045 ± 0.008	0.02 ± 0.01	−0.026 ± 0.014
xIso	0.978 ± 0.017	1.03 ± 0.02	0.96 ± 0.04

Notes. See Equation (6) for a definition of the model parameters. The numbering scheme for the ${x}_{{\rm{H}}\;{\rm{I}},{\imath }}$ coefficients is: (1) low-velocity, i.e., local gas, (2) and (3) IVCs and HVCs.

Asymmetric error bands are computed using Minos when the best-fit value is at the boundary imposed during the likelihood optimization process.

^a10²⁰ cm⁻² (K km s⁻¹)⁻¹. ^b10³² sr W⁻¹.

Download table as:  ASCIITypeset image

For each target (HVC or IVC) in Table 1 we calculate the test statistic TS, defined as

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{ln}\left(\displaystyle \frac{{\mathcal{L}}}{{{\mathcal{L}}}_{0}}\right),\end{eqnarray} \tag{ 7 }$$

where ${\mathcal{L}}$ is the maximum likelihood for the model that includes the source of interest, while ${{\mathcal{L}}}_{0}$ is the maximum likelihood of the minimal model where the contribution from the source is set to null. TS is expected to be distributed as a ${\chi }^{2}$ with a number of degrees of freedom given by the difference of the number of free parameters in the two models (e.g., Mattox et al. 1996). We consider a target cloud as detected if $\mathrm{TS}\gt 25$, which formally corresponds to a 5σ confidence level that the likelihood improvement is not due to statistical fluctuations. As detailed in Protassov et al. (2002), the distribution of TS may deviate from a ${\chi }^{2}$ if the minimal model corresponds to a point on the topological border of the phase space of the parameters of the model including the source, which applies to our situation. However, the TS values for our detected targets all greatly exceed the threshold.

The TS values are reported in Table 3. For target HVCs and IVCs that are not detected we calculate a 95% confidence level (c.l.) upper limit on x_{H i} based on a Bayesian integration of the likelihood profile as applied for measurements in the 2FGL catalog (Nolan et al. 2012, 3.5).

In Figure 10 we show the γ-ray emission components associated with detected IVCs, evaluated by subtracting from the total γ-ray counts the components attributed to individual sources and foreground/background interstellar emission by the likelihood fit. Contours from the $N({\rm{H}}\;{\scriptsize{\rm I}})$ maps inferred from the 21 cm observations in Figures 3–5 are overlaid. The morphologies of the γ-ray detected signals follow closely the known ${\rm{H}}\;{\scriptsize{\rm I}}$ structures in the low-latitude IV Arch and lower IV Arch. In the case of the IV Spur a cluster of γ-rays is associated with the densest part of the cloud, but the visual comparison is hampered by limited statistics in the γ-ray data set.

Zoom In Zoom Out Reset image size

A key assumption of our model is that the emission of the gaseous components has a spectrum described by the emissivity spectrum of local gas q_LIS from Casandjian (2015), modified only by a global scaling factor over the whole energy range from 300 MeV to 10 GeV. While this assumption is expected to hold for the local (low-velocity) gas, there may be spectral variations for clouds at great distance from the disk.

For clouds detected in our analysis we can use the data to verify how well the q_LIS spectrum reproduces the observed one. To do so we perform the γ-ray analysis independently in six broad energy bins logarithmically spaced from 300 MeV to 10 GeV. In these fits the spectral parameters are frozen to their determination over the entire energy range except for the normalizations. Figure 11 shows as a reference the scaling factors of local ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity in the three regions studied. Figure 12 shows the scaling factors for the emissivity of all the detected IVCs with respect to local.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The scaling factors are all consistent within statistical uncertainties with the value determined from the fit over the entire energy range. In region A the residuals show a negative trend that indicates a spectrum softer than the average in the local ISM. Similar trends may be present for both the local ISM and IV gas in region B. It is not clear if this is real or due to any small imperfections in the modeling, but further investigation on this aspect is beyond the scope of this work and left for future studies. For the scope of our analysis it is justified to use the global scaling factors over the whole energy range because they are compatible within 1σ with the weighted average of the values for the separate energy bins.

5.4.1. Uncertainties of Inputs to the γ-ray Interstellar Emission Model

Some inputs to the model in Equation (6) are subject to important uncertainties, most notably the DNM column densities, IC emission, and the spectrum of undetected clouds. We explored the associated systematic uncertainties by changing the modeling assumptions and evaluating the impact on the best-fit parameters for the γ-ray analysis.

For each region we modified the baseline DNM template by using the alternative templates from Section 4.2 based on different assumptions for the errors in the determination of the DNM template, and replacing dust radiance with optical depth.

We also changed the IC model by using alternative models as described in de Palma et al. (2013), which included assuming maximum heights of the CR confinement halo of 4 and 10 kpc, and replacing the CR source distribution inferred from pulsar observations (Lorimer 2004) with that by Case & Bhattacharya (1998) based on supernova remnant observations.¹⁴ These are two of the most important input parameters to GALPROP to determine the morphology and spectrum of IC emission (Ackermann et al. 2012a).

Finally, we have shown in Section 5.3 that for detected targets the assumption that the emission spectrum follows the emissivity in the local ISM is appropriate, but we cannot do the same for targets that are not significantly detected, i.e., Complex A and the Upper IV Arch. We explored possible variations with respect to the local emissivity spectrum based on the set of propagation models in Ackermann et al. (2012a). We calculated the ratios of the emissivity spectra in those models at all locations in the outer Galaxy (where our clouds are located) over the emissivity spectrum at the solar circle. We parametrized the extrema of these ratios with two functions of energy

$$\begin{eqnarray}&&{f}_{1}^{-}={\left(\displaystyle \frac{E}{1\ \mathrm{GeV}}\right)}^{-0.1},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{f}_{2}^{-}={\left(\displaystyle \frac{E}{1\ \mathrm{GeV}}\right)}^{-[0.08+0.07\mathrm{ln}(E/1\ \mathrm{GeV})]}.\end{eqnarray} \tag{ 9 }$$

Equation (8) captures the softening due to energy-dependent escape that is expected near the boundaries of the confinement region in the outer Galaxy for sources located mostly in the inner Galaxy. In addition, Equation (9) captures a low-energy hardening predicted by high-reacceleration models in Ackermann et al. (2012a) at larger Galactocentric radii.¹⁵ Note that both functions indicate a softer spectrum. In order to be conservative we explore the extreme assumptions that the emissivity spectrum may actually be harder at the locations of the undetected targets by inverting the signs of the exponents

$$\begin{eqnarray}&&{f}_{1}^{+}={\left(\displaystyle \frac{E}{1\ \mathrm{GeV}}\right)}^{+0.1},\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&{f}_{2}^{+}={\left(\displaystyle \frac{E}{1\ \mathrm{GeV}}\right)}^{+[0.08+0.07\mathrm{ln}(E/1\ \mathrm{GeV})]}.\end{eqnarray} \tag{ 11 }$$

We use these functions to modify the emissivity spectra of the undetected clouds, by multiplying the local emissivity spectrum q_LIS in Equation (6) by the correction function and scaling it so that the emissivity integrated in the 300 MeV to 10 GeV range is the same as for q_LIS. In this way, we can still interpret the upper limits on x_{H i} for the undetected components as scaling factors with respect to the local emissivity, and, at the same time, explore the impacts on the other fit parameters.

We repeat the baseline analysis described in Section 5.2 for all the combinations of the model variations described above, and, for each parameter, we take the extreme variations with respect to the baseline fit as our estimate of the systematic uncertainty related to the model. In the cases of upper limits for undetected targets we use the worst, i.e., largest upper limits in the discussion that follows.

We summarize the scaling factors for ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity in the different regions and complexes in Table 4, including the spread due to inputs to the γ-ray interstellar emission model as determined through the methodology described in this subsection. In general, the emissivity spread for the detected complexes is driven by the choice of the DNM template, and, to a lesser extent, the choice of IC model for local gas (which is less structured, hence more degenerate with IC emission). For non-detected complexes the systematic uncertainties are dominated by the choice of the emissivity spectrum. The results in ROI B show a stronger dependence on the assumed IC model, which affects also the determination of the lower IV Arch emissivity and the upper limit on the emissivity of the upper IV Arch at 2 statistical σ level.

Table 4.  ${\rm{H}}\;{\scriptsize{\rm I}}$ Emissivity Scaling Factors

Region	Complex	Scaling Factor	Stata	Sys (Model)b	Sys (Jackknife)c
A	Local	0.99	0.02	${}_{-0.14}^{+0.04}$	0.08
	Low-Latitude IV Arch	0.94	0.15	${}_{-0.01}^{+0.26}$	0.08
	Complex A	$\lt 0.21$		$\lt 0.22$	$\lt 0.001$
B	Local	0.94	0.07	${}_{-0.22}^{+0.13}$	0.08
	Lower IV Arch	1.08	0.04	${}_{-0.08}^{+0.10}$	0.09
	Upper IV Arch	$\lt 0.23$		$\lt 0.45$	$\lt 0.01$
C	Local	1.08	0.04	${}_{-0.09}^{+0.05}$	0.09
	IV Spur	0.68	0.06	${}_{-0.03}^{+0.09}$	0.06

Notes. Upper limits are provided at 95% confidence level.

^a $1\sigma$ statistical uncertainty from the likelihood analysis. ^bSystematic spread from varying some inputs to the γ-ray interstellar emission model (see Section 5.4.1). ^c $1\sigma$ uncertainty or 95% c.l. upper limit from the distributions obtained with the jackknife test (see Section 5.4.2).

Download table as:  ASCIITypeset image

5.4.2. Jackknife Tests

So far we have characterized the statistical uncertainties in the fit parameters in Section 5.2 and systematic uncertainties related to some inputs to the model in Section 5.4.1. However, other potential sources of systematic uncertainties are related to the assumption that the emissivity scaling coefficients for each component are constant across a given ROI, and to the presence of regions with large deviations or mismodeled individual sources driving the fits.

To characterize the magnitude of these uncertainties we perform jackknife tests where we repeat the analysis described in Section 5.2 150 times for each ROI, each time masking a circular region with random center within the ROI and covering 20% of the ROI area. The mask size was chosen to be much larger than the 68% event containment radius of the LAT for front-converting photons at the lowest end of the energy range at 300 MeV ($\sim 1\buildrel{\circ}\over{.} 5$) and large enough to probe possible mild variations of the CR densities across the complexes studied.

Figures 13–15 show the distributions of the ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity scaling coefficients for the three ROIs. The distributions all peak within $1\sigma$ from the best-fit values obtained from the baseline analysis and have standard deviations smaller than or comparable to the $1\sigma$ statistical uncertainties from the γ-ray fits. This shows that from a statistical point of view the obtained parameter values are representative of the entire ROIs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In ROI A, Figure 13, the distribution of the scaling coefficients for Complex A, ${x}_{{\rm{H}}\;{\rm{I}},3}$, narrowly peaks at zero, showing that this parameter is always at the boundary of the allowed values. In the rest of the work, we will make use only of the upper limit on the Complex A emissivity that, on the contrary, is a robustly determined quantity. Similar considerations hold for the upper IV Arch in ROI B. In all ROIs, most notably in ROI A, the distribution of the scaling coefficients ${x}_{{\rm{H}}\;{\rm{I}},2}$ for the IVCs show tails at low or high values, which suggests that emissivities lower or higher than average are possible within the gas encompassed by the respective templates.

To take into account these uncertainties in our estimates of the emissivities we compute, based on the distributions of the fit parameters in Figures 13–15, the $1\sigma$ confidence intervals for detected complexes and 95% c.l. upper limits for complexes that are not significantly detected, as reported in Table 4 (details of the derivation are given in Appendix B). For the detected complexes the confidence intervals are of the same magnitude as the statistical errors. Hence we will combine them with the other uncertainties for the discussion that follows. In the case of the complexes that are not significantly detected the sources of uncertainty probed through the jackknife tests (such as possible variations of the CR densities across the complexes studied) have an impact on the emissivity upper limits much smaller than those due to uncertainties in the analysis model inputs and those encoded by the statistical error from likelihood fit, so they will be neglected.

The ${\rm{H}}\;{\scriptsize{\rm I}}$ emissivity scaling coefficients in Table 4 can be converted into emissivity values. To do so we multiply them by the integrated emissivity between 300 MeV and 10 GeV from Casandjian (2015), i.e., $7.1\times {10}^{-27}$ s⁻¹ sr⁻¹ H⁻¹. Note that this rests on the measurement described in Section 5.3 that the shape of the emissivity spectrum is the same as in the local ISM for detected complexes, and is valid within the extreme variations parametrized in Equations (8)–(11) for those that are not significantly detected.

The conversion into emissivities requires taking into account one additional source of systematic uncertainty, i.e., the characterization of the LAT instrument response, and in particular its effective area. Based on the estimate of the systematic uncertainties of the LAT effective area in Ackermann et al. (2012c), the systematic error on the local emissivity by Casandjian (2015) integrated between 300 MeV and 10 GeV is $0.6\times {10}^{-27}$ s⁻¹ sr⁻¹ H⁻¹.

Note also that in many HVCs and IVCs the ratio of the ionized gas mass to the neutral atomic gas mass is much larger than in the local ISM and close to unity (e.g., Wakker et al. 2008). If the column densities of ionized gas are correlated with those of neutral gas the scaling factors in the γ-ray fits will be overestimated to encompass the emission from ionized gas. In this respect our estimates of scaling factors and emissivities should be regarded as an upper limit to the real values. The same consideration applies to the potential presence of CO-dark molecular gas in HVCs/IVCs. Aside from the emissivity scaling coefficient slightly larger than 1 for the lower IV Arch there is no clear indication of such effect from our analysis.

The final results for the emissivities of the targets are summarized in Table 5. For the purpose of interpretation we will rely on those findings and on the emissivity scaling factors in Table 4, which to a good approximation are not subject to uncertainties due to the LAT instrument effective area.

Table 5.  Emissivities of HVCs and IVCs

Region	Complex	z	γ-ray emissivity	Stata	Sys (LAT)b	Sys (model)c	Sys (Jackknife)d
		(kpc)	(10−27 s−1 sr−1 H−1)
A	Low-latitude IV Arch	0.6–1.2	6.7	1.1	0.6	${}_{-0.1}^{+1.8}$	0.6
	Complex A	2.6–6.8			$\lt 1.7$
B	Lower IV Arch	0.4–1.7	7.7	0.3	0.7	${}_{-0.6}^{+0.7}$	0.6
	Upper IV Arch	0.7–1.7			$\lt 3.5$
C	IV Spur	0.3–2.1	4.8	0.4	0.4	${}_{-0.2}^{+0.6}$	0.4

Notes. z brackets are replicated from Table 1 for the reader's convenience. Data sources are listed there.

γ-ray emissivities are integrated between 300 MeV and 10 GeV.

We provide 95% confidence level upper limits for the worst-case scenario from the model variations considered and taking into account uncertainties due to the LAT effective area.

^a $1\sigma$ statistical uncertainty from the likelihood analysis. ^bError from uncertainties in the LAT effective area. ^cSystematic spread from varying some inputs to the γ-ray interstellar emission model (see Section 5.4). ^d $1\sigma$ uncertainty or 95% c.l. upper limit from the distributions obtained with the jackknife test (see Section 5.4.2).

Download table as:  ASCIITypeset image

The γ-ray emissivities are a proxy for the CR densities in the complexes studied. Owing to the cross sections for γ-ray production in nucleon–nucleon inelastic collisions (e.g., Kamae et al. 2006), the γ-ray energies studied between 300 MeV and 10 GeV constrain the CR densities from energies of ∼3 GeV/nucleon up to ∼200 GeV/nucleon. Note that while we do not expect significant changes in the He/H ratio because both elements are primordial, the abundances of metals, i.e., elements heavier than He, are known to vary in HVCs and IVCs and depending on the metallicity of target gas the emissivity per H atom can vary up to $\lt 5\%$ of the total (e.g., Mori 2009). Furthermore, the total emissivity depends also on the CR composition and spectra of the different elements at a level up to 20%–30% of the total emissivity per H atom (Kachelriess et al. 2014). Part of the emissivity variations (Table 5) or, equivalently, differences in the emissivity scaling factors (Table 4) therefore can be attributed to variations in the CR elemental composition and of their spectra, as well as to the uncertain metallicity of the targets (Wakker 2001).

Our measurement of the local gas emissivity in each ROI agrees well with unity, indicating the robustness of the analysis procedure and assumptions. The emissivity of the IV Arch indicates a $\sim 50\%$ decrease of the CR densities with respect to the local value within 2 kpc from the Galactic plane, while the HVC Complex A provides the strongest limit on the CR densities at a few kpc above the disk. The implications of these results are discussed in the following sections.

If CRs originate in the disk of the Milky Way and diffusively propagate in the surrounding halo we expect the CR densities, and hence the γ-ray emissivities, to decrease with distance z from the disk itself.¹⁶ We test this hypothesis against our data in Table 4 by computing the Kendall τ rank correlation coefficient. Using the emissivity scaling coefficients in Table 4, and adopting for z the mean within the bracket for targets in Table 1, and z = 0 kpc for the local gas in each ROI, we obtain $\tau =-0.54$. The negative value indicates, indeed, a negative correlation of emissivity with $| z|$.

We calculate the significance of this trend by comparing the τ of the actual data with a distribution of the τ coefficients obtained from the null hypothesis of no correlation between emissivity and z. We generate null hypothesis data sets starting from the real data set using the following procedure.

• 1.  
  We draw a set of emissivity scaling coefficients with normal distribution centered at the measured values and with σ equal to the sum in quadrature of statistical uncertainties and systematic uncertainties from the jackknife tests in Table 4 for significantly detected targets and local foreground gas in each ROI.
• 2.  
  We add upper limits at the values measured in Table 4 for targets which are not significantly detected.
• 3.  
  We randomly shift both scaling factors and upper limits with uniform probability distribution within the systematic uncertainty brackets from the model input variations in Table 4 (assuming the bracket to be symmetric for upper limits).
• 4.  
  We draw a set of z values with uniform probability distribution within the brackets in Table 1, and add three elements with z = 0 to represent the local foreground gas in each ROI.
• 5.  
  We randomly shuffle in an independent fashion the emissivity scaling factors and z sets so obtained to be consistent with the null hypothesis of no correlation between the two quantities.

The distribution of 30,000 null hypothesis data sets shows that there is a $2.5\%$ chance probability to obtain a τ coefficient $\lt -0.54$.

Therefore, our measurements provide evidence at 97.5% confidence level that the γ-ray emissivity per H atom decreases with distance from the disk of the Galaxy. This corroborates the notion that CRs at the relevant energies from GeV to TeV originate in the Galactic disk by directly tracing the distribution of CRs in the halo for the first time. Our result goes far beyond the test of the Galactic origin of CRs proposed by Ginzburg (1973), i.e., measurement of the γ-ray emissivity of Magellanic Clouds, that was performed by Sreekumar et al. (1993).

We note that we also expect a decrease in gas metallicity, hence emissivity per H atom, with distance from the disk if the most distant targets like Complex A are primeval gas falling into the Milky Way (Wakker 2001). However, the expected effect of at most $\sim 5\%$ on the γ-ray emissivities (Mori 2009) is not sufficient to explain the magnitude of the variations observed between the scaling coefficients in Table 4. Changes in the spectra of the different CR elements could also cause emissivity variations. However, the magnitude of this effect is estimated to be at most 20%–30% of the total emissivity per H atom (taking into account He in both CRs and the ISM), lower than the decrease of emissivity of $\gt 50\%$ for the Upper IV Arch and of $\gt 75\%$ for Complex A that we measured.

Additionally, in the region of the outer Galaxy where our targets are located the distance to the center of the Milky Way increases with increasing distance from the disk as well. While some propagation models advocate a sizable decrease of CR densities with increasing Galactocentric radius, observations (Abdo et al. 2010; Ackermann et al. 2011) indicate that the gradient of CR densities is marginal up to 15 kpc from the Galactic center. We conclude that the observed emissivity decrease as a function of z is more likely due to diffusion in the halo. A quantitative comparison with propagation models follows in the next section.

We compare our results with predictions from GALPROP. We will rely on the set of GALPROP models in Ackermann et al. (2012a) that have been compared extensively to direct CR measurements and LAT data, and which include a number of variations of the most relevant and uncertain parameters of the calculations. In Ackermann et al. (2012a) the CR transport equations were solved for the case of cylindrical symmetry, based on the assumption that CRs are injected in the disk of the Milky Way and propagate in a cylindrical volume with boundaries at maximum radius from the center of the Galaxy R_max and maximum height from the disk z_max where the CR particle number densities are required to go to zero.

Figure 16 shows the emissivity scaling factors from Table 4 against predictions from the models in Ackermann et al. (2012a). The model predictions are calculated for distances along the line of sight approximately corresponding to the peak gas column densities for each ROI.¹⁷ We normalized the model emissivities to the value in the disk at the solar circle. The model emissivities at the solar circle integrated from 300 MeV to 10 GeV in γ-ray energy were found to be within $\sim 7\%$ from the measurement in Casandjian (2015). Note that the calculation of the emissivities in Ackermann et al. (2012a) only accounts for interactions between CR and interstellar gas nuclei with atomic number $\leqslant 2$. Rescaling to the locally measured emissivity takes into account the contribution from metals in the targets and heavier CR species. Variations of the target gas metallicities are not taken into account in these models, and could produce a further decrease of emissivity up to at most 5% with increasing z (Mori 2009).

Zoom In Zoom Out Reset image size

The model parameter with the largest impact on the vertical gradient of CR densities, or, equivalently, γ-ray emissivities, is the maximum height of the confinement halo z_max. Models in Ackermann et al. (2012a) assume values of z_max between 4 and 10 kpc. While there is a broad agreement between models and measurements, the rapid decrease of the emissivities within 2 kpc from the disk inferred from the upper limit for the upper IV Arch seems to favor smaller values of z_max. The lower and upper IV Arch have similar altitude brackets, 0.4–1.7 and 0.7–1.7 kpc, respectively, but their emissivities are consistent with local for the first and $\lt 50\%$ of local for the latter. This is not necessarily contradictory because the distance brackets are pairs of upper and lower limits on the distances; hence the two clouds may be separated by a physical distance up to $\sim 1$ kpc.

To obtain a crude estimate of z_max from our results we fit to the emissivity scaling factors of IVCs and HVCs from Table 4 a function

$$\begin{eqnarray}&&\varepsilon (z)=1-{\left(\displaystyle \frac{z}{{z}_{\mathrm{max}}}\right)}^{\beta }\end{eqnarray} \tag{ 12 }$$

where is the emissivity expressed as a fraction of the emissivity at z = 0, z_max is the maximum height of the confinement halo and β is an index which modifies the steepness of the vertical gradient. Equation (12) with $\beta =1$ is an accurate z-dependent part of the solution of the diffusion equation in the plane-parallel geometry (infinitely thin Galactic plane) with a uniform source distribution when only ionization losses are assumed (e.g., Jones et al. 2001). The GALPROP models in Ackermann et al. (2012a) have vertical emissivity gradients at given Galactocentric radius characterized by $1\lt \beta \lt 1.5$. Values of β greater than 1 result from energy dependent escape at the radial boundary of the propagation model and, to a lesser extent, from energy losses of nuclei. In any case, we stress that the functional form in Equation (12) holds only in diffusion-dominated scenarios where it is required that the CR densities go to zero at $z={z}_{\mathrm{max}}$.

We perform the fit by using a maximum likelihood method where we neglect the different Galactocentric distances of the targets. We assume for the emissivity scaling factors a Gaussian probability distribution with σ given by the quadratic sum of statistical uncertainties, jackknife uncertainties, and the largest brackets from model input variations from Table 4. Upper limits including systematic uncertainties from Table 4 are assumed to be the 95% containment value of a Gaussian probability distribution. Finally, we assume that the probability distribution of target distances is uniform within the brackets in Table 1. Owing to the paucity of the data points, rather than fitting the index β to the data we assume the extreme values $\beta =1$ and $\beta =1.5$ as found from the fits to GALPROP predictions. The best fit to the data is obtained for ${z}_{\mathrm{max}}=(2.2\pm 1.9)$ kpc, and ${z}_{\mathrm{max}}=(1.8\pm 2.1)$ kpc for the two values of β. Note that the profile would be steeper, hence z_max smaller, if the emissivities were biased by the presence of unaccounted ionized gas in HVCs/IVCs.

This method provides an estimate of ${z}_{\mathrm{max}}\lt 6$ kpc ($2\sigma$ upper containment value from the fits above). Most of the constraining power on the z_max parameter, however, at present comes from the upper limit on the emissivity of the upper IV Arch. The latter hints at a possible tension with indications from the modeling of other observables related to CRs. Differences with respect to studies of synchrotron emission (Orlando & Strong 2013) can be understood if the densities of CR electrons beyond the boundaries of the confinement halo are sufficient to produce the emission observed in the domain from radio to microwaves. However, the interpretation of the flat radial profile of γ-ray emission produced by interactions of CR nuclei in the outer disk of the Milky Way in terms of a large z_max of the order of ∼10 kpc in the context of GALPROP models similar to those considered here (Abdo et al. 2010; Ackermann et al. 2011) seems to be disfavored.

The observed emissivities of the HVC and IVCs were compared to predictions from the limited set of diffusive reacceleration models for a simplified axisymmetric model of the Galaxy in Ackermann et al. (2012a). Other models proposed in the literature include CR-driven Galactic winds and anisotropic diffusion (e.g., Breitschwerdt et al. 2002), and a non-uniform diffusion coefficient that increases with the Galactocentric radius and the distance from the Galactic plane (e.g., Shibata et al. 2007), as well as more sophisticated ways of modeling the Milky Way structure (e.g., Werner et al. 2015), but further comparison to theoretical models is beyond the scope of this article.

Finally, we note that the CR propagation equations are often solved assuming that the CR particle number densities go to zero at the boundaries of the propagation volume. This assumption has little effect on the predicted CR fluxes for the Galaxy as a whole and at a considerable distance from the boundaries. However, close to the halo boundaries one may expect significant deviations from the model predictions and tracing the CR distribution in the halo can provide unique model-independent information about its structure and extent. Studies of the CR outflow from the Milky Way into the intergalactic space that are still in their infancy (e.g., Everett et al. 2012) could also benefit from the method illustrated in this work.

In this section we take as an example ROI A to illustrate some details of the iterative fitting procedure applied to the dust maps in Section 4.2. We follow the fit of the radiance map, R, for the assumption that the errors are proportional to the model M from Equation (4) in calculating the ${\chi }^{2}$. The results from all iterations are shown in Figures 17–19.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 17 and 18 show how the scaling coefficients of the ${\rm{H}}\;{\scriptsize{\rm I}}$ and CO maps change significantly during the iterative procedure and how they stabilize at a final value at the end. This can be interpreted as removing the bias due to the missing DNM component in the initial iteration, which is alleviated as the missing component is recovered based on the data themselves.

Figure 20 shows the difference between the filtered dust residual maps from the last iteration and the initial one. It is evident that borders exist between the ${\rm{H}}\;{\scriptsize{\rm I}}$-dominated and DNM-dominated, and between the DNM-dominated and CO-dominated regions, where the missing component is not recovered from the residuals in the initial iteration. This can be understood because the ${\rm{H}}\;{\scriptsize{\rm I}}$ and CO components in the initial fit compensating for the missing DNM component.

Zoom In Zoom Out Reset image size

Finally, we note that while in the initial iterations of the fits there are often both positive and negative significant residuals, only positive ones are found at the end of the iterative procedure. This is consistent with our working hypothesis that the missing DNM component traced by the dust maps is recovered by the template determined from the iterative procedure.

The final fit coefficients, reduced ${\chi }^{2}$ values, and error scaling factors for all the fits performed are given in Tables 6–11. In ROI B there are no significant positive residuals when fitting ${\tau }_{353}$, so the iterative procedure was not performed. The number of iterations varies depending on the region and the number of free parameters considered, ranging from at least five in ROI B up to at most 30 in ROI A.

Table 6.  Results from R Fit in ROI A

	$\sigma \propto R$	$\sigma \propto M$
σ Scaling factor	15%	12%
Number of iterations	12	11
${\chi }^{2}$/n.d.f.	1.006	1.000
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	2.160 ± 0.008	2.025 ± 0.007
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	3.19 ± 0.04	2.41 ± 0.03
${y}_{{\rm{H}}\;{\rm{I}},3}$ a	0.00 ± 0.01	0.000 ± 0.006
yCOb	3.84 ± 0.18	4.48 ± 0.10
yres	1.001 ± 0.002	1.000 ± 0.007
Risoc	0.731 ± 0.019	0.811 ± 0.016

Note.

^a10⁻³² W sr⁻¹ H⁻¹. ^b10⁻¹² W cm⁻² sr⁻¹ (K km s⁻¹)⁻¹. ^c10⁻¹⁰ W m⁻² sr⁻¹.

Download table as:  ASCIITypeset image

Table 7.  Results from ${\tau }_{353}$ Fit in ROI A

	$\sigma \propto {\sigma }_{\tau }$	$\sigma \propto {\tau }_{353}$	$\sigma \propto M$
σ Scaling factor	4.9	23%	24%
Number of iterations	11	9	5
${\chi }^{2}$/n.d.f.	0.993	1.003	1.000
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	0.738 ± 0.002	0.728 ± 0.002	0.806 ± 0.003
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	0.879 ± 0.013	0.859 ± 0.016	0.899 ± 0.021
${y}_{{\rm{H}}\;{\rm{I}},3}$ a	0.01 ± 0.03	0.00 ± 0.01	0.000 ± 0.006
yCOb	2.30 ± 0.08	2.31 ± 0.09	3.2 ± 0.2
yres	0.995 ± 0.014	0.997 ± 0.013	1.001 ± 0.001
${\tau }_{\mathrm{iso}}$ c	0.00 ± 0.01	0.000 ± 0.012	0.000 ± 0.005

Note.

^a10⁻²⁶ cm² H⁻¹. ^b10⁻⁶ (K km s⁻¹)⁻¹. ^c10⁻⁹ mag.

Download table as:  ASCIITypeset image

Table 8.  Results from R Fit in ROI B

	$\sigma \propto R$	$\sigma \propto M$
σ Scaling factor	14%	17%
Number of iterations	6	5
${\chi }^{2}$/n.d.f.	0.996	1.008
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	2.187 ± 0.010	2.246 ± 0.013
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	1.982 ± 0.011	2.168 ± 0.014
${y}_{{\rm{H}}\;{\rm{I}},3}$ a	0.46 ± 0.03	1.07 ± 0.04
yres	0.986 ± 0.015	0.98 ± 0.06
Risob	0.935 ± 0.012	1.039 ± 0.016

Note.

^a10⁻³² W sr⁻¹ H⁻¹. ^b10⁻¹⁰ W m⁻² sr⁻¹.

Download table as:  ASCIITypeset image

Table 9.  Results from ${\tau }_{353}$ Fit in ROI B

	$\sigma \propto {\sigma }_{\tau }$	$\sigma \propto {\tau }_{353}$	$\sigma \propto M$
σ Scaling factor	7.0	36%	32%
Number of iterations	5	6	8
${\chi }^{2}$/n.d.f.	0.994	1.003	0.994
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	0.729 ± 0.005	0.680 ± 0.003	0.708 ± 0.006
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	0.390 ± 0.005	0.351 ± 0.003	0.401 ± 0.007
${y}_{{\rm{H}}\;{\rm{I}},3}$ a	0.146 ± 0.012	0.152 ± 0.012	0.221 ± 0.018
yres	1.03 ± 0.02	0.98 ± 0.02	1.07 ± 0.06
${\tau }_{\mathrm{iso}}$ b	0.50 ± 0.06	0.00 ± 0.01	2.60 ± 0.08

Note.

^a10⁻²⁶ cm² H⁻¹. ^b10⁻⁹ mag.

Download table as:  ASCIITypeset image

Table 10.  Results from R Fit in ROI C

	$\sigma \propto R$	$\sigma \propto M$
σ Scaling factor	11%	13%
Number of iterations	6	7
${\chi }^{2}$/n.d.f.	1.058	1.027
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	2.045 ± 0.011	2.259 ± 0.013
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	2.124 ± 0.019	2.32 ± 0.02
yres	1.000 ± 0.014	1.00 ± 0.04
Risob	1.38 ± 0.02	1.25 ± 0.02

Note.

^a10⁻³² W sr⁻¹ H⁻¹. ^b10⁻¹⁰ W m⁻² sr⁻¹.

Download table as:  ASCIITypeset image

Table 11.  Results from ${\tau }_{353}$ Fit in ROI C

	$\sigma \propto {\sigma }_{\tau }$	$\sigma \propto {\tau }_{353}$	$\sigma \propto M$
σ Scaling factor	5.4	23%	26%
Number of iterations	8	7	6
${\chi }^{2}$/n.d.f.	1.018	1.018	0.994
${y}_{{\rm{H}}\;{\rm{I}},1}$ a	0.745 ± 0.003	0.682 ± 0.003	0.865 ± 0.004
${y}_{{\rm{H}}\;{\rm{I}},2}$ a	0.320 ± 0.006	0.300 ± 0.005	0.377 ± 0.007
yres	1.001 ± 0.004	0.992 ± 0.011	1.033 ± 0.039
${\tau }_{\mathrm{iso}}$ b	0.000 ± 0.005	0.000 ± 0.003	0.000 ± 0.006

Note.

^a10⁻²⁶ cm² H⁻¹. ^b10⁻⁹ mag.

Download table as:  ASCIITypeset image

Table 12.  Local ISM Properties from the γ-ray Analysis

	ROI A	ROI C
${X}_{\mathrm{CO}}$ a	$1.13\pm {0.08}_{-0.73}^{+0.20}\pm 0.10$	⋯
$\bar{R/N({\rm{H}})}$ b	$4.1\pm {0.3}_{-1.0}^{+0}\pm 0.3$	$3.2\pm {0.3}_{-1.0}^{+0.7}\pm 0.3$
$\bar{{\tau }_{353}/N({\rm{H}})}$ c	$2.06\pm {0.06}_{-0.65}^{+0.11}$	⋯

Note. No CO emission is detected in ROI C. The γ-ray scaling coefficients for the DNM maps derived from ${\tau }_{353}$ in ROI C were always consistent with 0.

The statistical uncertainties take into account the covariances of the fit parameters.

^a10²⁰ cm⁻² (K km s⁻¹)⁻¹. ^b10⁻³² W sr⁻¹ H⁻¹. ^c10⁻²⁶ cm² H⁻¹.

Download table as:  ASCIITypeset image

The scaling factors for the errors on ${\tau }_{353}$ given by the Planck Collaboration XI (2014) are ∼5. Possible reasons for this are discussed in Section 4.2. The scaling factors obtained when assuming $\sigma \propto R$ are 10%–15%; when assuming $\sigma \propto {\tau }_{353}$ they are $\sim 25\%$. It can be expected that the dispersion on R is smaller because this is an integrated quantity rather than referring to a specific frequency.

We note that the ${y}_{{\rm{H}}\;{\rm{I}},{\imath }}$ coefficients can have a physical interpretation in terms of average specific power, i.e., power radiated by dust per H atom, per unit solid angle, $\bar{R/N({\rm{H}})}$, and average specific opacity, i.e., optical depth of the dust per H atom, $\bar{{\tau }_{353}/N({\rm{H}})}$. Non-zero coefficients for IVCs are consistent with previous measurements of dust thermal emission from this class of objects (e.g., van Woerden et al. 2005, Chapter 9, and references therein), and, for ROIs A and B, with the measurements of metallicities close to solar (Table 1). On the other hand, the ${y}_{{\rm{H}}\;{\rm{I}},3}$ coefficient for Complex A being compatible with 0 in Tables 6 and 7 is consistent with the low metallicity measured in this cloud (Wakker 2001).

In Figure 21 we show the relative difference between different determinations of the DNM column density templates for ROI A, namely the difference between the baseline determination of the DNM column densities based on R assuming $\sigma \propto M$ (see Equation (4) for a definition of the quantities), and from two alternative determinations based on R assuming $\sigma \propto R$ (left), and based on ${\tau }_{353}$ assuming $\sigma \propto M$ (right). For the second case we calculated the average $R/{\tau }_{353}$ ratio in non-empty pixels of the two maps and used it to convert ${\tau }_{353}$ in equivalent R.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The figure illustrates how the different assumptions influence the final DNM template, for example how some regions have DNM only for one of the cases and not the other, and how the contrast within the map can change. For these reasons we considered the different determinations of the DNM template in the assessment of the systematic uncertainties for the γ-ray analysis in Section 5.4.

The results in Section 5 can constrain some properties of the local ISM seen in the foreground of the HVCs and IVCs, in particular the ${X}_{\mathrm{CO}}=N({{\rm{H}}}_{2})/{W}_{\mathrm{CO}}$ ratio in the Ursa Major molecular clouds, and the average specific power radiated by dust per H atom and unit solid angle, $\bar{R/N({\rm{H}})}$, and the average specific opacity per H atom, $\bar{{\tau }_{353}/N({\rm{H}})}$ in the DNM for the other regions.

These quantities can be extracted from the parameters in Equation (6) under the assumption that the same CR fluxes illuminate the atomic gas traced by ${\rm{H}}\;{\scriptsize{\rm I}}$ and the molecular gas traced by CO or the gas in the DNM. Then, ${X}_{\mathrm{CO}}={x}_{\mathrm{CO}}/2{x}_{{\rm{H}}\;{\rm{I}}}$, and $\bar{R/N({\rm{H}})}={x}_{{\rm{H}}\;{\rm{I}}}/{x}_{\mathrm{DNM}}$ (when the DNM is traced using R), or $\bar{{\tau }_{353}/N({\rm{H}})}={x}_{{\rm{H}}\;{\rm{I}}}/{x}_{\mathrm{DNM}}$ (when the DNM is traced using ${\tau }_{353}$). In ROI B we do not detect any significant γ-ray emission associated with the sparsely populated DNM templates, therefore we exclude it from the following discussion.

First, we check the assumption that the spectrum of CRs illuminating the different phases in the local gas complexes is the same in each ROI. Figure 22 shows the x_CO and x_DNM coefficients as a function of energy from the analysis in energy bins in Section 5.3. Within the large statistical uncertainties the spectra of the emissivities in each ROI are consistent with those of ${\rm{H}}\;{\scriptsize{\rm I}}$ in Figure 11, all of them being consistent with the spectrum of gas in the local interstellar medium from Casandjian (2015). There is no clear analog of the softening in Figure 11 for ROI A for either CO-traced gas or the DNM. However, the softening may be present and hidden by the large statistical uncertainties.

Given the fact that the spectra are consistent within the statistical uncertainties, we extract from the fit coefficients the properties of the local ISM discussed above. We take into account the systematic uncertainties related to inputs to the interstellar emission model using the results of the analysis in Section 5.4.1 and from the jackknife tests presented in Section 5.4.2. The ISM properties are reported in Table 12. They are similar to those inferred from LAT data for other nearby interstellar complexes (e.g., Planck and Fermi Collaborations 2014), demonstrating that the modeling of the foregrounds in our analysis is robust. The ${X}_{\mathrm{CO}}$ ratio in the Ursa Major molecular clouds that we obtain is consistent with other estimates from alternative methods in the literature (e.g., de Vries et al. 1987).