Department of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

Abstract
We describe an algorithm to optimally extract individual spectra of blended sources from a long slit spectrum. A semi-analytic model for the spatial profile is used: a Voigt profile for the undersampled core with a numerical correction applied to the wings. This is derived from an isolated source which must also be placed on the slit. $\chi ^{ 2}$ fitting is used to separate the blended components with maximum weight given to the best data. We demonstrate the successful application of the algorithm to data on two X-ray binaries in crowded fields: XTE J2012+381 and V404 Cyg.

1 Introduction

A common situation in observational astronomy is to seek accurate photometric or spectroscopic information about sources within crowded fields, i.e. where the nearest neighbours are blended with the sources of interest. This is a particularly common situation for X-ray binaries, most of which are faint Galactic plane objects. A solution to the photometric problem has long been available in various implementations of Stetson's DAOPHOT program (Stetson 1987). For crowded field spectroscopy, however, no widely available solution exists. A possible approach to the problem is described here and has much in common with DAOPHOT. This spectroscopic solution also seeks to be optimal, in the sense that it should yield the highest signal-to-noise spectra available from the data. It thus derives from optimal extraction methods developed for unblended spectra (Horne 1986; Marsh 1989). We developed the method initially in the context of an observation of the X-ray binary XTE J2012+381 (Hynes et al. 1999). It has also provided useful for another X-ray binary, V404 Cyg (Hynes et al. 2002). We illustrate the application of the method to both these sources in Sect. 6.

Since developing the method, another approach has been presented by Buie & Grundy (2000). This is also an optimal extraction, but differs in that no template source is used; instead the spatial profile is reconstructed iteratively from the blended sources. These authors also use a purely numerical profile, rather than the semi-analytic model we adopt. Their approach was developed for HST/NICMOS observations of a Pluto and Charon; hence these requirements were necessary. For ground-based observations of point sources in crowded fields, using a template source and a semi-analytic model profile should allow our method to be applied to poorer quality data. We should emphasise that the assumptions made here do only apply to point sources; this method cannot be used for sources with extended or resolved structure.

2 Preliminary image processing

Before applying the method, the images should have been debiased and flat-fielded using standard techniques. It is assumed hereafter that the image is oriented with the dispersion direction roughly horizontal and the slit is aligned with columns.

Sky subtraction should also have been performed and to ensure correct weighting across the profile, this should be done in a way which yields a image of the fitted sky background. Horne (1986) discusses these issues and we essentially follow his methods. Our implementation fits a low-order polynomial to each column with Horne's iterative cosmic ray rejection. Regions containing sources are masked out. We then retain both the subtracted image and an image containing the fitted sky.

It is also necessary to do wavelength calibration separately using normal procedures. The method employed for the spectra presented in Sect. 6 was to extract and wavelength calibrate a spectrum of the template source using normal IRAF methods, and then assume the same wavelength-pixel correspondence for all the objects. This is acceptable for these data, since the spectra considered do have the image of the slit precisely aligned with the CCD columns, but in general a more rigorous procedure may be needed.

3 The spatial profile

3.1 The profile model

The first step in deblending spatial profiles is to define how the profile of a single point source should look. Several methods have been proposed for this for a single source. Horne (1986) normalises the observed count rates for each column, then fits a low order polynomial to each row of the data. This has the advantage of making no a priori assumptions about the spatial profile, but is inappropriate for spectrum deblending because his spatial profile is only defined in pixel steps; it cannot easily be transferred to another source which samples the spatial profile differently. In principle, one could interpolate between pixels, but this is dangerous when the spatial profile is undersampled, as is often the case (even where the seeing is poor enough that the spatial profile would become oversampled, it is common to bin to reduce readout noise). Horne's method is also only appropriate for spectra with small distortion, i.e. nearly parallel to CCD rows. Marsh (1989) describes an extension of this empirical approach which works even for very distorted spectra. His method also produces well sampled spatial profiles which can be transferred to another source on the slit, taking advantage of the fact that strong distortion of the spectrum will lead to different columns sampling the spatial profile differently. This could be used as the basis for a deblending algorithm. In some cases, however, the distortion is very slight so this method will not work either; even by combining columns, the spatial profile is not well sampled. This was the case for the two datasets discussed in Sect. 6, and so led us to develop a different method.

