INNOVATIONS IN THE ANALYSIS OF CHANDRA-ACIS OBSERVATIONS

Typically, the process of source detection is a distinct precursor to source extraction and analysis. Source detection algorithms must estimate some sort of significance for a putative source's signal with respect to the background expected to contaminate that signal. In other words, source detection algorithms must perform some type of extraction of a proposed source, estimate the local background, and compute some type of significance statistic from those two quantities.

Much of this paper describes the algorithms we have developed to perform those same three tasks—extraction (Section 5), background estimation (Section 5.4), and characterizing source significance(Section 7)—on complex ACIS campaigns involving multiple misaligned observations of crowded fields of point sources. We expect these careful algorithms to produce higher quality results for point sources than can possibly be produced by typical source detection schemes that operate on binned images, do not have knowledge of the Chandra PSF, and are not designed for multiple misaligned observations.

Thus, we have adopted and we recommend a source detection strategy that combines the convenience and efficiency of traditional source detection tools with the accuracy achieved only by detailed extraction and analysis of an entire catalog. First, a liberal catalog of candidate point sources is obtained from a variety of traditional source detection methods, run with aggressive thresholds that seek to find weak sources but are expected to produce significant numbers of spurious detections. We extract and characterize the significance of the candidates, and then prune those found not to be significant. Nandra et al. (2005) independently developed a similar strategy.

The trimmed catalog is then extracted again, re-pruned, and so on until all the remaining sources are deemed significant. This iterative pruning procedure, depicted in Figure 1 as a loop through the steps "point source catalog," "point source extraction," and "prune/reposition sources," is necessary because construction of both source apertures (Section 5.1) and background regions (Section 5.4) in crowded fields can be more appropriately performed when the existence of neighboring sources is taken into account. In contrast, typical detection algorithms estimate background for a putative source using the nearby pixels in an image, without explicit regard for possible contamination by nearby sources. Typically, only one iteration is needed for situations where source crowding is unimportant.

We test the significance of candidate sources in three energy bands (0.5–2 keV, 2–7 keV, and 0.5–7 keV) and adopt the highest value. We prefer the upper limit of 7 keV instead of the more typical value of 8 keV because there is a line in the instrumental background spectrum at 7.47 keV. When a source has been observed multiple times, we adopt the highest significance among all combinations of observations (Section 6.2).

We construct the list of candidate sources using a combination of several methods. Many sources can be identified by a wavelet-based detection algorithm for Poisson images that has been widely used on Chandra data—the CIAO tool wavdetect⁴³ (Freeman et al. 2002)—run with the liberal threshold of 1 × 10⁻⁵ rather than the commonly used level of 1 × 10⁻⁶ in order to nominate as many actual sources as possible without nominating an excessive number of candidates that will later be found to be insignificant. For a simple field with one pointing, 12 images are searched—corresponding to three energy bands (0.5–2 keV, 2–7 keV, and 0.5–7 keV) and four image pixel sizes (025, 05, 1'', and 2'') to sample appropriately the variable Chandra PSF. All the wavdetect catalogs are merged—pairs of catalogs are matched using the algorithm described in Section 8, and the most accurate position from each matching pair is retained. The resulting catalog is visually examined to remove any obvious duplicates remaining.

Wavelet decomposition of images is often not effective in resolving closely spaced sources, as the method is designed to find structures on a range of spatial scales. For such situations, image reconstruction algorithms that remove the blurring effects of the PSF can significantly improve source identification. We search for faint and crowded sources by locating peaks in reconstructed images obtained with the Lucy–Richardson algorithm (Lucy 1974), implemented in the IDL Astronomy User's Library.⁴⁴ Similar image reconstructions have been invaluable for achieving effective spatial resolution as good as ∼03 FWHM in observations of Chandra targets, for example, a jet in the Chandra first-light target (Chartas et al. 2000) and SN 1987A (Burrows et al. 2000). As the PSF varies strongly across the Chandra field due to the telescope optics, the image reconstruction is performed on many small overlapping tiles (Figure 3) using local PSFs. The candidate sources identified in these tiles are merged with those obtained with the wavdetect algorithm.

Zoom In Zoom Out Reset image size

Figure 4 shows an example of the effectiveness of image reconstruction in regions with crowded, faint sources. The wavdetect procedures located 50 sources in this sub-image of the center of the young stellar cluster Trumpler 14; an additional 50 are found as peaks in the reconstructed images (Townsley 2006; L. K. Townsley et al. 2010, in preparation). The reliability of most of these sources is confirmed; 89 of the 100 reconstruction sources coincide with stars detected in deep, high-resolution near-infrared exposures (T. Preibisch 2009, private communication; Ascenso et al. 2007), some with separations <1''. Of the 11 unconfirmed X-ray sources, nine appear in close pairs where the other member of the pair is confirmed by an IR counterpart. Testing the validity of these X-ray sources will require very high resolution, sensitive IR data, as it may be difficult to find faint, close companions in IR observations of such crowded regions suffused by diffuse IR emission. Since the X-ray flux of a young star does not correlate closely with its IR brightness, it is not unreasonable for us to find close X-ray pairs that are not seen in IR images, and vice versa.

Zoom In Zoom Out Reset image size

When a project's scientific goals include the search for faint X-ray emission from a previously defined catalog of interesting objects found in other wavebands, the candidate X-ray source catalog can be supplemented with source positions that are not derived from the ACIS image. The AE procedures can then judge whether significant emission is present from each single source, and the positions can be "stacked" for high-sensitivity examination of the collective properties of the user-provided catalog (Section 7.7).

Finally, we often supplement the candidate source catalog with the coordinates of suspected point sources identified from visual review (e.g., with the DS9 visualization program) of the X-ray data. This visual review also allows removal of obviously spurious candidate sources that occasionally emerge (e.g., from detector artifacts such as bright source readout streaks) and duplicate sources.

In typical source detection schemes, the existence of a candidate source is evaluated using the signal-to-noise ratio (S/N), which is a photometry value divided by its uncertainty. When photometry values have Gaussian distributions, then an S/N threshold directly corresponds to a level of significance in a statistical test of the null hypothesis that there is no source signal present, only background.

