IS THE ULTRA-HIGH ENERGY COSMIC-RAY EXCESS OBSERVED BY THE TELESCOPE ARRAY CORRELATED WITH ICECUBE NEUTRINOS?

The Telescope Array (TA) experiment is the largest hybrid detector in the Northern Hemisphere advancing the effort to understand the sources of ultra-high-energy cosmic rays (UHECRs). The TA Collaboration recently reported an analysis of five years of data measured by the surface detectors in which they find a cosmic-ray hotspot with 5.1σ significance within a 20° radius circle centered at R.A. = 146.7, decl. = 43.2 (a post-trial chance probability of 3.4σ; Abbasi et al. 2014). If this signal is real, it may signal the first detection of a UHECR point source in the universe.

Additionally, the IceCube Collaboration recently reported the first detections of extraterrestrial neutrinos with energies above 30 TeV (IceCube Collaboration 2013). These 28 events are distributed across the sky and are currently consistent with an isotropic distribution, although an excess of marginal statistical significance is detected in the Southern Hemisphere. The effective area of the southern sky is much larger for the IceCube detector, due to the fact that the Earth is opaque to very-high-energy neutrinos. The effective area of the southern sky is roughly double that of the north at an energy of 100 TeV, and roughly three times as large at an energy of 1 PeV. Additionally there are significant atmospheric backgrounds in the Southern Hemisphere due to muon showers, which are not present in the Northern Hemisphere. Due to these facts, 24 of the 28 IceCube neutrinos have best-fit positions in the Southern Hemisphere, and the effective area of IceCube is relatively small in the region surrounding the TA excess. Interestingly, out of only four neutrinos observed throughout the Northern Hemisphere, two high-energy neutrinos are found in a similar region as the TA source (Figure 1).

Zoom In Zoom Out Reset image size

High-energy neutrinos and UHECRs are known to be tightly connected. In addition to EeV cosmogenic neutrinos generated during the UHECR propagation from the source to the Earth, TeV–PeV high-energy neutrinos can be produced when cosmic rays interact inside the source, if the source has a dense photon background (e.g., active galactic nuclei, AGNs: Stecker et al. 1991; Murase et al. 2014, gamma-ray bursts, GRBs: Murase & Ioka 2013) or hadron background (e.g., newborn pulsars, Fang et al. 2012); and (or) in the source's local environment (e.g., galaxy clusters/groups and star-forming galaxies, Kotera et al. 2009; Murase et al. 2013).

Here we test the correlation between the excess observed by TA and the two spatially coincident events observed by IceCube. We employ Monte Carlo statistics to demonstrate that the overlap between these two excesses is statistically significant at the 2σ level. We then examine the spectrum of the neutrino signal and find that while huge uncertainties are present, the signal is consistent with that observed in UHE. Finally, we note that the single source that hosts both the UHECR and neutrino excesses could be a star-forming galaxy.

The hotspot observed by TA can be modeled as an ellipse which stretches approximately from an R.A. of 125°–170° and a decl. of 20°–65°(Abbasi et al. 2014; indicated as a red circle in Figure 1). This covers an area on the sky of 0.38 sr. We can compute the probability that two IceCube neutrinos land within this target region. We note that the IceCube Collaboration published instrumental effective areas for both Northern Hemisphere and Southern Hemisphere searches. They find the effective area in Northern Hemisphere searches to fall approximately a factor of 2–3 below that of Southern Hemisphere searches, a fact consistent with the observation of 24 neutrinos in the Southern Hemisphere and only 4 in the Northern Hemisphere.

While the effective area of IceCube likely has a complicated spatial dependence, we assume that it is constant within the Northern Hemisphere. While this assumption is naive, the primary contribution to the decreasing effective area in the Northern Hemisphere is the amount of matter a neutrino must travel through before reaching the IceCube detector. Thus, the first order correction is that the Northern Hemisphere effective area will be highest along the equatorial plane, and smallest at high declination. Since the TA excess is centered around decl. = 432, the first order correction to the effective area in this region would be approximately 0, or perhaps slightly negative, making our results conservative. Given a flat effective area, we compute the probability for two IceCube events to have central positions within the TA hotspot to be 0.017, and we can thus reject the coincidence of the IceCube events with the TA excess at 2.1σ. Including a fifth neutrino (corresponding to neutrino 5, which is centered at a declination of only −04) dilutes this result somewhat to 1.9σ. We note that additional modeling of the IceCube effective area (as well as the observation of additional neutrinos) will further clarify or rule out this enticing coincidence.

