Optimal Extraction of Fibre Optic Spectroscopy

Keywords: instrumentation: spectrograph — methods: data analysis

Fibre-optic multi-object spectroscopy is a powerful tool with which to segment the focal plane of a wide-field spectroscopic telescope. The technique, pioneered in the 1970s (see Hill (1988) for a historical review), has been repeatedly put to use for major undertakings in astronomy. The high multiplex advantage that can be achieved, typically over wide fields-of-view, using fibre optic feeds easily outweighs the expense of implementation and inherent limitations of a fibre feed system for many classes of astronomical observation. Notable recent examples of include the 2dF Galaxy Redshift Survey (Colless et al. 2001), 2dF Quasar Survey (Boyle et al. 2000) and the Sloan Digital Sky Survey (York et al. 2000).

The use of coherent fibre bundles to record two-dimensional spatially resolved information at the spectrograph slit is an obvious extension of the technology, and indeed a number of current-generation Integral Field Spectrograph (IFS) systems use fibre feeds to reformat the focal plane (e.g. AAT-SPIRAL, Gemini-GMOS, ESO-VIMOS & ARGUS).

However, increased multiplex (or field of view for IFS instruments) comes at a price. With CCD (and IR array) pixels still at a premium in any astronomical instrument, there is a need to tightly pack the fibres into the spectrograph slit, minimising the dead space between adjacent fibres in order to maximise the number of fibres available for a given instrument format. But tight fibre packing means poor sampling of individual fibre Point Spread Functions (PSF) and can result in significant overlap of the profiles of adjacent fibres leading to strong cross-contamination, or cross-talk, between fibres. Recovering the spectral information for spectra observed in tightly packed systems is therefore complex. Simple summation of the pixel values surrounding the peak of a fibres’ trace is doomed to failure due to contamination from adjacent fibres (particularly with high contrast observations as are common in IFS applications). Even if individual fibres are reasonably well resolved, the extraction will often suffer increased noise from contamination by the wings of adjacent fibre profiles, or poor weighting of pixel values introducing increased CCD read-noise in the limit of low background observations.

What is required is an extraction procedure which takes account of the interaction between adjacent fibres and provides a statistically optimal estimate of the true intensity for each fibre spectrum. We have recently implemented such an extraction mechanism for the AAOmega spectrograph system at the Anglo-Australia Telescope. The code has been developed in the idl programing language and implementation within the 2dfdr data reduction environment (constructed primarily in fortran). In this work we quantify the need for such an extraction algorithm and demonstrate that while there is limited gain in using the routine for the well resolved fibre profiles of the default AAOmega instrument Multi-Object Spectroscopy (MOS) mode, it is essential for high purity data taken with the tight fibre packing utilised for the recently implemented AAOmega mini-shuffling nod-and-shuffle observing mode and for the SPIRAL IFS feed to the AAOmega spectrograph. We describe the methodology used and demonstrate its application to observations taken with the AAOmega-SPIRAL system.

The paper is organised as follows. In Sections 2 & 3 we introduce the three alternative spectra extraction methodologies currently implemented with in the 2dfdr data reduction software and four modes of instrument operation for the AAOmega spectrograph system at the AAT. In Section 4 we demonstrate the need for an extraction algorithm which accounts for fibre-to-fibre cross-talk. The two non-trivial extraction algorithms from Section 2 are then described in detail in Section 5. The procedure for applying the the optimal extraction algorithm is given in Sections 6 & 7. Scattered light effects are addressed in Section 8 and example data from the AAOmega-SPIRAL system is presented Section 9.

Three spectral extraction methodologies are currently implemented within the 2dfdr data reduction environment used with data from the AAOmega spectrograph. In each case the starting point is a map of the fibre profile centroids on the CCD. Since this fibre map is a series of near parallel tracks it has historically been labeled the tramline map. The generation of the tramline map is conducted separately from the spectral extraction and is not discussed in detail here. Briefly, the fibre profile centroids are determined at intervals across the CCD (every ∼50 pixels along the dispersion axis) and a low order polynomial model is fitted to the centroids for each fibre, guided by an optical model for the expected spectrograph camera distortions.

The three spectral extraction methods currently available within 2dfdr are:

• Tramline Summation This most basic of extractions is obtained by simple summation of all pixel values associated with a given fibre. The spectra for the fibres are dispersed broadly along rows of CCD pixels and each column of the CCD is treated independently. The summation range for each fibre typically runs over the pixels bounded by the mid-points between the two adjacent fibre profiles. While quick to compute, the simple tramline extraction propagates the maximum CCD readout-noise into the final extraction since it gives equal weight to all pixels in the summation regardless of the flux level of the fibre profile in a given pixel. Tramline summation also suffers an aperture loss effect if the inter-fibre gap over which pixels are summed is not significantly larger than the fibre profile width. This aperture correction can of course be accounted for if the true fibre profile is known.

• Gaussian Weighted Summation An algorithm for performing a weighted Gaussian summation over a single isolated fibre profile (via a Least Squares Fit) is discussed in Section 5.1. This weighted summation minimises the contribution for CCD readout noise and does not suffer an aperture effect. This mode is the default for AAOmega-MOS spectroscopy and we show in Section 4 that for well separated spectra it delivers acceptable results.

• Multi-Fibre Deconvolution The ideal solution for the minimisation of cross-talk between fibres is an extraction algorithm which performs a multi-fibre deconvolution of the data given an underlying model assumption for the fibre profiles on the CCD.

The extent to which fibre-to-fibre cross talk is present under each of the extraction models described above is tested in Section 4. The algorithms adopted for the Gaussian and multi-fibre deconvolution solutions are then discussed in Sections 5.1 & 5.2. All three approaches have been implemented for the AAOmega spectrograph at the Anglo-Australia Telescope. The multiple operating modes of AAOmega are outlined in the next section.

Table 1 outlines four instrument modes for the AAOmega facility at the AAT (Saunders et al. 2004; Sharp et al. 2006). Data from the AAOmega system is primarily processed using the 2dfdr data reduction software and the three extraction algorithms presented above have been developed from the code base of the previous 2dF spectrographs.

Table 1.  Pertinent parameters for the four model instrument modes of the AAOmega system

The first instrument mode is the default multi-object spectroscopy (MOS) mode of operation. In this mode 392 science fibres are deployed on astronomical targets across the π deg² field-of-view of the 2dF prime focus corrector. The fibres feed the dual beam AAOmega spectrograph, each arm of which is equipped with a 2k × 4k E2V CCD. This results in a fibre-to-fibre separation, the fibre pitch, of ∼10 pixels. The cameras deliver largely Gaussian profiles of FWHM ∼3.4 CCD pixels.

The second mode under consideration is the newly implemented mini-shuffling mode. The high quality AAOmega PSF makes it possible to undertake fibre nod-and-shuffle observations for high quality sky subtraction (Glazebrook & Bland-Hawthorn 2001) by interleaving multiple on-sky exposures on a single CCD frame using the inter-fibre gaps to provide the required CCD storage areas during charge shuffling. This mode effectively doubles the number of AAOmega fibres on a CCD by halving the fibre pitch. This introduces significant fibre-to-fibre crosstalk as will be demonstrated in Section 4. A subsequent paper will discuss the effectiveness of mini-shuffling in detail.

The optimal extraction algorithm presented in Section 5.2 was developed for the AAOmega-SPIRAL Integral Field Spectrograph system used with AAOmega at the AAT. It is based on an implementation developed for the CIRPASS IR spectrograph (Parry et al. 2004). Both CIRPASS and SPIRAL use a tight fibre packing, and with good reason. CIRPASS initially utilised a Rockwell Hawaii-1K IR array and hence required a tight fibre packing to maintain the wide field of view in IFS mode (Krajnović, Sharp & Thatte 2006). The SPIRAL system is designed to allow nod-and-shuffle (Glazebrook & Bland-Hawthorn 2001) observations with the full 22.4 × 11.2 arcsec FoV of the SPIRAL system, while using the same detector real-estate as the AAOmega multi-object fibre feed from the 2dF fibre positioner (Sharp et al. 2006). Hence for the SPIRAL system CCD pixels are at a premium as each fibre essentially requires twice as many pixels on the CCD as dictated by its intrinsic foot print.

