TheAstrophysicalJournal,832:109(30pp),2016December1 doi:10.3847/0004-637X/832/2/109 ©2016.TheAmericanAstronomicalSociety.Allrightsreserved. BOW TIES IN THE SKY. I. THE ANGULAR STRUCTURE OF INVERSE COMPTON GAMMA RAY HALOS IN THE FERMI SKY Avery E. Broderick1,2, Paul Tiede2, Mohamad Shalaby1,2,3, Christoph Pfrommer4, Ewald Puchwein5, Philip Chang6, and Astrid Lamberts7 1DepartmentofPhysicsandAstronomy,UniversityofWaterloo,200UniversityAvenueWest,Waterloo,ON,N2L3G1,Canada 2PerimeterInstituteforTheoreticalPhysics,31CarolineStreetNorth,Waterloo,ON,N2L2Y5,Canada 3DepartmentofPhysics,FacultyofScience,CairoUniversity,Giza12613,Egypt 4HeidelbergInstituteforTheoreticalStudies,Schloss-Wolfsbrunnenweg35,D-69118Heidelberg,Germany 5InstituteofAstronomyandKavliInstituteforCosmology,UniversityofCambridge,MadingleyRoad,Cambridge,CB30HA,UK 6DepartmentofPhysics,UniversityofWisconsin-Milwaukee,1900E.KenwoodBoulevard,Milwaukee,WI53211,USA 7TheoreticalAstrophysics,CaliforniaInstituteofTechnology,Pasadena,CA91125,USA Received2016July15;revised2016August31;accepted2016September1;published2016November21 ABSTRACT Extended inverse Compton halos are generally anticipated around extragalactic sources of gamma rays with energiesabove100GeV.TheseresultfrominverseComptonscatteredcosmicmicrowavebackgroundphotonsby a population of high-energy electron/positron pairs produced by the annihilation of the high-energy gamma rays on the infrared background. Despite the observed attenuation of the high-energy gamma rays, the halo emission has yet to be directly detected. Here, we demonstrate that in most cases these halos are expected to be highly anisotropic, distributing the upscattered gamma rays along axes defined either by the radio jets of the sources or oriented perpendicular to a global magnetic field. We present a pedagogical derivation of the angular structure in the inverse Compton halo and provide an analytic formalism that facilitates the generation of mock images. We discuss exploiting this fact for the purpose of detecting gamma ray halos in a set of companion papers. Key words: BL Lacertae objects: general – gamma-rays: diffuse background – gamma-rays: general – infrared: diffuse background – plasmas – radiation mechanisms: nonthermal 1. INTRODUCTION which peaks in comoving units near z≈1 due to the peak in the cosmological star formation rate around that time. As The extragalactic gamma ray sky above 100GeV is such,observationsoftheabsorbedVHEGRspectraofnearby dominated by the unresolved emission of a subset of active galacticnuclei(AGNs)(Ackermannetal.2016a).Ofthese,the sources have resulted in direct measurements of the EBL vast majority are blazars—objects with relativistic jets pointed (Biteau & Williams 2015). From these, it is clear that even in our direction (see, e.g., Table5 of Ackermann et al. 2011). today the universe is effectively optically thick to VHEGRs, with While the mechanisms by which these very high-energy gamma rays (VHEGRs) are produced remain unclear (Man- nheim 1993; Ghisellini et al. 1998; Böttcher 2007), their D (E , z) = D (z)E0 = 35⎜⎛1 + z⎟⎞-z⎛⎜ Eg ⎞⎟-1 Mpc, propagation through the cosmos has provided an invaluable pp g pp,0 E ⎝ 2 ⎠ ⎝1TeV⎠ g means by which to probe the intervening universe (Gould & (1) Schréder 1967; Stecker et al. 1992; de Jager et al. 1994; Salamon & Stecker 1998; Domínguez et al. 2011; Gilmore whereE isafiducialenergyandζ=4.5for z<1andζ=0 et al. 2012; Vovk et al. 2012). 0 for z…1 (Kneiske et al. 2004; Neronov & Semikoz 2009). Extragalactic VHEGR sources are observed to be strongly More recently, the evolution of the electron–positron pairs biased toward low redshifts, with the number of known sourcespeakingataredshiftof0.1–0.2(see,e.g.,theredshift has provided a means to probe the magnetization of the interveningcosmos.ThehomogeneityoftheEBLcoupledwith distribution of high-syncrotron-peak sources in Ackermann et al. 2011, 2015). This is a natural consequence of the the large Dpp places many of these pairs within intergalactic annihilationofVHEGRonthenearlyhomogeneousinfrared– voids, where their propagation can be affected by the ultraviolet extragalactic background light (EBL) which intergalactic magnetic field (IGMF). These pairs thus osten- sibly permit the only measurement of the large-scale IGMF permeates the universe, generated by previous generations ofstarsandquasars(Gould&Schréder1967),andwhichcan located in the mean density regions far removed from galactic thus be probed by propagating VHEGRs (Ackermann activity. et al. 2012c; Dwek & Krennrich 2013). The center of The high energy of VHEGRs imply similarly high-energy pairs, which correspond to Lorentz factors of momentum energy of the VHEGRs and EBL photon exceeds the pair-creation threshold, i.e., E E  4m2c4; and thus g IR e E E VHEGRs can annihilate as they propagate through the EBL. g » g » 106 g . (2) Inpractice,themeanfreepath,Dpp,ofVHEGRstoabsorption 2mec2 1 TeV ontheEBLisbothenergy-andredshift-dependent,depending on the evolving density and spectrum of the EBL. Higher If nothing else happens, these pairs will cool on the cosmic EBL densities at larger redshifts correspond to shorter D , microwavebackground(CMB),producinganinverseCompton pp 1 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. cascade (ICC) of photons with typical energies of ⎛ E ⎞2 E E » 2g2E = 2⎜ g ⎟ CMB GeV. (3) IC CMB ⎝ ⎠ 1TeV 1 meV That is, the ICC effectively reprocesses an initial TeV gamma ray into many GeV gamma rays. It is the non-observation of this ICC component in known VHEGR sources that has provided the strongest lower limits on the IGMF to date (see, e.g., Neronov & Semikoz 2009). In a number of extragalactic VHEGR sources, the intrinsic gamma ray spectrum can now be constructed after making weak assumptions either about the intrinsic spectrum or the absorption on the EBL, and thus the resulting ICC emission is estimated. This is then limited directly by observations by the Fermi gamma ray space telescope, which has ruled out the presence of the ICC component with extraordinary confidence (e.g.,Neronov&Semikoz2009).ThisisnaturalifanIGMFis Figure1.AngularsizeofVHEGRmeanfreepathasafunctionofredshiftand energy. present—within the IGMF, the electron–positron pairs deflect away from the line of sight and, therefore, the upscattered gamma rays are beamed away from us. Based The situation would be immediately clarified by the direct upon this scenario, typical estimates for the IGMF range from detection of the ICC component, lying squarely within the 10−17–10−15G,dependingonassumptionsondutycycles(see energy range probed by the Large Area Telescope (LAT) on e.g., Neronov & Semikoz 2009; Neronov & Vovk 2010; Fermi.Notonlywoulddoingsoobviatetheaboveassumptions Tavecchio et al. 2010, 2011; Dermer et al. 2011; Dolag et al. (importantly including the third), it would settle questions 2011; Taylor et al. 2011; Takahashi et al. 2012; H. E. S. S. regardingthedurationofVHEGRoutburstsandtherebyreduce Collaboration et al. 2014; Prokhorov & Moraghan 2016). theuncertaintyontheIGMFlowerlimitssubstantially(Dermer This argument for non-zero IGMFs is predicated on three et al. 2011). Even in the presence of an IGMF, the ICCs are key assumptions: deflectedawayfromtheaxisalongwhichtheoriginalVHEGR emission is beamed. Thus, the ICC component should be 1. The intrinsic TeV emission is narrowly beamed. visible for observers who either do not see the VHEGR 2. The intrinsic TeV emission spectrum can be reasonably emission, or see only weak VHEGR emission. approximated, usually by an exponential cutoff Efforts to directly detect the ICC component are fundamen- power law. 3. No other processes control theevolution ofthe electron– tallycomplicatedbythelargemeanfreepathsoftheVHEGRs, typicallyresultinginhalosthatextendovermanydegreesand, positron pairs. therefore, have low surface brightness (see, e.g., Figure 1). The first is well supported by the prevalence of blazars Therefore, all efforts to date to detect this emission have among VHEGR-bright AGN specifically, and gamma ray stacked multiple Fermi gamma ray images to increase the bright AGN, generally (Wakely & Horan 2008; Ackermann significance with which the halo emission can be separated et al. 2011, 2015), which immediately implies that VHEGR from that due to the central source and background. These emission is localized near the axis of the radio jet. The second efforts are further complicated operationally by the uncertain- is reasonably well supported by the gamma ray spectra of ties in the point-spread function (PSF) of Fermi and the nearby VHEGR-bright AGN, the systematic softening of spatially varying gamma ray background (e.g., Neronov observed gamma ray spectra with increasing redshift (Ack- etal.2011;Ackermannetal.2013).Asaresult,thisprocedure ermann et al. 2012c), and the underlying assumption that the has led to now disproven detections of an excess (Ando & intrinsic spectra only weakly evolve. Kusenko 2010). A more recent attempt that utilizes the most The third remains unclear. Should any alternate cooling recentPSF,reportedinChenetal.(2015),neverthelessexhibits mechanism dominate inverse Compton cooling, it would similar sensitivies to the uncertain instrument response. preempt the generation of ICCs directly. A recently suggested All such efforts have ignored the possibility of structure examplewouldbecoolingmediatedbylarge-scalebeam-plasma withinthegammarayhalo.However,suchstructureisanatural instabilities driven by the bulk motion of the relativistic consequence of either the original beaming of the VHEGRs electron–positronpairsthroughtheionizedintergalacticmedium responsible for the generation of the pairs or the orientation of (Broderick et al. 2012; Schlickeiser et al. 2012, 2013; the IGMF where the pairs are created (see, e.g., Long & Chang et al. 2014). While the nonlinear development of these Vachaspati 2015). Thus, for a wide range of parameters for an instabilities is uncertain, there are a variety of lines of IGMF, we expect a highly anisotropic ICC halo. Here, we astronomical evidence that suggest a cooling mechanism with present semi-analytical computations of the halo structure, verysimilarpropertiesisatwork(Changetal.2012;Pfrommer explicitlydemonstratingthepresenceofthestructure,identify- et al. 2012; Puchwein et al. 2012; Broderick et al. 2014a; ingitsorigin,andcreatingthefacilitytogeneratemockimages Lambertsetal.2015).Regardlessoftheoriginoftheadditional of VHEGR sources with realistic ICC halo structures. In cooling, however, should the ICCs be preempted, the resulting principle,knowledgeofthehalostructurecanaidsubstantially gamma ray spectra would necessarily be consistent with the in efforts to directly detect the ICC component. We report on current lack of a detection of an ICC in known extragalactic anexplicitimplementationofamethodtodosoinacompanion VHEGR sources, independent of an IGMF. paper (P. Tiede et al. 2016a, in preparation). In a companion 2 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. Figure2.CartoonsofthemechanismsbywhichanisotropyintheICChalosisgenerated,distinguishedbythestructureoftheunderlyingIGMF.Left:foranIGMF tangledonsmallscales(correlationlengthsoflB3Mpc)theanisotropyisduetothestructureofthegammarayjet.Right:foranIGMFthatisuniformacrossthe gammarayjet(correlationlengthlB30Mpc),theanisotropyisduetothegeometryofthegyrating,relativisticpairs.Inthelatter,inverseComptongammarays fromelectronsandpositronsareshowninredandgreen,respectively.Alldistancesarenottoscaleandopeninganglesareexaggeratedforvisualpurposes. letter, we apply thisformalism to Fermi-LAT data anddiscuss generates a resultant structure in the GeV image. This is the consequences of this measurement for the IGMF (P. Tiede shown explicitly in the left-hand panel of Figure 2, along with et al. 2016b, in preparation). the associated gamma ray image of the ICC halo. In Section 2, we describe qualitatively the origin of the Alternatively, the process of gyration in the IGMF can also anisotropy and how it relates to the structure of the IGMF. impartstructureontheimage,ifweconsiderablazarwherethe Generalexpressionsdescribingthegenerationandevolutionof VHEGRs are beamed toward us. In the presence of an IGMF the energetic pairs are presented in Section 3. Applications to that is homogeneous on scales comparable to D , electrons pp cases of highly tangled and ordered fields are described in and positrons will gyrate on fixed trajectories that emit toward Sections 4 and 5, with typical applications shown. The anobserveronlyforasubsetofinitialinjectionpositions.This construction of mock Fermi images, including various is still superimposed on the jet structure, resulting in an componentsandinstrumentaleffects,isdiscussedinSection6. asymmetric image structure if the line of sight does not The dependence of the mock ICC halos for a typical bright coincide with the boresight of the jet, as shown in the right- Fermi source on intrinsic source parameters and IGMF hand panel of Figure 2. Gamma rays on opposite sides of the structure is explored in Section 7. Finally, concluding remarks original AGN are produced predominantly by different lepton aremadeinSection8.Inordertostreamlinethepaper,mostof species, i.e., positrons on one side and electrons on the other. the demonstrations are left to the appendix. Hence, different AGN populations can be used to probe the existence of small- or large-scale tangled IGMFs. While non- alignednonthermallydominatedAGNsaresuitableforprobing 2. QUALITATIVE ORIGIN OF STRUCTURE small-scale fields, blazar geometries are ideal for exploring large-scaleIGMFs,becauseinbothcasestheassociatedstrong Before describing the creation of physically realistic ICC anisotropy of the pair halos enables efficient stacking of the halos,webeginwithasummaryofthekeyideasunderlyingthe structures we will find. This is predicated on the standard angular power spectra. In the following sections, we make theseinstancesexplicit,computingphysicallyrealistichaloflux picture of the ICC halo formation described in Section 1: distributions, which connect the energy-dependent flux dis- VHEGRs are emitted from AGNs and travel cosmological distances prior to generating energetic electron–positron pairs tributions to theunderlying physical properties ofthe VHEGR emission and IGMF geometry. on the EBL via photon annihilation. Those pairs then inverse Compton upscatter of the CMB to GeV energies over a comparatively short distance. However, for two independent 3.FORMATIONANDEVOLUTIONOFINTERGALACTIC reasons, these ICC halos are not isotropic. PAIRS First, we consider the case of a small-scale tangled IGMF, The generation of the halos involves two critical steps that which would isotropize the generated electron–positron pairs, imprint anisotropy upon the resulting gamma ray sky: the in combination with a non-aligned nonthermally dominated generationoftheinitialpairsbytheVHEGRemissionandtheir AGN, i.e., the line of sight is not intersecting the jet opening subsequent evolution during the inverse Compton upscattering angle. The VHEGRs are originally beamed along the jet axis. of the CMB. To simplify the computation, we will exploit the This is evidenced by the overwhelming dominance of blazars disparityinscalesbetweentheinverseComptoncoolinglength intheextragalacticgammarayAGNsample(Ackermannetal. (0.01–0.1Mpc) and the VHEGR mean free path, and thus 2011, 2015). Because the VHEGR mean free path is long in assumethattheICCemissionisgeneratedinsitufollowingpair comparison to the inverse Compton cooling time of the production. Therefore, the necessary elements are: resulting pairs, this implies that the emission is essentially local,and,therefore,arisesfromabiconicalregionindicatedby 1. The spatial and energy distribution of VHEGR photons, theradiojetofthesourceAGN.IftheinverseComptongamma and the corresponding distribution of the injected pairs. rays are isotropically emitted, arising, e.g., from a highly 2. Models for the consequences of pair energy loss and tangled IGMF, the spatial structure in the gamma rays deflection in the IGMF. 3 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. 3. The ICC spectrum due to upscattering the CMB. Forconcreteness,wewillassumethatthespectralandspatial distributions of the VHEGRs are separable. This may not be a We first treat each element generally, and then specialize to good approximation at large angles if the emission varies specific limits of interest in Sections 4 and 5. throughout the jet. However, for the applications of primary interest, where the largest effects occur at viewing angles 3.1. Spatial and Energy Distribution of Injected Pairs withintherelativisticbeamingpatternofthejet,i.e.,G-1,thisis jet The origin of the ultrarelativistic pairs is the VHEGR well motivated. emission from the AGN. That nearly all of the bright For the spectrum, we assume a broken power law, extragalacticTeVsourcesareblazars,impliesthattheVHEGR characterized by different photon spectral indexes above (Γ ) h emission from the TeV-luminous AGN is strongly beamed andbelow(Γ)somepivotenergy(E );anexamplespectrumis l p along the jet. Thus, we begin with a description of the shown in Figure 3. The assumed values of the photon spectral anisotropic VHEGR flux from the AGN itself, before indexes are set by the typical values for the high synchrotron annihilationontheEBLmodifiesit.Specifically,westartwith peak sources (HSPs) in the Fermi sample: Γ≈1.8 and l the number of VHEGR photons passing through a solid angle Γ ≈2.5(Ackermannetal.2011,2015,2016b).The1σrange h dW¢ inthedirection xˆ¢,withenergiesbetweenE′andE′+dE′, about these values for HSPs is given by the dotted lines in over a time interval dt′: Figure 3. For the angular flux distribution, we assume a Gaussian jet profile with opening angle θ. That is, we set dN j AGN (E¢, xˆ¢). (4) ⎧ dt¢dE¢dW¢ dNAGN (E¢, xˆ¢) = f G(q¢)⎪⎨(E¢ Ep)-Gl E¢< Ep (9) Foroursystemofcoordinates,primedquantitiescorrespondto dt¢dE¢dW¢ 0 ⎪⎩(E¢ E )-Gh E¢ E , p p the frame of the AGN (including the sites of pair production), assumed to be at a fixed redshift, and centered upon the AGN where itself. Unprimed quantities correspond to the observer frame, centered upon the observer. Note that x¢ and Ω′ correspond to G(q¢) = e(cosq¢-1) q2j + e-(cosq¢+1) q2j. (10) the position in three-dimensional space and the solid angle The two terms correspond to the two oppositely oriented jets relative to the AGN, respectively. The former is related to withangularstructuressetbyG(q),whichforq  1rad,asis j angular position on the sky,a, and distance along the line of typicallythecasefortheAGNconsideredhere(seeSection7), sight, in the direction ℓˆ, by correspond to Gaussian jet profiles. Note that despite the fact x¢ = D (z)a + ℓ¢ℓˆ. (5) that we are considering a version of the Gaussian jets, the A approximate form simplifies significantly as a consequence of whereDAistheangulardiameterdistance.Hence,forexample, writing this in terms of cosq¢. the radial distance from the AGN is given by Thisemissionisnaturallynormalizedbythe1–100GeVflux r¢2 = D (z)a2 + ℓ¢2 and the angle relative to the jet axis (that observed by Fermi. This necessarily depends on Θ′, although A is oriented along zˆ¢), located at an angle Θ′ relative to the line forFermiobjects,inpracticedoessoonlyweakly;thefactthat of sight and Φ′ relative to a fiducial direction on the sky (e.g., thehardFermiAGNsareblazarsimpliesthattheviewingangle north) defining the x-axis in the image plane, is given by is smaller than the relativistic beaming angle of the jet (i.e., Q¢  G-1). Thus, in the absence of a dominant inverse jet cosq¢= zˆ¢ · xˆ¢ Compton halo, assuming a roughly fixed effective area from ℓ¢cosQ¢ + D sinQ¢(a cosF¢ + a sinF¢) 1–100GeV, thegamma rayflux observed by Fermidue to the A x y = . (6) intrinsic emission of the source between these energies is r¢ Note that along the line of sight, i.e., at a = 0, r′=ℓ′ F =ò100 GeV dEò dW¢ dt¢dEV¢HEGR dNAGN and θ′=Θ′. 35 1 GeV W¢Fermi dt dE dt¢dE¢dW¢ As VHEGRs propagate, they annihilate upon the EBL, and 1 100 GeV dN thus both the flux and spectrum of the VHEGRs deviate from = ò dE¢ò dW¢ AGN those emitted by the AGN. Due to the homogeneity of the 1 + z 1 GeV W¢Fermi dt¢dE¢dW¢ A EBL, the optical depth to annihilation is given by =(1 + z)-3 efff E N(G, z)G(Q¢), D2 0 p l r¢ A t(E¢, z, r¢) = , (7) (11) D (E¢, z) pp where the source-frame solid angle subtended by Fermi is wherethemeanfreepathisgiveninEquation(1)anddepends A D2(1 + z)2 and N(G, G, E , z)isanormalizationfactor uponbothVHEGRenergyandredshift(z).8Theresultingflux eff A l h p depending only on the spectral shape: distribution of VHEGRs far from the AGN is then N(G, G, E , z) dN dN l h p VHEGR (E¢, x¢) = AGN (E¢, xˆ¢)e-t(E¢,z,r¢). (8) dt¢dE¢dW¢ dt¢dE¢dW¢ =ò(1001+z) GeV Ep dx[x-GlQ(1 - x) + x-GhQ(x - 1)], (1+z) GeV Ep (12) 8 In principle, the optical depth should be integrated over r′ (and thus z′). in which Θ(x) is the Heaviside function. For our However, in practice, in the cases of interest D is much shorter than the pp Hubblelength,justifyingoursimplerexpression. purposes, we will presume that Ep > 100GeV generally, 4 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. determination of the evolved electron distribution. This evolution occurs under the action of three processes: injection by the VHEGRs, possible gyration in an IGMF, and inverse Compton scattering. Different assumptions regarding the scale and strength of the IGMF will impact the evolution and, therefore, the image. Throughout, we will make use of the ultrarelativistic approximation everywhere permissible, and thus pc and E are used interchangeably. Generally, the evolution of the distribution function of the pairs, f , is described by the Boltzmann equation, e f˙e + v · fe + p˙ · pfe = f˙scat + f˙inj, (15) where f˙ describes the impact of inverse Compton scattering scat and f˙ is the injected distribution of pairs. p and v are the inj momentum and velocity of the pairs, respectively. Ingeneral,thescatteringtermmustbewrittenintermsofan Figure3.TypicalassumedintrinsicspectrumforabrightFermiblazar,witha integral over initial and final particle states, describing the spectralbreakat1TeV.Dottedlinesshowthe1σvariationsinthelow-energy multitudeofwaysinwhichparticlesmayscatterintoandoutof and high-energy photon spectral indexes Γl and Γh, respectively, among the thestateofinterest. However,inthesoft-seedphotonlimitthe hardgammarayblazars.Thedark-grayshadedregionindicatesthe1–100GeV scatteringtermcanbesubstantiallysimplifiedasaconsequence region for which we generate ICC halo realizations; the light-gray shaded region indicates the 1–10TeV VHEGR band primarily responsible for the of the small momentum changes during each scattering (see ICChalos. Appendix B). That is, in this limit f˙scat » -p · (Wfe), (16) resulting in the simplified expression N(Gl, z) = (1 + z)1-Gl where [(1GeV Ep)1-Gl - (100GeV Ep)1-Gl] (Gl - 1). Therefore, the normalization is related to the observed flux via W = -4 sTus pp º - pp , (17) 3 m2c2 t m c F D2 (1 + z)3 e IC e f = 35 A . (13) 0 E A N(G, z)G(Q¢) p eff l in which t = 3m c 4s u » 2.4 ´ 1012(1 + z)-4 years is IC e T s The rate at which pairs are created by the VHEGR is set by the asymptotic inverse Compton cooling time for a nonrelati- therateatwhichVHEGRsareannihilated.Eachgeneratedpair vistic lepton. As a result, the Boltzmann equation simplifies to has an energy ¢ = E¢ 2.9 Thus, the production rate of electrons/positrons with energies between ¢ and ¢ + d¢ at a strictly partial differential equation. Strong variability has been observed in many extragalactic position xˆ¢ is VHEGR sources, indicating that variable emission is char- dNe (¢, x) acteristic. However, when this occurs over timescales much shorter than the typical cooling time, the impact on the pair dt¢d¢d3x¢ distribution is small. That is, even in the presence of a rapidly = dW¢dr¢dE¢ dNVHEGR (2¢, x¢) dt (2¢, z) varying f˙ ,thepairdistributioncanbewellapproximatedbya d3x¢ d¢ dt¢dE¢dW¢ dr¢ stationaryinjsolution. Thus, we set f˙ = 0. = 2 e-t(2¢,z) dNAGN (2¢, xˆ¢), (14) The large disparity between theeinverse Compton cooling r¢2 D (2¢, z)dt¢dE¢dW¢ lengthofthepairs(700kpcat0.5TeV)andthepairproduction pp mean free path of the VHEGR photon (currently 800Mpc at 1TeV), implies that this is a fundamentally local and nearly where the additional factor of dτ/dr′ is the rate at which pairs homogeneous process. The latter immediately implies that are being produced locally. It is straightforward to show that v · f issmall,andmaybeneglectedhenceforth.Theformer powerandnumberfluxareconservedintheaboverelation,and impliesethat the energy distribution of f˙ is given by inj is shown explicitly in Appendix A.1.1. Equation(14);wedeferadiscussionoftheangulardescription Atthispoint,wehaveonlytheinjectionspectrumofpairsas to Sections 4 and 5. a function of position. To progress, we need a model for how Therefore, in practice, we seek to solve the pairs evolve and their subsequent inverse Compton emission. p˙ · pfe + p · (Wfe) = f˙inj. (18) 3.2. Evolution of the Pair Distribution Functions For non-pathological p˙(p) and W(p), this is immediately We have greatly simplified the computation of the gamma solvable via the method of characteristics (see Appendix C). ray flux at the expense of pushing the difficulty onto the Fordeflections in a locally fixedmagnetic fielddescribed by a local spherical coordinate system (q¢, J¢, j¢), it is convenient 9 Toavoidconfusion,wewillrefertotheenergiesoftheemittedgammarays, toexpressthisintermsofasetofpolarmomentumcoordinates includingtheVHEGRs,asE′,leptonenergiesas¢,andobservedgammaray energiesasE. aligned with the field (p¢, J¢, j¢), yielding p p 5 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. fe(x¢, p¢) = mepc¢4tIC òp¥¢ dq¢q¢2f˙inj⎡⎣⎢⎢x¢, q¢, J¢p, j¢p  wBtIC2sinJ¢p⎛⎝⎜mqe2¢2c2 - mpe2¢c22⎞⎠⎟⎤⎦⎥⎥, (19) where w = eB m c is the cyclotron frequency. distribution. That is, the seed photon distribution function is B e ⎛ ⎞ u E¢ f (x¢, q¢) = s d⎜q¢- CMB⎟, (20) 3.3. Inverse Compton Emission s s 4pq¢3c ⎝ s c ⎠ s While we will consider a variety of potential models for the where E¢ = 0.7(1 + z) meV and u are the typical energy evolution of the pairs after production, their inverse Compton CMB s and energy density of CMB photons, respectively. Relaxing emissionmaybedescribedbyasingleframeworkinwhichthis thismakeslittledifferencetotheresultingemissionspectrumat evolution enters solely as an unknown electron and positron distribution function, f (x¢, p¢). We construct this framework the expense of complicating the analysis substantially. e here, paying particular attention to identifying a number of Thehighenergiesofthepairsincomparisontothesoft-seed photonsresponsiblefortheircooling,impliesthattheresulting simplifying assumptions. Thefirstoftheseisstationarity,whichinturncorrespondsto inverseComptonemissionishighlybeamedinthedirectionof theleptonpropagation.Wewillassumethatthisisexclusively a presumption regarding the duty cycle of VHEGR sources. thecase,i.e.,theinverseComptonphotonspropagatealongthe Short-timevariabilitywillnecessarilyimpartsimilarvariability direction of the scattering electrons and positrons.10 In the on the energy dependence and spatial distribution of the pair isotropic approximation, this implies that the differential population. Both of these will be smoothed, however, if the sourceisactiveforasufficientperiod.Inthecaseoftheenergy scattering cross-section for upscattering a seed photon with momentum q¢ to a gamma ray with momentum q¢ by an dependence, this is set by the inverse Compton cooling s timescale, roughly 106yr. For the spatial distribution this is electron with initial and final momentum p¢ and pf¢ is determinedbythepropagation-dependenttimedelay,whichfor ⎛ ⎞ 2° halos ranges from 102yr to 4 ´ 103yr, depending on ds = s d3⎜q¢ - 2qsp¢p¢⎟d3(p¢ - p¢ + q¢), redshift and gamma ray energy. Thus, if VHEGR activity d3q¢d3q¢ d3p¢d3p¢ T ⎝ m2c2 ⎠ f s f e persists for longer than 106yr we may assume that the (21) underlying pair population has reached steady state in both terms. Note that this is comfortably short in comparison to where the first δ-function encodes both the forward-propaga- typical radio duty cycles of a few times107yr to a few times tionapproximationandthesoft-seedlimit,whilethesecondδ- 108yr (Alexander & Leahy 1987; McNamara et al. 2005; functionenforces momentumconservation. Notethatsincethe Nulsen etal.2005;Shabalaetal.2008), suggestingthatthisis pair distribution functions will evolve, this need not be in the well justified. As such, we will not consider non-stationary direction of the original VHEGR. transients (though see Menzler & Schlickeiser 2015). Giventheseassumptions,therateatwhichICCgammarays The inverse Compton process is simplified by three with momenta q¢ are produced is additional assumptions, relaxing any of which will make at most small changes to the results. The first is that a single ⎛ ⎞3 2 ⎛ ⎞ gsacneandtetegrareetnidoenrgaatmeofmsuapbasriaeryqssueidsnotcngroeetantteehrdea,mtiosine.els.v,eosfthaaandtndiiihtniivolaentraesleopnaCitrhosem.EpTtBhoiLns dtd¢dN3IxC¢,de3q¢ = s2Tqcs¢ucs⎝⎜⎜2mqe2s¢cq2¢⎠⎟⎟ fe⎝⎜⎜x¢, m2eqcs¢qq¢¢ ⎠⎟⎟ qs¢=EC¢MB c, correspondstoajointconstraintontheenergyoftheVHEGRs (22) considered and the seed photon population that is inverse Compton scattered, effectively limiting it to the CMB. The thederivationofwhichcanbefoundinAppendixD.Fromthis, former constitutes a conservative assumption, limiting our wemayconstructthefluxasseenonEarthfromasourcealong attentiontoVHEGRswithenergieslessthan10TeV.Thelatter a line of sight in the direction ℓˆ¢, and thus corresponding to issuppressedbytheratioofthenumberdensitiesofCMBand momenta q¢ = E¢ℓˆ¢ c, EBL photons—the EBL provides the only substantial alter- (nDatwiveekp&opKurlaetnionnri—cho2f0a1t3l)e.aIsntc1lu0d4i,nagndadtdhiutisonisalragreeneinraptiroancsticoef dt¢dd2Nx¢gd3q¢ = å òdℓ¢dtd¢dN3IxC¢,de3q¢⎛⎝⎜x¢, Ec¢ℓˆ¢⎞⎠⎟, (23) pairs will modify the GeV signal by a comparable amount, e+,e- justifying its neglect. obtained by directly integrating the Vlasov equation (see, e.g., TheinverseComptoncoolingofthepairsisprimarilydueto Broderick 2006). In the above we have neglected any the CMB, with typical energies ¢ » 10-3(1 + z) eV. For subsequent scattering or absorption of the ICC gamma rays, VHEGR energies E¢  m2c4 ¢ ~ 3 ´ 102 (1 + z) TeV, e which is well justified by their typical energy, given by thisissufficientlylowthattheenergyoftheupscatteredphoton Equation (3) and generally100GeV. is well approximated by Equation (3). Typically, we will assume that only VHEGRs with energies less than 10TeV are relevant, hence this is well justified. We will further assume 10 This is an excellent approximation for a single scattering due to the high energiesoftheleptons.Thatthisremainstrueaftermanyscatteringsisshown that the CMB is characterized by a mono-energetic photon inAppendixE. 6 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. The above may immediately be converted into a surface 4.1. Pair Distribution Function brightness as measured by Fermi: Within this context, we approximate the effect of rapid isotropization as an instantaneous redefinition of the pair dNg = dt¢dE¢ ò q¢2dW¢q d2x¢ dNg injection model: dtdEd2a = ådt dEAefWfF¢ermi 3me4cc5 d2a dEt¢d2x¢d3q¢ f˙inj,e(x¢, p¢) = 4pcp¢2 dt¢ddN¢ed3x¢(p¢c, x¢). (25) (1 + z)416t E2 2E e+,e- IC CMB CMB Theisotropyrenderstheazimuthalevolutionofthedistribution ⎛ ⎞ ´ ò dℓ¢fe⎝⎜x¢, mec 2EE ℓˆ¢⎠⎟, (24) function during inverse Compton cooling moot, and thus CMB fe(x¢, p¢) = tI4Cpmpe¢c42 òp¥¢ dp˜¢dt¢ddN¢ed3x¢(p˜¢c, x¢). (26) where E¢ = (1 + z)E and E¢ E¢ = E E were CMB CMB CMB CMB used. At this point, it is necessary to explicitly define the For the expression for the pair injection spectrum in evolution model for the electrons and positrons and the Equation (14), this may be integrated explicitly at each r′, corresponding f (x, p). Thus, we now turn to explicitly producing e tchoensciodnetreixntgotfhtewionjleicmtiiotinngmeovdoellutdioesncarriybemdobdyelEs.qWuaetidoenfe(1r4u)ntiinl fe(x¢, p¢)= tI4CpmEec45dt¢ddN¢ed3x¢⎛⎝⎜E2p, x¢⎞⎠⎟et(Ep,z) Sections 6 and 7 the construction of mock images, which p includeadditionalcontributions(centralsourceandgammaray ⎛¢⎞-4 ¥ ⎛2¢⎞1-Gh background), instrumental effects (i.e., the PSF), and draw ´⎜ ⎟ ò d¢⎜ ⎟ e-t(2¢,z) explicit subsamples with the appropriate statistics. ⎝Ep⎠ ¢ ⎝ Ep⎠ The physical picture of ICCs reprocessing the VHEGR = tICmec5 dNe ⎛⎜Ep, x¢⎞⎟et(Ep,z) emission from the central source to lower-energy gamma rays 8pE3 dt¢d¢d3x¢⎝ 2 ⎠ p permits a variety of integral relationships between the two. These areusedtocheck both theresults ofthis section aswell ⎛¢⎞-(Gh+2)⎛¢⎞Gh-2 ⎛ 2¢⎞ ´⎜ ⎟ ⎜ ⎟ G⎜2 - G, ⎟, (27) as the following two, where the energy-dependent halo ⎝E ⎠ ⎝E ⎠ ⎝ h E ⎠ p D D structures differing in the assumed structure of the IGMF, are dealtwithinAppendixA.Atthesametime,acleardiscussion whereE º E t(E , z) andΓ(k,x)istheincompleteGamma D p p of the anticipated redshift dependence may be found in function of order k (Abramowitz & Stegun 1972). Appendix A.3. 4.2. Gamma ray Halo Structure Withthedistributionfunctioninhand,wemaynowcompute 4. GAMMA RAY HALOS FROM TANGLED FIELDS the gamma ray halo flux distribution from Equation (24) The first limit we consider assumes rapid isotropization of explicitly. This is simplified in this case by the fact that the pairs, i.e., the pair momenta become isotropic on a fe(x¢, p¢)isseparableinx¢and p¢.Uptoanenergy-dependent timescale short in comparison to the inverse Compton cooling coefficient,thefluxdistributioncanbeexpressedasanenergy- time. This naturally occurs if the IGMF is sufficiently strong dependentrescalingofasinglespatialfunction,i.e.,theimages and tangled. This occurs in two steps, first isotropizing in the at different energies form a homologous class of images azimuthal direction about the magnetic field due to gyration, described by and second isotropizing in the poloidal direction due to ¥ variations in the local magnetic field orientation. In practice, Yk(R¢) º ò dℓ¢r¢-k-2G(k, r¢)G(q¢), (28) -¥ this may occur only statistically, should the gyration radius be smallincomparisontothecorrelationlength,i.e.,aftercoarse- where R¢ = DA(z)a, r¢ º R¢2 + ℓ¢2, and θ′ is given by graining overthe coolingscale, or via diffusion if the gyration Equation (6). Hence, in practice, for a given source structure radius is large in comparison to the correlation length, λ . and geometry (i.e., θ, Θ, Φ, and Γ ) it is sufficient to B j h Quantitatively, the IGMF strength and value of λB at which numericallycomputeYk(R).ThisfunctionisshowninFigure4 the pairs effectively isotropize depends on the gamma ray for various viewing angles. energyandsourceredshift.Pairswithhighenergieswillgyrate IntermsofΨ ,thefullenergy-dependent,gammaraysurface k the slowest, placing the strongest limit on IGMF strength; at brightness distribution is given by 10TeV a lepton will make a full gyration within the inverse ComptoncoolingtimewhenB  10-15G.Ifthecurrentlower dNg =F mec2 3(1 + z)-1 ⎛⎜mec2⎞⎟1-Gh limitsonthestrengthoftheIGMFaretakenatfacevalue,they dtdEd2a 35d2E2 4pN(G, z)G(Q¢)⎝ E ⎠ imply this is likely to be the case (however, cf. Dermer A CMB l p et al. 2011). The statistical isotropization requires tangling of ⎛ 2E ⎞-Gh 2 ⎛ 2E a⎞ ´⎜ ⎟ Y ⎜ ⎟, (29) the IGMF on scales comparable to the jet opening angle at ⎝E ⎠ 2-Gh⎝ E d ⎠ CMB CMB A distances of D . In turn, this requires l  q D , which is pp B j pp smallest at 10TeV and z=0: l  3Mpc. Therefore, where we have summed over the two species and B nominally, the applicability of the rapid isotropization limit d º (E m c2)D (E , z) D (z) is a measure of the angular A p e pp p A depends solely on the structure of the IGMF. size of the mean free path of VHEGRs, modulo the large 7 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. Figure4. Y(R)forvariousviewingangles,Θ,andk=−0.6(i.e.,Γ=2.6).Atobliqueviewingangles(Θ10°),theanisotropicstructurefollowscloselythatof k h thejet.Atacuteviewingangles(Θ10°),geometricforeshorteningcurtailsthisstructure,andbyΘ≈2°itisnearlyisotropic. energy ratio appearing in the prefactor. Note that the normal- simplicityincomparisontoMonteCarlopropagationschemes, izationisfullydeterminedbytheobservedF andtheassumed this remains nontrivial, requiring a handful of practical 35 intrinsicspectralshape.Thatthisreproducestheintegratedflux optimizations. inTthheepsapreantitalVsHtrEuGctRurseisofvetrhiefieIdCiCn AhaplpoenfodlilxowAs.2t.h1a.t of the tabWuleatinbgegthine rebsyultcloongsatrriuthctminigcallYy2i-nGhR(R,)rannguinmgerfircoamlly10−an8d– 104. This is then integrated over R and E between 1 and underlying gamma ray jet. That is, since the pairs are 100GeV to obtain the ratio between the number of gamma isotropized, and, therefore, emit via inverse Compton iso- tropically,theirgammarayemissionmapsoutthepairinjection rays in the halo to that contained within the direct-source sites.SincetheVHEGRemissionisconfinedtojets,soarethe component,i.e.,thehalo-corefluxratio.Wethenconstructthe following additional marginalized probability distributions for pairs, and consequently so is their emission. As a result, the a halo photon on the image: orientation of the ICC halos in this limit follows that of the gamma ray jet and presumably the much smaller-scale dN radio jets. P(E)µò d2a g Generally,theICChalospectrumwillvaryacrosstheimage dtdEd2a due to the different mean free paths of the underlying P(a∣E)µò dq dNg VHEGRs. However, it is possible to characterize the large- adtdEd2a scale ICC halo spectrum. Yk(R) peaks near a fixed R, falling dNg rapidlythereafter,correspondingtotheexponentialsuppression P(qa∣E, a)µ , (30) dtdEd2a oftheinjectionofpairsbytheVHEGRafterameanfreepath. As a result, for sufficiently extended images the internal wherea = ∣a∣,q is the polar angle on the sky relative to the a structureofY (R)isunimportantandtheICChalospectrumis k projectedjetaxis,andthenormalizationsaresetappropriately. proportional to E-1-Gh 2. The resolved spectra are generally Thefirst isalreadyknown analytically, P(E) µ E-1-Gh 2. The harder due to the smaller D . pp lattertwoareconstructednumericallyintermsofintegralsover For the nearly flat spectrum sources under consideration Y (R) and tabulated. To generate a gamma ray event the here, typical values for the low- and high-energy photon 2-Gh spectral indexes are Γ≈1.8 and Γ ≈2.5. Given the latter, procedure is then: l h thetypicalICChalophotonspectralindexis2.25,intermediate 1. Drawarandomnumbertodetermineifthegammarayis to the two, and importantly, intermediate between that of the in the halo or the source based upon the halo-core flux source and the soft background (which has a typical photon ratio. If it is in the halo: spectral index of 2.4). This has two relevant consequences. 2. Draw an energy from P(E). First,itisdifficulttospectrallyseparatetheICChaloemission 3. Draw a value for α from P(a∣E). fromthesourceorthebackground.Second,sincetheICChalos 4. Draw a value for θα from P(qa∣a, E). are marginally harder than the background, the marginal value oflow-energygammarays(1GeV)issmall,typicallyadding more noise than ICC halo signal. For techniques that leverage 4.4. Example Halo Realizations the anisotropic structure of the halos, this is compounded Example realizations of ICC halos associated with a small- by thegrowth of the PSF at low energies. Hence, our decision scale,tangledIGMFareshowninFigure5forvariousviewing to restrict our attention to gamma rays with energies angles. As anticipated implicitly by the energy dependence of above 1GeV. D and explicitly by Equation (29), the energy of the gamma pp rays is strongly correlated with the distance from the central source. High-energy halo photons are found at the smallest 4.3. Numerical Implementation angular offsets, while low-energy halo photons compose the The creation of mock Fermi images will ultimately require bulk of the distant halo. the rapid generation of samples drawn from the energy- At large viewing angles, the expected angular structure of dependent flux distribution in Equation (29). Despite its thehaloiseasilyvisible,bothbeforeandafterconvolutionwith 8 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. Figure5.ExamplerealizationofICChalointhepresenceofasmall-scaletangledmagneticfieldviewedatajetinclinationof60°(left),30°(middle),and2°(right) before(top)andafter(bottom)convolvingwiththeFermiPass8R2_V6-frontPSF.Inallcases,theon-axisnumberofgammarayswasassumedtobe5000andthe sourceplacedatz=0.3,θ=3°,andΓ=2.6.Theenergyofthegammaraysisindicatedbycolor:with1–3GeVindarkred,3–10GeVinorange,10–30GeV j h ingreen. the Fermi Pass P8R2_V6 PSF for front-converted events (see 130Mpc.Thus,typically,theIGMFmaybetreatedasuniform Section 6.1). Small viewing angles, commensurate with those ifl  102Mpc. B implied by observations of gamma ray blazars, result in an This presents a natural counterpoint to the isotropized case. extremeforeshorteningthateliminatesmuchoftheanisotropy. Italsopermitsasimplifiedgeometricpicture,sinceonlyasmall Inallcases,thetotalnumberofintrinsicsourcegammarays subset of the gyrating pairs will gyrate into our line of sight, seenbyanobserverlookingdownthegammarayjetaxis,i.e., and therefore, generate gamma rays in our direction. the on-axis source flux, was held fixed, corresponding to a The characteristic radial extent of halos in this case is fluence of 5000ph, chosen to match a typical, bright, hard determinedbyboththedeflectionangleoftypicalpairsandthe Fermi source. The comparatively small halo fluence is a result openingangleofthejet.ForobservedgammarayenergyE,the of the isotropy of the ICC component, which follows from the former typically limits the angular extent of the halos to rapid isotropization of the pairs, and dilutes the energy flux from the initially highly beamed VHEGR jet by a factor of Da » wBtIC » 1.5⎜⎛ B ⎟⎞⎜⎛ E ⎟⎞-1(1 + z)-4, 2q-j2. For the parameters employed here, this corresponds to a def g2 ⎝10-16G⎠⎝10GeV⎠ reduction by a factor of nearly 7×102. (31) which for z1 is generally substantial for fields consistent 5. GAMMA RAY HALOS FROM ORDERED FIELDS with the arguments of Neronov & Vovk (2010). This is not The second limit we consider assumes that the distribution surprising given that those require similar size deflections to function in minimally isotropized, evolving over the entire explain the absence of an ICC component in the spectra of cooling time under the action of a uniform IGMF. We further VHEGR sources. The latter sets a typical angular scale for assumethatthisoccurscoherentlyacrossthejet,requiringlarge blazar jets independent of IGMF strength of IGMFcoherencelengths.Boththemagnitudeandnatureofthe limit on λ depends on redshift and energy. For sufficiently qjDpp Dpp B Da » » 2 , (32) nearby objects λ is limited by the proper distance, D —for jet D - D D - D B P A pp A pp Mkn421thisis130Mpc.Formoredistantobjectstherelevant limitissetbyD .ForatypicalICCgammaray(10GeV)from whichiscomparabletothelimitfromthemagneticdeflections pp atypicalmoderate-redshiftsource(z=0.2)thisisalsoroughly fornearbysources(z0.15)anddominatestheICChalosizes 9 TheAstrophysicalJournal,832:109(30pp),2016December1 Brodericketal. fordistancesources.Asaresult,onceagainthestructureofthe Foraleptoninjectedatinfiniteenergywithazimuthalanglej′ IGMFdominatesthegrossfeaturesoftheICChaloforalarge- the e are the energies at which it will have gyrated k times k scale, homogeneous IGMF that is consistent with the current through an azimuthal angle j¢. Thus, these represent the p lower limits on its strength. energiesatwhichthemultiplicityofleptonsorientedinagiven azimuthal direction increases by one.Note that thee depend 5.1. Pair Distribution Function k on species through the sign of ω , corresponding to the B In this case, the pairs retain the memory of their injection direction of the gyration. momentum, and thus, f˙ (x¢, p¢)= c dNe (p¢c, x¢) 5.2. Gamma ray Halo Structure inj,e p¢2sinJ¢ dt¢d¢d3x¢ p To constructthespatialdistribution ofICCgammarays, we ´d(J¢ - J¢)d[j¢ - j¢]. (33) againemployEquation(24),nowwiththedistributionfunction p p in Equation (36). Unlike the rapidly isotropized limit, a large- Perhaps counter-intuitively the j (E′) depend on E′ as a result scale homogeneous IGMF evolves the lepton distribution of the definition of p adoptped in Appendix C, which functionsonaconeoffixedopeningangleabouttheIGMF, B. 0 As a result, the integration along the line of sight receives effectively describes the orientation of the injected leptons contributionsonlywherethelineofsight,ℓ¢,liesonthiscone. not at injection, but by that of a fictitious lepton injected at This admits a simple geometric identification of the gross infinite energy and inverse Compton cooled until E′ (which geometryofthecontributingregions,whichmaymosteasilybe necessarily occurs over a finite time). The corresponding imagined by time-reversing the inverse Compton process. The stationary distribution function during inverse Compton cool- observed ICC gamma rays are produced by leptons directed ing is along ℓ¢, and thus leptons that are gyrating about B on a cone fe(x¢, p¢) = ptI¢C4msinecJ2¢pd(J¢p - J¢)òp¥¢ dq¢dt¢ddN¢ed3x¢(q¢c, x¢)d[j¢p(q¢) - j¢], (34) where we have subsumed all of the azimuthal argument in defined by ℓ¢. Since the initial VHEGRs are assumed to Equation (19) intoj¢ (q¢), i.e., propagate radially from the source, this implies that r¢ must p also lie on this cone. Thus, the ICC halo will receive j¢ (q¢) = j¢ + wBtICsinJ¢p⎛⎜me2c2 - me2c2⎞⎟. (35) contributionsonlyfromleptonslocatedonaconewithopening p p 2 ⎝ q¢2 p¢2 ⎠ anglecos-1(B · ℓ¢ Bℓ¢)andcenteredattheVHEGRsource.In practice, this is enforced via the remaining δ-function in f (x¢, p¢) in Equation (36). e Again the integration may be performed explicitly for the Performing the integration along ℓ¢ is then reduced to the injectionspectruminEquation(14),nowwiththeaidoftheδ- determinationoftherelevantposition onthecontributingcone function. Some care must be taken regarding the infinite atagivenangularpositionandperformingtheintegrationover number of zeros in the argument of the δ-function, corresp- the δ-function. Once again, it is convenient to define a onding to the infinite number of times the momentum of the coordinate system aligned with the line of sight as in electron/positron gyrates into azimuthal (but not necessarily Equation (5) and following Equation (28). Within this basis, poloidal) alignment with a given direction, yielding we also define the three projections B º B · ℓˆ, ℓ B º B · R¢ R¢, and B º B · (R¢ ´ ℓˆ¢) R¢. Then, the con- ⎛ ⎞-4 R X fe(x¢, p¢)= 8w mc3c2E d(sJipn2-JJ)⎝⎜E¢⎠⎟ tributing cone has B e p p p ´ dt¢ddN¢ed3x¢⎛⎝⎜E2p, x¢⎞⎠⎟et(Ep,z) r¢=R¢ B2ℓ2B+ℓBBRR2 , sinJ¢ℓ = B2B- Bℓ2 , ´ å¥ ⎛⎝⎜⎜2Eek⎞⎠⎟⎟4-Ghe-2ek ED Q(ek- ¢), tan(j¢ℓ - j¢)= B2(Bℓ2 -∣2BBR2)Bℓ-BRBBℓ2X(∣Bℓ2 + BR2). k=-¥ p (38) (36) where(J¢, j¢) is the orientation of ℓ¢ as seen from the source. where Θ(x) is the Heaviside function and The anglℓe θ′ℓ relative to the jet axis is then determined via Equation (6). Note that while ℓ¢, and thus Bℓ, is fixed, R¢, and eº m c2 wBtICsinJ¢ 2 . thus both B and B depend on a, i.e., they vary across the k e 2pk - (j¢p - j¢)  wBtICsinJ¢me2c4 2¢2 image.TheiRntegratioXnovertheδ-functionproducesaJacobian (37) in addition to restricting the contributing region to the above 10