Furthermore, we note that the statistical correlation between these IceCube neutrinos is also of interest. The IceCube Collaboration notes that the morphology of the 28 observed extraterrestrial neutrinos is consistent with isotropy, although a slight excess is observed in the Southern Hemisphere. However, this finding of isotropy does not necessarily hold for any particular set of neutrinos, once a specific region of interest (ROI) is picked out a priori. Using the IceCube Collaborations published error ellipses for each of its recorded neutrinos, the probability of two Northern Hemisphere neutrinos being located within their 1σ error ellipses is 0.21 (0.8σ), while the probability of two neutrinos being located closer than the observed offset of 13° is 0.14 (1.1σ). Though neither of these measures are statistically significant, the fact that the two neutrinos also overlap the TA hotspot is statistically interesting. Further data from the IceCube experiment may be able to soon conclusively determine whether there is an anisotropy in the IceCube neutrino flux which is consistent with the TA excess.

In addition to checking the spatial consistency of the IceCube and TA data, it is worth checking whether the spectrum and intensity of the neutrino events are consistent with the TA excess. While constraints on the spectrum are extremely loose due to small-number statistics, it is worth asking whether the intensities of each excess are consistent, given theoretical models for the spectra of each messenger. In order to calculate the neutrino intensity, we partition the neutrino events into two energy bins: one which stretches from 30 TeV to 3 PeV, and contains two events against a Northern Hemisphere astrophysical background of approximately 0 events and a diffuse background of 0.22 events, which is calculated based on the probability that an event from an isotropic direction in the northern sky happens to land inside the TA excess region. Additionally, we list a second energy bin from 3 PeV to 10 PeV, which contains no events and no background. Any flux yielding more than 2.3 events in the high energy bin is thus excluded at more than 90% confidence (Anchordoqui et al. 2014a).

Adopting the IceCube Northern Hemisphere effective area and 662 days live time of detection, we find the neutrino flux in the hotspot region to be (1.3 ± 0.9) × 10⁻⁸ GeV cm⁻² s⁻¹, as indicated by the blue cross in Figure 2 (left). The null bin poses a loose limit because IceCube has a relatively small effective area in Northern Hemisphere above 1 PeV. This upper limit can be enhanced with the IceCube point source sensitivity at the hotspot declinations, for muon neutrinos with energy from 1 PeV to 1 EeV in E⁻² spectrum, presented as the black dash line in Figure 2 (left; Aartsen et al. 2013). We note that these bins were picked in order to encapsulate the entire extraterrestrial IceCube neutrino flux into a single energy bin, and that it is difficult to calculate any spectrum from this process as the choices in binning parameters greatly affects the calculated best fit spectrum.

Zoom In Zoom Out Reset image size

Due to larger statistics, the flux of the TA excess can be determined with greater statistical significance. The TA hotspot contains 19 ± 4.49 out of the 72 TA highest energy cosmic ray events, corresponding to a flux above 57 EeV of:

$$\begin{eqnarray} E^2J_{H} & = & E^2J(E)_{4\pi }\left(\frac{N_{H}}{N_{4\pi }}\right)\left(\frac{4\pi }{\Omega _{20^\circ }}\right) \left(\frac{A_{4\pi }}{A_{H}}\right) \nonumber\\ & = & (4.4\pm 1.0)\times 10^{-8}\left(\!\frac{E}{10^{19.5}\,\rm eV}\!\right)^{-2.6} \,\rm \frac{GeV}{cm^{2}\,s\,sr} \end{eqnarray} \tag{ 1 }$$

where E³J_4π = 7.9 × 10⁻⁹ GeV cm⁻² s⁻¹ sr⁻¹ is the all-sky UHECR flux at 10^19.5 eV (Kistler et al. 2013), A_H/A_4π = 1/2 (The Telescope Array & Pierre Auger Collaboration 2013) and N_H/N_4π are the ratios of effective areas and event numbers between the hotspot region and all sky.