The final mode of operation considered in Table 1 is that proposed for the HERMES high resolution (R ∼30 000) multi-object stellar spectrograph under development for use with the 2dF positioner at the AAT. HERMES will employ 4k × 4k CCDs with 15 μm pixels. However at f1.5, the HERMES camera design is slower than the f1.3 cameras of AAOmega, resulting in a broader PSF (5 pixels) than that of the currently implemented AAOmega modes. This results in an increased profile overlap compared to the basic AAOmega mode.

Before presenting our new extraction algorithm which reduces fibre-to-fibre spectral cross-contamination (Section 5.2) we will motivate its development by consider the level of this cross-talk in models of the four observing systems of Table 1. The pertinent parameters for the models of each observing mode are:

• Fibre Pitch The spacing between fibres on the CCD.

• Fibre Profile e.g. the projected FWHM of a fibre profile on the CCD.

• Relative Intensity The intensity ratio between spectra to be extracted.

Figure 1 presents a visual representation of the instrument modes of Table 1. A pair of isolated fibres are considered. This represents the limiting case of cross-talk between two adjacent fibres only. Two model fibre profiles are shown with an input Δmag = 3, a flux ratio of ratio ∼15.8. The fibre centre, and ±3σ range are marked for each fibre (solid bars) as are the inter-fibre ranges¹ associated with each fibre (dashed bars). This later range is used for the Tramline summation extraction, while ±3σ is used for the Gaussian summation unless it is larger than the inter-fibre range, in which case it is truncated to this range. A visual inspection of Figure 1 indicates that for all but the default AAOmega mode there may be large fibre-to-fibre cross-contamination, due to the significant overlap between fibre profiles, if a simple pixel intensity summation extraction is used. For the later instrument modes an extraction process is needed that accounts for the flux of both fibres simultaneously and attempts to correctly distribute the intensity information between the two fibre profiles in the overlap region. A quantitative analysis follows.

Figure 1  Model fibre profiles are shown for a pair of fibres with each of the instrument modes given in Table 1. The logarithmic scaling heightens the visibility of the region of overlapping profiles. The profile intensity ratio is ∼15.8 giving Δmag = 3. The fibre centres, and ±3σ range are marked for each fibre (solid bars) as are the free inter-fibre ranges associated with each fibre (dashed bars). For all but the first instrument profile, the 3σ range extends beyond the inter-fibre range.

Examination of observational data from the AAOmega system indicates that a Gaussian fibre profile is an acceptable first approximation. We will use the Full Width at Half Maximum, FWHM,² as the measure of profile width in what follows. We express relative spectral intensities in terms of the magnitude difference of the integrated profile fluxes, Δmag, for two astronomical objects illuminating the fibres. Conventional wisdom for the AAOmega system has been to limit Δmag <3.

Figure 2 demonstrates that for the default AAOmega-MOS mode of operation (the two top figures in each column) cross-contamination of the fainter spectrum is limited to <1.6% of its integrated profile flux for Δmag<3. Improved accuracy requires a reduced Δmag. However the three remaining instrument modes detailed in Section 3 all suffer significant cross-contamination. In these cases not only is cross-contamination of the weaker spectrum increased significantly, but double counting of flux leads to >10% errors in flux normalisation when Δmag< 0.3. Clearly an improved extraction methodology is needed which can account for the overlapping fibre profiles in these cases. No figure is given for the Multi-profile deconvolution extraction since, under the assumption of this idealised initial analysis, it returns perfect results for all four instrument modes.

Figure 2  Using two model Gaussian fibre profiles and the instrument modes of Table 1 the ratio of extracted flux to input model profile flux is presented as a function of the relative intensity difference (expressed as a magnitude difference, Δmag) between the two input profiles. The left column shows results using the Tramline summation extraction, the right the Gaussian least squares fit. In all cases the second fibre is Δmag fainter than the brighter first fibre. An accurate extraction would follow the 1:1 locus. For all instrument modes the extracted flux for fibre 2, the fainter profile, is increasingly in error as Δmag increases due to cross-contamination from the brighter source fibre. No traces are shown for the multi-profile deconvolution extraction as the process is found to be 100% accurate (the correct flux is recovered for both fibres) in the noise free model limit for all four instrument modes.