However, most ACIS sources are quite weak with very few counts extracted from the source aperture, and in crowded fields sometimes very few counts available to estimate the background. Gaussian approximations to photometric confidence intervals in such cases can be quite poor. We prefer to test directly the null hypothesis that a candidate source does not exist using the method described by Weisskopf et al. (2007, Appendix A2) based on the Poisson distribution. Computation of this significance statistic, P_B, is described in Appendix B.

The observer must set, either a priori or after careful examination of the image and possible multi-wavelength counterparts, a threshold value of P_B that strikes a reasonable balance between the competing goals of low false detection rates and high sensitivity. Analytical methods to estimate false detection rates and sensitivity have been developed (e.g., by Nandra et al. 2005; Georgakakis et al. 2008), however, extending these methods to studies that involve crowded source apertures (Section 5.1) and multiple overlapping pointings (Section 6.2) is not straightforward. Monte Carlo simulation of the detection process (e.g., by Cappelluti et al. 2007; Kim et al. 2007) is a powerful technique to study false detection rates and sensitivity, however, simulation of our complex detection process would require unreasonable quantities of both human and computing resources. Thus, we do not currently have an objective or authoritative way to set detection thresholds. As with most source detection tasks, the scientist's subjective judgment is critical for setting criteria for source existence, and high-quality observations in other wavebands are invaluable for informing that judgment.

Many Chandra observers obtain extraction apertures for their point sources either from the CIAO tool wavdetect⁴⁵ (Freeman et al. 2002; the most popular source detection tool in the Chandra community), or by drawing the apertures by eye.⁴⁶ Because wavdetect does not have explicit knowledge of the Chandra PSF, and because it seeks to detect both point-like and extended sources, it sometimes produces aperture shapes that are radically different from the PSF (e.g., ellipses with extreme eccentricity). An additional concern with apertures from wavdetect or visual inspection is the risk of introducing an upward bias into photometry because both tend to "follow the light"—including only regions of the image where, by chance, the source produced an excess of events over the background. In our opinion, a region derived from the local PSF is the most objective and appropriate extraction aperture for point sources.

When a source is "extended" (resolved by Chandra) the spatial distribution of its observed events does not follow the local PSF, and application of the point source extraction techniques in this section (Section 5) will underestimate the flux of the source. Such sources are more properly designated as "diffuse" and should be extracted by procedures that do not make the point-like assumption (see Section 9). No automatic procedure for defining the extraction aperture of such a source seems feasible—one must either define an aperture that follows the observed light (as wavdetect does), or apply a priori information about the structure of the source to define a suitable aperture. Observers often struggle to decide if a source is extended. The Chandra Source Catalog provides a sophisticated analysis of source extent⁴⁷, and the CXC has recently introduced a science thread⁴⁸ for measuring source extent. AE does not currently address this task.

Several observers employ elliptical approximations to the PSF as extraction apertures, e.g., Nandra et al. (2005) and the Chandra Source Catalog⁴⁹ use ellipses that enclose 70% and 90% of the PSF power, respectively. AE extraction apertures are built from contours of the local PSF at 1.5 keV, as shown in Figures 5 and 6. By default, AE apertures enclose ∼90% of the PSF power, which we have found is a reasonable trade-off between maximizing the source's signal and minimizing the background's signal in typical fields. AE does not currently attempt to optimize the aperture size based on the source and background levels (e.g., use a large aperture for bright sources and a small one for sources just above the background).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For many Chandra observations, some of these default apertures will overlap significantly due to source crowding. AE iteratively reduces the aperture sizes of crowded sources until the apertures no longer overlap (Figure 6). The brighter member of a pair of crowded sources maintains its default aperture until the aperture for the weaker member has been driven down to a minimum allowed size (enclosing ∼40% of the PSF power). Of course, some light from the wings of nearby sources may contaminate even a reduced source aperture; this light constitutes an additional background component for the source. AE provides a sophisticated background algorithm that models and subtracts this component (Section 5.4).

In a particular observation, a source may be so highly crowded that assigning minimal apertures to it and its neighbor cannot prevent overlap. When multiple observations of the source are available, such overlapping extractions are discarded by AE (Section 6.2) in recognition that there are limits to our ability to estimate backgrounds (Section 5.4) under conditions of extreme crowding. If all extractions of a pair of candidate sources are overlapping, then our policy is to eliminate the weaker candidate in recognition that either the detection process has mistakenly split a single physical source into two candidates, or that two physical sources are not resolved by our observation (i.e., meaningful source properties cannot be determined for each source).

A PSF model and an extraction region must be constructed independently for each observation of a source. A pair of sources may be well separated with nominally sized apertures in an on-axis observation, but may suffer severe crowding with reduced apertures in an off-axis observation. Even a source observed twice at similar off-axis angles may have apertures with very different shapes, since the azimuthal shape of the Chandra PSF varies across the focal plane.⁵⁰

Once apertures are determined for each observation of a source, a mostly routine extraction process is followed for each observation using standard CIAO tools—dmextract selects the events within the aperture and constructs a Type I HEASARC/OGIP-compatible⁵¹ "source" spectrum; mkarf queries the Chandra Calibration Database and constructs an ancillary reference file (ARF); mkacisrmf queries the Chandra Calibration Database and constructs a response matrix file (RMF). When an aperture dithers across multiple CCDs, standard extraction methods⁵² will overestimate the source flux (by as much as 50%) because the response at the source position of only one CCD is computed (using a single call to mkarf). AE avoids this miscalibration of the extraction by computing and summing the response at the source position of all the CCDs (using multiple calls to mkarf).

Perhaps the most important calibration facilitated by the PSF model is accounting for the point source light falling outside the aperture and not extracted—a standard part of optical and infrared data analysis commonly referred to as "aperture correction." Such correction is necessary because the Chandra Calibration Database effective area data assume an infinitely large detector and extraction aperture.⁵³