Over a distance D from the source to our Galaxy, UHECRs are expected to experience a deflection angle δθ_EG ≈ 3°–5° Z (D/100 Mpc)^1/2 (E/100 EeV)⁻¹ in a turbulent extragalactic magnetic field (EGMF) with field strength B_EG ∼ 2 nG and coherence length λ_EG ∼ 300 kpc (Miralda-Escude & Waxman 1996). The Galactic magnetic field (GMF) adds another δθ_G ≈ 1°–5° Z depending on the arrival direction of the events (Giacinti et al. 2010; Farrar et al. 2013). This is somewhat smaller than the observed extent of the TA hotspot region. However, crossing the intergalactic magnetic field can further enhance the deflection between the arrival direction of super-GZK protons and the sky position of their actual sources to be ∼15° (Ryu et al. 2010). Additionally, if the UHECR flux is dominated by higher-Z particles, as suggested by data from the Pierre Auger Observatory (Abraham et al. 2010; Letessier-Selvon 2014), then the deflection angle would be enhanced by a factor of Z. Thus, the large diffuse region of the TA hotspot is entirely consistent with emission from a single bright point source. Interestingly, we note that the existence of a single bright UHECR source in the Northern Hemisphere could help elevate the tension between the chemical composition measurements (Xmax and RMS-Xmax) of the TA and Auger. Specifically, while TA observed a UHECR flux compatible with 100% protons, observations by Auger indicate the presence of heavy nuclei in the UHECR flux. If this bright UHE source is a proton accelerator, it would strongly contribute to the TA Xmax measurements, compared to the chemical composition of the isotropic background. Alternatively, if the source emits both light and heavy nuclei and if the EGMF and GMF are intense, the protons from this source will be less deflected by the magnetic fields, while the heavy nuclei will be greatly deflected. This would produce an overabundance in light nuclei in the TA ROI, and an overabundance of heavy nuclei in the Auger ROI. Future gamma-ray observations of the hotspot region can help constrain the strength of the magnetic fields in order to constrain the scenarios.

If the UHECR source powering the TA hotspot is at large distance, cosmic rays above the GZK energy would lose energy via photo-pion interaction on cosmic microwave background in their extragalactic propagation, and produce cosmogenic neutrino flux up to $E_\nu ^2\Phi _\nu =3/8\,E^2J(E)_{\rm H}\,f_\pi =(1.0\pm 0.2)\times 10^{-9}\,\rm GeV\,cm^{-2}\,s^{-1}$, assuming on average each pion takes 20% proton energy and a maximum pion production rate f_π = 1. As the red mark in Figure 2 (right) shows, this flux level is well below the high energy IceCube upper limits, suggesting that it is consistent with both UHECRs and neutrinos emitted by the same source.

In Figure 2 (right), on the other hand, we attempt to constrain the spectrum of the neutrino source in a bin-independent way. It is worth noting that the neutrino spectrum can be constrained not only from the two observed events at 63 TeV and 210 TeV, but also by the lack of detections at any other energy. We adopt a flat prior on the spectral index and then produce a Monte Carlo population of neutrinos for each spectral index. Weighting the spectrum by the energy dependent effective area of the IceCube experiment in the Northern Hemisphere (IceCube Collaboration 2013), we calculate the relative likelihood that emission spectra with varying spectral indices would produce the observed data. Given the existence of only two events, the constraints on the spectral index are extremely loose, with value of 3.0, and a 1σ (2σ) confidence interval of 2.35–3.95 (1.77–4.71).

Using this wide range of spectral indices, we can attempt to determine whether the neutrinos observed at TeV energies may be correlated with the UHECR observed by TA. If we assume that all cosmic rays are accelerated with the same mechanism with a single power-law, to a zero-order approximation, the neutrinos would be expected to follow the same spectral index until the cosmic-ray energy hits the pion production threshold.

In Figure 2 (left) we show blue dashed lines indicating the possible neutrino spectral indices for fits at our 1σ and 2σ confidence intervals, alongside orange-dashed lines indicating the corresponding cosmic-ray spectrum. We note that while our neutrino data are centered around a spectral index α = 3.0, cosmic-ray all-sky observations by the TA (Kistler et al. 2013) and KASCADE (Apel et al. 2013) disfavor indices softer than 3.0. A spectral index between 2.0 and 3.0 is highly consistent with neutrino and UHECR observations and produces a coherent picture for all data.