Two parameters of the assumed fibre profiles dominate the extraction error, errors in the assumed fibre centroid and errors in the profile widths. Figure 3 demonstrates the effects of introducing random modifications to the fibre centroids and widths (in isolation). Fifty thousand realisations are computed with a random error introduced to each fibre profile drawn from a normal distribution with the indicated width. The ratio of the mean extracted flux to the input model value (and associated ±1σ scatter) is given for each error distribution. The flux ratio between the fibres was set to unity (Δmag = 0). The simulations show that while errors in the fibre centroid value may be tolerated at the 0.2 pixel level (6% of the FWHM here), an accurate knowledge (<0.1 pixel, ∼3% of FWHM) of the fibre profile is required to control errors under all extraction methodologies.

Figure 3  The ratio of the mean measured flux to the input model intensity (and associated RMS scatter) is shown for the three extraction methods of Section 2 after 5000 independent realisations of the error model. For each realisation the model fibre profile centroid and FWHM parameters are drawn from normal distributions with the indicated widths. The summation aperture correction factor has not been applied to the Tramline extraction, leading to its offset below the other two extraction methods.

A final parameter to consider when testing the accuracy of the extraction of fibre profiles is the signal-to-noise ratio (SN) of the data. To investigate this, the two fibre profile model is retained, but a constant readout noise per pixel and a shot-noise component is added to each realisation of the model data based on the intrinsic strengths of the input model profiles. Profiles are generated with a range of normalisations while the input magnitude ratio (Δmag) is held fixed for each realisation. Figure 4 shows the results obtained for the three extraction methodologies, for a range of SN ratios, assuming the default AAOmega-MOS mode as described in Table 1.

Figure 4  Variations in the extracted flux are considered as a function of the signal-to-noise ratio intrinsic to the model date realisation. Signal-to-noise is quoted for the integrated profile flux. Two sets of simulations are shown, both for the AAOmega-MOS instrument mode. Fifty thousand independent trials of the noise model are run with the fibre profile flux ratio set to Δmag = 0 and 3 (top and bottom rows respectively). Fibre 2 is the fainter fibre in the Δmag = 3 case. The left column presents the ratio of the extracted flux to that intrinsic to the model. The right column compares the intrinsic SN to that recovered for the extracted profiles. In the case of Δmag = 0 the results for the two fibres are indistinguishable. The optimal extraction is seen to return a smaller scatter in the recovered flux ratio, the scatter remaining largely constant as SN decreases. For the Δmag = 3 case both the Tramline and Gaussian extraction return an increased scatter in the ratio for the fainter fibre (fibre 2) of the pair. The RMS error bars for fibre two in the lower right plot are similar for all three extraction algorithms and so for clarity only those associated with the OptEx algorithm are given.

The figure shows the intrinsic signal-to-noise ratio used when generating the model profiles and the signal-to-noise ratio of the extracted profile flux estimate. This recovered signal-to-noise is defined as the ratio of the extracted count rate to the noise level as estimated from the model data by the extraction process. When reviewing the results of the simulations, one should note that the SN ratio recovered for each spectrum is a statistical measure derived from the data, and does not directly address the issue of the double counting of flux due to fibre-to-fibre crosstalk.

Figure 4 indicates that for the default AAOmega-MOS mode the choice of extraction methodology has only limited effect on the final SN ratio. In the case of Δmag = 0 there is some indication that the multi-component decomposition (optimal extraction) reduces the scatter in the ratio between recovered and intrinsic profile flux. This increase becomes more pronounced for lower SN ratio data. The situation is similar for Δmag = 3 with the error in the flux value derived for the fainter fibre suffering a greater inflation for the tramline and Gaussian extraction in comparison to the optimal extraction.

The process is repeated for each of the instrument configurations from Table 1. Figure 5 shows the results for Δmag = 0 and Figure 6 for Δmag = 3. It is immediately apparent that in the case of tightly packed fibres fibre-to-fibre cross-talk plays a significant role in modifying the extracted flux values for the tramline and Gaussian extraction methods.

Figure 5  Figure 4 is repeated for each instrument configuration presented in Table 1 and with Δmag = 0. For the compact fibre systems a significant systematic flux error is seen (for both fibres) in the extracted flux due to double counting in the overlap regions by the Tramline and Gaussian extraction methods.