The alternative is simpler to visualise, since it involves fitting a analytical profile in the spatial direction rather than a polynomial in the dispersion direction. Analytic spatial profiles have been applied by Urry (1988) to IUE data using a Gaussian profile. For ground based data a Gaussian model will often give a good fit to the seeing dominated cores of spatial profiles, but will usually underpredict the extended wings present due to instrumental imperfections. An effective refinement to this is to use a Voigt profile, a convolution of a Gaussian profile with a dispersion profile (e.g. Gray 1992):

$\begin{displaymath}P_i(u_i,a) \propto \int^{\infty}_{-\infty} \frac{\exp(-u_1^2)}{(u_i-u_1)^2+a^2}{\rm d}u_1 \end{displaymath}$

(1)

3.2 Fitting the profile

The first step in the fitting process is to perform an unconstrained fit to each column of the template data, with the centre, Gaussian width, Voigt damping parameter and profile scaling as free parameters. We have used a downhill simplex (amoeba) algorithm to obtain the best-fit parameters (e.g. Press et al. 1992). No attempt is made to reject bad pixels or cosmic rays at this stage, as these will show up as anomalous parameter values. Of these parameter values, it is the profile centre, width and damping that define a normalised profile. All of these parameters can be expected to vary smoothly in wavelength, so after obtaining fitted values, a low order polynomial in wavelength is fitted to each parameter, rejecting anomalous values that may have been affected by cosmic rays. If the wavelength range is small then a zero-order fit (i.e. a constant) may suffice. If the data quality is not sufficient to constrain the fits to a single row of the image then several rows can be binned, since the wavelength dependence is likely to be slow, and the derived polynomial resampled back to the original resolution. These polynomial fits, together with the normalisation such that $\sum_i P_i = 1$ define the spatial profile as a function of wavelength.

If the difference in brightness of the two sources is very large then it becomes crucial that the wing of the brighter source be very well fitted to avoid contamination of the fainter source. This can be achieved using a numerical correction to the (Voigt) model profile. The error involved in assuming a Voigt profile is usually only important in the extreme wings of the profile: the seeing-dominated core is typically well fitted by this model. We can therefore define a numerical correction factor (observed profile divided by model) at the sampled points and interpolate this to obtain a general correction function for any pixel sampling. Since this involves the extreme wings of the line, where signal-to-noise will be poor, such a correction typically cannot be defined meaningfully as a function of wavelength; instead only a single, averaged, wavelength independent correction can be determined. A symmetric profile is also assumed. These assumptions appear adequate for the spectra we have considered; hence the correction is simply a function of distance from the profile centre. Given higher quality data, an asymmetric, wavelength dependent function could readily be used, and this generalisation is straightforward.

Assuming that there are no changes in focus along the slit (or at least the restricted part of it which is used) and that the image scale (and hence the separation between stellar images) is independent of wavelength, the only additional information required to define the profiles of the blended sources is the separation of each from the template source. For spectra with minimal distortion, this is easy to obtain by performing an initial fit to the average of several rows of data, leaving the centres of the two blended profiles as free parameters. Distorted spectra would require that the distortion (known from the fits to the template profile) be removed by resampling the data for this initial fit only.

4 Spectral extraction