If the IceCube neutrinos and TA hotspot are indeed related, then it is worth considering possible sources of UHECRs which could produce both signals. While it is possible that these sources are produced by an over-density in extragalactic UHECR production in this direction, the high fractional intensity of the TA excess coupled with the relative smoothness of the Fermi-LAT extragalactic γ-ray sky disfavor this scenario and instead indicate emission from either a single bright UHECR source, or possibly a handful of such sources (Abdo et al. 2010; Ackermann et al. 2012a).

Any potential source of both the TA and IceCube excesses should be able to power both UHECRs and TeV–PeV neutrinos. In addition, it should not overproduce PeV–EeV neutrinos compared to the point-source limits produced by IceCube (blue dash line in Figure 3). Below we investigate a few possible scenarios in detail and summarize them in Table 1.

Zoom In Zoom Out Reset image size

Gamma-ray bursts (GRBs) have long been suggested as potential UHECR sources (see Mészáros 2006 for a review), and high-energy neutrino emission peaked around PeV can happen when cosmic rays interact with fireball photons (Waxman & Bahcall 1997; Murase & Nagataki 2006; Cholis & Hooper 2013; Liu & Wang 2013; Murase & Ioka 2013). The non-thermal photon spectra of GRBs can be described by a broken power law, n()∝^−α below the break energy _b and n()∝^−β above, with α ∼ 1 and β ∼ 2. The photopion production efficiency is

$$\begin{eqnarray} f_{p\gamma }(\epsilon) &=& 0.2\,\frac{L_{\gamma,51}}{\epsilon _{b,\rm MeV}\,\Gamma _{2.5}^4\,\delta t_{\rm ms}} \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\hspace{57pt}\times \left\lbrace \begin{array}{@{}ll}(\epsilon /\epsilon _{b})^{\beta -1} & \epsilon \le \epsilon _b \\ (\epsilon /\epsilon _{b})^{\alpha -1}& \epsilon > \epsilon _b\end{array}\right. \end{eqnarray} \tag{ 3 }$$

for a GRB with luminosity L_{γ, 51} = L_γ/10⁵¹ erg s⁻¹, bulk Lorentz factor Γ_2.5 = Γ/10^2.5 and variability time δt_ms = δt/1 ms (corresponding to R = 2Γ² c δ t ≈ 10¹³ cm). A photon break energy in observer frame at $\epsilon _{b, \rm MeV} = \epsilon _b/1\,\rm MeV$ corresponds to a break in neutrino spectrum at $E_\nu ^b = 7\times 10^{14}\, \Gamma _{2.5}^2\,\epsilon _{b,\rm MeV}^{-1}\,\rm eV$.

The neutrino spectrum of a GRB can be approximated by

$$\begin{eqnarray} E_\nu \,\Phi _\nu &=& \frac{3}{8}\,f_{p\gamma }\,E_{\rm CR}\,J_{\rm CR} \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &\propto & \left\lbrace \begin{array}{@{}ll} \left(E_\nu /E_\nu ^b \right)^{\beta +1 -s} & E_\nu \le E_\nu ^b \\ \ms \left(E_\nu /E_\nu ^b \right)^{\alpha +1 -s} & E_\nu ^b < E_\nu < E_\nu ^c \\ \ms \left(E_\nu /E_\nu ^c \right)^{\alpha -1 -s} &E_\nu > E_\nu ^c \end{array}\right. \end{eqnarray} \tag{ 5 }$$

Here $E_\nu ^c = 2\times 10^{16}\, L_{\gamma,51}^{-1/2}\,\Gamma _{2.5}^4\,\delta t_{\rm ms}\,\rm eV$. This additional break in the neutrino spectrum is caused by the synchrotron cooling of charged pions in the magnetic field of the dissipative fireball, and the ratio of synchrotron cooling time to pion decay time gives a suppression factor $\sim E_\nu ^{-2}$. Assuming an intrinsic cosmic-ray spectrum E^−2.2 normalized by the TA excess, a GRB with typical parameters is capable to produce the neutrino excess (as shown by the red dashed line in Figure 3).