Figure 6  Figure 5 is repeated with Δmag = 3. Even with a large flux ratio between the fibres the Optimal extraction method recovers the correct flux for both the bright (fibre 1) and faint (fibre 2) profiles. The apparent increase in the SN ratio of the faint fibre with the tramline and Gaussian extraction is due to the erroneous double counting of flux from the adjacent fibres, which would represent a degradation of the SN ratio actually achieved.

In each of the three spectral extraction methodologies introduced in Section 2 it is assumed that each fibre projects a spectrum onto the detector such that the light is dispersed along rows of pixels. At each spectral pixel the fibre profile is assumed to be well modelled by an analytic function aligned precisely along columns of the detector, orthogonal to the dispersion axis. These assumptions are approximately true at the pixel-to-pixel level, the macroscopic spectral curvature on the CCD being on a much larger scale. Such an approximation is made in almost all approaches to the extraction of spectra from 2D CCD data.

In the discussion of the two non-trivial spectral extraction algorithms presented below we use the terms defined in Table 2, in which the desired products of the fitting process are the integrated fibre intensities η_k for each fibre k, and the associated variance in η_k expressed as var_k. The process includes all pixels, i, along the spatial (slit) axis of the CCD which is assumed orthogonal to the dispersion axis. It is then repeated for each pixel of the dispersion axis.

Table 2.  Definition of terms used in the discussion of the extraction algorithms

An extraction method which minimises the read-noise contribution from pixels in the wings of the fibre profiles, while maintaining a full summation of all of the observed flux, is in principle rather trivial for the case of a single simple spectral profile (although careful consideration of the system is required for a truly optimal result, see Horne (1986) for a comprehensive discussion).

We consider each column of the dispersed data in isolation, and construct a model for the intensity at a given spatial pixel assuming an appropriate fibre profile φ_i, and integrated intensity η.

In the case of non-overlapping spectra, and given a recorded CCD pixel intensity D_i with associated error estimate σ_i,³ an estimate of the integrated profile η is given by the minimisation of

This yields

and, on minimisation,

with the summation running over a sufficiently large range of pixels, i, to cover the full fibre profile, without overlapping with adjacent profiles. A range of ±3σ (or half the inter-fibre separation for tightly packed fibres) is typically adopted. The resulting value of the integrated profile flux η is calculated for each fibre in isolation and for each element along the dispersion axis in tern.

The core of the fibre profile achieved by the AAOmega system is observed to be well matched to a Gaussian profile; and for speed and simplicity, while still maintaining sufficient accuracy, the current 2dfdr data reduction environment used with AAOmega data implements a simple Gaussian profile.

Furthermore, we assume that in the first instance the precise position of each fibre profile as a function of wavelength (i.e. detector column) is known from the previous trace of an appropriately illuminated flat field frame as discussed in Section 2.

We then define a normalised spatial profile for fibre k, as a function of the CCD pixel i along the spatial axis (the fibre slit axis, orthogonal to the dispersion axis), such that

Hence the contribution of fibre k to the observed intensity value at CCD pixel i is given by the product of the integrated profile intensity for the fibre, η_k, and the fibre profile at pixel i. For any given pixel of the detector a model value for the recorded count rate at pixel i is then simply the sum over all contributing fibre profiles,

Given D_i and σ_i, the recorded count rate and an estimate of the statistical error at each pixel (along a single column of the detector), we wish to evaluate η_k for each spectrum. This is achieved via the minimisation of the residual,

and on setting δR/δη_k = 0 and substituting for M_i we find,

Letting

and

we find

Equation (11) can be solved using any of the multitude of methods for solving coupled linear equations. Numerical stability considerations have guided our implementation towards a Singular Value Decomposition (SVD) approach. SVD allows terms of low significance to be removed from the matrices during fitting, providing control over numerical stability for high contrast data sets. One notes that, even for very tightly packed fibre data, the overlap between fibres separated by more than a few FWHM will be close to zero in most instances (with the exception of very high contrast ratio data, or for a model PSF with wide scattering wings). This makes c_kj at least band diagonal, and tri-diagonal in many cases, allowing trivial solutions to Equation (11) in these cases.

Propagation of error information is critical for all observational science, but is woefully missing in many spectral extraction routines. It can be achieved as follows. For each pixel i in the input image we calculate the total of all fibre profiles, φ_ki, that contribute to that pixel,

Then for each spectrum k, we calculate the fractional contribution of this spectrum’s profile to pixel i,

