text:
# The Gamma-Ray Blazar Content of the Northern Sky

## 1 Introduction

The EGRET telescope on the CGRO satellite has detected 271 sources in a survey of the Gamma-ray ($`100`$MeV to 10 GeV) sky. Of these roughly a quarter have been identified as blazars and along the Galactic plane there are a half dozen objects confirmed as young pulsars through their pulsed $`\gamma `$-ray emission. Thus the bulk of the sources remain to be identified. There is a young Galactic population along the plane clearly correlated with high mass stars (Kaaret & Cottam, 1996; Yadigaroglu & Romani, 1997). An intermediate latitude excess, especially in the direction of the Galactic bulge, suggests an older Galactic population whose nature is uncertain. Finally there remain a substantial number of high latitude sources with no AGN identification.

The selection of blazar candidates has largely proceeded by correlation with existing radio surveys (Hartman et al., 1999; Mattox et al., 2001). In contrast, the Galactic plane sources and individual intermediate latitude sources have been the subject of targeted multi-wavelength campaigns (Roberts, Romani & Kawai, 2001; Halpern et al., 2001; Wallace et al., 2002). In this project we attempt to obtain a more complete census of plausible blazar counterparts, sifting sources with extant radio survey data and then conducting a multi-wavelength follow-up. We are obtaining spectroscopic confirmation of the candidate AGN with Hobby$``$Eberly Telescope (HET) Marcario LRS spectroscopy. This survey has already discovered a number of new likely $`\gamma `$-ray blazars, including good candidates for the most distant persistent $`\gamma `$-ray sources known.

### 1.1 Blazar Properties

The ‘blazar’ label is somewhat heterogeneous, but in the context of the unified AGN model, these sources are believed to be viewed close to the axis of a powerful relativistic jet. As such they are compact flat spectrum radio sources, with apparent superluminal motion at VLBI scales. The optical counterparts exhibit significant polarization and OVV (optically violently variable) behavior (Urry & Padovani, 1995). Optical spectroscopy often shows large equivalent width emission lines (flat spectrum radio quasars) while a significant fraction are continuum-dominated BL Lac-type objects. The broad-band spectral energy distribution (SED) is sometimes used to divide these into two classes with ‘red’ blazars showing a synchrotron peak in the IR-optical with a Compton peak in the $`\gamma `$-ray while ‘blue’ blazars have a synchrotron component extending into the X-ray with a Compton peak inferred to extend to the TeV range (Urry 1999, and references therein). There are surveys underway to substantially increase the number of known blazars (e.g. DXRBS; Landt et al. 2001) which will eventually help in systematizing the broad-band properties of these sources.

Of these properties bright, flat spectrum radio emission seems best correlated with $`\gamma `$-ray activity, but this may be primarily a selection effect of present counterpart lists, which focused on the brightest radio sources as the principal candidates. In an effort to make a less biased census of the counterparts, we have re-examined the correlation with X-ray and radio properties in selecting candidates.

## 2 Counterpart Selection and ‘Figure of Merit’ Ranking

The existing $`\gamma `$-ray blazar lists were largely selected from the Green Bank 6cm and 21cm single dish surveys (Condon, et al., 1991; White & Becker, 1992). With relatively poor resolution, confusion and extended jet/host galaxy emission were a serious impediment to selecting flat spectrum cores as blazar counterparts. This was especially true at lower radio fluxes and toward the Galactic plane. Accordingly few sub-Jy counterparts have been identified. Real progress can now be made since interferometric surveys covering much of the Northern sky are available at 21cm (NVSS; Condon, et al. 1988) and 3.5cm (CLASS; Myers, et al. 2002 – the source list was kindly supplied in advance of publication by Ian Browne). The CLASS survey, targeting compact gravitational lens candidates, preselected flat spectrum sources by comparing NVSS 21cm and Green Bank single dish 6cm fluxes, mapping sources with spectral index $`\alpha 0.5`$ ($`S_\nu \nu ^\alpha `$) with the VLA A array, to resolve structure at the 0.2<sup>′′</sup> scale. Coverage was complete for Dec$`0^{}`$ and $`|b|10^{}`$; there was also partial coverage beyond these limits. We are primarily interested here in the isolated unresolved cores that are not lens candidates.

By comparing the 3.5cm and 21cm interferometric fluxes we have a greatly improved estimate of the spectral index of the compact cores, which are well identified by sub-arcsecond positional matches. We find that a number of sources previously selected as ‘flat’ have extended high frequency emission. Of course with such a variable population, non-simultaneous observations can have erroneous spectral indices, but our estimates are certainly much better than extrapolations from the confused GB6 fluxes.

When CLASS coverage is not available, we do need to use GB6 fluxes. Following the CLASS prescription, we subdivide the 6cm flux between the coincident 21cm NVSS sources to make spectral index estimates. We then extrapolate with the 4.85GHz/1.4GHz spectral index to 8.4GHz to estimate a flux in the CLASS range. In general, this extension was needed in the Galactic plane where a plethora of 21cm sources make spectral indices unreliable. However, in a number of cases we are able to show that no target meets the CLASS survey criterion in the EGRET error box. This will be used to rule out typical blazar-type counterparts. All sources classified as blazars ended up having a CLASS detection, except 3EG J2016+3657.

We augment this radio selection by using the convenient, relatively uniform RASS bright and faint source catalogs. Counterpart candidates are evaluated for suitability of optical follow-up with USNO-B1.0/POSS photometry. Our analysis will use the positions of these counterparts in the ‘Test Statistic’ (TS) maps of the Third EGRET (3EG) catalog, which plot the estimated likelihood of the observed $`\gamma `$-rays as a function of point source location.

### 2.1 Radio/X-ray properties

High frequency radio emission is clearly a good discriminant. However to make the least biased selection of counterparts, we compute the over-density of sources near high latitude ($`|b|>20^{}`$) gamma-ray detections in bins of radio flux, spectral index and x-ray flux. To do this, we compare the number of sources detected within a given TS probability contour with many random realizations of the sky. The latter were obtained by shifting the EGRET TS maps on the sky in 2 degree increments, computing the random source TS values and correcting for variation in sky area from the survey cuts. The source counts within a given confidence contour as a function of, e.g., 8.4GHz flux, were then compared to the random realizations after normalization by the effective sky coverage. We define the excess fractional source density associated with the 3EG source as

$$n=\frac{N_{3EG}N_{Random}}{N_{3EG}}$$
(1)
for each flux bin. The distribution for 8.4 GHz flux within the 95% confidence region is shown in figure 1. A simple least-squares power law fit to the binned data then defines our excess source density function for evaluation of individual candidates. This exercise shows significant excess at fluxes well below the $``$Jy limit considered in earlier analyses. Other authors have also recently concluded that fainter radio sources should be considered as plausible candidates (Wallace et al., 2002).

We have similarly measured the excess of flat spectrum sources as a function of spectral index. Here the correlation is strongest for $`\alpha 0.6`$, becoming negligible for $`\alpha 0.53`$. Correlation with the RASS source catalog also shows an excess. However the maximum excess for the brightest sources is only 0.5 and about half of the strong radio counterpart candidates were not detected in the relatively shallow all-sky survey (some are present in deeper pointed observations). Clearly we do not wish X-ray non-detection to exclude a candidate, although X-ray bright sources are somewhat more likely to be correlated with the EGRET detections. We thus add 0.5 to the best fit X-ray $`n`$ so that it runs from 0.5 to 1.0. In this way X-ray detection can increase the source’s selection by up to a factor of 2 (for $`>`$1PSPC cps), but X-ray non-detection is not seriously constraining. Thus the X-ray component of our analysis is not crucial – if it is excluded no sources are lost from our final source list, although their detailed ranking does change somewhat. The excess probability functions exhibited no strong dependence on the gamma-ray source localization. Therefore, the area within the 95% contours were taken to be representative of the gamma-ray sky for the derivation of the merit functions.

### 2.2 Source Position Weighting

The combination of 8.4GHz flux, spectral index and X-ray flux provides a useful tool for identifying EGRET -like blazars. Indeed, we are studying a set of objects selected by these cuts, without any 3EG coincidence, as possible GLAST-detectable blazars. Given the strong $`\gamma `$-ray variability and the limited duty cycle in the flaring state, it would not be very surprising for relatively bright sources to have eluded detection during the limited EGRET exposure. However, to select individual EGRET counterparts, we will also use the radio source position in the likelihood analysis. Previous blazar identification efforts have used the circular or elliptical fits to the TS maps as an approximation to the allowed sky region. However, many TS maps contain open or complex contours at the 95% level, so we employ instead the precise TS at the radio source position.

The value at our counterpart candidate position $`(\alpha ,\delta )`$ is compared to the map maximum, giving $`\mathrm{\Delta }TS=TS_{max}TS(\alpha ,\delta )`$. $`\mathrm{\Delta }TS`$ has been calibrated by Monte Carlo simulations (Mattox et al., 1996) to give the statistical probability $`L[\mathrm{\Delta }TS(\alpha ,\delta )]`$ that the source lies within the contour of constant $`\mathrm{\Delta }TS`$ for 50%, 68% 95% and 99% confidence levels. We linearly interpolate between these values to obtain a smooth $`L[\mathrm{\Delta }TS(\alpha ,\delta )]`$; sources at large $`\mathrm{\Delta }TS`$ have little probability of association with the 3EG source, even if very blazar-like. The statistical error estimate for $`\mathrm{\Delta }TS`$ does not, however, account for position uncertainties associated with nearby unmodeled $`\gamma `$-ray sources and with other systematic and instrumental biases. This is confirmed by plotting the distribution of bright 3EG blazars (Mattox et al., 2001) as a function of $`\mathrm{\Delta }TS`$. Considerably fewer than half of the sources lie within the 50% contour, and several lie outside the 99% contour. We find that the additional spread can be well modeled by dividing the map value of $`\mathrm{\Delta }TS`$ by 1.5 and then computing the nominal statistical positional probability. The probabilities renormalized in this way are a good match to the observed source distribution.

### 2.3 Combined Figure of merit

We combine the fractional excesses for each of the counterpart candidate’s properties with the estimate of the positional probability to compute a total ‘Figure of Merit’ for the source as a radio/X-ray counterpart:

$$FoM=n_{8.4GHz}\times n_\alpha \times n_{xray}\times L(\alpha ,\delta )$$
(2)
Notice that this product of excess source fraction and positional uncertainty is not a normalized probability of source identification. It is however an unbiased ranking of the counterpart likelihood that can be compared across the entire (high latitude) 3EG population. To avoid confusion, we multiply by $`100`$ in the quoted FoM values.

To evaluate the significance of this FoM statistic, we produced random realizations of the radio/X-ray sky by random, independent RA and DEC draws from the true CLASS position list, after excluding the actual positions of the $`70`$ highest-FoM coincident sources \[which would otherwise imprint the observed correlation on the random realizations\]. We compute the average distribution of FoM for $`10^3`$ realizations of such random skies, and compare with the true FoM distribution in Figure 2. The hashed region shows the $`\pm 1\sigma `$ range for our estimate of the fraction of sources in a given FoM bin that are ‘real’, i.e. in excess of random coincidence. For FoM$`<0.1`$ there is little excess correlation. There is a puzzling deficit of true sources with FoM$`2`$, but we have not been able to trace this to any one of the FoM source properties and so conclude that this is a statistical fluctuation. We select sources with FoM$`0.25`$ (i.e. $`<20\%`$ false positive, even for the lowest FoM) as good counterparts; this is about a factor of two above the lowest bins showing significant correlation. To facilitate comparison with earlier ID lists, we divide these sources in half with FoM$`>1`$ designated as ‘likely’ counterparts and $`0.25<\mathrm{FoM}<1.0`$ as ‘plausible’ counterparts (note however that our ‘plausible’ sources have a relatively high fraction of true associations). The line in figure 2 shows a simple linear fit to the source probability with a strong decrease below FoM=0.1. Integrating through the distribution of measured FoM, we find that of the 35 ‘likely’ sources we expect less than 3 false positives ($`92\%`$ good IDs) and of the 32 ‘plausible’ sources $`6`$ may be false positives ($`82\%`$ good IDs). We can compare the FoM distributions of our IDs with those designated in the 3EG catalog and in Mattox et al. (2001). Many of the highest FoM sources are common to all lists. However, both the 3EG and Mattox lists claim high confidence for several sources that we only assign low probability. In addition each list has a number of sources that are not flat spectrum or compact and do not meet our interferometric selection criteria. These sources, 11 in the 3EG and 7 in Mattox et al. (2001) are mostly lower confidence (‘a’ and ‘plausible’ designations in these catalogs). However, in addition to selecting previous IDs, we also find $`50\%`$ more high confidence sources than either, and virtually all of our ‘plausible’ sources are new. Noting that Mattox et al. (2001) retained ‘plausible’ sources with estimated likelihoods as low as a few percent, while our lowest FoM sources have a likelihood $`>80\%`$, we believe that our list is less biased, more complete and more reliable than earlier efforts. In particular, we can more easily identify multiple associations with a given EGRET source, rather than taking only the ‘best’ blazar ID.