Since the size of the Chandra PSF varies significantly with photon energy,⁵⁴ the aperture correction is a function of energy. AE builds images of the local PSF at five monochromatic energies (Appendix C), and calculates for each PSF image the fraction of the power that falls within the aperture. Those five "PSF fractions" are interpolated to estimate a PSF fraction at every energy. The product of this PSF fraction curve and the nominal ARF produced by CIAO represents the true observatory effective area that is supplying photons within the aperture; this corrected ARF is carried forward in the analysis. The energy and off-axis dependence of this aperture correction is illustrated in Figure 7.

Zoom In Zoom Out Reset image size

Omitting an aperture correction has two scientific effects. First, flux estimates are biased downward (by ∼10% for the examples in Figure 7 for a typical observed spectrum peaking near 1.5 keV). Second, the inferred shape of the source spectrum is somewhat distorted, especially for off-axis sources, because the PSF fraction varies with energy (Figure 7, lower panel).

Assessment of the detection significance and the spectral properties of a weak source depend critically on estimation of unbiased background spectra for each extraction (each observation) of the source. Often, background estimates must be performed "locally" for each source to account for spatial variations due to diffuse emission, or wings of nearby point sources. AE supports three methods of constructing local background spectra.

• 1.  
  When all sources are well separated, a traditional approach using annular background regions is adequate. AE implements this approach by first masking (removing) virtually all the point source light from the data set using circular mask regions with a nominal radius that is 1.1 times a radius that encloses 99% of the PSF. Circular background regions are then defined independently for each source to encompass a number of background events specified by the observer.
• 2.  
  AE also provides an alternative masking algorithm that first models the surface brightness expected from the point sources (using the PSF models and rough estimates of source fluxes), and then iteratively seeks to mask (remove) all portions of the field where the ratio of that point source glow to the local background level exceeds a threshold. This algorithm allows the large wings of bright sources to be removed with large masks, while avoiding excessive loss of background area around weak sources, as shown in Figure 2. Circular background regions are then defined independently for each source to encompass a number of background events specified by the observer.
• 3.  
  When sources are crowded, the concept of creating background spectra from a data set cleaned of all point source light is not appropriate. By definition, a significant component of the background for a crowded source arises from the wings of its neighbors; those wings must be modeled.The ultimate analysis strategy would conceptually separate a source's background into a point source wing component and a traditional smooth component not associated with detected point sources. The latter would be estimated via the masking techniques described above, and the former would be modeled via a complex multi-source spatio-spectral modeling process. Lacking the resources to implement that strategy, we have implemented a less complex strategy that seeks to estimate both components within each source's aperture by carefully sampling the observed event data.AE first constructs a spatial model for every source (using the PSF models and rough estimates of source fluxes), and then iteratively constructs a background region for each source that seeks to sample fairly the light that each neighbor's wing is depositing into the source aperture. In the simple case where the source has one neighbor that is contributing background to the source aperture, the algorithm will tend to build a background region that is an annulus around the neighbor, as shown in Figure 8. As the background region grows during the iteration, it will seek to avoid the wings of more distant sources that are not contributing light to the source aperture; thus for uncrowded sources this algorithm produces similar background regions as algorithm No. 2. AE tries to balance the competing goals of acquiring the desired number of background events, fairly representing each background component (wings from each neighboring point source, diffuse emission, and instrumental background), and sampling the background locally (acquiring background events in a compact region around the target source).

Zoom In Zoom Out Reset image size

For either masking-based approach, the observer can manually specify masks to remove portions of the field that are not suitable background estimators, such as CCD readout streaks associated with bright sources. For the spatial modeling-based approach, the observer can supply a spatial model for background structures that are not point sources, such as CCD readout streaks; such structures will be avoided or included in individual background regions to the extent that the source aperture is contaminated by them.

Regardless of how the background region is defined, we believe that the scaling traditionally applied to the background spectrum—the ratio of the geometric areas of the source aperture and background region—is not appropriate for ACIS data because the ACIS exposure map contains significant small-scale features caused by the dithering of CCD edges and bad columns. For example, a source lying in a so-called chip gap may have an effective exposure time that is less than half of nominal, while the majority of its background region may have nominal exposure. AE accounts for these exposure map features by adopting a background scaling equal to the ratio of the integrals of the exposure map over the source aperture and background region.

In principle, X-ray properties for a source could be estimated directly from the multiple extractions of the source (multiple source spectra, background spectra, ARFs, and RMFs). For example, photometric quantities could be estimated by maximizing a likelihood function involving terms for all of the observations (source and background extractions). Similarly, a spectral model for the source could be derived by simultaneously fitting all of the extractions. However, this approach is impractical for the majority of faint Chandra sources, where a few counts are spread over a number of separate observations.

AE adopts an alternate strategy—multiple extractions of a source are merged into "composite" data products (PSF model, source spectrum, background spectrum, ARF, and RMF), which are carried forward for analysis. Source spectra and exposure times are summed. Weighted averages of the PSFs, ARFs, and RMFs are computed using IDL and the FTOOLS addarf and addrmf. In order to support "grouping" of low-count spectra at later stages of the analysis, the composite background spectrum is represented in the same way as the source spectrum, as an integer "counts" vector (the sum of the extracted background spectra) with an associated scaling value (to account for the different sizes of the source and background regions), rather than as a real-valued "rate" vector with uncertainties on each rate value.⁵⁶ Extraction data products from each observation are retained so that the observer can analyze them separately (or simultaneously) if desired.

An important constraint on the design of the individual extractions arises from the decision to sum the integer-valued background spectra, namely, that all the individual background regions must be designed to have similar scaling values. This constraint arises because each individual background spectrum must be fairly represented in the composite background in order for Poisson $(\sqrt{N})$ uncertainty estimates to apply to the composite background. This can be understood by considering an extreme example. Imagine a background channel in which observation No. 1 produced 25 counts with a scaling of 25 and observation No. 2 produced 10,000 counts with a scaling of 10,000. The scaled composite background rate would be (25/25) + (1 × 10⁴/1 × 10⁴) = 2.0. However, the Poisson uncertainty estimated (by a fitting package) for that background rate would be far too small, since it is based on the notion that 10,025 counts were observed. In fact, the 25-count extraction totally dominates the uncertainty on the scaled composite background, which actually has an approximate uncertainty of $\sqrt{(\frac{\sqrt{25}}{25})^2 + (\frac{\sqrt{10,000}}{10,000})^2} \simeq 0.2$. AE includes code and workflow recommendations to ensure that all extractions of a source result in similar scaling of their background spectra.