An estimate of the total variance in the extracted intensity η_k is then given by

The variance estimate associated with any given detector pixel, var(D_i) is determined directly from the observed pixel value and the predetermined read-noise and gain parameters for the CCD system. Given Equation (6) the value is given by the sum of the contributions from each fibre at that pixel plus the additional detector read-noise. Under the assumption of shot-noise in the photon arrival rate at a pixel, the contribution of a particular fibre to the variance estimate for that pixel is given by the produce of the total count rate for the fibre and the normalised profile intensity at the pixel.

When multiple fibre profiles contribute to a pixel, a model for the observed variance in that pixel is then

We then wish to minimise the residual between this model and the variance in the observed data, var(D_i).

By analogy with Section 5.2, and solving for δR/δvar_k = 0, we therefore arrive at an estimate for the variance, var_k, associated with each integrated fibre profile intensity η_k by solving

with c_kj′ = ∑_i(φ_kiφ_ji) and b_j′ = ∑_i [var(D_i) – n_rd²]φ_ji. This can be solved directly alongside Equation (11), with minimal additional computational overhead.

As demonstrated in Figure 3, variation in the fibre profile FWHM with wavelength is of particular concern.⁴ Poor correction for background scattered light will also prevent a satisfactory solution being achieved. Results can be significantly improved via an iterative solution to the problem of determining these free parameters.

We adopted an iterative solution with a number of free external fitting parameters. The model is derived from a high signal-to-noise flat field frame and then the parameters are locked at the resulting values when fitting for science frame intensities. The free parameters are:

• Background We include a low order polynomial fit to the broadband scattered light distribution across the CCD (pedestal plus gradient). In this model, scattered light is effectively treated as a pedestal correction to the observed data. For the AAOmega system, the scattering is well modelled as a DC pedestal offset across the CCD. This background level is a free parameter in the iterative solution of Equation (11), but see Section 7 for a further refinement.

• Fibre Profile The normalised fibre profile under the assumption of a Gaussian PSF has only a single free parameter, σ, the Gaussian profile width. A more accurate accounting for the SPIRAL fibre PSF is considered in Section 8.2.

• 2D Parameter Distribution A full 2D polynomial fit to the variation of these free parameters as a function of CCD position should also be performed. This can be achieved by fitting a subset of data columns at regular positions across the dispersion axis. A low order model would then be fitted for the full data set guided by instrument design considerations. This has not yet been implemented within the 2dfdr software.

With these external parameters in place we iteratively solve Equation (11), for all fibres simultaneously, for each element of the dispersion axis in turn. The idl prototyping code used mpfit⁵ (Markwardt 2009), adopting the χ² statistic as the figure of merit applied to the difference between the data and model.

The multi-fibre deconvolution algorithm described above is computationally intensive and to be of value the time taken to process observational data must be considered. The AAOmega system (Table 1) requires the solution for ∼400 simultaneous fibre profiles (∼800 for the mini-shuffle mode) across 4000 spatial pixels for each of the 2000 spectral columns of a 2k × 4k E2V CCD. The algorithm has been successfully implemented in the 2dfdr software environment (written primarily in fortran95) and the extraction is completed for a single frame in ∼1–2 minutes on a modern desktop PC.

In its current format, fifty percent of the calculation time for each instance of the fitting process is concerned with the calculation of fibre profiles and the population of the c_jk and b_j matrices from Equation (11). The current implementation does not distribute calculations across multiple PCUs and so this minor code extension should allow for a near linear speed increase, for these calculations, proportional to the a number of processors available on a modern multi-core computer.

For data taken with the SPIRAL-AAOmega system we undertake the following procedure:

1. A fibre-flat-field frame is taken in which all the fibres are illuminated with a uniform continuum source in order to trace the centroid of the fibre profiles across the CCD. Fibre-flat-field frames can also be used to correct for variations between fibres in the relative response function with wavelength. They do not provide an accurate correction for pixel-to-pixel CCD response variations due to the high degree of spatial structure across the fibre profiles. If such a correction is required, e.g. to correct for interference fringing in the CCD, a more uniform long-slit-flat field frame would be required.

2. All CCD frames are processed for overscan/bias correction and population of the variance array information.