Assuming that the fit to the spatial profile for the template source is satisfactory, the centre position and normalised profile as a function of wavelength are known, and the spectral extraction itself can be performed. This is simply a case of finding the scaling of each component, although complicated by the need to reject cosmic ray events. Horne (1986) describes his method in terms of an optimally weighted sum of counts over an aperture containing all of the stellar light. He notes the equivalence between this weighted sum method and fitting a known profile to data of known variance, and it can be shown that his optimally weighted sum yields the solution of minimum $\chi ^{ 2}$ . Using this equivalence we can generalise the method to blended spectra, by minimising the $\chi ^{ 2}$ of the fit with respect to the two profile scalings. This is straightforward and for the simple problem of scaling two profiles, the solution of minimum $\chi ^{ 2}$ can be calculated analytically; see Appendix A. This is done iteratively, to allow rejection of cosmic rays. Given that there are two components to be fitted, however, the combined profile is unknown, and so more care is needed in cosmic ray rejection than for a single source.

5 Cosmic ray rejection

One of the great strengths of optimal extraction algorithms is that the spatial profile is known. Cosmic rays or bad pixels, even when on top of the spectrum, will distort the spatial profile, and hence can be recognised and rejected. This means that instead of removing cosmic rays from the extracted one-dimensional spectra, and rejecting a whole pixel step in wavelength, they can be removed from the two-dimensional spectra, retaining the information in the uncontaminated pixels. While this process is easy to do by eye, it proves difficult to train a computer to recognise cosmic rays automatically, especially when multiple profiles are being fit. Given a first approximation to the profile, bad values can be rejected by an iterative sigma-clipping algorithm (Horne 1986). The strategies we have used to obtain the crucial first approximation are described here.

Beginning with the simplest case of recognising contamination of a single profile, we rely on the fact that the scaling of the profile can be estimated from a single pixel value; hence each pixel gives an independent estimate of the integrated flux. A contaminated pixel will then give an anomalous estimate. If this estimate is severely wrong, as is often the case with cosmic rays, then a straight average may be very much in error, leading to subsequent rejection of good data. A better method is to take the median of the single pixel estimates as the first estimate for the integrated flux. This then gives a first approximation to the fitted profile which forms the basis for subsequent refinement with iterative sigma-clipping.

Where two independently fitted profiles are involved more care is needed. The only case considered here is that one source is brighter than the other at all wavelengths in the spectrum (true for all data the method has so far been applied to.) For the brightest source, a profile region is selected from the midpoint between the two sources to an arbitrary cutoff on the other side. Within this region, the profile of the brighter source is expected to dominate, so the median ratio method described above for single sources can be applied. This gives an approximate flux scaling for the brighter stellar profile, which can then be subtracted. The median ratio method is then applied to the remaining profile of the fainter source to estimate a flux scaling for this. With an approximate scaling for both profiles, iterative sigma-clipping can be used to refine the list of rejected pixels.

Clearly if three or more blended sources are involved then the situation is more complex, and there are more possible permutations of relative brightness to consider. This makes it harder to treat generally and specific cases will need to be considered individually.

6 Examples

6.1 XTE J2012+381

$\begin{figure} \par\includegraphics[angle=90,width=8.8cm,clip]{MS1993f1.eps}\end{figure}$

Figure 1: Fit to spatial profiles in the red part of the WHT spectrum of XTE J2012+381. Data are shown by points, the model profile is dashed and the dotted lines indicate the De-blended components. Both data and fits have been summed over 20 pixels in the dispersion direction to reduce noise and illustrate the quality of fit achieved. On the left is the fit to the spatial profile of the template star. On the right is the two-profile fit to the blended stars.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{MS1993f2.eps}\end{figure}$	Figure 2: Wavelength dependence of the Voigt profile fitting parameters, the Gaussian width and the Voigt damping parameter. Points mark fits to individual pixels in the dispersion direction, with no binning. The solid lines are low-order polynomial fits to the parameters.
Open with DEXTER

XTE J2012+381 presents an especially difficult challenge for spectroscopy. The X-ray source was originally identified with what appeared to be a normal F star. It was only upon closer examination of images of the field taken in good seeing that it became clear a second fainter star was present (Hynes et al. 1999). This was heavily blended with the F star, but lay closer to the precise radio counterpart and now appears the most likely optical counterpart to the X-ray source. We obtained spectroscopy of these two stars with the William Herschel Telescope (WHT) on 1998 July 20, with the slit aligned to pass through both, as well as through a third unblended star. These data provided the original motivation for development of this method. A fuller description of the data is provided by Hynes et al. (1999). We focus here on aspects relevant to the deblending algorithm.

