Fermi-LAT Observations of High- and Intermediate-Velocity Clouds: Tracing Cosmic Rays in the Halo of the Milky Way
It is widely accepted that cosmic rays (CRs) up to at least PeV energies are Galactic in origin. Accelerated particles are injected into the interstellar medium where they propagate to the farthest reaches of the Milky Way, including a surrounding halo. The composition of CRs coming to the solar system can be measured directly and has been used to infer the details of CR propagation that are extrapolated to the whole Galaxy. In contrast, indirect methods, such as observations of -ray emission from CR interactions with interstellar gas, have been employed to directly probe the CR densities in distant locations throughout the Galactic plane. In this article we use 73 months of data from the Fermi Large Area Telescope in the energy range between 300 MeV and 10 GeV to search for -ray emission produced by CR interactions in several high- and intermediate-velocity clouds located at up to kpc above the Galactic plane. We achieve the first detection of intermediate-velocity clouds in rays and set upper limits on the emission from the remaining targets, thereby tracing the distribution of CR nuclei in the halo for the first time. We find that the -ray emissivity per H atom decreases with increasing distance from the plane at 97.5% confidence level. This corroborates the notion that CRs at the relevant energies originate in the Galactic disk. The emissivity of the upper intermediate-velocity Arch hints at a 50% decline of CR densities within 2 kpc from the plane. We compare our results to predictions of CR propagation models.
Cosmic rays (CRs) are a wide-spread phenomenon, pervading galaxies and intergalactic space. They are directly detected within the solar system with spacecraft, balloon-borne, and ground-based instruments. Other methods to study CRs include observations of radio synchrotron emission and rays produced in interactions of CRs with interstellar gas, radiation and magnetic fields. CRs consist of all known stable and long-lived isotopes ranging from the most abundant hydrogen and helium nuclei through the traces of very rare Actinides, and particles such as electrons, positrons, and antiprotons. The energy density of CRs in the interstellar medium (ISM) is comparable to those of the Galactic interstellar radiation field, magnetic field, and turbulent motions of the interstellar gas. CRs are one of the essential components of the Galactic ecosystem (e.g., Ferrière, 2001), affecting the thermal, chemical and magnetohydrodynamical state of the ISM (e.g., Ptuskin et al., 2006; Indriolo et al., 2009). Emission from CRs interacting within the Milky Way is also a source of foregrounds/backgrounds for many observations spanning from radio to rays.
Observations of X-ray and -ray emission from Galactic objects such as supernova remnants, pulsars and massive-star clusters reveal the presence of energetic particles, suggesting that efficient acceleration processes occur in their vicinities (for recent reviews see, e.g., Reynolds, 2008; Bykov, 2014). The Galactic origin of CRs is also supported by -ray observations of external galaxies, which indicate that the CR densities are not universal (e.g., Sreekumar et al., 1993) and are correlated with the star-formation properties of the galaxies (e.g., Ackermann et al., 2012b). Also, the radial sizes of radio halos measured in external, edge-on spiral galaxies appear to be correlated with active star formation in their disks (Dahlem et al., 1995).
CR particles accelerated in sources in the Galactic disk are injected into the ISM where they propagate through a surrounding halo before escaping into the intergalactic space (e.g., Strong et al., 2010). This paradigm is substantiated in models that can reproduce reasonably well all the related observables (for a recent review see, e.g., Strong et al., 2007). The propagation of particles in the ISM leads to the destruction of primary nuclei via spallation and gives rise to secondary nuclei and isotopes which are rare in nature, and other particles. The observed abundances of stable secondary CR nuclei (e.g., boron) and radioactive isotopes (Be, Al, Cl, Mn) allow the separate determination of the confinement halo size and diffusion coefficient (e.g., Strong & Moskalenko, 1998) which is somewhat model dependent. Typical constraints on the halo height arising from such method are 4–6 kpc (e.g., Moskalenko et al., 2001; Trotta et al., 2011).
Observations of edge-on external galaxies show directly the existence of radio halos extending to kpc distances perpendicularly to their disks (e.g., Dahlem et al., 1994), and sometimes up to kpc (Irwin et al., 1999). Radio observations of large-scale synchrotron emission from the Milky Way were also interpreted in the context of diffusion models as indicating a vertical extent of the halo of kpc to kpc (e.g., Strong et al., 2000), with values of kpc favored based on a recent analysis including WMAP and 408 MHz observations (Orlando & Strong, 2013). The discrepancies between halo sizes determined with different methods are not necessarily surprising. The method based on the direct detection of CR nuclei provides the size of the effective volume filled with CRs that are still confined to the Milky Way. In contrast, observations of synchrotron emission show larger halos because part of the emission may be generated by particles that are leaking into the intergalactic space, i.e., that are already beyond the “point of no return”.
The Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope mission (Atwood et al., 2009) has been surveying the sky since 2008 in a wide energy range of 20 MeV to GeV. These observations provide information about CR particles (mostly protons and electrons) in distant locations that can be used to study CR propagation (e.g., Ackermann et al., 2012a). Measurements with the LAT confirmed earlier observations (e.g., Strong et al., 1988) that the radial gradient of the -ray emissivity (-ray emission rate per hydrogen atom) in the ISM outside the solar circle is much smaller than that in the density of putative CR sources (Abdo et al., 2010; Ackermann et al., 2011). This can be interpreted as another hint of a kpc CR confinement halo, because it facilitates radial diffusion of CRs and compensates for the paucity of CR sources in the outer Galaxy.
In this work we propose a more direct method to constrain the propagation of CRs in the halo of the Milky Way with LAT data. Although the interstellar gas in the Milky Way is largely confined to the Galactic plane with a scale height for atomic hydrogen of pc across much of the disk (Kalberla & Kerp, 2009), certain populations of interstellar clouds are known to exist at much greater distances from the plane. So-called high-velocity clouds (HVCs) and intermediate-velocity clouds (IVCs) are interstellar clouds of atomic hydrogen, , with line-of-sight velocities inconsistent with Galactic rotation (see, e.g., van Woerden et al., 2005). The demarcation between HVCs and IVCs is conventionally taken to be Doppler shifted velocity with absolute value between km s and km s with respect to the local standard of rest.
The distances and metallicities of HVCs/IVCs can be inferred from spectroscopic observations of foreground and background stars. HVCs typically have metallicities less than Solar and are generally considered to be remnants of the formation of the Milky Way, having existed in the halo of the Galaxy since the disk formed, or old tidal streams extracted from nearby dwarf galaxies (Wakker, 2002). Some IVCs have metallicities much closer to Solar and these may be gas ejected from the plane, e.g., by Galactic fountains (Wakker, 2002). No molecular gas, e.g., CO emission, has been detected in HVCs (Wakker et al., 1997), and IVCs are also typically free of CO emission as well.
These clouds have low column densities, cm, but can subtend hundreds of square degrees. The exposure of the LAT sky survey is deep enough for HVCs/IVCs to conceivably be detected in diffuse -ray emission and used to estimate the CR densities in the Galactic halo. These clouds are not actively forming stars and are believed to be free of internal sources of CRs and hence any -ray emission from them samples the large-scale distribution of Galactic CRs.
The article is structured as follows. In Section 2 we describe our criteria for selecting a sample of HVCs/IVCs for study. In Section 3 we discuss the -ray and multiwavelength survey data we use. Section 4 provides the details of how we derive the maps, in particular how we fit spectral components with different line-of-sight velocities to achieve kinematic separation of interstellar gas by distance along the line of sight, and other maps of ISM tracers to accurately account for foreground gas in the local neighborhood. In Section 5 we introduce the -ray analysis procedure to measure or constrain the CR densities in the targets, and the procedure we implemented to estimate systematic uncertainties. In Section 6 we discuss the results for the individual regions that we studied and the implications for the propagation of CRs in the halo of the Milky Way. The conclusions are in Section 7. Appendix A provides more technical information about the derivation of the maps of the ISM, including an iterative approach to unbiasing the maps of the column density in the dark neutral medium (DNM). Appendix B describes how we computed uncertainties and upper limits with the jackknife method for evaluating the systematic uncertainties.
Wakker (2001) reports the most complete census to date of IVCs and HVCs, including lower and upper limits on their distances (which we refer to as distance brackets), from spectrophotometric distances of foreground and background stars. In order to be able to constrain the CR densities in well-defined regions of the Milky Way halo, we select as targets for our analysis clouds in Wakker (2001) for which distance (altitude) brackets have been determined. We summarize the selected targets in Table 1, grouped in three regions of interest (ROIs) as used for the analysis later. In each case the ROIs enclose the parts of the clouds associated with the foreground and background stars used for determining the distance brackets.
|Center||Pixels||(kpc)||(kpc)||(relative to solar)||(km s)|
|A||Low-Latitude IV Arch||0.9–1.8||0.6–1.2|
|B||Lower IV Arch||0.4–1.9||0.4–1.7|
|Upper IV Arch||0.8–1.8||0.7–1.7|
Note. – The ROIs are defined as rectangular regions in a plate carrée projection in Galactic coordinates, centered at each ROI center position and with the number of pixels of width given above. See text for details. Note. – The brackets in and , and the metallicities are taken from Wakker (2001). The metallicity of the IV Spur is unknown. Note. – The ranges in velocities with respect to the local standard of rest, , are the preliminary boundaries that we use in Section 4.1 to construct the maps.
Region A encompasses an IVC known as the low-latitude intermediate velocity (IV) Arch, and an HVC known as Complex A. The low-latitude IV Arch has a mass of at a distance of kpc from the Sun, corresponding to a distance above the plane of kpc. According to Wakker (2001) the distance (comparable to that of the Perseus arm of the Galaxy), the metallicity (close to solar), and negative velocities in the range of km s to km s (higher than the Perseus arm itself in this direction) are all compatible with the characteristics of a high-interarm cloud that is part of the return flow of gas injected from the disk into the halo by a fountain process. Complex A is a distant HVC with the narrowest distance bracket111HVC Complex C has a broader distance bracket of kpc (Wakker et al., 2007). We leave it for future studies. of kpc ( kpc), which implies a mass of . Its low metallicity suggests that this HVC is a relic of the ancient Universe, a clump of gas that never underwent star formation and is now falling into the Milky Way.
Region B encompasses two IVCs that are together known as the IV Arch. The gas in the IV Arch is seen in two distinct velocity ranges: km s (the lower IV Arch) and km s (the upper IV Arch222The nomenclature refers to the absolute value of the velocity.). We consider in this analysis only the portion of the IV Arch seen at Galactic longitudes , for which there is a distance bracket for the lower part of the Arch of kpc ( kpc). At lower longitudes the distance is less reliably determined, and part of the lower IV Arch may be associated with a CO core with distance kpc from the solar system (Wakker, 2001). On the other hand, the upper IV Arch has a distance bracket of kpc ( kpc). Therefore, whether the two velocity components of the IV Arch are also physically separated entities is uncertain. The total mass of the IV Arch is of the order of several , and its metallicity is not well determined, but possibly close to solar (Wakker, 2001).
Region C encompasses the IV Spur, which is an extension of the IV Arch at lower velocities and higher Galactic latitudes. The distance bracket is kpc ( kpc), which implies a mass of . The metallicity is unknown.
All the target objects discussed reside in a relatively small sector of the Milky Way with and . This is in fact one of the regions richest in IVCs and HVCs, and for which Wakker (2001) could get access to stellar probes observations.
Furthermore, in recent years distance brackets have been determined for additional, lower-mass IVCs and HVCs (e.g., Wakker et al., 2008). Such clouds have ratios of mass over distance square of the order of kpc. Therefore, for normal -ray emission rates per H atom, as measured in the Galactic disk, we expect total fluxes of cm s at energies MeV. These fluxes are one order of magnitude or more below the LAT sensitivity level (Ackermann et al., 2012c). Therefore, we do not expect to be able to use LAT observations of these clouds to constrain their CR densities, and we do not include them among our targets. Conversely, the clouds selected for the analysis all have mass over distance square ratios of the order of kpc, and hence expected -ray fluxes cm s, in the range detectable by the LAT.
We have not included among our targets Complex GCP, also known as the Smith Cloud, because this HVC is seen at Galactic latitudes toward the inner Galaxy, where the diffuse -ray foregrounds from the Galactic disk are brighter and more challenging to model. We leave Complex GCP for future studies. We note that Drlica-Wagner et al. (2014) and Nichols et al. (2014) searched for -ray emission associated with Complex GCP due to exotic processes related to hypothetical dark matter particles, and did not report on any significant detections.
3.1 -Ray Data
The LAT is a pair-conversion -ray telescope for the energy range 20 MeV to 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, MET333Fermi Mission Elapsed Time, measured in seconds since 1 January 2000. 239500801 (2008 August 4), up to MET 431481603 (2014 September 4), for a total of 73 months. We use the P7REP dataset (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 dataset used:
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).
We retain only candidate rays detected at inclinations from the Earth zenith , 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.
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 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 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 v09r34p01444The Science Tools are available along with LAT data from the Fermi Science Support Center, http://fermi.gsfc.nasa.gov/ssc..
3.2 Interstellar Medium Tracer Data
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 emission line by Kalberla et al. (2005). The survey provides the line brightness temperatures over a grid of (with an effective angular resolution of ), and covers velocities with respect to the local standard of rest from km s to km s, at a resolution of 1.3 km s. Data were corrected for stray radiation reaching a residual noise contamination K for most directions in the sky. We use data from Kalberla et al. (2005) to derive 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, , 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 map on a grid. The intensities are assumed to be proportional to the molecular hydrogen column densities .
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 hundred parsecs 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.
Various observations indicate that a linear combination of and 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 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 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 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, , and the map of total dust radiance, . The angular resolution of these maps is 5′. They are used to construct maps of the DNM as described in 4.2.
4 Construction of the Interstellar Medium Tracer Maps
For each ROI we construct maps of the column densities for different complexes along the line of sight, and maps of the DNM column densities. The maps are constructed over regions that include at least 5° boundaries around the ROI in order to properly take into account the point-spread function (PSF) of the LAT in fits of the models to the -ray data.
4.1 Construction of the Maps
The column densities of atomic gas, , 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.
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).
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 associated with each component is proportional to the integral of the brightness temperature over velocity within the adjusted velocity boundaries.
For each line of sight on the 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).
In region A we set the preliminary velocity boundaries at km s to separate local gas from the low-latitude IV Arch, and at km s to separate the low-latitude IV Arch from Complex A. At latitudes 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 point-spread function.
In region B we set the preliminary velocity boundaries at km s to separate local gas from the lower IV Arch, and at km s to separate the lower from upper IV Arch. In the case of region C we set the preliminary velocity boundaries at km s to separate local gas from the IV Spur, and at 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 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 lines555Out of 11421, 34001, 32421 analyzed in region A, B, and C, respectively. of sight with differences up to . The magnitude of the correction for the spillover is on average of the uncorrected column densities. (The average is the largest for ROI B, where it is .) 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 of the total column densities. In a small subset of the lines of sight () the fits do not converge properly, mostly due to the faintness of emission in these high-latitude regions. In such cases we do not apply any correction for spillover.
We have calculated the column densities under the approximation of small optical depth. An approximate correction for 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 and maps as described in 4.2, we are effectively accounting also for densities underestimated due to the approximation of small optical depth.
We replace the column density seen in ROI A at , with the average of neighboring lines of sight owing to an anomalous temperature profile indicative of an artifact in the LAB data.
Figure 4 (a) and (b) show some level of anticorrelation between higher/lower intensity regions. This is due to some genuine anticorrelation between the intensities of emission in velocity space. At the same time, the lower IV Arch shows a bright feature at km s that blends into the local components, which makes the separation particularly challenging. This is addressed by the spillover correction.
4.2 Construction of the DNM Maps
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 and maps, have been used to trace the DNM, neutral interstellar matter not properly traced by the linear combination of and . 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; 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 dataset, 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 (i.e., the map of optical depth at 353 GHz, , or radiance, ) by assuming that either quantity is proportional to the atomic and molecular column densities as traced by and , respectively, plus an isotropic term representing any zero-level shift, for example, due to residual contamination from the cosmic infrared background.
Here and in the following equations we indicate with the dust specific power (opacity) per hydrogen atom, and with the specific power (opacity) per unit in the molecular phase. We use indices to label the complexes separated kinematically along the line of sight as described in Section 4.1. The model is fit to the dust map using the minimum method. The is defined, for a generic model , as
Following the work of the Planck and Fermi Collaborations (2014), we consider different possible assumptions on the uncertainties . 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 , we assume the uncertainties to be . We empirically consider two alternative hypotheses to assess the impact of this choice:
, which tests how the agreement between and affects the results of the fits;
, where is the map of statistical uncertainties on provided by the Planck Collaboration XI (2014) (for fits of only).
We initially assume that the proportionality coefficients are 30% for the model or the dust map itself, and 1 when using . These are adjusted later as a part of the iterative procedure to obtain a reduced . The is minimized to determine the best-fit coefficients, and .
From the results we derive a dust residual map
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 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 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 :
The model is defined:
We fit the model to the dust map minimizing the as defined in Eq. 2;
The dust residual map is built:
where and are the best-fit coefficients;
We filter to extract the significant positive residuals 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 are compatible with 0 within statistical uncertainties from the fit; 2) the difference between minimum in iteration and iteration is .
After the first pass just described, the scaling coefficients for the uncertainties are re-defined, so that the reduced in the latest iteration. Requiring the reduced 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 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 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 as the best dust column-density tracer at high latitudes, we take as our baseline DNM templates those obtained from 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 North Celestial Pole (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.
5 -Ray Analysis
5.1 Model Building
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, 2014) 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 (Sec. 5.4). For each tracer map we allow a free scaling coefficient in the model, to define the overall CR flux in that component666HVCs and IVCs span volumes much larger than those of interstellar complexes studied so far using LAT data. Were there any variations of CR densities within the clouds, the scaling coefficient would represent an average CR flux over it..
The IC component is modeled based on the work by de Palma et al. (2013). They considered 8 models of IC emission calculated using the GALPROP777http://galprop.stanford.edu 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 dataset. 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 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 of the source in 3FGL is , just the global normalization free if , and all the parameters fixed otherwise. The thresholds and are chosen for each ROI so that the number of free parameters associated with sources is to prevent convergence problems in the fitting procedure. The pairs of thresholds chosen are , , and , for regions A, B, and C, respectively. We verified that none of the 3FGL sources are positionally coincident with structures in the maps in Figures 3, 4, and 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.22914 in ROI C. We will show later in Figures 7, 8, and 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.
To summarize, for each ROI the model for the -ray intensities is described by Equation 6.
There, stands for the tabulated emissivity spectrum of in the local interstellar space from Casandjian (2014). are the column density maps that we constructed in Section 4.1, so that, for each complex 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. is the CO intensity map described in Section 3.2 (used only for ROI A); hence can be interpreted as the scaling factor that transforms the local emissivity into -ray emissivity per intensity unit in the region studied. Similarly, is one of the templates of the DNM constructed in Section 4.2, and is the scaling factor that transforms the local emissivity into -ray emissivity per dust radiance (or optical depth) unit. IC represents the IC emission model: is an overall global scaling factor, and are the parameters of the log-parabolic spectral correction, while is a reference energy which was set to 700 MeV (the choice is arbitrary because the parameter is degenerate with the others). Finally, is the tabulated isotropic background spectrum (different for front- and back-converting events), with as common scaling factor, and SRC represents the model for individual sources described above.
5.2 Analysis Procedure and Baseline Results
Using the dataset described in Section 3.1, for each ROI in Table 1 we build binned count cubes with spacing in angle and 15 logarithmically spaced energy bins from 300 MeV to 10 GeV for front- and back-converting events. The model in Eq. 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, 8, and 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.
|Parameter||ROI A||ROI B||ROI C|
|aa10 cm (K km s)|
|bb10 sr W|
Note. – See Eq. 6 for a definition of the model parameters. The numbering scheme for the coefficients is: 1) low-velocity, i.e., local gas, 2) and 3) IVCs and HVCs. Note. – Asymmetric error bands are computed using Minos when the best-fit value is at the boundary imposed during the likelihood optimization process.
For each target (HVC or IVC) in Table 1 we calculate the test statistic TS, defined as
where is the maximum likelihood for the model that includes the source of interest, while 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 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 , which formally corresponds to a 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 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 based on a Bayesian integration of the likelihood profile as applied for measurements in the 2FGL catalog (Nolan et al., 2012, 3.5).
|Region||Complex||TS||95% c.l. upper limit|
|A||Low-Latitude IV Arch||127|
|B||Lower IV Arch||750|
|Upper IV Arch||0||0.23|
Note. – See Eq. 7 for the definition of TS. Note. – The upper limits are expressed as a fraction of the emissivity of local .
In Fig. 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 maps inferred from the 21-cm observations in Figures 3, 4, and 5 are overlaid. The morphologies of the -ray detected signals follow closely the known 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 dataset.
5.3 Spectra of the Detected Clouds
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 from Casandjian (2014), 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 spectrum reproduces the observed one. To do so we perform the -ray analysis independently in 6 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 emissivity in the three regions studied.
Figure 12 shows the scaling factors for the emissivity of all the detected IVCs with respect to local.
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 Assessment of Systematic Uncertainties
5.4.1 Uncertainties of inputs to the -ray interstellar emission model
Some inputs to the model in Eq. 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 kpc 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 observations888We note that the derivation of the Galactic supernova remnant distribution in Case & Bhattacharya (1998) is subject to uncertainties and discordant with some later works (e.g., Green, 2014). In our study, though, we use it only as a way to probe the uncertainties due to the modeling of CR propagation relying on the previous work by Ackermann et al. (2012a); de Palma et al. (2013).. 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
Eq. 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, Eq. 9 captures a low-energy hardening predicted by high-reacceleration models in Ackermann et al. (2012a) at larger Galactocentric radii999 Our choice to include this effect in the evaluation of systematic uncertainties is conservative, because the low-energy hardening is a feature that appears only in models with high reacceleration, while models with convection or pure diffusion models are also viable (Ackermann et al., 2012a).. 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
We use these functions to modify the emissivity spectra of the undetected clouds, by multiplying the local emissivity spectrum in Eq. 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 . In this way, we can still interpret the upper limits on 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 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.
|Region||Complex||Scaling Factor||Stataa statistical uncertainty from the likelihood analysis.||Sys (Model)bbSystematic spread from varying some inputs to the -ray interstellar emission model (see Sec. 5.4.1).||Sys (Jackknife)cc uncertainty or 95% c.l. upper limit from the distributions obtained with the jackknife test (see Sec. 5.4.2).|
|Low-Latitude IV Arch||0.94||0.15||0.08|
|Lower IV Arch||1.08||0.04||0.09|
|Upper IV Arch|
Note. – Upper limits are provided at 95% confidence level.
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 () and large enough to probe possible mild variations of the CR densities across the complexes studied.
Figures 13, 14, and 15 show the distributions of the emissivity scaling coefficients for the three ROIs. The distributions all peak within from the best-fit values obtained from the baseline analysis and have standard deviations smaller than or comparable to the 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.
In ROI A, Figure 13, the distribution of the scaling coefficients for Complex A, , 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 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, 14, and 15, the 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.
6 Results and Discussion
6.1 The -Ray Emissivities of HVCs and IVCs
The 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 (2014), i.e., 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 Eq. 8, 9, 10, 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 (2014) integrated between 300 MeV and 10 GeV is 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.
|Region||Complex||-ray emissivity||Stataa statistical uncertainty from the likelihood analysis.||Sys (LAT)bbError from uncertainties in the LAT effective area.||Sys (model)ccSystematic spread from varying some inputs to the -ray interstellar emission model (see Sec. 5.4).||Sys (Jackknife)dd uncertainty or 95% c.l. upper limit from the distributions obtained with the jackknife test (see Sec. 5.4.2).|
|(kpc)||( s sr H)|
|A||Low-Latitude IV Arch||0.6–1.2||6.7||1.1||0.6||0.6|
|B||Lower IV Arch||0.4–1.7||7.7||0.3||0.7||0.6|
|Upper IV Arch||0.7–1.7|
Note. – brackets are replicated from Table 1 for the reader’s convenience. Data sources are listed there. Note. – -ray emissivities are integrated between 300 MeV and 10 GeV. Note. – 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.
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 GeV/nucleon up to 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 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 to 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 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.
6.2 Testing the Origin of Cosmic Rays in the Galactic Disk
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 from the disk itself101010Owing to the distances to other massive star-forming galaxies, any contributions from CRs escaping from even the closest systems would be too small to affect this conclusion and also smaller that the uncertainty of our emissivity measurements.. 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 the mean within the bracket for targets in Table 1, and kpc for the local gas in each ROI, we obtain . The negative value indicates, indeed, a negative correlation of emissivity with .
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 . We generate null hypothesis datasets starting from the real dataset using the following procedure:
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;
we add upper limits at the values measured in Table 4 for targets which are not significantly detected;
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);
we draw a set of values with uniform probability distribution within the brackets in Table 1, and add three elements with to represent the local foreground gas in each ROI;
we randomly shuffle in an independent fashion the emissivity scaling factors and sets so obtained to be consistent with the null hypothesis of no correlation between the two quantities.
The distribution of 30000 null hypothesis datasets shows that there is a chance probability to obtain a coefficient .
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 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 to of the total emissivity per H atom (taking into account He in both CRs and the ISM), lower than the decrease of emissivity of for the Upper IV Arch and of 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 is more likely due to diffusion in the halo. A quantitative comparison with propagation models follows in the next section.
6.3 Comparison with CR Propagation Models
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 and maximum height from the disk where the CR particle number densities are required to go to zero.
The model predictions are calculated for distances along the line of sight approximately corresponding to the peak gas column densities for each ROI111111This accounts for the different ranges of Galactocentric radius in each line of sight. The effect is clearly visible for ROI A in Figure 16 where the profiles are concave, and the curves for kpc bifurcate into two families for kpc and kpc due to particle escape at the radial boundary.. 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 from the measurement in Casandjian (2014). 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 . 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 with increasing (Mori, 2009).
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 . Models in Ackermann et al. (2012a) assume values of between 4 kpc 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 . The lower and upper IV Arch have similar altitude brackets, kpc and kpc, respectively, but their emissivities are consistent with local for the first and 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 kpc.
To obtain a crude estimate of from our results we fit to the emissivity scaling factors of IVCs and HVCs from Table 4 a function
where is the emissivity expressed as a fraction of the emissivity at , is the maximum height of the confinement halo and is an index which modifies the steepness of the vertical gradient. Eq. 12 with is an accurate -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 . 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 Eq. 12 holds only in diffusion-dominated scenarios where it is required that the CR densities go to zero at .
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 and as found from the fits to GALPROP predictions. The best fit to the data is obtained for kpc, and kpc for the two values of . Note that the profile would be steeper, hence smaller, if the emissivities were biased by the presence of unaccounted ionized gas in HVCs/IVCs.
This method provides an estimate of kpc ( upper containment value from the fits above). Most of the constraining power on the 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 of the order of 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.
We have searched for high-energy -ray emission from a sample of HVCs and IVCs in the halo of the Milky Way at different distances from the Galactic disk using 73 months of data from the LAT in the energy range between 300 MeV and 10 GeV.
We have achieved the first detections of IVCs in rays, the low-latitude IV Arch, the lower IV Arch, and the IV Spur, and set constraints on the emission from the remaining targets, the upper IV Arch and the HVC Complex A (Table 4). The spectra of the detected complexes were all consistent within statistical uncertainties with that of gas in the local interstellar space.
We find evidence at 97.5% confidence level that the -ray emissivity per H atom of the clouds decreases as a function of distance from the disk. This corroborates the notion that CRs are accelerated in the Galactic disk and then propagate in a surrounding halo, for the first time in a direct way.
The -ray emissivity per H atom of the upper IV Arch hints at a 50% decline of the CR densities within 2 kpc from the disk. The upper limit on the emissivity of Complex A at 22% of the local value gives the most stringent constraint to date on the fluxes of CRs at few kpc from the Milky Way disk.
Our results can be compared to CR propagation models and inform further development to more realistically take into account the escape of CRs from the halo of the Milky Way and the outflow of particles merging into the local intergalactic medium. At the same time, the observational constraints can be improved if distance brackets for more HVCs and IVCs in the mass-distance range detectable by the LAT are measured or existing brackets are made more constraining, e.g., by the ongoing Gaia survey (e.g., de Bruijne, 2012).
Appendix A Additional Results on the Interstellar Medium
The analysis described in this paper is aimed at determining the CR content of IVCs and HVCs. Nevertheless, we obtained results relevant to the properties of interstellar gas and dust in such clouds, as well as in the local complexes seen in their foregrounds. This appendix summarizes these results, including some methodological aspects.
a.1 Example of Results from the Iterative Dust Fitting Procedure
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, , for the assumption that the errors are proportional to the model from Eq. 4 in calculating the . The results from all iterations are shown in Figures 17, 18, and 19.
Figures 17 and 18 show how the scaling coefficients of the 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 -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 and CO components in the initial fit compensating for the missing DNM component.
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.
a.2 Fitting Coefficients for the Dust Maps
The final fit coefficients, reduced values, and error scaling factors for all the fits performed are given in Tables 6, 7, 8, 9, 10, and 11. In ROI B there are no significant positive residuals when fitting , 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.
|Number of iterations||12||11|
|aa10 W sr H|
|aa10 W sr H|
|aa10 W sr H|
|bb10 W cm sr (K km s)|
|cc10 W m sr|
|Number of iterations||11||9||5|
|aa10 cm H|
|aa10 cm H|
|aa10 cm H|
|bb10 (K km s)|
|Number of iterations||6||5|
|aa10 W sr H|
|aa10 W sr H|
|aa10 W sr H|
|bb10 W m sr|
|Number of iterations||5||6||8|
|aa10 cm H|
|aa10 cm H|
|aa10 cm H|
|Number of iterations||6||7|
|aa10 W sr H|
|aa10 W sr H|
|bb10 W m sr|
|Number of iterations||8||7||6|
|aa10 cm H|
|aa10 cm H|
The scaling factors for the errors on given by the Planck Collaboration XI (2014) are . Possible reasons for this are discussed in Section 4.2. The scaling factors obtained when assuming are 10% to 15%; when assuming they are . It can be expected that the dispersion on is smaller because this is an integrated quantity rather than referring to a specific frequency.
We note that the 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, , and average specific opacity, i.e., optical depth of the dust per H atom, . 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 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).
a.3 Differences between Alternative Determinations of the DNM Column Density Templates
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 assuming (see Eq. 4 for a definition of the quantities), and from two alternative determinations based on assuming (left), and based on assuming (right). For the second case we calculated the average ratio in non-empty pixels of the two maps and used it to convert in equivalent .
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.
a.4 , Dust Specific Opacity and Power in the Local DNM from the -Ray Analysis
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 ratio in the Ursa Major molecular clouds, and the average specific power radiated by dust per H atom and unit solid angle, , and the average specific opacity per H atom, in the DNM for the other regions.
These quantities can be extracted from the parameters in Eq. 6 under the assumption that the same CR fluxes illuminate the atomic gas traced by and the molecular gas traced by CO or the gas in the DNM. Then, , and (when the DNM is traced using ), or (when the DNM is traced using ). 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 and 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 in Figure 11, all of them being consistent with the spectrum of gas in the local interstellar medium from Casandjian (2014). 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 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).
|ROI A||ROI C|
|aa10 cm (K km s)|
|bb10 W sr H|
|cc10 cm H|
Note. – No CO emission is detected in ROI C. The -ray scaling coefficients for the DNM maps derived from in ROI C were always consistent with 0. Note. – The statistical uncertainties take into account the covariances of the fit parameters.
Appendix B Evaluation of Errors and Upper Limits With the Jackknife Method
In this appendix we summarize the method used to evaluate uncertainties and upper limits from the jackknife tests in Section 5.4.2. Following Dudewicz & Mishra (1988), given the jackknife estimates with let
be their arithmetic mean. A quantity known as is defined as
If we define
the uncertainty of is