3. An iterative fitting algorithm is applied to extract the integrated fibre profile intensities for the flat field. Experience with spiral data shows that re-adjustment of the fibre centroids measured from the fibre-flat-field is not required for AAOmega.

4. The iterative fitting algorithm is applied to the science data, holding the centroid and PSF values fixed (determined from the fibre-flat-field above) and fitting only for the pedestal offset as a free external parameter when iteratively solving for the integrated fibre intensities using Equation (11). Using the alternative background subtraction method outlined in section Section 7 below, the extraction of science data can be made non-iterative once the extraction parameters have been determined from the fibre-flat-field data.

If a simple pedestal/polynomial model is adopted for the background component then it is possible to fit the background without resorting to an iterative solution. For example, a simple pedestal can be thought of as an additional fibre profile with a uniform profile intensity at each pixel.

Higher order polynomial terms of the form i^N could also be introduced, up to some number of terms N_BG. One would then solve Equation (11) for η_k with 1 < k < N_fib for each fibre with normalised fibre profile φ_ki and also for N_fib + 1 < k < N_fib + N_BG governing the background model intensity terms.

This more elegant formalism allows the solution of Equation (11) without iteration in the case that all other external variables are predefined.

Two basic assumptions of the spectral extraction processes discussed in this work are:

1. The spatial (slit) axis on the CCD can be considered orthogonal to the dispersion (spectral) axis at each pixel along the dispersion axis.

2. Each spectral pixel is treated as if it where independent of the adjacent spectral pixels, i.e. no accounting is take of the 2D nature of the instrument PSF.

Essentially we are replacing the real-world spectrograph Point Spread Function (PSF) with a pixel-by-pixel Line Spread Function perpendicular to the dispersion axis.

Such assumptions are made primarily for computational expedience. For a compact and well-sampled PSF they are valid. This is fortunate since a full 2D modeling, which amounts to simultaneous deconvolution involving every pixel on the CCD (and a number of virtual pixels extending beyond the physical device), will likely remain beyond the computing power available for routine multi-fibre observations for some time due to the inherent scale of such a problem. As an example of this scale consider that the 512 fibres of the AAOmega-SPIRAL system each have ∼2000 spectral elements (∼830 independent resolution elements) covering ∼8 million CCD pixels.

The model breaks down if the fibre PSF possesses a measurable scattering wing which extends significantly beyond the core of the fibre profiles. If such a broad scattering profile is present then the assumption that data can be accurately modeled as a 1D line spread function rather than the 2D PSF breaks down since scattered light present at any given position along the dispersion axis will be largely dominated by the intensity of the local average spectral intensity of the surrounding spectrum and unrelated to the specific spectral intensity of the current spectral pixel.

In high signal-to-noise data (including the fibre-flat-field frames for science data that will ultimately be in the low signal-to-noise regime), even a scattering fibre profile with low overall percentage light level in the scattering component can have significant side effects on the ultimate data quality if the profile has a broad PSF. The scattered light from an increasing number of more remote fibre profiles will become more and more important.

Some scattered light is inevitable in all optical systems. Low level background structure which is largely devoid of the high spatial frequency signature of the fibres themselves can be removed relatively easily. One may fit directly to the background level observed in regions free from contamination by the fibre profiles. Alternatively a low order model can be included as an additional fibre profile as demonstrated in Section 7. This model profile would ideally be informed by an independent modeling of the background scattered light component.

For scattering more directly associated with the extended wings of the individual fibre PSFs, a more complex consideration is required. As previously mentioned, a full 2D deconvolution across the entire CCD array is likely beyond the capabilities of modern systems, at least at a rate comparable with the data rate from modern multi-fibre spectrographs.

A solution is to maintain the 1D orthogonal extraction assumption, but using fibre profiles with extended wings. This approach is flawed however. For an extreme example of why consider the case of a strong stellar absorption line within a strong continuum source. The scattered light at line centre would be dominated by the strong local continuum average, due to the underlying 2D nature of the scattered light PSF. The absorption line core profile would however be rather weak in comparison. This local variation in the spectral shape essentially induces a local variation in the assumed line spread function. One cannot merely fit for a core profile with a scattering wing modeled as a fixed percentage of the core profile flux.