$\begin{figure} \par\includegraphics[angle=90,width=8.8cm,clip]{MS1993f3.eps}\end{figure}$

Figure 3: WHT spectra of the F star and the faint red companion believed to be the optical counterpart of XTE J2012+381. Both have been binned $\times 4$ in wavelength for clarity. Atmospheric absorption features (marked $\oplus$ ) have been corrected for in both spectra. The only prominent line in either spectrum is H $\alpha$ : absorption in the brighter star and possible weak emission in the suggested counterpart to the X-ray source.

Open with DEXTER

The severity of the blending can be clearly seen in Fig. 1. The suggested counterpart to the X-ray source is the fainter of the two blended stars. Clearly to successfully extract its spectrum requires a good fit to the wing of the brighter star. Fortunately, the signal-to-noise ratio of the two-dimensional image was quite good, so it was possible to apply the full model described above, i.e. a Voigt profile with a numerical correction to the wings. The latter significantly improved the fit and with this no significant discrepancy was seen between the model profile and the data. The fits to the Voigt profile parameters are shown in Fig. 2. Both the width and damping parameter are relatively well constrained as a function of wavelength.

The extracted spectra are shown in Fig. 3. There are a shortage of spectral features in both spectra, but the extraction does clearly distinguish the much redder spectrum of the proposed X-ray source and the contrast between H $\alpha$ absorption in the F star and possible weak emission in the X-ray source is dramatic. The reality of this emission is discussed by Hynes et al. (1999); suffice here to note that no obvious artifacts are present, but that these data are perhaps not the best test of the method. Fortunately our second target, V404 Cyg, does provide a better test, albeit in a less heavily blended situation.

6.2 V404 Cyg

Another X-ray binary for which this deblending method is appropriate is V404 Cyg, as this has a fainter companion star 1.5arcsec away. There is no physical association between the two. Observations were obtained of V404 Cyg in quiescence using the WHT on 1999 July 6-8. These are described in more detail in Hynes et al. (2002). To summarise, the goal was to perform spectroscopy near H $\alpha$ with as high a time resolution as possible ( $\sim$ 240s) to examine spectral variability; hence we used a low spectral resolution to minimise readout noise and a wide slit to minimise slit losses. Figure 4 shows two examples of averaged (in wavelength) spatial cuts with respect to the centre of the template star. The images with the best and worst seeing from the night of July 7/8 were chosen. The fainter star is clearly blended with V404 Cyg even in the best image (seeing $\sim$ 0.8arcsec).

$\begin{figure} \par\includegraphics[angle=90,width=8.8cm,clip]{MS1993f4.eps}\end{figure}$	Figure 4: Best (solid) and worst (dashed) average spatial profiles of V404 Cyg from 1999 July 7/8. These are the observed profiles not models. Columns containing cosmic rays have been excluded from the average. Peaks are from left to right, the blended companion star, V404 Cyg and the profile template star. The profiles have been scaled such that the peak flux from the template star is unity.
Open with DEXTER

The optimal deblending algorithm was applied with a simplification. Because the goal of this project was to obtain the highest time-resolution practical, the signal-to-noise ratio of individual frames was not high. Because of this, it was not possible to perform a wavelength dependent fit and so a single average profile (and a profile correction) were calculated. Given the limited wavelength coverage of these data this is an acceptable approximation.

$\begin{figure} \par\includegraphics[angle=90,width=8.8cm,clip]{MS1993f5.eps}\end{figure}$

Figure 5: Deblended spectra of V404 Cyg and the blended star. All the spectra from 1999 July 7/8 have been averaged. The deblending algorithm has appears to have separated the spectra cleanly with the H $\alpha$ emission from V404 Cyg clearly distinguished from absorption in the blended star. The absorption feature in the blended star is not exactly aligned with the dip in the line profile of V404 Cyg.

