The BMW Detection Algorithm applied to the Chandra Deep Field south: deeper and deeper A. Moretti1, D. Lazzati2, S. Campana1, G. Tagliaferri1 [email protected] 2 0 ABSTRACT 0 2 n a Chandra deep fields represent the deepest lookat theX–ray sky. We analyzed J the Chandra Deep Field South (CDFS) with the aid of a dedicated wavelet- 7 1 based algorithm. Here we present a detailed description of the procedures used to analyze this field, tested and verified by means of extensive simulations. We 1 v show that we can safely reconstruct the Log N–Log S source distribution of the 0 CDFS down to limiting fluxes of 2.4 × 10−17 and 2.1 × 10−16 erg s−1 cm−2 in 9 2 the soft (0.5–2 keV) and hard (2–10 keV) bands, respectively, fainter by a factor 1 0 ∼> 2 than current estimates. At these levels we can account for ∼> 90% of the 2 1–2 keV and 2–10 keV X–ray background. 0 / h p Subject headings: X–ray background – number-counts – detection algorithm - o r t s a 1. Introduction : v i X The Chandra observatory is providing the astronomical community with the deepest X– r a ray look at the sky (Mushotzky et al. 2000; Hornschemeier et al. 2000; Giacconi et al. 2001). Two 1 Ms observations have been recently carried out: one in the northern hemisphere on the Hubble Deep Field North (HDFN) and the other in the southern hemisphere on a field, named Chandra Deep Field South (CDFS), selected for its low column density (the southern twin of the Lockman hole) and for the lack of bright X–ray and optical sources. The data are available from the Chandra public archive (http://asc.harvard.edu/udocs/ao2-cdf-download.html) A further 1 Ms data set on the HDFN will be available in the next future. 1 Osservatorio Astronomico di Brera, Via E. Bianchi 46, Merate (LC), 23807,Italy. 2 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK. – 2 – The main goal of these observations is to look at the X–ray sky at the deepest level and togaininsight inthepopulationofemittingsourcescomprisingthecosmicX–raybackground (XRB). The details of the data reduction and analysis procedures that have been applied to manage the CDFS data-set are discussed in Section 2. In order to fully exploit the potential of these data, refined detection algorithm have to be used. We developed a wavelet-based source detection algorithm (Lazzati et al. 1999; for its main characteristics see also Section 3), that we applied to the full sample of ROSAT HRI fields (Campana et al. 1999; Panzera et al. 2002, in preparation). We modified this detection algorithm, called Brera Multi-scale Wavelet (BMW), to account for the specific characteristics of the Chandra ACIS Imaging and Spectroscopic instruments. The algorithm (BMW-Chandra) has been extensively tested in the extreme conditions provided by the CDFS (Section 4). Our final goal is to obtain a source detection and a source Log N–Log S at the faintest limits in order to resolve as much as possible of the XRB in point sources. In order to compare our results with previous investigations, we carry out the analysis in the soft (0.5–2 keV) band and in the hard (2–10 keV) band. We are able to safely reconstruct the Log N–Log S source flux distribution down to 2.4×10−17 and 2.1×10−16 erg s−1 cm−2 in the soft and hard bands, respectively (Section 5). A first account of these results have been given in Campana et al. (2001), here we extend these results discussing in more details the resolved background and source characteristics. Conclusions are reported in Section 6. 2. The data All exposures were taken with the Chandra X–ray Observatory (Weisskopf et al. 2000) AdvancedCCDImagingSpectrometer(ACIS-I)detector(Bautzetal. 1998). ACIS-Iconsists × of four CCDs arrangedin a 2 2 arrayassembled in an inverse shallow pyramid configuration to better follow the curved focal surface of the mirrors. The full ACIS-I has a field of view of 16.9′ × 16.9′. The on-axis image quality is ∼ 0.5′′ FWHM increasing to ∼ 3.0′′ FWHM at ∼ 4.0′ off-axis. The CDFS has been obtained performing eleven exposures of the same area of the sky with slight position and orientation offsets. A total exposure time of 966 ks is obtained by summing them together. In Tab. 1 we report the main properties of each exposure. First we have to match the position of the different exposures to add them into a single image. To this aim we performed a preliminary source detection on the 11 sub-frames at full resolution. We used the grid of the positions of the brightest sources to calculate the ten roto–translation matrices with respect to the first image. In all ten matrices we find that the rotation terms are negligible. The translation factors of each frame with their uncertainties – 3 – arereportedinTab. 1. Thereporteduncertainties arethestandarddeviationsoftheposition distributions of the grid sources after the translations. Starting with the level 2 event files, the data were filtered to include only the standard event grades (corresponding to the ASCA grades 0, 2, 3, 4, 6). In order to compare our results with previous studies and existing estimates of the cosmic X–ray background, we selected a 0.5–2 keV soft energy band and a 2–7 keV hard energy band to carry out the scientific analysis. The limit at 7 keV has been selected because at higher energies the effective area decreases and the background increases, resulting in a lower signal to noise of celestial sources. The fluxes are then extrapolated to the 2–10 keV energy band. Chandra observations are affected by the ‘space weather’, resulting in periods with high background. In order to eliminate these periods we calculated the total counts light curve ∼ for each exposure using a time resolution such that in each bin we had 400 events (i.e. ∼ 1,000 s in the soft band and ∼ 600 s in the hard band). We excluded all time intervals above 3σ of the mean. From a total of 966 ks, we excluded 24 ks in the soft band (28 ks in the hard band), giving a final effective exposure time of 942 ks (939 ks in the hard). The ∼ lost fractions amount to 3% of the total time, but allowed us to reduce considerably the background in each band. In particular, we registered a very strong flare in the observation ID1431 which contained about 90% of the total events both in the soft and hard band in about 10% of the exposure time. CCDs are usually affected by defects such as bad columns, bad and hot/flickering pixels and cosmic ray events. The great majority of these defects are already eliminated by the standard Chandra processing pipeline. Random flickering pixels (probably excited by en- ergetic cosmic rays) however do occur and have to be removed manually. These are active pixels for only 2–6 frames and might be easily detected as X–ray sources. Following Tozzi et al. (2001), we defined a pixel as “flickering” when it registers two events within the time of two consecutive frames and in each observation we eliminated all the events registered from that pixel (this can be done safely thanks to the lack of bright sources). In this way we excluded about 600 events in the soft band whereas in the hard band we did not find any flickering pixels. We also eliminated the data from chip 3 of the ID581 observation, because it suffered from Good Time Intervals inconsistences with other chips, resulting in an anomalous high background level. The count-rate to flux conversion factors in the 0.5–2 keV and in the 2–10 keV bands were taken from Tozzi et al. (2001). These were computed using the response matrices at the aim point and amount to (4.6 ± 0.1) × 10−12 erg s−1 cm−2 per count s−1 (0.5–2 keV), and (2.9 ± 0.3) × 10−11 erg s−1 cm−2 (2–10 keV) per count s−1 in the hard band, assuming a Galactic absorbing column of 8 × 1019 cm−2 and a photon index Γ = 1.4, i.e. – 4 – the spectrum of the XRB. Flux uncertainties are derived considering power law indices in the range Γ = 1.1−1.7. Fluxes are corrected for vignetting by using the exposure map (see below). The analysis has been carried out on images rebinned by a factor of 2 (i.e. 1 pixel corresponds to 0.98 arcsec), allowing us to deal directly with the entire image. We found that working on sub-images at the natural scale improves slightly the positional accuracy at the price of a four times longer computational time and with the problem of matching the four sub-images. In the analysis of the CDFS we have restricted our analysis to the central 8 arcmin radius circle, assuming as the center the aim point of ID581 observation. This is a good compromise between the simplicity of the geometry and the efficiency of the data analysis (see below). The minimum value of the exposure map in this region is about 23% of the maximum value. This occurs at the border of the circle where the exposure map has a very steep decrease. More than 90% of the exposure map within the 8 arcmin central circle have values larger than 80% of the maximum value, corresponding to an effective exposure time greater then 753 ks. The average background levels in this region are 0.18 and 0.29 counts per binned pixel in the soft and hard band, respectively. These correspond to 0.07 and 0.11 counts s−1 per chip and are in very good agreement with the expected value reported in the Chandra Observatory Guide. 3. The BMW algorithm The wavelet transform (WT) is a mathematical tool able of decomposing an image in a setofsub-images, eachofthemcarryingtheinformationoftheoriginalimageatagivenscale. These features make the WT well suited for the analysis of X–ray images, where the scale of sources is not constant over the field of view. In addition, the background is automatically subtracted since it is not characterized by any scale. The use of WT as X–ray detection algorithms was pioneered by Rosati et al. (1995; 1998) for the detection of extended sources in ROSAT PSPC fields and subsequently adopted by many groups (Grebenev et al. 1995; Damiani et al. 1997a; Pislar et al. 1997; Vikhlinin et al. 1998; Lazzati et al. 1998; Freeman et al. 2001). The BMW detection algorithm is a WT-based algorithm for the automatic detection and characterization of sources in X–ray images. It is fully described in Lazzati et al. (1999) and has been developed to analyze ROSAT HRI images producing the BMW-HRI catalog (Campana et al. 1999; Panzera et al. 2002). We have recently updated this algorithm to – 5 – support the analysis of Chandra ACIS Imaging and Spectroscopic images. Here we sum- marize the basic steps of the algorithm. First, the WT of the input image is performed. The BMW WT computation is based on the discrete multi-resolution theory and on the “´a trous” algorithm (Bijaoui et al. 1991). This is different with respect to algorithms based on a continuous WT, which can sample more scales at the cost of longer computing time (Rosati et al. 1995; Grebenev et al. 1995; Damiani et al. 1997a). We used a Mexican hat mother, which can be analytically approximated by the difference of two Gaussians (Slezak et al. 1994). We performed the Chandra ACIS-I data analysis with a rebin factor of 2 and we used the scales a = 1, 2, 4, 8, 16 and 32 pixels, where a is the scale of the transform (see Lazzati et al. 1999). Candidate sources are identified as local maxima above the threshold in the wavelet space at each scale. A catalog for each scale is obtained. Different scale catalogs are then cross–correlated and merged in a single catalog. For each scale the threshold in the wavelet space is calculated by means of Monte Carlo simulations: for a grid of values of Poissonian backgrounds in the range between 10−4 −101 counts per pixel we estimated the number of expected spurious detections as a function of the threshold value. The background value is measured as an average on the whole image by means of a σ-clipping algorithm. The number of spurious sources per field can be fixed arbitrarily by the user (typical values are in the range 0.1 − 10). At the end of this detection step we have a first estimate of the source position, an estimate of the counts from the local maximum value of the WT and a guess of the size from the scale where the WT is maximized. The final catalog is obtained through the characterization of the sources, which is per- formedbymeansofaχ2 minimizationwithrespect toaGaussianmodelsourceinthewavelet space. The WT of a detected source is used at three different scales, with the central being the scale at which the source has been found (corrections are introduced when the best scale is the first or the last). In order to fit the model on a set of independent data, the WT coef- ficients are decimated according to the scheme described in Lazzati et al. (1999). Neighbor sources can be fitted simultaneously in the deblending process, allowing the characterization of faint objects located near bright sources. 3.1. Total background map The WT detection is carried out on top of a map representing at best the image back- ground (named total background map, see Lazzati et al. 1999). The ACIS-I background consists of two different components, the cosmic X–ray background and cosmic ray-induced events (also called particle background). The first is made by focalized X–ray photons – 6 – from not resolved sources and therefore suffers from mirror vignetting. Its spectrum can be modelled as a power law with photon index Γ = 1.4. Following the CIAO 2.1 Science Threads (http://asc.harvard.edu/ciao/threads) we built the exposure maps for each of the eleven observations (both in the soft and in the hard band) assuming a photon index Γ = 1.4 as input spectrum. The second contribution has a spatial pattern which is basically due to the electronics of the system and depends on its temperature. Based on the Chandra ACIS background web page (http://asc.harvard.edu/cal/Links/Acis/acis/Cal prods/bkgrnd/current/index.html) we foundthatfor theACIS-I chips inthe energy rangeandinthe spatialregionofinterest for the only particle background we can assume a spatial flat pattern with an inaccuracy of less than 10%. We assumed that the particle background in the soft band is 50% of the total and the 66% in the hard band with a flat energy spectrum (Baganoff 1999). Summing the two components of the background we obtain a total background map and taking a Poissonian realization of it we can accurately account for the data in both bands as shown in Fig. 1. This procedure is slightly different from the one adopted in Campana et al. (2001), where a map based only on the cosmic X–ray background was considered. 4. Simulations We tested the detection algorithm, trying to reproduce the real data in terms of back- ground level and source flux distribution, matching the characteristics of the CDFS fields. The analysis is carried out in two energy bands: the soft (0.5–2 keV) and the hard (2–7 keV) band. From the detection algorithm point of view the two bands differ for two main reasons: the shape of the exposure maps and the different background level. This forced us to carry out two different sets of simulations. ′ × ′ X–ray sources were generated at random positions in a 30 30 square region around the center of the image (larger than the real field of view) and with fluxes distributed in the range 10−17 −10−13 erg s−1 cm−2 for the soft band and 10−16 − 10−13 erg s−1 cm−2 in the hard band. The number counts have been generated according to the Log N–Log S integral distributions reported in Tozzi et al. (2001), i.e. power laws with slope of 0.66 in the soft band and 0.92 in the hard band. The CDFS is made of eleven observations. For this reason we added a further element to our simulations: for each source we calculated its position in each observation and derived the corresponding off-axis angle (each source has different positions in the detector in each different observation, depending on the coordinate of the pointing, which are slightly differ- ent, and on the roll angle). Then we calculated the expected number of photons for each – 7 – observation. These depend on the conversion factor and on the local value of the exposure map. Finally, for each single observation we multiplied the expected number of photons with the corresponding point spread function image (from the Chandra CALDB) and we summed the eleven images to the background image. This procedure is aimed at reproducing at best the image characteristics. We used these routines instead of the MARX simulator essentially because of computa- tional reasons: we found that dealing with the input/output of data was easier for us using our IDL routines. We repeated the procedure for 100 times, for a total sample of about 150,000 input sources (here and in the following we refer to simulations in the soft band, unless otherwise stated; in the hard band results are similar). We recovered 25,000 of them within the central 8 arcmin radius, from which we were able to verify the performance of the procedure. 4.1. Spurious detection and position determination In the BMW detection algorithm the expected number of spurious sources is one of the fundamental input parameters. In our analysis we fixed the number of spurious sources to a total of 4.3 in the 8 arcmin radius circle (which corresponds to 6 expected spurious in a × 1024 1024 pixels image). In 100 simulated frames we found a total of 425 spurious sources (431 expected, see left panel of Fig. 2). This shows that the contamination of our sample is well understood and under control. This corresponds to sources detected at a significance level larger than 4σ. Source positions are recovered accurately. In the right panel of Fig. 2 we plot the differences between the input and output position of simulated sources. The difference distributions on both axes are well approximated by Gaussian functions centered on zero and with a r.m.s. of 0.6′′ (i.e. ∼ half a pixel at rebin 2, at which the analysis has been carried out). 4.2. Sky–coverage The distribution of the detected sources is usually depleted at the faint end. In order to compute the Log N–Log S it must be corrected for (i) detection probability lower than one close to the sensitivity limit (completeness function; CF) and (ii) the different areas effectively surveyed at different fluxes (sky coverage; SC). The CF is due to the fact that sources are preferentially detected if they sit on a positive background fluctuation, while – 8 – they are missed in the opposite case. The proper SC takes into account that the sensitivity limit can change in the same observation if the background or spatial resolution are not uniform or in different observations when different exposure times are considered. In single observations the SC depends on the off–axis angle, since the detector sensitivity decreases (lower effective exposure time) and the spatial resolution worsen. In the CDFS the geometry is complicated by the sum of observations with slightly different aim points and different roll angles. We have verified that this fact does not affect significantly the symmetry of the SC. More importantly, the SC has strong variations in those locations of the sky that were not imaged in all the eleven observations. In order to avoid large SC fluctuations, we restricted ′ our analysis to the inner 8 of the field, where the background and exposure are roughly constant (see above). Using WT based detection algorithms the SC can be computed either analytically or by means of simulations. The first approach is recommended when the intrinsic size of sources is poorly known (like in surveys of cluster of galaxies, see Rosati et al. 1995). In the CDFS the presence of extended sources is a minor effect for the computation of the Log N–Log S. For this reason we computed the SC by comparing the input and output number counts of ′− ′ ′− ′ simulated fields. Wedivided theselected field ofview into four circular regions (0 3, 3 5, 5′ − 6.5′ and 6.5′ −8′) and derived the SC and CF from simulations as a single correction (which we hereafter will call sky coverage). In the upper panels of Fig. 3 the derived CF as a function of flux and off-axis angle is shown for the soft and hard band; in the lower panels the total sky-coverages are shown in unit of square degrees. 4.3. Fluxes and bias correction In Fig. 4 we show the comparison between the input and output counts of simulated sources (output fluxes are corrected for the flux lost during the WT characterization, i.e. the so-called Point Spread Function (PSF) correction, evaluated on bright isolated sources for every arcmin of the field of view, and vignetting, evaluated directly from the exposure map). For high signal to noise sources the counts are measured correctly, while at lower signal to noise ratio the count distribution of the detected sources suffers from the well known Eddington bias (e.g. Hasinger et al. 1993). In our case the bias starts affecting the data below ∼ 30 input counts in the soft band (corresponding to ∼ 10−16 erg s−1 cm−2). This bias is due to the superposition of faint sources (below the detection threshold) on positive background fluctuations, making them detectable. Given the background level we expect to detect sources with less than 10 counts only in correspondence of background peaks (or when two faint sources are merged together into a brighter one). It follows that the counts of – 9 – the very faint sources are necessarily overestimated (otherwise they could not be detected). We corrected for this bias following the approach of Vikhlinin et al. (1995). We fitted the output source counts as a function of the input counts with polynomial functions for a grid of circular coronae, since the bias is a function of the source width (see Fig. 4). By inverting the fitting functions, we get an unbiased estimate of the source counts in the faint tail of the simulated sample. Input counts versus corrected output counts are shown in the right panel of Fig. 4. To assess the goodness of our correction procedure we perform a maximum likelihood fit to the differential (unbinned) flux distributions of the 100 simulated fields individually and we compare the results of the fits with the expected values from the input distributions. In Fig. 5 we compare the results before and after the bias correction. As explained below, in the final test section, for the corrected fluxes we set the flux limit at 5 and 7 counts in the innermost region in the soft and hard band, respectively, whereas for the uncorrected fluxes we are forced to use limiting flux of 10 and 12 counts, respectively. From the results of the maximum likelihood fit it is evident that after the flux correction results are in very good agreement with the input value, whereas the distribution of the uncorrected fluxes suffers for a bias. This effect is more evident in the hard band where the background level is higher. A possible alternative to the correction of the measured fluxes is the use of a complete- ness function allowing for a probability higher than unity for those fluxes where faint sources are wrongly detected. Such a probability function can be obtained by comparing the number of input and output sources in simulations as a function of their output fluxes instead of the (usually considered) input ones. In any case, this probability function will depend on the input Log N–Log S and it will be necessary to perform simulations with a number count distribution as close as possible to the unknown distribution of the sky sources (an iterative approach would be required). On top of that, the sensitivity limit of a survey adopting ∼ this procedure to correct for the Eddington bias will be a factor of 2 brighter than the correction of fluxes discussed above. Thus we followed the first procedure. 4.4. Final tests and flux limits In order to investigate any systematic error in our analysis we verified how the entire procedure works on different and independent simulated samples. Thus, we have generated different samples of sources with different slopes and normalizations of the integral flux distribution, comparing the input sample with the output of the detection procedure, after the correction for the flux bias and the sky coverage. Again, to assess the goodness of our analysis and to derive the flux limit, we perform a maximum likelihood fit to the differential – 10 – (unbinned) flux distributions. In Fig. 6 we compare the input integral flux distribution (dashed line) with the output one before (dotted line) and after (solid line) the correction for the Eddington bias. As stated above, without the flux correction for this bias, a reliable flux limit is about a factor 2 larger. This happens because all the detected sources with input fluxes in the range 5–10 counts range are detected typically with 10–15 counts. As a consequence the slope of the faint tail of the distribution of uncorrected fluxes is overestimated. This effect depends on the background level and on the slope of the distribution. It is more important in the hard band where the background is higher, and where the source distribution has a larger slope. As illustrated in Fig. 6 in all cases we are able to recover the input distributions with high accuracy and with different input slopes down to 5 and 7 counts in the soft and hard band, respectively (corresponding to fluxes of 2.4 × 10−17 and 2.1 × 10−16 erg s−1 cm−2). These numbers refer to the innermost regions and rising to 5, 6 and 8 counts in the other regions (7, 7 and 10 in the hard band). 5. Log N–Log S distribution ′ Within the 8 radius region in which our analysis has been carried out we detected 244 and 177 sources in the soft and hard band, respectively (Fig. 7). These sources are detected independently in the two images for a total number of 278 sources. This is at variance with e.g. Rosati et al. (2001) who made a detection over the entire 0.5–7 keV energy band and then pick up sources with signal to noise ratio larger than 2.1 in the 0.5–2 keV and 2–7 keV energy bands. Given a typical source spectrum and taking into account the different background levels in the two bands this latter approach is in general more efficient in the source detection, because it uses all the signal available, but it is not well suitable for our purposes. Infact, our final aim is drawing the flux distributions in the two bands down to the faintest limits, where the detection is highly incomplete and to do this we based on the results of the Monte Carlo simulations: if we perform the detection and the photometry in the two bands indipendently we can simulate all the procedures directly, whereas if we perform the detection in the full band we have to assume an input spectrum for the sources to simulate the photometry in the two bands making the evaluation of the sky coverage at the faint end tricky. To further characterize the detected sources we computed the hardness ratio HR = (H−S)/(H+S) where H and S are the net counts in the hard and the soft band corrected ∼ for PSF and vignetting losses, respectively. There are 34 sources ( 12% of the total number of detected sources) that are revealed only in the hard band (HR = 1), and 101 sources