When multiple observations of a field are available, astronomers commonly choose to disregard observations that are judged to be unhelpful to whatever astrophysical measurement is desired. For example, in ground-based optical studies, exposures with unusually bad seeing or unusually high background might be ignored. In all ACIS studies, observers are encouraged to discard data obtained during periods of very high instrumental background⁵⁷ (due to solar activity).⁵⁸ In ACIS studies involving multiple misaligned pointings, observers often choose to ignore data taken far off-axis (where the Chandra PSF is significantly degraded⁵⁹) if on-axis coverage is available (e.g., Plucinsky et al. 2008). ACIS observers also commonly choose to search for sources both within each observation separately and within observer-defined combinations of observations (e.g., Muno et al. 2009).

AE tries to formalize this common practice of ignoring unhelpful observations by implementing three data-selection strategies that are applied independently to each source.

• 1.  
  When the observer is interested in the validity of proposed sources using AE's P_B statistic (Section 4.3), then we recommend an AE option that selects whatever subset of the available extractions that optimizes (minimizes) P_B for each source. Under this option, an extraction will be included only when the particular signal and background it contributes lowers P_B. This option is appropriate when the observer wishes to adopt the scientific policy that a source is deemed to exist if it is significant in any observation, or in any combination of observations. For highly variable astrophysical sources, such as young stars, this is a reasonable strategy even when multiple observations have similar PSFs. The price paid for increased sensitivity to variable sources is an increased false detection rate from the additional number of random "trials" that can produce spurious detections.
• 2.  
  When the observer is interested in the position of sources, then we recommend an AE option that selects whatever subset of the available extractions optimizes (minimizes) the expected position uncertainty (Section 7.1) for each source. An extraction much farther off-axis than its peers will be included only when the disadvantage of its larger PSF⁵⁹ is outweighed by the advantage of averaging over more counts. Since we register each of our observation to an absolute astrometric reference frame, all observations of a source are well aligned and it is appropriate to estimate the position using only the best data we have.
• 3.  
  When the observer is interested in time-averaged photometric properties (e.g., fluxes, spectra) then selecting the extractions to merge becomes more problematic. Anytime extractions are discarded in order to optimize a photometric quantity (e.g., S/N) a bias can be introduced into all photometric properties because the discarded extraction may have, by chance, lower observed flux than the long-term average. Thus, the observer must balance two undesirable outcomes: sources whose photometric accuracy is degraded by including very poor quality extractions, and sources whose photometry suffers the suspicion of bias because some extractions were discarded after examining their data. AE offers an option that strikes this balance by discarding extractions only when retaining them would drive the S/N of the merged data set significantly below the optimal S/N. The observer specifies the minimum acceptable ratio between the S/N achieved by the merge and the optimal S/N achievable by discarding more extractions. This strategy tolerates limited deterioration to the S/N in order to avoid photometric bias arising from data selection. Bright sources will tend to incorporate all observations since background is unimportant, whereas weak sources will tend to reject observations that have backgrounds much higher than their peers (e.g., those with much larger apertures).

While initial source positions are provided by a source detection process, these may not be optimal. AE provides three procedures for improving source positions. One position estimate is constructed by calculating the centroid of the extracted events. A second estimate is determined from the peak of the spatial correlation between the Chandra-ACIS PSF (Appendix C) and the events in a neighborhood around the source. A third position estimate is obtained from the location of the peak in a maximum-likelihood (Lucy 1974) image reconstruction of the source neighborhood. In a multi-observation reduction, the correlation and reconstruction operations are performed using a composite (multi-observation) event image and a composite PSF image because a given source may be undetected in an individual observation.