Open with DEXTER

Figure 5 shows the average extracted spectra of the two blended stars from July 7/8. The spectra appear to have been separated cleanly and there is no obvious crosstalk around H $\alpha$ between the double-peaked disc emission in V404 Cyg and the narrower (stellar) absorption line in the blended star. There is also a lot of structure in the spectrum of V404 Cyg (e.g. 6600-6900Å) which is not present in the spectrum of the blended star which appears largely featureless. As V404 Cyg is brighter, the signal-to-noise ratio of its spectrum should be higher and hence these features must be real; they likely arise from the K0IV companion star which dominates in quiescence. The blended star, in contrast, appears to be of earlier spectral type, probably F (Hynes et al. 2002).

Figure 6 shows continuum light curves from both stars for the night of July 7/8. Spectra were scaled relative to the template star to approximately remove variable slit loss and transparency effects. Clearly the light curve of the blended star shows little or no contamination from the large flare occurring at the beginning of the night in V404 Cyg, indicating again that there is no detectable cross-talk between the two spectra.

7 Conclusions

We have demonstrated an algorithm for extracting separate spectra of sources whose spatial profiles are blended in a long slit spectrum given an isolated template source with which to define the spatial profile. We implement this method using $\chi ^{ 2}$ fitting with rejection of cosmic rays, resulting in a generalisation of the optimal extraction method to the multi-object case. Finally we have demonstrated the application of the method to two X-ray binaries in crowded fields that require separation from nearby contaminating sources.

Appendix A: Analytic $\chi ^{ 2}$ fitting of multiple profiles

For the simple case considered here, the spatial profile of blended sources is known a priori from a template source and it is assumed that the offsets of these blended sources relative to the template are known. The extraction then reduces to finding the flux scaling of each profile. The optimal solution is equivalent to minimising $\chi ^{ 2}$ , with care to reject any cosmic rays on the profile.

Consider N model profiles P_ij ( $j=1 \ldots N$ ) to be fitted to an observed profile D_i (variance V_i), yielding Nscalings F_j. The badness of fit statistic is then:

$\begin{displaymath}A = \sum_i p_i^2 / V_i \hspace{5mm} B = \sum_i q_i^2 / V_i \hspace{5mm} C = \sum_i p_i q_i / V_i \end{displaymath}$

$\begin{displaymath}f = \frac{D/C - E/B}{A/C - C/B} \hspace{5mm} g = \frac{E/C - D/A}{B/C - C/A}\cdot \end{displaymath}$

(A.4)

		$\displaystyle \sigma_f^2 = \sum_i \left(\frac{\partial f}{\partial D_i}\right)^2V_i = \frac{B}{AB-C^2}$
		$\displaystyle \sigma_g^2 = \sum_i \left(\frac{\partial g}{\partial D_i}\right)^2V_i = \frac{A}{AB-C^2}\cdot$	(A.5)

An optimal extraction of spatially blended spectra

1 Introduction

2 Preliminary image processing

3 The spatial profile

3.1 The profile model

3.2 Fitting the profile

4 Spectral extraction

5 Cosmic ray rejection

6 Examples

6.1 XTE J2012+381

6.2 V404 Cyg

7 Conclusions

Appendix A: Analytic $\chi ^{ 2}$ fitting of multiple profiles

References

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{MS1993f6.eps}\end{figure}$	Figure 6: Continuum light curves of V404 Cyg and its blended companion.
Open with DEXTER

An optimal extraction of spatially blended spectra

1 Introduction

2 Preliminary image processing

3 The spatial profile

3.1 The profile model

3.2 Fitting the profile

4 Spectral extraction

5 Cosmic ray rejection

6 Examples

6.1 XTE J2012+381

6.2 V404 Cyg

7 Conclusions

Appendix A: Analytic fitting of multiple profiles

References

Appendix A: Analytic $\chi ^{ 2}$ fitting of multiple profiles