### 2.4 Radio Counterpart IDs

We find at least one blazar candidate counterpart for 66 of the Northern 3EG sources. In many cases, the previously claimed 3EG ‘high confidence’ and ‘plausible’ blazars are recovered with the largest FoM value for a given error ellipse. However, in a number of cases, our criteria select a different blazar as much more likely than the claimed catalog association. In addition, we find some 2 dozen new radio associations in the Northern sky, whose FoM confidence are at least as good as that of the lowest of the previously identified ‘high confidence’ sources. In several cases there are multiple radio blazar IDs in a single EGRET error ellipse. One rather complex example is 3EG J0118+0248. In the 3EG catalog this was possibly identified with 0119+041; our analysis does select it as a (low FoM) possible association, but it is well outside the 99% error contour. In the Mattox et al. (2001) analysis, the 3EG source was attributed to 3C37 ($`z=0.672`$); this AGN however has a quite steep spectrum for the compact component and does not meet our blazar identification criteria. Instead, our analysis selects two compact flat spectrum sources as the most likely counterparts. The first is a relatively low luminosity radio galaxy core at z=0.047, the second is a newly discovered flat spectrum radio QSO at z=4.0. The 3EG error contours are quite elongated with the major axis spanning these two sources (Figure 3), although 0119+041 is also a plausible source of this extension. A second example is 3EG J0808+5114, which hosts a flat spectrum radio quasar at z=1.14 and a BL Lac at z=0.13, again separated along the major axis of an extended $`\gamma `$-ray uncertainty region. Given that the 3EG survey is strongly flux limited, that the $`\gamma `$-ray blazar luminosity function is quite steep, and that the poor $`\gamma `$-ray resolution causes substantial source overlap and confusion, it is not surprising that in some cases a combination of fainter sources can push a location above the detection threshold. In several other cases support for multiple IDs comes from complex or elongated likelihood (‘TS’) maps.

### 2.5 Radio non-IDs

Since the GB6/NVSS/RASS provide a uniform survey (abetted by CLASS follow-up), the absence of plausible IDs for some EGRET sources is significant. For example, we find no IDs for two high latitude sources believed to be likely pulsars (3EG J0010+7300= CTA 1 and 3EG J1835+5918, see Reimer et al. 2001). Our technique allows us to select individual 3EG sources as non-blazar IDs. This is particularly useful at intermediate latitude, as it lets us separate the blazar fraction of the population, leaving a cleaner non-blazar sample. North of Dec=0, 28 previously unidentified sources are selected as not like EGRET blazars. There is some concentration to the Galactic plane, but much of this sample is at high latitude, suggesting either a new class of extragalactic $`\gamma `$-ray sources or an old high scale height Galactic populations of pulsars or similar objects.

### 2.6 Summary Results

We summarize our source classification in Table 1 and Figure 4. In summary, of the 116 Northern 3EG sources (excluding the Solar flare), 66 have at least one plausible blazar-like radio counterpart. In the Table our high confidence associations are listed in bold face and our lower confidence associations in plain text. Previously claimed AGN associations not supported by our analysis are given in italics in parentheses, with FoM values computed when possible. The Table also notes 10 sources associated with pulsars or known plerion/young SNR (indented). These are neutron star/pulsar/SNR associations from the 3EG catalog, with a few additions from new discoveries which are noted in the individual source comments. Of course many additional associations have been claimed for Galactic objects; we do not evaluate or exhaustively review these here, but we believe that beyond the few pulsed detections, the identifications rapidly become quite speculative.

An additional 28 sources are unidentified but according to our analysis are definitively unlike the known 3EG blazars. The remaining sources are at low $`|b|5^{}`$ latitude. Only faint flat spectrum sources are allowed in the error boxes of 7 of these and, given the high density of Galactic sources, these are almost certainly chance coincidence. The final 5 show bright extended HII regions which could, in principal, hide a bright flat spectrum compact blazar, but again in practice we expect the sources to be Galactic. Further radio imaging study could better rule out AGN counterparts in these cases.

Figure 4 summarizes the present 3EG source classifications in Galactic coordinates. Above the Dec=0 line, we use our new classifications from Table 1, with blazars as filled circles. In the South, the same symbols are used for the 3EG ‘A’ and ‘a’ blazar designations. Pulsar IDs and strong pulsar/plerion candidates are shown by filled and open stars. The definitive non-blazars are shown by open circles. The improved completeness of the Northern classification is evident. Clearly pushing south into the bulge population will be very interesting.

Since AGN are expected to be isotropic, it is of interest to examine the distribution of the northern blazar IDs. Detailed comparison requires folding through the EGRET exposure and sensitivity maps for a model blazar luminosity function; we defer this to a future communication. However since exposure is relatively uniform for $`|b|>20^{}`$, we report simple number counts. From $`30^{}<|b|<90^{}`$ we find 36 blazar identifications, or 11.6/sr within our survey boundary. Extrapolating toward the plane, we find for $`20^{}<|b|<30^{}`$ 12 sources against 11 expected. At $`10^{}<|b|<20^{}`$ we find 8 sources vs. 12 expected, and a similar fraction within $`10^{}`$ of the plane. Given the higher background and detection threshold near the plane, this does not however mean that our identification of blazars at low $`|b|`$ is incomplete.

## 3 Optical Follow-up

Of our radio-selected counterparts, 50 are new sources, not selected in previous EGRET counterpart lists. Further, several additional selected radio sources, while flagged as possible counterparts in earlier studies, have not had spectroscopic confirmation of their AGN nature. We first checked for published spectroscopy of these radio sources, comparing with the Tenth Edition of the Quasar Catalog (Veron-Cetty et al., 2001) and cross checking by querying NED. We find that 16 (5 high probability IDs) of our newly selected counterparts have known redshifts. However over thirty sources in the Northern sky were found to have previously unstudied blazar-like radio/X-ray candidate counterparts, many with quite high FoM. These are the target of our HET spectroscopy, along with a number of similar sources having no EGRET association.

Just before submission, we learned that Halpern, Eracleous & Mattox (2003) have recently completed follow-up spectroscopy of a set of previously claimed radio counterparts of 3EG (and 2EG, GEV) sources with lower 4.8GHz fluxes. Their list includes 8 sources also identified in our exercises. For three sources we have new spectroscopy in common; for two their new redshift agrees with that which we obtain (below). We were also able to obtain a redshift for the third source (J1826+0149). For three other sources Halpern, Eracleous & Mattox (2003) obtain no $`z`$, for two the $`z`$ they measure was already found in the QSO catalog, and for the last (J1239+0443) they obtain a new $`z`$. This value is included in our table as the only entry not from cataloged data or our new spectroscopy.

### 3.1 HET LRS Spectroscopy

Spectroscopy was obtained using the Marcario Low Resolution Spectrograph (LRS)(Hill et al., 1998) on the 9.2m HET (Ramsey et al., 1998). This novel telescope, of tilted optical Arecibo design is still under development, but with a range of magnitudes the blazar candidates were suitable targets for the early operations phase. These targets were observed in regular queue observations from 3/02 -2/03. We obtained 2$`\times `$450s exposures for most targets; a few of the most interesting fainter targets were observed with $`2\times 900`$s exposure. Observations were made employing a 300 l/mm grating and a 2<sup>′′</sup> slit for a dispersion of 4$`\AA `$ per (binned) pixel and an effective resolution of 16$`\AA `$ covering $`\lambda \lambda 420010500`$. The new spectra are displayed in Figure 5.

Standard IRAF CCD reductions, optimal extraction and calibration were performed. Unfortunately, the varying pupil during the HET tracks is not yet fully monitored and so small (few %) flat-field variations were not corrected, as with untracked flat field frames the varying illumination worsened the (already substantial) fringe features and sky subtraction errors in the near IR. Redshifts were measured by cross correlation analysis with AGN and galactic templates, using the IRAF RVSAO package. A few identifications are based only on the MgII $`\lambda `$2798 line, but in each case the proposed $`z`$ is the only plausible value, given the absence of other expected strong lines. Several of our objects are continuum dominated with very faint line features, for these the substantial light grasp of the HET was important. All of the observed targets were found to be AGN, although some spectra were taken under rather poor conditions. Some of these have relatively uncertain $`z`$ estimates, these are noted by (:) in Table 2.

As of this publication we have obtained 29 new spectroscopic identifications and an additional 16 new associations with previous spectroscopy. The new and archival redshifts for our new blazar IDs are given in Tables 1 and 2. The later table also contains precise positions and sky survey flux estimates; as expected for blazars, the flux at our observation epoch was often substantially different. Most of the AGN were found to be flat spectrum radio quasars. We denote sources with emission line equivalent width $`5\AA `$ as BL Lacs. Several sources have the narrow lines (FWHM$`<1000`$km/s) and lower excitation states characteristic of radio galaxies, we denote these as ‘R’ in Table 1. It should be noted that within the classical AGN paradigm, narrow-line radio galaxies would not be expected to show blazar activity. The redshift distribution of our identified sources is shown in Figure 6. For comparison, we also plot the redshift distributions of the sources selected by the 3EG and Mattox et al. (2001) criteria. The $`z`$ distributions are rather similar, largely displaying spectroscopic ID selection effects, but it is clear that we have substantially higher sensitivity to high redshift and have nearly doubled the maximum $`z`$ value.

## 4 Notes on Individual Sources

3EG J0118+0248 As noted above, the Mattox et al. (2001) identification with the bright steep spectrum source 3C37 is not supported by our analysis and the 3EG association with 0119+041, is apparently superseded (or at least augmented) by our discovery of two flat spectrum candidates along the error ellipse major axis. If these contribute similarly to the $`\gamma `$-ray flux, they would represent extrema of the $`\gamma `$-ray blazar population, with a factor of $`10^4`$ difference in luminosity. If the association is confirmed, J0122+0310 at z=4.0 would be the highest redshift non-GRB $`\gamma `$-ray source known. This 3EG source will be particularly interesting for study in the GLAST era, since with higher angular resolution the contributions of the several posited associations can be disentangled and GLAST’s good sensitivity above 10GeV should allow intergalactic absorption of J0122+0310 against a cosmic background of starlight to be detected at $`20`$GeV.

3EG J0204+1458 In addition to the likely association J0204+1514 this source may also incorporate flux from the flatter and more centrally placed radio source J0205+1444.

3EG J0215+1123 with a reasonably large FoM J0213+1213 is an interesting low $`z`$ BL Lac candidate for TeV studies.