Another possible channel of producing TeV–PeV neutrinos is via hadronuclear interaction (pp or Np) in baryon-rich sources, including newborn pulsars (Fang et al. 2012, 2013) and magnetars (Arons 2003; Murase et al. 2009), galaxy clusters (Inoue et al. 2005; Kotera et al. 2009; Murase et al. 2013) and star-forming galaxies (Loeb & Waxman 2006; Katz et al. 2013; Murase et al. 2013; Tamborra et al. 2014; Anchordoqui et al. 2014b). The neutrino spectrum from hadronuclear interaction inherits the cosmic-ray spectrum, $E_\nu ^2 \Phi _\nu = 3/8 \, f_{\rm pp} \, E_p^2\,J_p$. In order to connect the flux of the neutrinos and the UHECRs of the excess, a pion production efficiency f_pp ∼ 0.2 is assumed for an intrinsic UHECR spectrum falling as E^−2.2. Specifically, f_pp = n_p σ_pp κ_pp c t_int, where n_p is the target baryon density, t_int is the duration of the interactions and σ_pp ∼ 10⁻²⁵ cm² and κ_pp ∼ 0.5 are the proton–proton interaction cross section and elasticity.

In pulsars and magnetars, the baryon target is the surrounding supernova ejecta. For an ejecta mass M₁ = M/10 M_☉, the expansion speed is β_−1.5 c = 10⁹ cm s⁻¹, $f_{\rm pp} = 0.2\,M_1\,\beta _{-1.5}^{-2}\,t_{\rm yr}^{-2}$ for a duration t_yr = t/1 yr. The neutrino production from a newborn pulsar with surface magnetic field 10¹² G and initial spin period 1 ms is indicated as the red solid line in Figure 3. Notice that the typical time scale for a magnetar is t ≈ 300 s, so while UHECRs can be accelerated, they are unlikely to escape from the source (Fang et al. 2012).

Alternatively in galaxy clusters, cosmic rays can be accelerated by the accretion shocks formed around large scale structures, and then subsequently interact with nucleons of the intergalactic medium. A typical intergalactic gas density of n ∼ 10⁻⁴ cm⁻³ (Voit 2005) yields a value of f_pp ranging from 10⁻⁴ to 10⁻² (Kotera et al. 2009). Hence the neutrino flux from a single galaxy cluster providing the TA UHECR flux is far below the IceCube detection threshold.

In star-forming galaxies, cosmic rays can be accelerated in the collisionless shocks produced by supernova explosions, and/or by the high energy accelerators contained in the galaxy. The starburst interstellar medium then serves as the baryon target with typical column density Σ_{g, −1} = Σ_g/0.1 g cm⁻², corresponding to a pion production time $t_{\rm pp} = 3\,\Sigma _{g,-1}^{-1}\,H_{\,\rm kpc}\,\rm Myr$. The time cosmic rays stay in the galaxy, t_esc, depends on energy. At low energy cosmic rays are confined in the starburst-driven wind and transport via advection, so $t_{\rm esc}\approx 10\,H_{\rm kpc}\,v_{g,7}^{-1}\,\rm Myr$ for a galaxy with scale height of ∼1 kpc and wind speed v_{g, 7} = v_g/100 km s⁻¹. These particles undergo strong interactions as f_pp = t_esc/t_pp ∼ 1. At higher energy diffusion dominates, then $t_{\rm esc}\approx H^2/4D = 1.6\,D_{0,26}^{-1}\,H_{\rm kpc}^2\,(E/100\,\rm PeV)^{-1/3}\,\rm Myr$, assuming Kolmogorov turbulence and diffusion coefficient D₀ ≈ 10²⁶ cm² s⁻¹ at GeV (Murase et al. 2013). A spectral break at $E_\nu ^b\approx 1\,\Sigma _{g,-1}^{3}\,H_{\rm kpc}^3\,D_{0,26}^{-3} \,\rm PeV$ is expected as cosmic-rays with higher energy are less confined for interactions. This case is shown as red dash-dotted line in Figure 3.

Another important class of high-energy neutrino sources is AGNs (Stecker et al. 1991; Winter 2013; Murase et al. 2014; Dermer et al. 2014). In BL Lac objects, non-thermal photon emission from the inner jet is the dominant interaction target for UHECRs. The photomeson production rate is $f_{\rm p\gamma } = 10^{-4}\, L_{\gamma,45} \, t_5^{-1} \, \Gamma _1^{-4}\, (\epsilon _b / 100\,\rm eV)$ at time t₅ = t/10⁵ s in a jet with Lorentz factor Γ₁ = Γ/10 and gamma-ray luminosity L_{γ, 45} = L_γ/10⁴⁵ erg s⁻¹ (Murase et al. 2014). The TA excess is located in the direction of Markarian 421, which is among the best candidates for a blazar producing UHECR. This is due to the fact that Mrk 421 is both among the brightest γ-ray blazars (Abdo et al. 2011) and is relatively close (∼134 Mpc Gaur et al. 2012), allowing UHECRs to reach Earth at energies above the GZK cutoff.