A solution is to model each fibre profile as a composite of two normalised profiles, one for the core profile, φ_i^c, and a second broader scattering profile φ_i^s. Each component will have an independent intensity, η^c and η^s. The model for the CCD pixel intensity at pixel i is then the sum, over all fibre profiles k which contribute to the pixel, of both the core and scattering profiles:

Equation (11) is now solved for 2 × N_fib fibre profile intensities η_k.

There are inherent degeneracies in this composite profile description for each fibre. At low signal levels one can encounter significant problems with numerical stability when solving for values of η_k in real data. We have found the best results are achieved via a Single Value Decomposition (SVD) approach. The opportunity to suppress poorly behaved elements of the intermediate solution matrices within the SVD solution (akin to discarding from the analysis the coupled linear equations which contain limited information, as discussed by Press et al. (1996)) recovers numerical stability and more than compensates for the increased computational burden presented by SVD over simpler solvers. The SPIRAL IFU data extraction shown in Figure 7 were achieved using a composite profile extraction and an SVD solver.

Figure 7  A reconstructed 2D image, created by collapsing the spectral data cube over wavelength, is shown for a SPIRAL-IFU observation of a standard star. Four dead fibres have been interpolated in the data, and are marked with small spots. Left – The reconstructed image obtained using a simple tramline extraction model. The vertical structure is not a telescope diffraction spike, but rather an artifact of fibre-to-fibre crosstalk. Centre – The Gaussian summation extraction suppresses noise but not the fibre-to-fibre cross contamination. Right – The reconstructed image obtained using the optimal extraction model. The vertical extraction artifact structure is removed and a rounded PSF is recovered.

An interesting side-effect is encountered as one approaches the edges of the CCD. For a spectrograph with a flat response function in wavelength, scattering into an individual fibre profile from shorter and longer wavelengths will be much the same as that scattered out of the element in question (provided the spectral shape does not possess strong spectral breaks). However, as one approaches the limit of the free spectral range of the spectrograph this symmetry no longer holds. The exigencies of spectrograph design will typically place sharp spectral features such as those from order sorting filters, dichroic beam splitters or the grating blaze just beyond the operational range of the spectrograph. This means that the scattered light arriving in pixels towards the edge of the detector array may manifest a significant reduction in intensity since little light is scattered back onto the detector from just beyond the observed wavelength range. This has been found to be particularly true for the AAOmega system since the Volume Phase Holographic (VPH) gratings it employs have strong gradients in the blaze profile at the extreme ends of the intended observational ranges.

A striking comparison of fibre-to-fibre crosstalk using the tramline, Gaussian and optimal extraction techniques is given in Figure 7.

The left panel of Figure 7 shows the reconstructed image obtained using the simple tramline extraction model. The image is created by collapsing a 3D SPIRAL data cube over the spectral dimension, to create a low resolution (0.7 arcsec pixels) image of a standard star observation. The vertical structure is not a telescope diffraction spike, but rather the manifestation of fibre-to-fibre crosstalk. The SPIRAL fibres are arranged in 16 banks (IFU short axis) of 32 fibres (IFU long axis). Hence all fibres in the central column of the IFU are subject to scattered light in the wings of the profiles of the fibres close to the bright star under observation. With a simple tramline extraction, which fails to account for this crosstalk, the central columns of the IFU exhibit an enhanced light level in all spectra in the central columns. For more exotic astronomical objects, any and all columns of the IFU with high contrast data (either in the form of intensity or spectral type variations) will suffer a similar fate.

For the central pane of Figure 7 the extraction was performed using the Gaussian weighting. While this methods optimises the signal-to-noise ratio for each spectrum, it fails to account for the cross-contamination between the spectra.

The right most image shows the same data after applying the optimal extraction technique. A round stellar image is finally achieved, with an order of magnitude reduction in any signature of fibre-to-fibre cross contamination.

We have presented an optimal extraction methodology for use with multi-fibre spectroscopy in the regime where fibres are tightly packed onto the detector array. We demonstrate that high accuracy can be achieved with this approach, minimising the impact of fibre-to-fibre cross-contamination. Additionally we show that the Gaussian least squares fitting extraction used by default for AAOmega-MOS data reduce within 2dfdr is adequate for the well-separated fibre profiles of this instrument mode provided the range of input target magnitudes is kept below Δmag <3.