3EG J0222+4253 This object has been argued to be a composite by Kuiper et al. (2000) with a possible pulsed detection of the millisecond pulsar PSR J0218+42 accounting for part of the sub-GeV flux. The remainder was attributed (also in the 3EG catalog) to the BL Lac 3C 66A. This source, with $`\alpha _{1.4/8.4}=0.62`$, lies just below the spectral cut-off for inclusion in our FoM. However, 2.3 GHz VLBA data are available (Fey & Charlot, 2000) which give a much flatter index for the compact component, which we adopt in Table 1. This does make 3C 66A highly likely. CLASS analysis also selects the 215mJy source J0223+4259 (3C 66B, a radio galaxy with a core+jet at $`z=0.02`$) as a plausible counterpart. Confusion in the NVSS image prevents detection of the core, so using the 15GHz core flux measured by Jackson, et al. (1993) we obtain $`\alpha _{8.4/15}=0.38`$ and a strong plausible association. The GeV localization of Kuiper et al. (2000) supports either 3C object as the source of the higher energy $`\gamma `$-rays. Certainly GLAST observations will be helpful sorting out this complex region. Interestingly, a TeV detection is claimed in this region Neshpor et al. (1998) which would be somewhat surprising from $`z`$ as large as 0.44. The apparent position is consistent with either 3C 66A or 3C 66B and future air Cerenkov observations may be able to select the preferred source, given their $`0.1^{}`$ separation.

3EG J0245+1758 The source closest to the TS maximum is a strong-lined FSRQ at z=3.59. This source shows a likely damped Ly$`\alpha `$ system at $`z3.15`$, raising the possibility that the source is gravitationally magnified by an intervening Galaxy.

3EG J0329+2149 Mattox et al. (2001) noted the radio source as a possible association. We confirm it as a likely blazar counterpart and have measured z=2.07, in agreement with Halpern, Eracleous & Mattox (2003).

3EG J0404+0700 The better of the two likely associations is our most puzzling case, spectroscopically. When first observed, the spectrum was highly continuum dominated. The two strongest lines were picked up as Ly$`\alpha `$ and CIV $`\lambda 1449`$ at $`z=3.13`$ in cross-correlation analysis. Both are (local) $`5\sigma `$ detections; they have rest equivalent widths of $`0.44`$Å and $`0.39`$Å, respectively. CIII $`\lambda 1909`$ is obscured due to the fringing and no other strong lines are expected. This redshift ID is supported by an apparent strong MgII absorption doublet at $`z=1.578`$ (REW 0.19/0.22$`\AA `$). With such low EW for the emission lines this source is (in this state) clearly a BL Lac (Marcha et al., 1996). However a second observation several months later caught the source in a low state for which broad MgII line and several Ne and O lines give a highly significant measurement of $`z=1.133`$. The MgII line is in fact also present in the high state at identical flux. The difference spectrum shows no MgII line, but does have the higher $`z`$ features noted above. We conservatively adopt the lower redshift, as $`z=3.13`$ more than doubles that of the highest $`z`$ BL Lac cataloged to date (PKS 1309-216 $`z=1.49`$ Blades, Murdoch & Hunstead 1980). Further spectroscopy in a variety of flux states is certainly warranted.

3EG J0450+1105 Halpern, Eracleous & Mattox (2003) have observed J0449+1121, identified in the 3EG catalog as the likely counterpart, and have not confirmed the catalog $`z`$. Their nearly featureless spectrum suggests a BL Lac.

3EG J0459+3352 This Galactic error region has at least two maxima. The blazar candidate J0503+3404 is well associated with one peak, but the second peak, nearly 2 degrees away is more likely Galactic.

3EG J0634+0521 This source has been plausibly associated with the Be X-ray binary SAX J0635+0533 (Kaaret, 2001), so it might be classified as ‘p’.

3EG J0808+5114 This source is probably composite with two likely counterparts (Figure 3).

3EG J0917+4427 The bright, 1.3Jy, flat spectrum association is only plausible, being just outside the 99% contour. However, this error region is highly elongated and source is offset in the short direction. If an unmodeled systematic of $`0.3^{}`$ induces this offset, this would be a likely association.

3EG J1133+0033 and 3EG J1329+1708 we have observed a plausible radio source in each, identifying them as featureless BL Lacs, but have not obtained redshifts.

3EG J1605+1553 Our new identification at z=0.11 is brighter at 8.4GHz and has a smaller $`\mathrm{\Delta }TS`$ than the association proposed by Mattox et al. (2001).

3EG J1621+8203 Mukherjee et al. (2002) have previously argued for identification with this Seyfert.

3EG J1835+5918 This source has been the subject of intense study (Reimer et al., 2001; Halpern et al., 2002), who have both argued that it is an isolated Geminga-like pulsar.

3EG J2016+3657 This crowded Galactic error region contains the SNR CTB 87, several bright HII regions in addition to the radio source B2013+370, claimed as a possible BL Lac (Halpern et al., 2001). Our spectrum of this source, taken under very poor conditions, is sufficient to support the continuum-dominated nature of the heavily absorbed counterpart, but is insufficient to determine a redshift. We tentatively adopt the blazar designation of these authors, but confirmation requires improved angular resolution and possibly variability correlation with GLAST.

3EG J2021+3651 Is likely associated with the newly discovered pulsar and plerion discussed in Roberts et al., (2002); Roberts, Romani & Kawai (2001).

3EG J2035+4441 May be associated with an X-ray/radio plerion (Roberts, Romani & Kawai, 2001).

3EG J2036+1132 Our analysis produces FoM=0.22 for B2032+107 (a $`z=0.6`$ BL Lac), the 3EG A designated counterpart of this source. This is just below our cutoff, but as it is significantly less than the FoM of J2034+1154, we propose re-assignment to this source. Our second source J2031+1219 is bright and flat, but lies significantly outside the 99% confidence contour. A re-evaluation of the likelihood contours in terms of multiple sources could affect the rankings of these plausible counterparts.

3EG J2206+6602 The 3EG catalog also selected this source with low probability; we confirm it as a $`z=1.12`$ FSRQ.

3EG J2209+2401 This source was claimed as an identification in (Mattox et al., 2001), but the listed redshift was evidently a typographical error, repeating that of another source. We have measured a single line $`z`$ that excludes the published result.

3EG J2227+6122 This source is likely identified with the recently discovered X-ray/radio pulsar PSR J2229+6114 and its associated wind nebula (Halpern et al., 2001b).

## 5 Conclusions and Implications

We believe that we have made a relatively unbiased assessment of the association of flat spectrum radio sources with EGRET AGN, identifying plausible counterparts down to 8.4GHz fluxes of $`100`$mJy. Assuming, as we have argued, that the bulk of our likely and plausible IDs are correct, we have substantially increased the completeness of identification of the Northern $`|b|10^{}`$ EGRET sources from $``$40% to $`70`$%. We have also argued that in a number of instances the 3EG sources are composites; including this, the increase in the number of proposed $`\gamma `$-ray blazars is even larger.

Our identification allows selection of sources with substantially smaller 8.4GHz radio flux. It should not be too surprising that the number of candidates tapers off smoothly below the previous typical 1Jy value. We suspect that with smaller GLAST error ellipses or improved multi-wavelength constraints, identification could continue well below the $`0.1`$Jy limit of our survey. Our identification of $`30`$ sources that have flat spectrum counterparts (if any) well below 0.1Jy does not, of course, preclude that some of these may be radio-faint or steep spectrum AGN. Indeed, excess coincidences do continue slightly below our FoM identification limit (Figure 2). Also a number of the 3EG and Mattox et al. (2001) identifications not duplicated here are steep spectrum, but very bright low $`z`$ AGN. In these cases (with the exception of 3EG J0118+0248/3C037, discussed above) the steep spectrum proposed counterpart is located at large $`\mathrm{\Delta }TS`$. Nevertheless, if such sources are truly associated, it is reasonable that these form a distinct subset of the 3EG population, not identified in this analysis. However, the fact that an appreciable fraction of the ‘non-blazar’ 3EG sources correlate with the Galactic disk suggests that many represent a new (old) Galactic population. Detailed assessment of the completeness of blazar IDs through the plane requires careful treatment of the 3EG exposure maps; we defer this and other population analysis to a later paper.

Along with fainter radio associations, we have pushed back the horizon of plausible blazar identifications (Figure 4). It will be interesting to compare the SED and VLBI $`\beta _{}`$ properties of these fainter blazar counterpart candidates. This survey has found several individual targets of particular interest, including good $`\gamma `$-ray source associations at $`z=34`$, several distant BL Lacs, and several low $`z`$ AGN that may be suitable targets for ground-based TeV observations. These targets should be the blazars brightest to GLAST and they will provide the most detailed spectra and light curves; as such they merit further study in preparation for the GLAST era. Since spectroscopic identifications are continuing, we defer detailed discussion of the red shift/ luminosity function distributions to a future paper.

The ‘non-blazars’ are also excellent targets for further study in preparation for GLAST, since this sample will most likely produce new classes of high energy emitters. While we do not discuss here the correlation of $`\gamma `$-ray spectra and variability with the lower energy SEDs, we do note that our ‘non-blazars’ correlate fairly well with the ‘steady’ sources of Gehrels, et al. (2000) and the ‘persistent’ sources of Grenier (2000), although exceptions occur. Clearly, pushing the improved identifications south to cover the Galactic bulge and beyond will be very important for characterization of these populations.

We thank Roger Blandford for helpful discussions on the nature of $`\gamma `$-ray blazar sources and for assistance with access to the CLASS source list. Ian Browne and Thomas York are also thanked for passing along an early copy of these important data. The referee, R.C. Hartman, and S. Digel provided careful readings and comments that substantially improved the manuscript. We are particularly grateful to the HET team, especially Gary Hill for bringing the Marcario LRS into operation and Matt Shetrone, Jeff Mader, Brian Roman and the rest of the telescope team for getting these early observations done. We dedicate these HET results to the late Stanford provost Gerald J. Lieberman, whose early support of Stanford’s participation in the HET made this project possible. The Hobby-Eberly Telescope is operated by McDonald Observatory on behalf of The University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. The Marcario Low Resolution Spectrograph is a joint project of the Hobby-Eberly Telescope partnership and the Instituto de Astronomia de la Universidad Nacional Autonoma de Mexico. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. DSE was supported by SLAC under DOE contract DE-AC03-76SF00515 and PFM acknowledges support from NASA contract NAS5-00147.


format: text

text:
# References

On the Difference of Energy between the Einstein and Møller Prescription
I-Ching Yang 111E-mail:icyang@cc.ntttc.edu.tw and Irina Radinschi 222E-mail:iradinsc@phys.tuiasi.ro

Department of Natural Science Education and

Advanced Science and Technology Research Center,

National Taitung Teachers College, Taitung, Taiwan 950, Republic of China

and

Department of Physics, “Gh. Asachi” Technical University,

Iasi, 6600, Romania

ABSTRACT

In some black hole solutions, these do not exist the same energy-momentum complexes associated with using definition of Einstein and Møller in given coordinates. Here, we consider the difference of energy between the Einstein and Møller prescription, and compare it with the energy density of those black hole solutions. We found out a special relation between the difference of energy between the Einstein and Møller prescription and the energy density for considered black hole solutions.

PACS No.:04.20.-q, 04.50.+h

In the theory of general relativity, many physicists, like Einstein , Landau and Lifshitz , Tolman , Papapetrou , Møller , and Weinberg , had given different definitions for the energy-momentum complex. Specifically, the Møller energy-momentum complex allows to compute the energy in any spatial coordinate system. Some results recently obtained  sustain that the Møller energy-momentum complex is a good tool for obtaining the energy distribution in a given space-time. Also, in his recent paper, Lessner  gave his opinion that the Møller definition is a powerful concept of energy and momentum in general relativity. In his paper Virbhadra  point out that several energy-momentum complexes (ELLPW) give the same result for a general non-static spherically symmetric space-time of the Kerr-Schild class.

In particular, whatever coordinates do not exist the same energy complexes associated with using definitions of Einstein and Møller in some space-time solutions . According to the definition, the Einstein energy complex is