However, given the low pion production efficiency of Mrk 421, it is difficult for this source to produce the bright neutrino flux observed by IceCube (Dimitrakoudis et al. 2014). One caveat to this analysis is the addition of thermal gas in the source region, which would make it more opaque to UHECRs. However, this is difficult to reconcile with X-ray/TeV observations that suggest the ratio of thermal and synchrotron photons to be ∼0.1 Fossati et al. 2008). One possible solution is an existence of high-energy hardening in the UHECR spectrum, significantly increasing the PeV cosmic-rays responsible for the IceCube neutrino flux.

In contrast to AGNs such as Mrk 421, flat spectrum radio quasars (FSRQs) provide a photopion production efficiency of 1%–10% due to interactions with photons of the broad-line region and the accretion disk, but the cosmic-ray spectrum must be softened above 10¹⁶ eV in order to avoid overproduction of EeV neutrinos. Therefore FSRQs cannot be significant UHECR sources (Dermer et al. 2014).

We summarize the above discussions in the first three columns of Table 1. In the fourth column we indicate the availability of such sources within 20° of the hotspot center. Notice that GRBs and pulsars are transient sources, and their neutrino emission time ($\delta t_\nu ^{\rm GRB} \sim \rm s$, $\delta t_\nu ^{\rm pulsar}\sim \rm yr$) is much shorter than the dispersion time in UHECR arrivals, $\delta t_{\rm UHECR} \approx 4\times 10^4\,(l/100\,\rm Mpc)\,(\delta \alpha /1^\circ)^2\,\rm yr$ for every 100 Mpc propagation and 1° deflection in the EGMF (Kotera & Lemoine 2008). Therefore time correlation between neutrino and UHECR arrivals is not expected for these types of sources, unless that sources are born periodically in certain star-rich region. On the other hand, for the steady types, we investigate the relevant sources in the available high-energy catalogs and surveys (Gao & Solomon 2004; Piffaretti et al. 2011; Ackermann et al. 2012b, 2014; The Fermi-LAT Collaboration 2013). We list a selection of results in Table 2. Future measurements will narrow down the candidates, with high energy neutrino observation limiting or confirming the source direction(s), and more statistics of UHECRs helping to constrain the source distance(s) and other properties.

The recent observation of a UHECR hotspot by the TA observatory has potentially opened the doors to a new era of UHECR point source detections. Interestingly, the TA hotspot also correlates with a marginally statistically significant overabundance in the IceCube neutrino flux. We find that the morphology and spectrum of the TA and IceCube observations are consistent with each other. Moreover, among possible scenarios, a single source with the type of star-forming galaxy can successfully host both excess at the same time.

Again, we caution that the correlations presented in this work are highly tenuous, given the extremely small statistics in both the IceCube and TA data sets. However, this fact conversely makes the possible detection of a single point source appealing, since upcoming data releases by IceCube and the TA should confirm or rule out this detection.

Shortly after we submitted this work, the IceCube Observatory updated the high-energy neutrino list obtained from unbinding the 2012–2013 data (Aartsen et al. 2014). Three additional events were found in the Northern Hemisphere (not including event 37 which is a track event with very low deposition energy and is highly likely to be a background event). With the three-year IceCube data the significance of events 9 and 26 overlapping with the TA excess drops to 1.6σ, but the conclusions of this paper remain unchanged.

We would like to thank Markus Ahlers, Jacob Feintzeig, Jeff Grube, Dan Hooper, Albrecht Karle, Kohta Murase, and Nathan Whitehorn for helpful conversations. We also acknowledge comments by the anonymous referees, which lead to significant improvements in this manuscript. K.F. and A.V.O. acknowledge financial support from NASA 11-APRA-0066 and NSF grant PHY-1068696 at the University of Chicago. T.F. is supported by a Grant-in-Aid for the Japan Society for the Promotion of Science Fellowship for Research Abroad. T.L., is supported by the National Aeronautics and Space Administration through Einstein Postdoctoral Fellowship Award Number PF3-140110. K.F., T.L., and A.V.O. acknowledge support from the Kavli Institute for Cosmological Physics through grant NSF PHY-1125897 and an endowment from the Kavli Foundation.