In our experience, all three methods give very similar positions for most on-axis sources; the centroid position is simplest to calculate and is often used. The PSF correlation position is best for sources far off-axis (θ ≳ 5'); centroid positions are biased estimates of the true locations due to the asymmetry of the off-axis PSFs. For very closely spaced sources with overlapping PSFs (Figure 6), the best positions are provided by the maximum-likelihood reconstructed images. Since the centroid position is calculated using the extracted events, it should be re-calculated after a source is repositioned to verify convergence.

Positional errors are estimated separately along the right ascension and declination axes according to

$$\begin{eqnarray} \sigma _{\rm R.A.} &=& \sigma _{\rm model,R.A.}/\sqrt{N}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \sigma _{\rm Decl.} &=& \sigma _{\rm model,Decl.}/\sqrt{N}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \hspace{10pt}\sigma _{r} &=& \sqrt{\sigma _{\rm R.A.}^2 + \sigma _{\rm Decl.}^2}, \end{eqnarray} \tag{ 3 }$$

where σ_model,R.A. and σ_model,Decl. are the standard deviations along the right ascension and declination axes of a model of the counts falling within the extraction region. That model consists of the PSF projected onto each axis plus a flat background component appropriately scaled, both truncated by the source's extraction region; the spatial distribution of the observed data is not considered.⁶⁰ N is the total number of extracted counts. The quantity σ_r, commonly known as the distance root-mean-square error, is generally reported as the "1σ" error, or 68% confidence interval, for the source position.

Standard ACIS spectral files divide the energy range covered by the instrument into many hundreds of pulse-invariant (PI) energy channels, many of which have very few counts for typical sources. Energy channels are therefore frequently grouped. AE offers only one of the several grouping criteria found in standard grouping tools—groups are constructed to have similar S/Ns—but improves on standard implementations in several important ways:

• 1.  
  Standard algorithms provide no convenient way for the observer to define precisely the energy range over which spectral models will be fit. AE allows the observer to specify the boundaries of the first and last groups. For example, if the energy range 0.5–8.0 keV is specified, then the first group will span all the channels below 0.5 keV and the last group will span all the channels above 8.0 keV; these terminal groups can be subsequently "ignored" within a fitting package.
• 2.  
  The criterion AE uses to define a group is a target S/N for the net counts in the group. Standard grouping tools do not consider the background spectrum, and thus produce background-subtracted grouped spectra with S/N values lower than the target.
• 3.  
  Standard algorithms are asymmetric, filling groups from low to high energies. For sources with fewer than several hundred counts, groups often begin with a string of empty channels and end on a channel that happens to have an observed event. This permits anomalies in the grouped spectra, such as narrow-width groups with overestimated flux adjacent to wide groups with underestimated flux. These distortions can bias the spectral fitting and inflate the best-fit χ² value. AE's grouping algorithm attempts to mitigate this problem by selecting group boundaries mid-way in the run of empty channels (if any) that lie between events belonging to adjacent groups.

X-ray spectra observed with ACIS energy resolution typically show some combination of continuum components, line emission, and absorption features from interstellar matter. Accurate characterization requires fitting with nonlinear astrophysical models (Section 7.5). However, this fitting procedure does not give meaningful spectral parameters for faint sources when a wide range of models are statistically consistent with the sparse data set. Nonparametric estimators of spectral shape—statistics that describe the shape of the observed spectrum without assuming any astrophysical model—are commonly used for weak sources.

The oldest and most widely used statistic characterizing the shapes of X-ray spectra is the "hardness ratio" involving various ratios of soft and hard band counts. These estimators (known as "Poisson proportions" in the statistical literature) are mathematically unstable for sources with extreme spectra where the numerator or denominator has few or zero counts (Brown et al. 2001). Hong et al. (2004) demonstrated that quartiles (three values of a random variable that divide a population into four equal groups) of the observed event energies are superior characterizations of spectral shape for low-count observations, avoiding the problem of empty or nearly empty pre-defined energy bands. Quartile–quartile plots can be mapped to astrophysical parameters for specified simple model families.

AE provides estimates of the 25% quartile, median (50% quartile), and 75% quartile of the observed event energies over a variety of bands. These statistics are background-corrected, meaning that they seek to characterize the observable spectrum of the astrophysical source as if background were not present. AE implements an intuitive and straightforward background correction for standard quartiles, based on the observed cumulative distribution of the net spectrum, as shown in Figure 9. This method appears to be equivalent to that described by Hong et al. (2004, Appendix C), which was developed independently. AE does not yet estimate individual confidence intervals for these background-corrected quartiles.

Zoom In Zoom Out Reset image size

Our own Chandra studies of young stars have chosen to use the median energy statistic, rather than multiple quartiles, for deriving intrinsic properties of absorbed thermal plasmas (Section 7.4). Many sources in these studies are so weak that dividing the few observed counts into multiple quartiles is not feasible.

The common procedure for obtaining broadband source fluxes and luminosities (given astronomical distances) is to fit the spectrum derived above with astrophysical models convolved with the ARFs and RMFs representing the observatory response. This path (which we follow in Section 7.5 for strong sources) is less reliable for faint sources with insufficient photons to constrain the nonlinear fitting process. We describe here a photometric procedure to estimate photon fluxes.

AE begins by performing standard aperture photometry, computing "net counts" quantities for various energy bands, S(E), by subtracting a scaled background from the counts observed near the source.

$$\begin{eqnarray} &&S(E) = C^s(E) - (A^s/A^b) C^b(E). \end{eqnarray} \tag{ 4 }$$

Here C^s(E) is the number of counts observed in the source aperture in energy band E and C^b(E) is the number of counts observed in the background region in energy band E. The source and background "areas," A^s and A^b, are derived by integrating exposure maps rather than computing geometric areas (Section 5.4). Note that these raw "net counts" photometric quantities are not adjusted to correct for finite apertures. Instead, an energy-dependent aperture correction (Section 5.3) is applied to the calibration of the extraction (the ARF file); thus, any calibrated photometric quantities (such as the flux estimates below) will include the aperture correction.

Asymmetric confidence intervals for both C^s and C^b are estimated using the common analytical approximations to confidence intervals of the Poissonian distribution constructed by Gehrels (1986, Equations (7) and (12)). These confidence intervals are propagated through Equation (4), using the method described by Lyons (1991, Equation (1.31)), to estimate an asymmetric confidence interval for S(E). This technique is not ideal, and we anticipate adopting the Bayesian technique for estimating confidence intervals used by the Chandra Source Catalog,⁶¹ which is implemented in the CIAO tool aprates.

These standard "net counts" quantities are used to construct two estimates of incident photon flux onto the Chandra telescope, designated here F^⋆ and F^#, in units of photon cm⁻² s⁻¹ (not the usual erg cm⁻² s⁻¹). For the first estimate, the net count rate is divided by the observatory response in each of many narrow channels, and these nearly monochromatic photon flux densities are summed over a chosen energy band:

$$\begin{eqnarray} F^\star _{\rm phot}(E) &=& \frac{S(E)}{{\rm EXPOSURE} \times {\rm ARF}(E)}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} F^\star _{\rm phot}(E_{\rm min}<E<E_{\rm max}) &=& \sum _{E_{\rm min}}^{E_{\rm max}} F^\star _{\rm phot}(E), \end{eqnarray} \tag{ 6 }$$

where S(E) is the net observed counts (Equation (4)), ARF(E) is the observatory effective area (Section 5.2), EXPOSURE is the source exposure time, and E_min and E_max are the energy range provided by the user.⁶² Commonly used energy bands are E_min = 0.5 keV and E_max = 2 keV for the Chandra "soft band" and E_min = 2 keV and E_max = 8 keV for the Chandra "hard band."

For the second photon flux estimate, the net count rate in a user-supplied band is divided by the observatory response averaged over the band:

$$\begin{eqnarray} {\rm ARF}(E_{\rm min}<E<E_{\rm max}) &=& \sum _{E_{\rm min}}^{E_{\rm max}} {\rm ARF}(E), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &&F^\#_{\rm phot}(E_{\rm min}<E<E_{\rm max}) \nonumber\\ &&\qquad =\frac{\sum _{E_{\rm min}}^{E_{\rm max}} S(E)}{{\rm EXPOSURE} \times {\rm ARF}(E_{\rm min}<E<E_{\rm max})}.\qquad \end{eqnarray} \tag{ 8 }$$

The F^⋆ estimator can suffer from large Poisson errors because events (either source or background) at energies where the ARF is very small have a large effect on the estimator, and is thus not recommended for weak sources. For example, an event at 8 keV where the ARF value is tiny makes a much larger contribution to F^⋆ than an event at 2 keV where the ARF value is large. The F^# estimator suffers from a systematic bias with respect to the true incident photon flux because the effective area averaged over the energy band (Equation (7)) is the correct calibration of the net counts photometry only for the non-physical case of a flat incident spectrum.

When a priori information provides constraints on the shape of a source's intrinsic spectrum, these photon flux estimates F^⋆ and F^# can be combined with the background-corrected median energy (Section 7.3), E_median, and with the astronomical distance d to estimate apparent source luminosities L_x in a chosen broad energy band:

$$\begin{equation} L_x \approx 4 \pi d^2 [F^\star {\rm or } F^\#] E_{\rm median}. \end{equation} \tag{ 9 }$$

Getman et al. (2010) have studied in detail the accuracy and precision of this estimate for both bright and faint sources, particularly when absorption from line-of-sight interstellar matter is present. E_median is an accurate, though nonlinear, predictor of interstellar column density N_H when the intrinsic family of spectral models (e.g., power law, thermal plasma) for the source is known (Feigelson et al. 2005; Getman et al. 2010, Figure 4). E_median thus enters the calculation twice, once to scale the photon flux to the observed luminosity in Equation (9), and again to scale the observed luminosity to the intrinsic luminosity (corrected for absorption; Getman et al. 2010, Figure 5).

Simulations (Getman et al. 2010) indicate that these nonparametric estimates of observed luminosities and absorption column densities are quite accurate, with only moderate biases and statistical errors. Systematic errors are larger for estimates of the absorption-corrected luminosities, and can be extremely large for heavily absorbed sources when the chosen spectral band includes the soft X-ray regime. In our work on faint sources, we choose to use the more statistically stable F^# photon flux estimator and prefer absorption-corrected luminosities in the hard band (2–8 keV, which is less vulnerable to errors in absorption correction) over luminosities in the total band (0.5–8 keV).

AE relies on the widely used XSPEC package⁶³ (Arnaud 1996) for spectral fitting. XSPEC provides a wide range of intrinsic spectral models (such as non-thermal power laws and thermal plasmas), absorption by intervening material, and convolution with the Chandra/ACIS instrumental response. Models derived from χ² minimization become inaccurate for sources with few counts, and the problem is exacerbated when a significant fraction of the extracted events are background that is subtracted. Spectral fitting of faint sources is thus usually pursued using the C-statistic applied to unbinned data, an application of the Likelihood Ratio Test under the assumption that the data follow the Poisson distribution (Cash 1979). Background subtraction is not possible for likelihood-based statistical procedures; instead, a background model must be included in the spectral model that is fit to the data (source plus background) extracted from the source aperture, and that background model should be simultaneously fit to the data extracted from the background region. This fitting procedure is depicted in Figure 10 as a data flow diagram.

Zoom In Zoom Out Reset image size

XSPEC implements such a background model, which was derived by Wachter et al. (1979), for use with the C-statistic.⁶⁴ This model consists of a free parameter representing the background flux in each of the hundreds of spectral channels used in an ACIS spectrum. Such a model raises two theoretical concerns. First, in most cases where the C-statistic would be used, there are many more free model parameters than background events. Second, such a model is completely unconstrained—extremely narrow and quite non-physical features (one channel wide) in the incident background spectrum are presumed to be possible. Indeed, observers have found that this algorithm sometimes provides very poor "best fits" to the observed spectrum, with incorrect normalizations and non-physical spectral parameters.⁶⁵

We found that substituting a more constrained model for the ACIS background stabilized the algorithm and lessened the problem of non-physical fits. For sources with relatively few (∼100) background counts, we use a continuous piecewise-linear function with up to 10 vertices, which we call the cplinear model;⁶⁶ an example is shown in Figure 11. Figure 12 gives an example where the background model of Wachter et al. (1979) gives a poor fit and the cplinear model gives a good fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The highly constrained cplinear model is inadequate to model high-quality background spectra which exhibit complex structure. More physically based background models would be useful to the ACIS community; our goal here is to emphasize that severely underconstrained background models should be avoided, rather than to claim that cplinear is optimal. The ultimate solution to modeling low-count X-ray spectra may lie in Bayesian approaches, such as the innovative technique of van Dyk et al. (2001). Those authors also suggest that background models with appropriate functional forms are useful, however, the background model presented in that published work employed a free parameter for each channel in the spectrum (van Dyk et al. 2001, Equation (22)), similar to the Wachter model.

Another algorithmic innovation to spectral fitting provided by AE is the improved method for grouping spectra (either for fitting with the χ² statistic or for plotting ungrouped spectra fit with the C-statistic) described earlier (Section 7.2). Finally, AE also offers the observer several procedural conveniences. XSPEC scripts that drive fits to three common astrophysical models (absorbed one- and two-temperature thermal plasmas, absorbed power law), using either χ² or the C-statistic with the cplinear background model, are provided. AE automates the execution of these fitting scripts; fitting multiple models to thousands of sources is feasible. When multiple models have been fit to a source, the observer can visualize those models and select the most appropriate one using the graphical interface in the AE tool Spectra Viewer; an example screen shot is shown in Figure 13. Numerical output of the fitting process (in the form of a FITS table) includes best-fit spectral parameters with their confidence intervals and model fluxes integrated over various spectral bands.

Zoom In Zoom Out Reset image size

For visualization of source variability, light curves for all observations of the source are depicted on a single plot, which is shown in calibrated units (photon ks⁻¹ cm⁻²) rather than observed units (count ks⁻¹) to avoid spurious apparent variability arising from differences in observatory response or PSF fraction between observations. Two light curves are constructed. One is a histogram with variable-width time bins, constructed by the same algorithm used to group spectra (Section 7.2). The second is adaptively smoothed using a box-kernel whose size is adjusted to encompass a specified number of events. The adaptive smoothing process also produces a time series reporting the median energy of the observed events falling within each kernel, providing a visualization of spectral variation over time. For each observation of a source, AE also produces a "photon arrival diagram"—a scatter plot of event arrival time and energy. Figure 14 shows an example multi-observation light curve and a photon arrival diagram for a flaring source.

Zoom In Zoom Out Reset image size

For quantification of variability, the null hypothesis of a uniform source flux is tested using a one-sample Kolmogorov–Smirnov (KS) statistic.⁶⁷ For multi-observation cases, variability is tested both within and between observations (accounting for variations in effective area among the observations).

The current timing analysis has some deficiencies. No tests for autocorrelation, periodicity, or other temporal behaviors are made. The KS statistic is insensitive to variations near the beginning or end of the observation. Variability in the background is not considered in the uniform flux model, and background is not accounted for in the light curves or median energy time series. Although the ACIS background is often quite stable within an observation, multiple observations of a source can easily suffer very different background count rates within the source apertures because the aperture sizes can be quite different. Thus, for example, a spurious indication of variability can be produced for a weak source observed at two very different off-axis angles. We plan to eliminate this problem by adding to the uniform light curve model independent background levels for each observation.

It is often desirable to "stack" a sample of similar sources to construct a summed spectrum with higher S/N than available from individual objects. This can be performed either for sources detected in the ACIS data or at locations of objects known only from independent observations. Stacking many similar objects can thus give extremely long effective exposures and sensitivities, however, stacking results can be dominated by the few brightest objects in the group. This technique has been widely used in Chandra deep extragalactic field studies (e.g., Hornschemeier et al. 2001; Brandt et al. 2001; Lehmer et al. 2007).

In AE, stacking is achieved by relabeling extraction subdirectories of interesting sources to resemble extractions of distinct observations of a single source. Then the extractions are merged (Section 6) and source properties are estimated normally (Section 7).

The properties of individual sources obtained by AE can be collated into a summary FITS table with >100 columns and a row for each source. The scientist can then study relationships among quantities using, for example, the fv or chips FITS viewers, IDL two- and three-dimensional plotting, and R⁶⁸ multivariate statistical analysis. AE also prepares standard quantities for publication as LaTeX tables, illustrated by Tables 1 and 2. This feature is particularly valuable for fields with hundreds or thousands of sources.

As described in Section 5.4 and Figure 2, AE provides two methods for constructing mask regions that remove point source photons that would contaminate an analysis of diffuse emission. The remaining event distribution, consisting of instrumental background and possible astrophysical diffuse emission, is often best studied in smoothed images. Although the CIAO task csmooth (Ebeling et al. 2006) is effective and popular for smoothing unmasked data that contain both point sources and diffuse emission, it is poorly suited to working with images that contain unobserved regions, either masked point sources or regions beyond the edges of the detector, because the tool does not allow a field mask to be supplied. The CIAO thread for diffuse emission⁷⁰ with embedded point sources replaces pixel values in source regions of an image with values interpolated from surrounding background regions using random Poisson numbers or local polynomial regression.

We prefer to apply an adaptive kernel smoothing tool in TARA that handles excised pixels and field edges, which is similar to the asmooth tool in the XMM-SAS⁷¹ package and the new dmimgadapt⁷² tool in CIAO 4.2. Three smoothing kernels are available: a top hat (or Heaviside function), a symmetric two-dimensional Gaussian, and a two-dimensional Epanechnikov kernel (Silverman 1992). The observer supplies a threshold S/N value; at each pixel position, the algorithm finds the smallest kernel that achieves this S/N target and estimates an incident photon flux (photon s⁻¹) by dividing the number of counts under the kernel by the integrated exposure map under the kernel. Maps of photon flux, flux error, and kernel size are produced. The observer can choose to retain or smooth over the holes in the data introduced when point sources are masked. The method is described in detail, with sample images, in Townsley et al. (2003, Appendix C).

As the AE package does not attempt to delineate the morphology of diffuse structures, the observer must define each diffuse source of interest in the form of a DS9 region file that is accepted by CIAO. We have found the WVT Binning algorithm⁷³ (Diehl & Statler 2006), which is a generalization of the Voronoi binning algorithm described by Cappellari & Copin (2003), to be a very effective method for tessellating a field with compact regions of similar S/N.⁷⁴ Figure 16 shows an example of its use.

The procedure for extracting each diffuse source with AE is analogous to the point source procedure. For each observation, events within the source region are extracted, calibration products are constructed, and background spectra are extracted. All the extractions are merged, photometric quantities are estimated, and spectra are fit via automated scripts. However, more complex procedures are required to calibrate diffuse extractions and to account for background, as described in the next two sections.

The ARF data product⁷⁵ constructed for point sources by the tool mkarf can be thought of more abstractly as an observatory response function, ARF(E, p), computed at a position on the sky, p, for a monochromatic energy E, with units of cm² count photon⁻¹.⁷⁶ For Chandra data, this function represents both variations in the observatory response across the focal plane and the exposure time variations across the sky caused by dithering over bad detector pixels and detector edges.⁷⁷ Essentially, ARF(E, p) describes the "depth" of the observation, since effective area and observing time are degenerate terms in a flux calculation. For example, ARF(E, p) trends downward with off-axis angle, due to vignetting of the mirrors. ARF(E, p) is reduced at any sky position that dithered across bad pixels or dithered off of ACIS. Conceptually, ARF(E, p) is zero in regions on the sky that have been masked to remove point sources (Section 5.4).

The ARF(E, p) function for any particular observation will vary within a diffuse extraction region, and in a multi-observation study may be zero over much of that region. Figure 15 shows an example diffuse region that is fully contained in one ACIS-I observation (left panel) but barely intersects another ACIS-I observation (right panel). Just as a diffuse extraction can be viewed as an integration of the observed counts over the extraction region, the proper calibration of that extraction can be viewed as an integration of ARF(E, p) over the region.

$$\begin{equation} {\rm ARF}_R(E) = \int _R {\rm ARF}(E,\mathbf {p})\,d\mathbf {p}. \end{equation} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Note that this integration over a region on the sky gives ARF_R(E) units that include a geometric area term, such as cm² count photon⁻¹ arcsec². As ARF_R(E) is carried forward into photometry calculations and spectral fitting, all "flux" quantities derived should be interpreted as surface flux quantities, with arcsec⁻² appended to the units. Flux integrated over the diffuse region is then estimated by multiplying the surface flux by the region's total geometric area on the sky; the area lost to point source masks, bad pixels, and detector edges is already accounted for in ARF_R(E).

Recasting the diffuse ARF of each extraction to include the notion of "area on the sky" is particularly convenient when multiple extractions are to be combined (i.e., merged by AE, just as is done with point sources). Just as each extraction of a point source may have a different PSF fraction (Section 5.3), each extraction of a diffuse source may have a different "area on the sky" (due to masking and detector boundaries). In both situations, corrections to the calibration are most intuitively applied to each extraction prior to merging.

ARF_R(E) cannot be directly calculated with existing tools. The CIAO point source tool mkarf implements most of the behavior of the integrand ARF(E, p) described above, but has no mechanism to set ARF_R(E, p) to zero at positions where the event data have been masked. The CIAO tool mkwarf ("make weighted ARF") implements an effective area averaging calculation that is related to Equation (10), however, for each observation mkwarf averages over only the detector area that intersects the diffuse region, not over the entire diffuse region, which may include areas with zero response due to masking or detector boundaries. The difference between the diffuse ARF we have described here, ARF_R(E), and the data product returned by mkwarf can be summarized with the aid of Figure 15. For both extractions shown there, the mkwarf data products will have similar normalizations, because responses are averaged over the detector areas intersecting the extraction region. In contrast, ARF_R(E) will be much larger for the left-hand extraction than for the right-hand extraction, because responses are averaged over the entire region.

Although the full ARF_R(E) function is not available as a standard data product, the value of ARF_R(E) at a single energy E₀ is easily calculated by integrating over the extraction region a standard CIAO exposure map built for a monochromatic energy E₀ and masked in the same way as the event data. In the context of an AE extraction, the exposure map is "unnormalized" with units of s cm² count photon⁻¹, built by supplying the "normalize=no" option to the CIAO tool mkexpmap. This exposure map is related to the integrand ARF(E, p) as

$$\begin{equation} {\rm emap}_{E_0}(\mathbf {p}) = {\rm EXPOSURE} \times {\rm ARF}(E_0,\mathbf {p}), \end{equation} \tag{ 11 }$$

where EXPOSURE is the "exposure time" of the observation. Thus, the diffuse ARF we seek is easily calculated at one energy as

$$\begin{equation} {\rm ARF}_R(E_0) = \frac{\int _R {\rm emap}_{E_0}(\mathbf {p})\,d\mathbf {p}}{{\rm EXPOSURE}}. \end{equation} \tag{ 12 }$$

With the normalization of ARF_R(E) established at energy E₀, AE then relies on the mkwarf result to establish the shape of ARF_R(E) as a function of energy. AE's complete calculation of ARF_R(E) is thus

$$\begin{equation} {\rm ARF}_R(E) = \frac{{\rm ARF}_{mkwarf}(E)}{{\rm ARF}_{mkwarf}(E_0)} \times \frac{\int _R {\rm emap}_{E_0}(\mathbf {p})\,d\mathbf {p}}{{\rm EXPOSURE}}, \end{equation} \tag{ 13 }$$

where ARF_mkwarf is constructed by mkwarf.

To summarize this section, we calibrate each diffuse extraction by constructing ARF_R(E), which has units of cm² count photon⁻¹ arcsec² and accounts for the area on the sky lost to point source masks, bad pixels, and detector edges. "Flux" quantities subsequently derived should be interpreted as surface flux quantities, with arcsec⁻² appended to the units. Flux integrated over the diffuse region is estimated by multiplying the surface flux by the region's total geometric area on the sky.

Diffuse sources are often heavily contaminated by background—both instrumental background and emission from foreground and background astrophysical sources that are not of interest. Several strategies to account for background in diffuse sources are described in the AE manual. We choose to subtract from each diffuse extraction a standard instrumental background spectrum obtained by scaling the so-called ACIS stowed event data⁷⁸ provided in the Chandra Calibration Database (Hickox & Markevitch 2006). After merging the multiple observations of a diffuse region (with the same methods that are used for point sources in Section 6, with no optimization options enabled), standard "source" and "background" spectra with calibration data products are available. Astrophysical background contaminating a diffuse region can be handled in either of two ways. First, if a nearby "sky" region thought to be mostly free from the emission under study is available, then its observations and "stowed backgrounds" are extracted, and two net spectra from the diffuse region and the sky region are simultaneously fit using a shared model for the astrophysical background in both regions plus a source model for the emission of interest in the diffuse region. Alternatively, if no suitable "sky" region is available, then the astrophysical background must be directly modeled (Snowden et al. 2008). If appropriate fitting scripts are constructed, AE can automate this spectral fitting in the same way that point source spectral fitting is automated.

When more than a few diffuse regions are defined, many observers have found that spectral fitting results can be best understood by creating maps that show various results, e.g., model parameters and fluxes. Standard grayscale or false-color maps effectively depict to the human eye spatial variations in a parameter. We are currently experimenting with a complementary map-coloring technique that seeks to depict both a parameter's value and its uncertainty, using hue to encode the parameter value and brightness to encode its uncertainty, as shown in Figure 16.

Zoom In Zoom Out Reset image size