$$E_{\mathrm{Ein}}=\frac{1}{16\pi }\frac{H_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}l}}{x^l}d^3x,$$
(1)
where

$$H_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}l}=\frac{g_{00}}{\sqrt{g}}\frac{}{x^m}[(g)g^{00}g^{lm}],$$
(2)
and the Møller energy complex is

$$E_{\mathrm{M}\mathrm{ø}\mathrm{l}}=\frac{1}{8\pi }\frac{\chi _0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}l}}{x^l}d^3x,$$
(3)
where

$$\chi _0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}l}=\sqrt{g}g^{0\beta }g^{l\alpha }\left(\frac{g_{0\alpha }}{x^\beta }\frac{g_{0\beta }}{x^\alpha }\right).$$
(4)
Where the Latin indices take values from 1 to 3, and the Greek indices run from 0 to 3. Let us look into the difference of energy between the Einstein and Møller prescription, which be defined as

$$\mathrm{\Delta }E=E_{\mathrm{Ein}}E_{\mathrm{M}\mathrm{ø}\mathrm{l}}$$
(5)
In this article, we would discuss the problem within the difference between Einstein and Møller energy-momentum complexes.

In the first case, we think of two solutions of Einstein vacuum field equation:
(i) Schwarzschild space-time
The metric form of Schwarzschild space-time is

$$ds^2=fdt^2f^1dr^2r^2d\theta ^2r^2\mathrm{sin}^2\theta d\phi ^2,$$
(6)
where $`f=12M/r`$. It is a well-known results that the energy complexes of Einstein and Møller of Schwarzschild space-time are

$`E_{\mathrm{Ein}}`$ $`=`$ $`M,`$ (7)
$`E_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $`=`$ $`M,`$ (8)
and the difference is

$$\mathrm{\Delta }E=0.$$
(9)
(ii) Kerr solution
The metric form of Kerr solution is considered as

$$ds^2=\alpha dt^2\beta dr^2\gamma d\theta ^2\delta d\varphi ^22\sigma dtd\phi ,$$
(10)
where $`\alpha =12Mr/\mathrm{\Sigma }`$ , $`\beta =\mathrm{\Sigma }/\mathrm{\Delta }`$ , $`\gamma =\mathrm{\Sigma }`$ , $`\delta =r^2+a^2+2Ma^2r\mathrm{sin}^2\theta /\mathrm{\Sigma }`$ and $`\sigma =2Mar\mathrm{sin}^2\theta /\mathrm{\Sigma }`$ . Here $`\mathrm{\Sigma }r^2+a^2\mathrm{cos}^2\theta `$ and $`\mathrm{\Delta }r^22Mr+a^2`$. To use the results in the Virbhadra articles  and to set these $`Q=0`$, we could obtain the enrgy-momentum complexes of Einstein and Møller of Kerr space-time are

$`E_{\mathrm{Ein}}`$ $`=`$ $`M,`$ (11)
$`E_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $`=`$ $`M,`$ (12)
and the difference is

$$\mathrm{\Delta }E=0.$$
(13)
For Einstein’s vacuum field equation, the energy density is

$$T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}=0.$$
(14)
We would find that $`\mathrm{\Delta }E`$ equal to the value of $`T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}`$.

Next, we consider two case of the coupled system of the Einstein field and electromagnetic field:
(iii) Reissner-Nordström space-time
The metric form of Reissner-Nordström space-time is

$$ds^2=fdt^2f^1dr^2r^2d\theta ^2r^2\mathrm{sin}^2\theta d\phi ^2,$$
(15)
where $`f=12M/r+Q^2/r^2`$. Previously, the energy-momentum complexes of Einstein and Møller of Reissner-Nordström space-time had been calculated with

$`E_{\mathrm{Ein}}`$ $`=`$ $`M{\displaystyle \frac{Q^2}{2r}},`$ (16)
$`E_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $`=`$ $`M{\displaystyle \frac{Q^2}{r}},`$ (17)
and the difference is

$$\mathrm{\Delta }E=\frac{Q^2}{2r}.$$
(18)
Notice that the energy density in the Einstein-Maxwell field equation is

$$T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}=\frac{Q^2}{r^4}.$$
(19)
(iv) charged regular black hole
The metric form of charged regular black hole is

$$ds^2=fdt^2f^1dr^2r^2d\theta ^2r^2\mathrm{sin}^2\theta d\phi ^2,$$
(20)
where $`f=1\frac{2M}{r}(1\mathrm{tanh}(\frac{Q^2}{2Mr}))`$. Using the results of Radinschi’s articles , the energy-momentum complexes of Einstein and Møller of charged regular black hole are

$`E_{\mathrm{Ein}}`$ $`=`$ $`M\left[1\mathrm{tanh}({\displaystyle \frac{Q^2}{2Mr}})\right],`$ (21)
$`E_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $`=`$ $`M\left[1\mathrm{tanh}({\displaystyle \frac{Q^2}{2Mr}})\right]{\displaystyle \frac{Q^2}{2r}}\left[1\mathrm{tanh}^2({\displaystyle \frac{Q^2}{2Mr}})\right],`$ (22)
and the difference is

$$\mathrm{\Delta }E=\frac{Q^2}{2r}\left[1\mathrm{tanh}^2(\frac{Q^2}{2Mr})\right].$$
(23)
However, the energy density of the coupled system of the Einstein field and nonlinear electrodynamics field is

$$T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}=\frac{Q^2}{r^4}\left[1\mathrm{tanh}^2(\frac{Q^2}{2Mr})\right].$$
(24)
Here the relation between $`\mathrm{\Delta }E`$ and the energy density is written as

$$\mathrm{\Delta }E=T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}\times (\frac{r^3}{2}).$$
(25)
Then, in the last case, we use a special black hole solution without singularity:
(v) static spherically symmetric nonsingular black hole
The metric form of the static spherically symmetric nonsingular black hole solution  is

$$ds^2=fdt^2f^1dr^2r^2d\theta ^2r^2\mathrm{sin}^2\theta d\varphi ^2,$$
(26)
where $`f=1R_g/r`$ and $`R_g=r_g(1\mathrm{exp}(r^3/r_{}^3))`$ , $`r_{}^3=r_gr_0^2`$ ,$`r_g=2M`$ , $`r_0^2=\frac{3}{8\pi \epsilon _0}`$ . According to the results of our articles , the energy-momentum complexes of Einstein and Møller of the static spherically symmetric nonsingular black hole are

