We have applied the maximum entropy method to a large data set based on Compton Gamma Ray Observatory (CGRO)OSSE, Solar Maximum Mission, Transient Gamma-Ray Spectrometer, Gamma-Ray Imaging Spectrometer, HEXAGONE, and FIGARO observational results of the Galactic center 511 keV radiation in order to produce a map of the 511 keV line emission. This map suggests two components: a central bulge with a flux of 5.7 ± 0.29 × 10-4 photons cm-2 s-1 and a Galactic plane (GP) component with a flux of 2.2 ± 0.4 × 10-4 photons cm-2 s-1. The central bulge is located at l = -0
47 ± 0
24, b = 0
1 ± 0
18, and FWHM = 5
28 ± 0
48. We note that the position of the GP component coincides with a strong hot spot in the COMPTEL map of 1.8 MeV 26Al line emission. A comparison between CGRO
OSSE and other instruments with larger field of view suggests an extended diffuse 511 keV emission. An interesting hot spot at l = -4°, b = 7° with a flux of
2 × 10-4 photons cm-2 s-1 is shown on our map. A bootstrap test indicates that the significance level of this feature is 3.5
.
Subject headings: Galaxy: generalgamma rays: observations
1 Department of Astronomy, University of Maryland, College Park, MD 20742.
2 Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-2900.
3 Laboratory for High Energy Astrophysics, NASAGSFC, Greenbelt, MD 20771.
4 Space Science, Morse Hall, University of New Hampshire, Durham, NH 03824.
5 University of California, Institute of Geophysics and Planetary Physics, Riverside, CA 92521.
6 Code 7650, Naval Research Laboratory, 4555 Overlook Avenue, SW, Washington, DC 20375-5352.
The 511 keV emission from the Galactic center (GC) region was detected initially in 1970 (Johnson & Haymes 1973) and was identified as annihilation radiation by the first high-resolution observation in 1977 (Leventhal, MacCallum, & Stang 1978). Many balloon and satellite experiments have since reported this radiation from the general direction of the GC (see Tueller 1993, Skibo, Ramaty, & Leventhal 1992, and Lingenfelter & Ramaty 1989 for reviews). Detectors with larger fields of view (FOVs) tend to detect higher fluxes, and detectors with smaller FOVs have reported variable signals, suggesting that the GC region 511 keV source consists of at least two components: a steady source extended over several tens of degrees and a variable point source located near the GC. The distributed 511 keV emission could be produced by the +-decay products from radioactive nuclei (e.g., 56Co, 44Sc, and 26Al) produced by supernovae, novae, or Wolf-Rayet stars (Woosley & Pinto 1988; Lingenfelter & Ramaty 1989; von Ballmoos 1991). However, it is not clear where the variable component comes from. It should be pointed out, however, that the most recent data do not indicate source variability (Purcell et al. 1993; Teegarden et al. 1996).
Knowledge of the spatial distribution of the GC 511 keV line emission is needed to address its origin. The Oriented Scintillation Spectrometer Experiment (OSSE) on the Compton Gamma Ray Observatory (CGRO) has completed numerous observations of the GC and Galactic plane (GP) regions. Purcell et al. (1994) fitted these data with various distribution models. The only model that was not rejected is a two-component distribution consisting of spheroidal and disk components. All other models could be rejected at greater than 4.5 confidence level.
We have constructed a large data set from the results of CGROOSSE and other experiments. The angular responses and measured 511 keV line flux for each of these observations are included in the data set. In this Letter, we present the GC region 511 keV map reconstructed by the maximum entropy method (MEM) as an early result of the 511 keV mapping campaign.
At the present time, focusing instruments at 511 keV have not been flown. These instruments have an intrinsically small FOV and are ill-suited to detect diffuse emission. For now, detecting the 511 keV emission from the GC region requires coded-aperture or collimator systems with a FOV of many square degrees, and the map must be deconvolved from differences between detector signals. This kind of transform imaging has been reviewed by Caroli et al. (1987) and Lei, Fraser-Mitchell, Yearworth (1991).
The image reconstruction problem normally can be summarized as follows:
where Fk is the kth data value, Mj is the number of image pixels, fj is the jth image value to be reconstructed, Rkj is the response of sample Fk to pixel j, and nk is the noise on sample k. Due to the noise in the observed data and ill-conditioning of the matrix R, the exact solution for the above equation is always in oscillation, and there exist many images that can fit the observed data. There are many approaches to choosing the most reasonable image, and MEM is one of the most successful (Gull & Daniell 1978).
We chose MEMSYS5, the most recent commercial software available for MEM analysis, as our image reconstruction tool. The algorithm of MEMSYS5 differs greatly from previous MEM algorithms in that it is based on a fully Bayesian analysis of the image reconstruction problem. Two basic concepts, the prior probability P(f) and the likelihood, are defined in MEMSYS5. P(f) gives the initial probability that image f represents the true
sky before any data relating to that part of the sky are taken into account. MEMSYS5 uses an entropy function for this prior probability. The likelihood is the probability that the observed data D could have been generated from a given image f. The response model can be used to generate the data that would have been observed (in the absence of noise) if the true sky looked exactly like the given image. The noise model can then be used to find the probability that all the deviations of this simulated data from the real data could be explained as noise. As an example, if F is the simulated data from some image f, and D is the corresponding real data, for Gaussian noise of standard deviation
, the likelihood P(D|f) can be calculated as
The probability of the image f given the data D, written as P(f|D), can be calculated via Bayes's theorem as
where P(f) is the prior probability of image f, P(D|f) is the likelihood of the data given f, and P(D) is a constant that normalizes P(f|D) to a sum of unity.
MEMSYS5 uses an iterative algorithm. It starts iterations by beginning with a map of the highest prior probability, which is usually a flat surface. After each successive iteration, the prior probability goes down a little while the likelihood goes up. MEMSYS5 continues to introduce more structure into the iterative process until P(f|D), which is proportional to the product of the prior probability and the likelihood, reaches a peak.
The data set we use is a part of the one described by Purcell et al. (1997). It consists of 323 511 keV observations. These observations include all available OSSE observations of the GP and the GC region, and many other previous observations by balloon and satellite.
The OSSE observations consist of over 400 days of data collected from 55 separate viewing periods in phase 14 of the CGRO mission. Each viewing period contains from one to several source pointings of the 3
8 × 11
4 FWHM OSSE FOV. An angular response array was calculated for each pointing. A single set of background pointings was taken for each viewing period by offset-pointing the detectors by 10°
20° away from the fields along the narrow direction of the collimator. The data were analyzed using the standard OSSE analysis technique; however, since the same background pointings were used for all of the source pointings in one viewing period, the 511 keV fluxes measured are correlated. The correlations between the source positions were calculated, and Singular Value Decomposition (Press et al. 1988) was used to decorrelate the OSSE data by modifying the 511 keV flux and angular response arrays (Purcell et al. 1997).
We have also included data from the following instruments into our data set: FIGARO (Niel et al. 1990), Gamma-Ray Imaging Spectrometer (Gehrels et al. 1991; Leventhal et al. 1993), HEXAGONE (Smith et al. 1993), Solar Maximum Mission (Share et al. 1990), and Transient Gamma-Ray Spectrometer (Teegarden et al. 1996). The image area was divided into pixels; the intensities in the pixels are the image values f required by MEMSYS. The instrument angular response R (i.e., the contributions from the pixels to the observed flux) for each observation was estimated from the available information. The 323 observed fluxes are the data (F), and the errors are the quoted statistical uncertainties.
As shown in Purcell et al. (1993), the measured 511 keV fluxes by the older, wide FOV instruments do not agree with each other, and this inconsistency could come from real source variability or unknown systematic errors. We removed these observations (by RICE, HEAO 3, BellSandia, GSFC, CESR, and UNH) from our data set for two reasons: (1) These observations will contribute disproportionately to the total
2 of the map. Since MEMSYS5 will keep on iterating until the reduced
2 is acceptable, the resultant map is likely to include many false features because of this overiteration. (2) We expect that 511 fluxes from wide FOV observations will provide a normalization of the total 511 keV flux in the GC region, while OSSE results will provide all the features in the map. The remaining data sets should be enough to constrain the total flux in the map.
The MEMSYS5 reconstruction iteration process will continue until the stopping condition set by the user has been satisfied. MEMSYS5 provides various possible stopping conditions. We chose our stopping condition by Monte Carlo simulation. We built up a 511 line flux model distribution based on the best-fit model from Purcell et al. (1993), and we convolved the model with the real response in order to generate a theoretical data set of 323 observed 511 line fluxes. Then the real errors were introduced to blur the theoretical fluxes, and a simulated data set was formed. We built many such data sets to explore the best stopping condition for our reconstruction. We found that the classical reduced = 1 stopping criteria is good enough for our case, i.e., all the expected features could be seen in the reconstructed maps, and no significant false features occur. The nonnegative restriction has been applied by MEMSYS5 to prevent any of the pixels in the map from being given a negative flux. Figure 1a is one of the reconstructed maps from a simulated data set. The functional form of the model is
where DC is the distance between the GC and the pixel, and DP is the Galactic latitude of the pixel. The values of c and
p were properly set to make the FWHM of the central bulge to be 5° and the FWHM of the disk to be 1
5. The longitude extension of the disk was restricted to -20°
20°.
We applied the same stopping criteria to the real data, and the resultant map is shown in Figure 1b. We used a bootstrap method to estimate the total flux and significance levels of the individual features. We generated 20 simulated data sets of 511 line flux through randomizing the real flux by incorporating the corresponding errors, and we ran MEMSYS for each of the data sets in order to generate a map. The flux of the three features were calculated for each map, and the average values and the standard deviations of the flux were used to estimate the flux and significance of the features. Figure 1c shows the distribution of the errors for the intensities in the pixels, plotted in the same scale as Figures 1a and 1b.
A total flux of 2.6 × 10-3 photons cm-2 s-1 is detected if we assume that all the 511 keV emission comes from a region of 80° × 40° (l × b) centered at the GC. Three significant features, i.e., the central bulge, the extension toward positive longitude along the GP, and the hot spot near l = -4°, b = 7°, could be seen in the map. The central bulge has a flux of 5.7 ± 0.29 × 10-4 photons cm-2 s-1. The extended feature along positive longitude has a flux of 2.2 ± 0.4 × 10-4 photons cm-2 s-1. The feature at l = -4°, b = 7° has a flux of 2.0 ± 0.58 × 10-4 photons cm-2 s-1. The centroid of the central bulge, which was determined by the same bootstrap method, is at l = -0.47 ± 0.24 and b = 0.11 ± 0.18 (where the errors are purely statistical).
Our map is consistent with Purcell et al. (1994). We have subtracted the three features from the total 511 keV line flux in order to estimate the more extended diffuse emission. However, this flux depends strongly on the distribution model. Our data do not contain sufficient information to map the extended diffuse distribution. If we assume that all the 511 keV flux comes from our mapping area, the residual flux is 1.6 ± 0.1 × 10-3 photons cm-2 s-1.
Figure 2 shows the total exposure for the OSSE observations included in the mapping analysis (Purcell et al. 1997). Obviously, the OSSE exposure is nonuniform, and this may lead to artificial features in the map. We understand that MEMSYS5 is a biased technique that tends to suppress the flux from regions with lower exposure, thus creating apparent features in regions with higher exposure. We suggest that the hot spot at l = -4°, b = 7° is not caused by nonuniform OSSE exposure for three reasons: (1) this feature is not located at a region of high exposure; (2) there is no similar feature in the reconstructed map from the simulated data (Fig. 1b); and (3) this feature is clearly present in each of the 20 maps reconstructed from the bootstrap data sets, which means it was not caused by the fluctuation in any single observation. However, the large error of the flux of this feature indicates a significance of detection of only 3.5 . It is not clear if this feature corresponds to any known X-ray or
-ray source.
The centroid of the central bulge lies between the GC and 1E 1740.7-2942 and deviates from the GC in longitude at the 2 level, which might suggest that part of the central bulge flux comes from 1E 1740.7-2942. The imaging
-ray spectrometer SIGMA on the Granat spacecraft detected a daylong burst of broad (
200 keV FWHM) annihilation line radiation from 1E 1740.7-2942 located in the GP but removed by 0
8 from the GC (Bouchet et al. 1991; Sunyaev et al. 1991). Two years later, the BATSE data showed no similar events (Smith et al. 1993). The identification of a radio counterpart exhibiting jetlike morphology (Mirabel et al. 1992) and the fact that 1E 1740.7-2942 lies along the line of sight of the dense molecular cloud G0.86-0.08 suggest an association between 1E 1740.7-2942 and the variable compact 511 keV source (Bally & Leventhal 1991; Mirabel et al. 1991; Ramaty et al. 1992).
The extension along the positive longitude direction could be the disk emission. We do not see significant negative longitude disk emission in the map, but the OSSE exposure is larger at positive longitudes than at negative longitudes (see Fig. 2). The COMPTEL map of 26Al line emission (Diehl et al. 1994) shows its strongest hot spot on the GP between l = 0° and 10°, which agrees with the GP extension on our map. We also note that this extension is near GRS 1758-258 (l 4
5, b
-1
4), which is another black hole candidate with a core and jet structure in the radio.
However, the possibility that these features are caused by nonuniform exposure effects could not be ruled out totally. A number of OSSE observations have been scheduled for the upcoming CGRO cycle 6. These observations will make the total OSSE exposure to the 25° × 25° region around the GC much more uniform. We expect improved mapping results after these observations.
There are many other deconvolution methods that have been applied successfully to astronomical problems. The Richardson-Lucy method is one of them. This method is often used in the image analysis of the Hubble Space Telescope data. However, since our data have much worse statistics than the optical data, and the number of observations is much less than the number of pixels, this method may not be the best approach for our analysis.
Our data set is susceptible to a wavelet analysis, which is very useful for separating point sources and small-scale structure from diffuse structure. MEMSYS works well when the solution is close to the former, which is often zero. However, for some cases, such as a solution with extended components, the solution in the wavelet space could be better than the one in the pixel space. We expect to compare different types of basis functions (pixels, different flavors of wavelets) in future work and see if a better map can be obtained.