$`E_{\mathrm{Ein}}`$ $`=`$ $`MM\mathrm{exp}({\displaystyle \frac{r^3}{r_{}^3}}),`$ (27)
$`E_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $`=`$ $`MM\mathrm{exp}({\displaystyle \frac{r^3}{r_{}^3}}){\displaystyle \frac{3r^3}{r_0^2}}\mathrm{exp}({\displaystyle \frac{r^3}{r_{}^3}}),`$ (28)
and the difference is

$$\mathrm{\Delta }E=\frac{3r^3}{r_0^2}\mathrm{exp}(\frac{r^3}{r_{}^3}).$$
(29)
Notice that the energy density of the static spherically symmetric nonsingular black hole be assumed as

$$T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}=\frac{3}{r_0^2}\mathrm{exp}(\frac{r^3}{r_{}^3}).$$
(30)
The relation between $`\mathrm{\Delta }E`$ and the energy density is written as

$$\mathrm{\Delta }E=T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}\times r^3.$$
(31)
Although, we could summarize that the general relation between $`\mathrm{\Delta }E`$ and the energy density $`T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}`$ be written as

$$\mathrm{\Delta }E=T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}\times (kr^3),$$
(32)
with $`k=1/2`$ and $`k=1`$. But, it is still an open question why the special relation has between $`\mathrm{\Delta }E`$ and the energy density $`T_0^{\mathrm{\hspace{0.33em}\hspace{0.17em}0}}`$. Further study is needed to understand the difference between the Einstein and Møller energy complexes of more varied black hole solutions.

Acknowledgements

I.-C. Yang thanks the National Science Council of the Republic of China for financial support under the contract number NSC 90-2112-M-143-003.


text:
# On particle acceleration around shocks. I

## 1 Introduction

The theory of particle acceleration at shocks is currently unsatisfactory. We cannot follow this process from first principles, i.e., the development of a collisionless shock, with the ensuing injection of non-thermal particles, nor do we understand the later shock evolution, which should include the interaction between the population of high–energy, non–thermal particles, and the shock structure. Even in the test–particle limit, to which this paper is totally confined, there are at least two major problems. First, the properties of the scattering agents are very poorly known; this problem will not be dealt with in this paper. Second, the general foundations of the theory have not been laid down. In particular, there is no general equation for the particle distribution function (DF) for arbitrary (i.e., even relativistic) shock speed, and no general treatment is available when the DF is anisotropic. Most results available to us either have been derived in the vanishing speed (thus fully isotropic) limit (Bell 1978, Blandford and Ostriker 1978), or descend from numerical (Bednarz and Ostrowski, 1998) or semi-numerical (Kirk and Schneider 1987) approaches. This is true even in the idealized situation of a plane shock expanding forever in a uniform medium,with particle scattering due only to pitch–angle scattering. It is the aim of this paper to derive a sufficiently general equation describing in a relativistically covariant way this process, to elucidate the conditions leading to the determination of the particle spectrum, and to apply these results to the extremely relativistic limit, where the DF anisotropy is expected to be maximal. In passing, a number of old results will be recovered: it will be shown that the probability of reaching infinity is independent of the particle momentum, that the spectrum is the superposition of many bumps, each corresponding to a set of particles which have crossed the shock $`0,1,2,\mathrm{},N,..`$ times, and that, for very large particle energy, the superposition of these bumps spectrum leads to a power law in the particle impulse: all these results will be shown to hold for all shock speeds, Newtonian, relativistic or intermediate ones. Also, the explicit dependence upon the injection spectrum will be presented. In a future paper, this formalism will be applied to the extremely relativistic limit, including the effects of small or large angle scatterings, and a mean magnetic field.

The plan of this paper is as follows. In Section 2, we derive the new equation for the particle DF, in an explicitly relativistic covariant way. In Section 3, we reformulate the boundary conditions to which this problem is subject, in a physically motivated way, finding eventually the dependence of the distribution function upon the injection spectrum. In Section 4, we concentrate on the distribution function at energies large compared with those of injection; here we show that it must be a power law, showing what fixes the power law index analytically, and deriving the relativistic generalization of Bell’s law. We also make contact with the Newtonian limit. The important probability distributions $`P_d`$ and $`P_u`$, whose existence is merely postulated in the previous sections, are explictly written out in terms of eigenfunctions of the angular part of the scattering equation in Section 5. In Section 6, we summarize our results.

## 2 An equation for the distribution function

We consider a shock propagating at arbitrary speed in an homogeneous fluid. We place ourselves in the shock reference frame, and choose coordinates such that $`z=0`$ identifies the shock position. The upstream fluid is located at $`z<0`$, so that the fluid speed is always $`>0`$, both upstream and downstream. Suitable jump conditions at the shock are assumed (Landau and Lifshitz 1987), but they will not be necessary in this section.

We assume now that we can break down the effects of the magnetic field into two parts, as is customarily done. The first part is due to a small–coherence length magnetic field, possibly self–excited by the particles; this provides for an effective scattering, which is included in the collisional term. Another part, however, is due to a long–coherence length field, which provides instead for a smooth deflection of particles in phase space; this effect is included in the convective term (see below). We are thus using a dichotomic description of the magnetic field; most of the results to be described here, which are wholly unavoidable within this description, ultimately depend on the correctness of this assumption on the magnetic field.

The distribution function is of course subject to the collisionless Vlasov equation, supplemented with suitable collisional and injection terms. In somewhat abstract form, we have

$$\frac{df}{ds}=\left(\frac{f}{s}\right)_{coll}+\varphi _{inj}.$$
(1)
The right hand side is the collisional term, with collisions per unit proper time, plus the rate of particles injected at the shock.

The convective term is

$$\frac{df}{ds}\frac{dx^\nu }{ds}\frac{f}{x^\nu }+\frac{dp^\nu }{ds}\frac{f}{p^\nu }.$$
(2)
In the above equation, the term containing the derivative with respect to the energy is not usually explicitly written out, imposing instead ab initio that particles be on their mass shell; here we do this instead (showing later on that it makes no difference) in order to emphasize the covariant character of the equation. In fact, both the first and the second term on the right hand side are invariant for Poincaré group transformations. We make use of this by computing each term in different frames of reference, and with different orientations for the coordinate axes. In the first term, we place ourselves in the shock frame and orient the $`z`$ axis along the shock normal, so that we can make use of $`/t=/x=/y=0`$, because we assume the shock to be in steady state, and to have planar symmetry. We also write the $`z`$ component of the particle four–speed in terms of quantities in the fluid frame,

$$\frac{dz}{ds}=\gamma _p\gamma (u+\mu ),$$
(3)
where $`\mu `$ is the cosine of the angle which the particle speed makes with the shock normal, $`z`$, in the fluid frame, $`u`$ and $`\gamma `$ are the shock speed (in units of $`c`$) and Lorentz factor speed with respect to the fluid, and $`\gamma _p`$ is the particle Lorentz factor also with respect to the fluid. All particles are here assumed relativistic in any reference frame (so that their three-speeds are always $`c`$), even those which are just injected. Strictly speaking, this is accurate only for relativistic shocks. The generalization to include Newtonian particles, though easy in principle to do, introduces some cumbersome formulae (with radicals) which complicate this work uselessly: it will not be considered any further.

The second term on the right hand side of Eq. 2 will be instead computed in the fluid frame, so that we may neglect any motional electric field. In this case, the particle energy is conserved, and the term including the derivative $`f/E`$ disappears. The remaining components can be written as (Berezinsky et al., 1990, Ch.IX, Eq. 9.16)

$$\frac{dp^\nu }{ds}\frac{f}{p^\nu }=\frac{dt}{ds}\frac{dp^\nu }{dt}\frac{f}{p^\nu }=\gamma _p\omega \frac{f}{\varphi }$$
(4)
where $`\omega =eB/E`$ is the particle Larmor frequency in the fluid frame, and $`\varphi `$ is the longitudinal angle around the direction of the magnetic field, in momentum space.

To estimate the collisional term, we first place ourselves in the fluid frame, and then remark that, in full generality, such term has the form

$$\left(\frac{f}{t}\right)_{coll}=d(\mu )f+w(\mu ,\mu ^{})f(\mu ^{})𝑑\mu ^{}.$$
(5)
Here, we have made use of the assumption of linearity: we assume, that is, that the magnetic field turbulence responsible for the scattering is independent of the particle density. We consider particle scattering when there may be a correlation between the directions before and after the scattering event, so that the scattering coefficient $`w(\mu ,\mu ^{})`$ depends on both. Also, we are not assuming the diffusion approximation: the above equation applies to arbitrarily few, and arbitrarily large–angle, scattering laws. It would reduce to the well–known Ginzburg–Syrovatskii (1964) form in the diffusive approximation. The first term on the right hand side represents scattering of particles away from their direction of motion, while the second term represents the total contribution to the distribution function in the direction of motion under consideration, by scattering from all other directions. Clearly, $`d`$ and $`w`$ are related through

$$d(\mu )=_1^{+1}w(\mu ^{},\mu )𝑑\mu ^{}$$
(6)
which simply expresses probability conservation. $`w`$ (and thus also $`d`$) may depend upon $`p`$. Often, one makes the simplifying assumption that the $`p`$–dependence in $`w`$ (and thus in $`d`$) factors out; we shall comment later on why this may come in handy, but will not assume this. We also remark that $`w,d>0`$.

Since we need $`f/s`$, while above (Eq. 5) we have evaluated the collisional term in the fluid frame, we correct it by means of $`dt/ds=\gamma _p`$. Also, we wish to include the injection term, $`\varphi _{inj}`$. We assume that injection takes place at the shock only, but we leave the energy and direction dependence of the term otherwise arbitrary: $`\varphi _{inj}=G_{inj}(p,\mu )\delta (z)`$ where $`\delta (x)`$ is Dirac’s delta.

Putting together all of the above, we find

$$\gamma (u+\mu )\frac{f}{z}=\left(d(\mu )f+w(\mu ,\mu ^{})f(\mu ^{})𝑑\mu ^{}\right)+\omega \frac{f}{\varphi }+G_{inj}(p,\mu )\delta (z),$$
(7)
which is the equation we were searching for. It is worth remarking that all terms in the above equation are computed in the fluid frame, except for the space dependence, and its attaining derivative.Also, the above equation is valid both upstream and downstream, but now care must be taken because the quantities $`p,\mu ,\mu ^{},\omega ,\mathrm{}`$ are all evaluated in the local fluid frame, and thus quantities bearing the same name in the two distinct reference frames are not identical: they are related instead by a Lorentz transformation, which we give here for future reference. Using subscript $`u`$ and $`d`$ to denote quantities in upstream and downstream frame respectively, he have

$$p_u=p_d\gamma _r(1u_r\mu _d);\mu _u=\frac{\mu _du_r}{1u_r\mu _d},$$
(8)
where $`u_r`$ and $`\gamma _r`$ are the modulus of the relative speed between the upstream and downstream sections, and its associated Lorentz factor.

Whenever the scattering in each individual event is small, it is possible to apply to the above equation the very same treatment which leads to a Fokker–Planck equation. In this case one obtains (Kirk and Schneider 1987)

$$\gamma (u+\mu )\frac{f}{z}=\frac{}{\mu }\left(D(\mu ,p)(1\mu ^2)\frac{f}{\mu }\right)+\omega \frac{f}{\varphi }+G_{inj}(p,\mu )\delta (z),$$
(9)
where the effective diffusion coefficient is related to the r.m.s. angle of diffusion. It is worth remarking that it may occur that the above equation holds for the downstream section, but it is not justified for the upstream section, where the full Eq. 7 has to be employed. The reason is that the diffusive approximation on which the above equation is based requires each scattering event to lead to a r.m.s. deflection $`𝒪(1/\gamma )`$, where $`\gamma `$ is the fluid Lorentz factor with respect to the shock. When the shock becomes extremely relativistic, thus, the correctness of the diffusive approximation breaks down: it follows that one may use Eq. 9 for the downstream section, but is forced to use Eq. 7 for the upstream one.

In the following, we drop the term due to the long–coherence length magnetic field, which will be considered in Paper II.

## 3 Obtaining the particle spectrum

We begin by remarking that in this problem there is a net mean flux of particles across any surface; if the surface is steady in the shock frame, then the flux is also time–independent. Since the particles’ speed is always assumed $`c`$ ($`1`$ in our units), the infinitesimal flux across a unit shock area, per unit time, is given by

$$dJ=\mu _sp_s^2fd\mu _sdp_s=\gamma _u(u+\mu _u)p_u^2fd\mu _udp_u=\gamma _d(u_d+\mu _d)p_d^2fd\mu _ddp_d$$
(10)
depending upon whether we choose to express $`f`$ in terms of variables in the shock, upstream, or downstream frames. The representation in terms of downstream variables is especially useful, because in the downstream section the downstream momentum $`p_d`$ is conserved, so that a net flux across the shock of particles of impulse $`p_d`$ can be due only to the fact that some fraction of these particles is advected to downstrean infinity, with probability given by Eq. 21.

From now on, in this section, we shall use coordinates in the downstream frame exclusively, which we thus indicate, for ease of notation, without subscripts.

The mathematical treatment of this problem must be physically motivated: we envision a system where particles are injected at low energies at the shock, from which they cannot reach upstream infinity, because they would have to swim upstream, but from which they do reach downstream infinity. In order to study how particles generated at the shock reach downstream infinity, we first concentrate on the downstream part of the problem.

Consider now a downstream section of finite length $`L`$, beginning at the shock; then place two observers at either end, and have them inject particles into the downstream section. Then they may collect outgoing particles at either end. The problem is obviously physically well defined. The speed of a particle (in downstream frame coordinates) with respect to the shock is

$$v_s=\frac{u_d+\mu }{1+u_d\mu }.$$
(11)
Thus, the observer located at the shock injects particles for $`\mu >u_d`$ (here $`u_d`$ is the downstream fluid speed with respect to the shock), and collects the outgoing ones (i.e., those with $`\mu <u_d`$). So reasonable boundary conditions at the shock can be posed for $`\mu >u_d`$ only. The observer at the other end of the downstream section disappears as $`L\mathrm{}`$, so that there can be no entering flux there; we can only impose the regularity condition that $`f`$ does not grow to infinity, nor goes to zero.

Similarly, for the upstream section we see that we can impose boundary conditions for $`\mu <u_d`$ at the shock; at upstream infinity, the only reasonable boundary condition is that $`f0`$, because we cannot expect our particles to swim against the fluid advection all the way to upstream infinity.

We now see that the boundary conditions for the upstream problem are provided by the outgoing particles of the downstream section, and vice versa the outgoing particles of the upstream section provide the boundary conditions of the downstream section, thusly setting up an obvious (fixed point) problem. In order to make this explicit, let us introduce the conditional probability $`P_d(\mu _{in},\mu _{out})d\mu _{out}`$ that a particle, given that it crossed the shock toward downstream along a direction $`\mu _{in}`$, will recross it along the direction $`\mu _{out}`$. Then, calling $`J_{in}=(u_d+\mu _{in})f(\mu _{in})`$ the entering flux, the outgoing one will be given by

$$J_{out}(\mu _{out})=_{u_d}^1𝑑\mu _{in}P_d(\mu _{in},\mu _{out})J_{in}(\mu _{in})+G(p,\mu _{out})P_dJ_{in}+G_{out}$$
(12)
where we have also included in the outgoing flux the part deriving from the injection $`G_{out}=G_{inj}(\mu _{out})`$ for $`\mu <u_d`$. The symbolic notation will be helpful in the following.

Before specifying the equivalent relationship for the upstream part of the fluid, we pause to establish the connection between $`P_d`$ and the solutions of Eq. 7 (or Eq. 9). Given the definition of $`P_d`$, it seems obvious enough that $`P_d`$ is the particle flux $`(u+\mu )f(\mu )`$ given by the solution of Eq. 7 (or of Eq. 9) with boundary condition $`(u_d+\mu )f=\delta (\mu \mu _{in})`$ (unit entering flux) for $`\mu >u_d`$ at the shock ($`z=0`$), and $`f`$ regular at downstream infinity. In fact, we can see that this solution represents a monochromatic (in energy) flux of particles, all moving initially in the same direction, subject to pitch angle scattering, which is in turn responsible for kicking them back toward the shock, occasionally, along a direction $`\mu _{out}`$, with nonuniform probability.

One way of looking at $`P_d`$ is by realizing that it is somewhat related to the Green’s function for a diffusive problem. As an analogy, consider the paradigmatic diffusion equation

$$\frac{n}{z}=\frac{^2n}{\mu ^2}$$
(13)
which is solved by convolving the given boundary conditions with its Green’s function

$$2(\pi z)^{1/2}\mathrm{exp}((\mu \mu _i)^2/4z)\delta (\mu \mu _i)$$
(14)
as $`z0`$. Here too, the Dirac’s delta appears in the boundary conditions, and is smoothed out as the time variable ($`z`$) evolves. The major difference is that suitable boundary conditions for the above equation include the whole range in $`\mu `$, at $`z=0`$, while in our problem we can specify boundary conditions only for $`\mu >u_d`$. The function $`P_d`$ is the projection of the solution of the problem on exactly that part of the boundary at $`z=0`$ where we cannot specify the boundary conditions, a part which does not exist for Eq. 13. Still, because of this relationship, we will refer to $`P_d`$ and to $`P_u`$, shortly to be determined, as the Green’s functions of the problem.

As an aside, we show that, under some physically realistic conditions, $`P_d`$ does not depend upon the particle impulse $`p`$. This may be at first surprising because both coefficients $`d(\mu )`$ and $`w(\mu ,\mu ^{})`$ in Eq. 7 are allowed to depend on $`p`$: the proof follows. Rewrite Eq. 7 outside the shock (which means that we can drop the injection term) without the term depending upon the magnetic field, using as a new variable $`yd(\mu )z`$; this can always be done because, as noticed above, $`d>0`$. We obtain

$$\gamma (u+\mu )\frac{f}{y}=f+g(\mu ,\mu ^{})f(\mu ^{})𝑑\mu ^{}.$$
(15)
Since we assumed that the $`p`$–dependence in $`w`$ and $`d`$ factors out, we see from Eq. 6 that $`g`$ is independent of $`p`$. Then from Eq. 15 we see that $`p`$ has disappeared from the problem altogether, and that any solution can depend on $`p`$ only through the term $`y=d(\mu )z`$. However, $`P_d`$ is the flux at $`z=0`$, the shock position, so that we see that when this is computed ($`y=0`$), all dependence on $`p`$ drops out: $`P_d`$ does not depend on $`p`$. This result was known to Bell (1978), who proved it for Newtonian shocks; the proof given here is however valid for all shock speeds.

Upstream, we shall also have a probability distribution $`P_u(\mu _{out},\mu _{in})`$ that a particle leaving downstream along a direction $`\mu _{out}`$ will recross the shock along a direction $`\mu _{in}`$. An argument entirely identical to the one given above shows that, often, $`P_u`$ does not depend upon the particle’s impulse. There is however a small complication, in that particles of impulse $`p_i`$ in the downstream frame, exiting it along a direction $`\mu _{out}`$ and re-entering it along a direction $`\mu _{in}`$, emerge with a different impulse, $`p`$ (as measured again in the downstream frame) given by

$$p=\frac{1u_r\mu _{out}}{1u_r\mu _{in}}p_iGp_i.$$
(16)
$`G`$ is the energy amplification a particle receives, as it cycles once around the shock. Since $`1\mu _{out}u_d`$, and $`u_d\mu _{in}1`$, we have

$$1G\frac{1+u_r}{1u_r}.$$
(17)
Thus, particles leaving the downstream section with impulse $`p_i`$ provide the reentering flux at a different impulse. A trivial computation yields the entering flux at impulse $`p`$ as

$`J_{in}(p,\mu _{in})={\displaystyle _1^{u_d}}𝑑\mu _{out}P_u(\mu _{out},\mu _{in})\left({\displaystyle \frac{1u_r\mu _{in}}{1u_r\mu _{out}}}\right)^3J_{out}(p_i,\mu _{out})`$
$`+G_{in}(p,\mu _{in})P_uJ_{out}^{}+G_{in}(p,\mu _{in})`$ (18)
where again we added to the incoming flux the injected one $`G_{in}=G_{inj}(\mu _{in})`$ for $`\mu _{in}>u_d`$, and the prime over $`J_{out}`$ reminds us that the outgoing flux is to be computed at particle energy $`p_i`$ given by Eq. 16.

A comment on normalization is in order. All particles that enter upstream will in due time be brought backwards across the shock by the fluid’s advection, so that we expect

$$_{u_d}^1𝑑\mu _{in}P_u(\mu _{out},\mu _{in})=1$$
(19)
while the analogous statement for the downstream Green’s function $`P_d`$ does not hold for the same reason: a fraction of all particles will be advected toward downstream infinity, and the returning flux will be smaller than the departing one. In general, thus, $`_1^{u_d}𝑑\mu _{out}P_d(\mu _{in},\mu _{out})<1`$. The average probability of coming back to cross the shock is given by averaging $`P_d`$ over the whole incoming flux:

$$P=\frac{_1^{u_d}𝑑\mu _{out}_{u_d}^1𝑑\mu _{in}P_d(\mu _{in},\mu _{out})J_{in}(\mu _{in})}{_{u_d}^1𝑑\mu _{in}J_{in}(\mu _{in})}.$$
(20)
When the contribution of injection can be neglected, we see from the above, and from Eq. 12 that

$$P=\frac{_1^{u_d}𝑑\mu _{out}J_{out}(\mu _{out})}{_{u_d}^1𝑑\mu _{in}J_{in}(\mu _{in})}.$$
(21)
The reason for this is quite simple: in a steady state, like that assumed here, there can be no accumulation of particles except at infinity; since there is a larger flux entering the downstream region than leaving it, there must be an accumulation of particles in it, which can thus be realized only at downstream infinity, i.e., by leaving the region of the shock.

Our cardinal equations, Eqs. 12,3, can now be simply combined to obtain

$$J_{in}=P_uP_dJ_{in}^{}+(P_uG_{out}^{}+G_{in})QJ_{in}^{}+X.$$
(22)
This equation can be solved iteratively, using the term in parentheses as the first guess:

$$J_{in,(0)}=X;J_{in,(1)}=QX^{}+X,J_{in,(N)}=QJ_{in,(N1)}^{}+X\mathrm{}$$
(23)
to arrive at the solution in the form

$$J_{in}=X+QX^{}+QQX^{\prime \prime }+QQQX^{\prime \prime \prime }+\mathrm{}.$$
(24)
The interpretation of this equation is simple: $`X`$ is the flux of injected particles, including both those that are directly moving in and those that were initially moving out, but have been turned back exactly once, The operator $`Q`$ transforms an ingoing flux into another ingoing flux, made of the fraction of all particles that have completed exactly one loop around the shock. Thus the terms $`QQX^{\prime \prime }`$ and so on represent the particles which have made $`2,3,\mathrm{}N,\mathrm{}`$ loops around the shock, after injection. The fact that this theory manages to recover the well–known fact that the spectrum is the superposition of infinite bumps, each corresponding to particles that have looped an integer number of times around the shock, can be considered a useful test which has been successfully cleared. It also has the advantage of yielding the distribution function at the shock for the first time both very close to the injection energy, and very far from it. Also, it substantiates physical intuition that the solution of the problem must depend only on the injected flux ($`X`$) and on the scattering properties of the medium. Lastly, the equation concerning $`J_{out}`$ carries no new information: once $`J_{in}`$ has been found, it can be plugged into Eq. 12 which yields a $`J_{out}`$ which automatically satisfies Eq. 3.

## 4 The large-$`p`$ limit

We see from Eq. 16 that after half a cycle (only the upstream section changes the particles’ energies), the energy of a particle just injected with energy $`p_{inj}`$ may span the range ($`1mu_{out}u_d\mu _{in}1`$)

$$p_{inj}<p<p_{inj}\frac{1+u_r}{1u_r}p_{max}.$$
(25)
This shows that the evolution of that part of the particle distribution function which describes particles making many loops around the shock, must extend toward larger and larger energies; in fact, since the minimum amplification is $`G=1`$, the process under discussion does not describe a spreading of the injection distribution around the initial value, but a drift towards larger and larger energies (plus a spreading, of course): Fermi acceleration of Type I. So the question naturally arises of what is the asymptotic spectrum of particles, as we consider asymptotically larger energies.

If injection provides particles of greatest energy $`p_{inj}`$, we may neglect the injection terms in the above equations for energies exceeding $`p_{max}`$, and we thus consider Eq. 22 in this case:

$$J_{in}(p,\mu )=_1^{u_d}𝑑\mu _{out}P_u(\mu _{out},\mu )\left(\frac{1u_r\mu }{1u_r\mu _{out}}\right)^3_{u_d}^1𝑑\xi P_d(\xi ,\mu _{out})J_{in}(\frac{1u_r\mu _{in}}{1u_r\mu _{out}}p,\xi ).$$
(26)
We show now that this equation, except for the null solution, has a solution which is a simple power law in $`p`$. Let us first consider two situations, where injection occurs at slightly different energies, $`p_0`$ and $`p_0+\delta p_0`$. At large energies, $`p>p_{max}`$, the spectra of these two situations will differ by very little, $`J_{in}/p_0\delta p_0`$. The difference between two solutions of the inhomogeneous problem (i.e., including injection) must be a solution of the associated homogeneous problem (the one without injection, Eq. 26), so that $`J_{in}/p_0`$ is a solution of Eq. 26. For obvious dimensional reasons, $`J_{in}=h(p/p_0,\mu )/p^3`$ so that (using $`\dot{h}(z,\mu )h/z`$)

$$\frac{J_{in}}{p}=\frac{3h}{p^4}+\frac{\dot{h}}{p^3p_0}$$
(27)
$$\frac{J_{in}}{p_0}=\frac{\dot{h}}{p^2p_0^2}$$
(28)
and, upon eliminating $`h`$,

$$p\frac{J_{in}}{p}=3J_{in}\frac{1}{p_0}\frac{J_{in}}{p_0}$$
(29)
which shows $`pJ_{in}/p`$ to be the linear combination of two solutions of the homogeneous equation 26 (for $`p>p_{max}`$), and thus a solution itself. If we now assume Eq. 26 to have a unique solution, we find that $`pJ_{in}/p`$ must be proportional to $`J_{in}`$, when injection can be neglected. We thus have

$$\frac{1}{J_{in}(p,\mu )}\frac{J_{in}}{p}(p,\mu )=\frac{s}{p}$$
(30)
where $`s`$ is a constant yet to be determined, and the obvious solution

$$J_{in}(p,\mu )=\frac{(u_d+\mu )g(\mu )}{p^s}$$
(31)
follows, where we have arbitrarily factored out $`u_d+\mu `$ to remind us that $`J`$ is a particle flux; thus, $`g(\mu )/p^s`$ is the particle distribution function. It must be remarked that this power–law dependence has been derived without assuming that $`P_d`$ and $`P_u`$ are independent of the particles’ momentum $`p`$.

The importance of the above is obvious: it shows that the particle distribution function is a power law in the particles’ momentum, for every shock speed. This result had been obtained before in all numerical solutions for arbitrary shock speeds, and in analytic form in the Newtonian limit, but in the full relativistic regime it was put in by hand (Kirk and Schneider 1987), and the boundary conditions could have equally well been satisfied by any other function with a free paramemter. The present derivation is universal, and it has the rather unexpected result (at least, for this author) that the power law solution is exact, not asymptotic, as soon as we consider energies exceeding the injection ones, i.e., for $`p>p_{max}`$, and not just $`pp_{max}`$.

With the result provided by Eq. 31, our principal equations, Eqs. 12 and 3, assume the following form, when injection may be neglected:

$$(u_d+\mu )g(\mu )=_{u_d}^1𝑑\xi P_d(\xi ,\mu )(u_d+\xi )g(\xi )$$
(32)
$$(u_d+\mu )g(\mu )=_1^{u_d}𝑑\mu _{out}P_u(\mu _{out},\mu )\left(\frac{1u_r\mu }{1u_r\mu _{out}}\right)^{3s}(u_d+\mu _{out})g(\mu _{out}).$$
(33)
The first one concerns the outgoing flux (with respect to the downstream section), and thus is valid for $`1\mu u_d`$, while the second one concerns the entering flux, and applies for $`u_d\mu 1`$. Also, use of Eq. 31 in Eq. 21 shows that the probability of reaching downstream infinity does not depend upon the particle’s momentum, $`p`$.

Again, they can be combined to obtain

$$(u_d+\mu )g(\mu )=_{u_d}^1𝑑\xi Q^T(\xi ,\mu )(u_d+\xi )g(\xi )$$
(34)
where

$$Q^T(\xi ,\mu )_1^{u_d}𝑑\nu P_u(\nu ,\mu )P_d(\xi ,\nu )\left(\frac{1u_r\mu }{1u_r\nu }\right)^{3s}.$$
(35)
which applies to the incoming flux only, $`u_d\mu 1`$. Once again, the outgoing flux carries no new information: if the above is solved, plugging the solution into Eq. 32 returns an outgoing flux which automatically solves Eq. 33. All of the problem’s information is in one of the subintervals only, even in the case of the homogeneous equation.

The above equation, for an arbitrary value of $`s`$, does not have a solution. The equation which does have a solution is

$$\lambda (u_d+\mu )g(\mu )=_{u_d}^1𝑑\xi Q^T(\xi ,\mu )(u_d+\xi )g(\xi ).$$
(36)
In fact, the above equation is an homogeneous equation of the Fredholm type, with smooth and bounded kernel over a compact interval (Courant and Hilbert 1953), which has at most a finite number of solutions, each belonging to an eigenvalue $`\lambda `$ which is to be determined simultaneously with the eigenvector $`g(\mu )`$. We thus see that what fixes $`s`$ is the requirement that $`\lambda =1`$: the physical value of $`s`$ is that which has a unitary eigenvalue. We are forced by our limited mathematical skills to assume that Eq. 34 has exactly one solution, neither more nor fewer, but, apart from this, the task of determining the particle distribution function is completed: Eq. 36 with $`\lambda =1`$ yields the particle distribution function at the shock, and the all–important particle spectral index $`s`$, even in the case of arbitrary particle anisotropy.

### 4.1 The Newtonian limit

In order to gain physical insight into the condition that fixes $`s`$ ($`\lambda =1`$), we rewrite Eq. 36 somewhat. First, we write Eq. 33 as

$$\lambda (u_d+\mu )g(\mu )=_1^{u_d}𝑑\mu _{out}P_u(\mu _{out},\mu )\left(\frac{1u_r\mu }{1u_r\mu _{out}}\right)^{3s}(u_d+\mu _{out})g(\mu _{out}).$$
(37)
This can now be integrated over the whole range in $`\mu `$, $`u_d\mu 1`$, and divided by the whole outgoing flux, to obtain

$$\lambda \frac{_{u_d}^1𝑑\mu (u_d+\mu )g(\mu )}{_1^{u_d}𝑑\mu _{out}(u_d+\mu _{out})g(\mu _{out})}=\frac{_{u_d}^1𝑑\mu _1^{u_d}𝑑\mu _{out}P_u(\mu _{out},\mu )\left(\frac{1u_r\mu }{1u_r\mu _{out}}\right)^{3s}(u_d+\mu _{out})g(\mu _{out})}{_1^{u_d}𝑑\mu _{out}(u_d+\mu _{out})g(\mu _{out})}$$
(38)
The term on the left hand side can be rewritten by means of Eq. 21 as $`\lambda /P`$. The term on the right hand side is clearly the average value of the $`(s3)`$–th power of the amplification $`G`$ (Eq. 16), over the whole particle distribution: $`<G^{s3}>`$. The above is thus

$$\lambda =P<G^{s3}>,$$
(39)
and demanding that $`\lambda =1`$ implies that

$$P<G^{s3}>=1.$$
(40)
In the Newtonian limit, in which $`G11`$, we obviously have $`<G^{s3}><G>^{s3}`$, and thus the energy spectral index $`k=s2`$ is given by

$$s2k=1\frac{\mathrm{log}P}{\mathrm{log}<G>}$$
(41)
which is exactly Bell’s equation for the spectral index. We have thus established that the condition $`\lambda =1`$ is equivalent to $`P<G^{s3}>=1`$, which appears to be the correct relativistic generalization of Bell’s Newtonian results.

In the Newtonian limit, where $`g(\mu )`$ constant, independent of $`\mu `$ (the isotropic limit), it is easy to derive from Eq. 21 that $`P1u`$, and from Eq. 38 that $`<G>1+u`$, from which $`s=4`$, as per Bell’s results for a strong shock.

There are some differences between the fully relativistic approach and the Newtonian one. First, in the relativistic approach in general $`<G^{s3}>`$ cannot be approximated by $`<G>^{s3}`$, contrary to Peacock’s (1981) arguments. Second, here $`G`$ and $`P`$ are functions of $`s`$, so that Eq. 40 sets up a transcendental equation for $`s`$. Third, Eq. 40 is only of symbolic value, since it requires knowledge of the anisotropic particle distribution function, which is not known a priori. But luckily, the full relativistic approach developed here returns simultaneously the anisotropic distribution function (Eq. 34).

### 4.2 Bits and pieces

A useful identity is obtained by rewriting Eq. 3 without the injection term as

$$\frac{J_{in}(p,\mu _{in})}{(1u_r\mu _{in})^3}=_1^{u_d}𝑑\mu _{out}P_u(\mu _{out},\mu _{in})\frac{J_{out}(p_i,\mu _{out})}{(1u_r\mu _{out})^3},$$
(42)
integrating over $`\mu _{in}`$, using Eq. 19, and then inserting Eq. 31, to obtain

$$_{u_d}^1𝑑\mu _{in}(1u_r\mu _{in})^{s3}(u_d+\mu _{in})g(\mu _{in})=_1^{u_d}𝑑\mu _{out}(1u_r\mu _{out})^{s3}(u_d+\mu _{out})g(\mu _{out})$$
(43)
which of course is not trivial, because the intervals of integration are different.

It is worth noticing that, in all of Section 3, no reference whatsoever has been made to the explicit form of the right hand side of Eq. 7; all that is necessary is to determine the equation’s Green’s functions $`P_d`$ and $`P_u`$, but otherwise all results apply to all right hand side members, and can then be seen as a property of the left hand side member $`(u+\mu )f/z`$. They are thus of very general validity.

Lastly, we have assumed that Eqs. 34 and 36 has one, and just one solution. If it had none, then the steady state problem would have no solution, and the problem would be intrinsically time–dependent. If it had several, then a rather interesting case would arise, whereby shocks moving at the same speed in the same media, but with different injection properties would have different asymptotic spectral indices $`s`$. The near universality of the index $`s`$ points of course in the opposite direction, but it is also true that a proof of the existence of a unique solution is beyond our modest abilities.

## 5 $`P_u`$ and $`P_d`$

In this subsection we show how to build $`P_u`$ and $`P_d`$. For sake of definiteness, we concentrate on $`P_u`$, but a wholly analogous treatment holds for $`P_d`$.

We now wish to find an explicit expression for $`P_u(\mu _{out},\mu _{in})`$. That this quantity exists, it is physically obvious: it represents the probability that a net flux of particles leaving the downstream section along a direction $`\mu _{out}`$ (in downstream variables), reenters along a direction $`\mu _{in}`$. We can imagine an experiment to determine this quantity: an experimenter located next to the shock with a small cannon shooting out particles along the direction $`\mu _{out}`$, may collect them on a screen, determining how many come back along $`\mu _{in}`$. It corresponds to a solution of Eq. 7 or 9 without the injection term, with boundary condition

$$dJ=\delta (\mu _d\mu _{out})d\mu _d=\delta (\mu _u\mu _{out}^{(u)})d\mu _u.$$
(44)
Here, $`\mu _{out}^{(u)}`$ is the value of $`\mu _{out}`$ in the upstream frame variables, and the last identity is a trivial property of Dirac’s deltas.

In order to find $`P_u`$, we now use separation of variables to solve Eq. 7 or Eq. 9. For instance, writing $`f(\mu _u,z)=A(z)B(\mu _u)`$, we find the solutions of

$$\gamma \frac{1}{A(z)}\frac{dA}{dz}=\lambda _n=\frac{1}{(u+\mu _u)B(\mu _u)}\frac{d}{d\mu _u}\left(D(\mu _u,p_u)(1\mu _u^2)\frac{dB}{d\mu _u}\right).$$
(45)
Calling $`B_n`$ the angular eigenvectors, we find a generic solution of the above as

$$f=\underset{n}{}a_n\mathrm{exp}(\lambda _nz/\gamma )B_n(\mu _u)$$
(46)
where the $`a_n`$’s are coefficients which we must choose so as to satisfy the initial conditions, Eq. 44. We must thus have, at $`z=0`$,

$$\gamma (u+\mu _u)f=\underset{n}{}a_n(\mu _{out}^{(u)})(u+\mu _u)\gamma B_n(\mu _u)=\delta (\mu _u\mu _{out}^{(u)}).$$
(47)
Obviously, the coefficients $`a_n`$’s will depend upon the point where the delta is located, $`\mu _{out}^{(u)}`$.

Here there is a tricky, but important point. The above eigenvector problem is well–known to have solutions belonging to both positive and negative eigenvalues (plus the $`\lambda =0`$ case). Eigenvectors belonging to negative values of $`\lambda _n`$, of which there is an infinite number, are physically ill–behaved at upstream infinity, as can be seen from Eq. 46, so that we surely have to restrict ourselves, in all sums above to the well–behaved eigenvectors, i.e., those with $`\lambda _n>0`$. But does the above equation then hold? The answer is yes, as we now show.

The fact that an infinite set of functions can be arranged to satisfy the above equation is called by mathematicians completeness, indicating that any function can be written as the superposition of $`B_n`$’s with suitable coefficients. If the above equation holds, in fact, multiplying both sides by any function $`F(\mu _{out}^{(u)})`$ and integrating over the whole range of $`\mu _{out}^{(u)}`$, we find that

$$F(\mu _u)=\underset{n}{}\left(𝑑\mu _{out}^{(u)}a_n(\mu _{out}^{(u)})F(\mu _{out}^{(u)})\right)\gamma (u+\mu _u)B_n(\mu _u),$$
(48)
proving that we can write any function as the superposition of functions from the set, with suitable coefficients, and thus that the set is, by definition, complete. Now, the set of all eigenfunctions $`B_n`$ is well–known to be complete, but here is the rub: the full set includes both the eigenvectors which have $`\lambda >0`$ (which thus are physically well behaved at upstream infinity) and those which have $`\lambda 0`$, which diverge at upstream infinity, and which we thus cannot use.

The half set of the well–behaved eigenvectors is not complete over the whole range in $`\mu `$, but here we are helped by an unusual property of the eigenvectors of this equation: it can be shown in fact (Freiling, Yurko and Vietri 2002) that the half–set of all well–behaved eigenvectors (i.e., those with $`\lambda >0`$) is complete in the restricted range $`1\mu _uu`$. This property (sometimes called half-range half-completeness) implies that, for $`\mu _uu`$, we can always write

$$\gamma (u+\mu _u)f=\underset{n}{\overset{}{}}a_n(u+\mu _u)\gamma B_n(\mu _u)=\delta (\mu _u\mu _{out}^{(u)})$$
(49)
where however now the prime over the sum reminds us that the summation is extended only to well–behaved eigenvectors. For $`\mu _u>u`$, the above equation is not verified, and the sum will give some function, which we can now identify with our required $`P_u`$. Thus, $`P_u`$ is the continuation to the rest of the interval ($`\mu _u>u`$) of the sum above, which inside the restricted range $`\mu u`$ satisifes the above equality.

One word of caution should be stated here: the property of half-range half-completeness has been proved for the eigenvectors of Eq. 9, but not for those of Eq. 7. Still, since they are known to be so similar (both are called of the third, or polar type, one within the realm of integral equations, the other one of differential equations, to emphasize several formal similarities), and since the above physical argument makes it extremely plausible that $`P_u`$ exists also for this case, we shall make no further distinction between the two cases, and proceed as if this argument went through for both Eq. 9 and Eq. 7.

We now find the coefficients $`a_n`$’s. The well–behaved eigenvectors are not orthormal over the restricted range $`\mu _uu`$. They can easily made so by means of a well–known procedure (Schmidt’s diagonalization, Courant and Hilbert 1953), which leaves the lowest order eigenvector unaffected (except for the normalization), and suitably modifies the others in a finite number of steps. Once this new basis $`P_n`$ has been found, since it is complete (from the half-range half completeness), we easily find, for $`\mu _uu`$,

$$(u+\mu _u)\underset{n}{}P_n(\mu _{out}^{(u)})P_n(\mu _u)=\delta (\mu _u\mu _{out}^{(u)})$$
(50)
and thus, by working backwards to Eq. 10, (for $`\mu _{out}<u_d`$, $`\mu _{in}>u_d`$)

$$P_u(\mu _{out},\mu _{in})d\mu _{in}=(u+\mu _u)\underset{n}{}P_n(\mu _{out}^{(u)})P_n(\mu _u)d\mu _u.$$
(51)
## 6 Summary

The major results presented in this paper are the following:

* We have presented a relativistically covariant equation (Eq. 7) for the particle distribution function, which includes pitch–angle scattering (without assuming that the diffusive approximation holds), the effect of smooth magnetic fields, and particle injection at the shock; this equation applies separately in the upstream and downstream frames.
* We have described the transport properties of the upstream and downstream media by means of two Green’s functions, $`P_d`$ and $`P_u`$, which are independent of the particle distribution functions, and which have been explicitly constructed in terms of the eigenfunctions of the angular part of the problem (Eq. 51).
* We have established two relations (Eq. 12 and 3) which provide the true boundary conditions of the problem, giving the flux entering the downstream section in terms of the one leaving the upstream section, and vice cersa, plus the injection terms. This naturally sets up an equation for the particle flux, Eq. 22.
* We have solved this equation, Eq. 24, in a general way, which yields the particle spectrum both close to and far from the injection energies.
* We have shown that under physically realistic conditions, the probabilities of returning to the shock, $`P_u`$ and $`P_d`$, do not depend upon the particles’ energy $`p`$.
* In the limit of large particle energies, we have shown that the spectrum is a pure power law in the particles’ momentum, even when the probabilities of returning at the shock, $`P_u`$ and $`P_d`$, depend on $`p`$ (Eq. 31).
* In this limit, we have then simplified our main equations, arriving at a system of equations (Eq. 32, 33, or their combination, Eq. 34), and further showing that the requirement that Eq. 34 has the eigenvalue $`\lambda =1`$ fixes the all–important asymptotic particle spectral index $`s`$.
* We have recovered the low speed, Newtonian limit of Bell (1978) and Blandford and Ostriker (1978), and provided a relativistic generalization (Eq. 40) which is however of little use because it requires knowledge of the particle anisotropic distribution function, which can be obtained only through the full solution of Eq. 34, which automatically incorporates Eq. 40.

Though the treatment presented in this apper may appear rather abstract, it is not without practical consequences. For instance, the fact that the asymptotic spectral index $`s`$ is not asymptotic, but it applies as soon as $`p>p_{max}=p_{inj}(1+u_r)/(1u_r)`$, where $`u_r`$ is the modulus of the relative speed between the upstream and downstream fluids, implies that numerical simulations trying to compute $`s`$ need not extend to very large particle energies, but can save precious computational time by sticking to energies exceeding $`p_{max}`$ by a small factor, dictated exclusively by numerical questions. Also, the treatment presented in Section 3 may be somewhat simplified by a judicious use of Fredholm’s formulae, resulting in a considerable simplification for the calculational side of the theory. And lastly, they allow the treatment of the hyperrelativistic limit, as will be shown in a future paper.


text:
# References

Efficient Hopfield pattern recognition on a scale-free neural network

Dietrich Stauffer1,2,4, Amnon Aharony1, Luciano da Fontoura Costa 3 and Joan Adler 4

1 School of Physics and Astronomy, Raymond and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel

2 Institute for Theoretical Physics, Cologne University, D-50923 Köln, Euroland

3 Cybernetic Vision Research Group, IFSC-USP, Caixa Postal 369, 13560-970 São Carlos, SP, Brazil

4 Department of Physics, Technion-IIT, Haifa 32000, Israel

e-mail: aharony@post.tau.ac.il, stauffer@thp.uni-koeln.de,

luciano@if.sc.usp.br, phr76ja@techunix.technion.ac.il

Neural networks are supposed to recognise blurred images (or patterns) of $`N`$ pixels (bits) each. Application of the network to an initial blurred version of one of $`P`$ pre-assigned patterns should converge to the correct pattern. In the “standard” Hopfield model, the $`N`$ “neurons” are connected to each other via $`N^2`$ bonds which contain the information on the stored patterns. Thus computer time and memory in general grow with $`N^2`$. The Hebb rule assigns synaptic coupling strengths proportional to the overlap of the stored patterns at the two coupled neurons. Here we simulate the Hopfield model on the Barabási-Albert scale-free network, in which each newly added neuron is connected to only $`m`$ other neurons, and at the end the number of neurons with $`q`$ neighbours decays as $`1/q^3`$. Although the quality of retrieval decreases for small $`m`$, we find good associative memory for $`1mN`$. Hence, these networks gain a factor $`N/m1`$ in the computer memory and time.

PACS: 05.40-a, 05.50+q, 87.18.Sn

Traditional neural network models have nodes $`i`$ (“neurons”) coupled to all other nodes $`k`$ with some coupling constant $`J_{ik}`$ (”synaptic strength”), similar to Sherrington-Kirkpatrick infinite-range spin glasses . Here we consider one of the simplest neural network models, due to Hopfield . This model was mostly applied to infinite range and was only rarely put onto a square lattice with short-range interactions . Real neural networks seem to have neither infinite nor only nearest-neighbour connections. The spatial structures of neural networks were investigated  and compared with small-world and scale-free networks . Now we present computer simulations of the Hopfield model  with Hebb couplings between neighbours restricted to a Barabási-Albert (BA) scale-free network .

In the Hopfield model, each of $`N`$ neurons or sites can be firing ($`S_i=+1`$) or not firing ($`S_i=1`$). Neurons are coupled through $`J_{ik}`$, and are sequentially updated according to

$$S_i\mathrm{sign}(\underset{k}{}J_{ik}S_k).$$
$`(1)`$
(We mostly ignore the diagonal terms $`i=k`$ in our sums.) This rule corresponds to a low-temperature Monte Carlo simulation of a spin glass. The model has stored $`P`$ different patterns $`\xi _i^\mu (\mu =1,2\mathrm{},P),`$ which we take as random strings of $`\pm 1`$. The couplings are given by the Hebb rule:

$$J_{ik}=\underset{\mu }{}\xi _i^\mu \xi _k^\mu .$$
$`(2)`$
The first of these patterns is presented to the network in a corrupted form $`S_i`$, with ten percent of the $`S_i`$ reversed in comparison to the correct $`\xi _i^1`$. The question is whether the iteration through Eq. (1) transforms the erroneous $`S_i`$ into the correct $`\xi _i^1`$. The quality of this pattern recognition is given by the overlap

$$\mathrm{\Psi }=\underset{i}{}S_i\xi _i^1/N,$$
$`(3)`$
which is related to the Hamming distance and equals 1 for complete recognition and $`\pm 1/\sqrt{N}`$ for only accidental agreement at random sites; it is $`0.8`$ at the beginning of the pattern recognition process, due to the ten percent reversal.

Now we restrict the synaptic connections $`J_{ik}`$ to neurons which are neighbours in the BA network, but we still use Eqs. (1-3). In these networks, we start from a small core of $`m`$ sites, all connected with each other. Then $`Nm`$ nodes are added, one after the other. Each new site $`i`$ selects exactly $`m`$ sites from the already existing network sites as its neighbours $`k`$, with a probability proportional to the number of neighbours which the existing site $`k`$ has already: The rich get richer. When the network has added $`N`$ sites with a total of $`N+m`$ sites, its growth is stopped and the neural process of Eqs. (1-3) starts. Synaptic connections $`J_{ik}`$ exist only between sites $`i`$ and $`k`$ which are neighbours.

Since no longer every neuron is connected to all other neurons, the memory-saving trick of Penna and Oliveira  to avoid storing the $`J_{ik}`$ no longer applies. 400 Megabyte were needed for $`N=10,000`$ nodes and $`P=20,000`$ patterns. To save computer time, the $`J_{ik}`$ should be determined after and not before the construction of the BA network.

When only one pattern is stored, it is recognised completely after two iterations. With $`P>1`$, however, no complete recognition takes place, the overlap $`\mathrm{\Psi }`$ is usually at the final fixed point (reached after about five iterations) lower than at the beginning, as shown in Fig. 1a . However, the model can still recognise the first pattern as the one presented to it, since the overlap $`\mathrm{\Psi }0.19`$ for $`P=N=10^4`$ is still appreciably larger than the overlap $`|\mathrm{\Psi }|<500`$ with the other $`(P1)`$ patterns.

Rather similar results are obtained if we work on a nearest-neighbour hypercubic lattice with $`N=L^d`$ sites, similar to the studies made in  in two dimensions. Fig.1b shows that only for small numbers $`P`$ of patterns an increased $`d`$ means an increased final overlap. For $`d=7,\mathrm{\hspace{0.33em}10}`$ and 15 the overlaps with up to 20 patterns did not differ appreciably from $`d=5`$. No significant size effects were seen for $`4L20`$ in five and $`4L13`$ in seven dimensions.

Fig. 1 is based on one sample only for each point; using instead 100 samples at $`m=3`$ and $`N=10,000`$, we see in Fig. 2 that the overlap varies roughly as $`\mathrm{\Psi }(P)0.19P^{0.6}`$, except for very small $`m`$. A similar power law $`P^{0.6}`$ is also found for hypercubic lattices (not shown.) It would be interesting to understand this power law from some analytical analysis.

A much better recovery of the corrupted pattern is obtained if we take a larger inner core of the BA network, that means if $`m`$ is no longer small. (The first $`m`$ network sites are all mutually connected, as in the traditional Hopfield model.) Using 100 patterns, Fig. 3a shows the overlap $`\mathrm{\Psi }`$ for $`N+m=10^4`$ total sites as a function of $`m`$. Already at $`m=200,N=9800`$ the final overlap is 88 percent; at $`m=2000,N=8000`$ we have complete recovery. For $`1mN`$ the number of connections (counting each bond twice) is $`mN`$ in our case and $`N^2`$ in the fully connected case; thus we saved connections by a factor $`m/N`$. If we include the diagonal terms $`J(i,i)=P`$ in Eq. (1), we prevent the overlap from becoming worse than the initial overlap 0.8 for small $`m`$ and still get overlap near unity for large $`m`$, Fig. 3b.

Of course, with a large $`m`$ the network is no longer scale free, as shown in Fig. 4: The simple power law $`1/q^3`$ for the number of sites with $`q`$ neighbours each  persists for $`10^2<q<10^3`$, but a Gaussian peak is added for large $`q`$. However, the additional bump concerns only a relatively small number of sites, and is probably negligible for any practical purposes.

For infinite range, $`m=N`$, the usual Hopfield model  gives an overlap $`\mathrm{\Psi }`$ close to $`1`$ for $`P/N<0.14`$ and a relatively small overlap $`\mathrm{\Psi }0.2`$ for $`P/N>0.14`$, with a sharp jump at $`P/N=0.14`$. Our simulations, in contrast, show a gradual deterioration as soon as more than one pattern is stored, but the value of $`\mathrm{\Psi }`$ is still of order 0.2 and distinctly larger than for the other $`(P1)`$ patterns. Using a medium-sized fully connected core, like $`m10^3`$ at $`P10^2`$, surrounded by a larger BA network with $`N10^4`$ sites, gives a good compromise between good recovery and not too many connections.

So far it is not clear if the efficient recovery simply results from the relatively large average coordination number $`m`$, or by some additional ingredients in the problem. It would also be interesting if Nature takes advantage of a similar efficiency. If it does, do natural neural networks share some geometrical features with the large-$`m`$ (but finite) scale-free networks ?

We thank D. Horn for a discussion and the German-Israeli Foundation for supporting our collaboration.