RVSAO 2.0: Digital Redshifts and Radial Velocities

Michael J. Kurtz and Douglas J. Mink  

Harvard‐Smithsonian Center for Astrophysics, Cambridge, MA 02138; ,

ABSTRACT

RVSAO is a set of programs to obtain redshifts and radial velocities from digital spectra. RVSAO operates in the IRAF environment. The heart of the system is xcsao, which implements the cross‐correlation method and is a direct descendant of the system built by Tonry & Davis. emsao uses intelligent heuristics to search for emission lines in spectra, and then fits them to obtain a redshift. sumspec shifts and sums spectra to build templates for cross‐correlation. linespec builds synthetic spectra given a list of spectral lines. bcvcorr corrects velocities for the motion of the Earth. We discuss in detail the parameters necessary to run xcsao and emsao properly. We discuss the reliability and error associated with xcsao‐derived redshifts. We develop an internal error estimator, and we show how large, stable surveys can be used to develop more accurate error estimators. We also develop a new methodology for building spectral templates for galaxy redshifts, using the new templates for the FAST spectrograph as an example. We show how to obtain correlation velocities using emission‐line templates. Emission‐line correlations are substantially more efficient than the previous standard technique, automated emission‐line fitting. Using this machinery, the blunder rate for redshift measurements can be kept near zero; the automation rate for FAST spectra is ∼95%. We use emsao to measure the instrumental zero‐point offset and instrumental stability of the Z‐Machine and FAST spectrographs. We compare the use of RVSAO with new methods, which use singular value decomposition and χ2 fitting techniques, and conclude that the methods we use are either equal or superior. We show that a two‐dimensional spectral classification of galaxy spectra can be developed using our emission‐ and absorption‐line templates as physically orthogonal basis vectors.

Received 1998 March 13; accepted 1998 April 21

Subject headings: instrumentation: spectrographs— methods: data analysis— techniques: radial velocities

1. INTRODUCTION

 

Radial velocities are, along with position and brightness, among the fundamental measured values of astronomy. Recent technical advances are substantially increasing our ability to acquire radial velocity data; in the decade of the 1990s, the rate at which radial velocity measurements are taken will increase by 2 or 3 orders of magnitude. Substantial effort is required for these data to be reduced and analyzed in an accurate and timely fashion; here we describe the current reduction methods that we have developed for use by the Center for Astrophysics radial velocity and redshift programs, as well as by others.

Doppler (1841) understood that radial velocities would affect the color of stars (by analogy with the pitch of sound); Fizeau (1848, 1870) first recognized that this would mean a shift in the position of the Fraunhofer lines. Huggins (1868) made the first (visual) attempt (in 1862) to observe the shifts. Vogel (1892) made the first accurate photographic measurements and established most of the procedures necessary to calibrate and reduce the measurements of line positions to radial velocities.

Correlation methods for obtaining radial velocities were first suggested by Fellgett (1953, 1955), who was influenced by radar studies during World War II. Griffin (1967) was the first to implement these techniques. Note that Griffin credits Evershed (1913) with inventing the basic technique; Griffin further notes that Babcock (1955) had already built a similar instrument. Griffin's instrument performed analog correlations by physically shifting a template spectrum in the focal plane of the spectrograph, a technique that is still in heavy use today (e.g., Baranne, Mayor, & Poncet 1979).

Digital power spectrum techniques for the estimation of lag have long been known (e.g., Blackman & Tukey 1958 and references therein). Their use first became practical with the advent of digital detectors, fast digital computers, and the fast Fourier transform (FFT) algorithm (Cooley & Tukey 1965).

Simkin (1974) first showed how Fourier techniques could be used to obtain radial velocities and velocity dispersions from digital spectra. Several groups used power spectrum techniques to obtain velocity dispersions (see Sargent et al. 1977 and references therein) and obtained velocities as a by‐product of their analysis, but apparently the first use of digital cross‐correlation specifically to obtain radial velocities was by Lacy (1977), who did not use Fourier techniques, but used direct convolution with a digital mask, emulating Griffin's (1967) analog technique.

Tonry & Davis (1979, hereafter TD79) studied the use of power spectrum techniques to obtain redshifts from digital spectra and demonstrated the effectiveness of the method. TD79 invented the r‐statistic, which can be calibrated to give both the confidence and error of a measurement. The techniques and software described here are directly descended from the TD79 system; in 1990 September, radial velocity reductions at the CfA were moved from the old Data General Nova computer where the TD79 system resided onto a Unix workstation and the IRAF (Tody 1986, 1993) environment.

We began with the IRAF task XCOR by G. Kriss and routines from TD79, translated into Fortran by J. Tonry (Tonry & Wyatt 1988); these were extensively modified and resulted in xcsao version 1.0 (Kurtz et al. 1992, hereafter Paper I). Additionally, we used algorithms from the REDUCE/INTERACT system (Maker, Kurtz, & La Sala 1982), as modified by J. Thorstensen. The emission‐line finding programs had a somewhat different history. When the radial velocity reductions were moved from the Novas onto Unix, the emission‐line programs were implemented as stand‐alone C programs translated from the FORTH of TD79 by Mink & Wyatt 1992; work began in 1991 on a new IRAF task, resulting in emsao (Mink & Wyatt 1995).

The software described here has been used extensively; examples include the redshift surveys of Huchra, Geller, & Corwin (1995), da Costa et al. (1994), Shectman et al. (1996: Las Campanas Redshift Survey [LCRS]), Vettolani et al. (1997), and Geller et al. (1997). Stellar use revolves around the CfA Digital Speedometry program (Latham 1985) and its many projects and collaborations. Latham (1992) lists several of these projects.

Nordström et al. (1994) have described the techniques used by the CfA stellar group in detail; we will concentrate on issues related to galaxy redshifts, especially with the FAST spectrograph (Fabricant et al. 1998) although some stellar data will be used in § 3 on error. We will conform to the convention that user‐settable parameters will be in italics (e.g., st_lambda), parameters in the spectrum header will be in capitals (e.g., DIVCONT), and IRAF task names will be in lowercase bold (e.g., xcsao). Extensive on‐line documentation, help files, and examples, as well as the source code and executables, can be found at http://tdc‐www.harvard.edu/.

2. PRACTICAL USE OF xcsao

 

xcsao is the heart of the RVSAO system. In this section we give an overview of the correlation system. We point out critical features that investigators need to consider when setting up the reductions for a new project to maximize the efficiency of measurement and minimize systematic errors. The RVSAO system consists of several IRAF tasks, including xcsao; the algorithmic details of each of them are described in the Appendix.

2.1. Spectrum Preparation

Tokarz & Roll (1997) discuss the steps we take to obtain one‐dimensional wavelength‐calibrated spectra suitable for redshift measurements. Once the one‐dimensional spectra are in hand, it is necessary to tell xcsao the wavelength range over which the data are good, using the parameters st_lambda and end_lambda. Depending on the details of the instrumental setup, substantial additional error can occur if these parameters are not set or are set incorrectly. For example, if a substantial portion of the spectrum is from a region of the detector with little or no sensitivity, or is from a region of the spectrum with poor sky subtraction and strong night sky lines, the final redshift will be compromised.

2.2. Continuum and Emission‐Line Suppression

For all spectra the continuum must be removed; contpars (see § A5) performs these tasks here and was only slightly modified from the IRAF ONEDSPEC continuum package, as implemented in the RV package (Fitzpatrick 1993).

The continuum can be subtracted or divided out. Subtraction preserves the correct relative amplitudes for the lines in the data and thus the correct signal‐to‐noise (S/N) behavior in the cross‐correlation. Division preserves the correct equivalent widths of the lines and thus the correct parameterization of the spectrum. For spectral classification studies, division is preferred (Kurtz 1982); for radial velocity measurements, subtraction is superior. Division by the continuum results in amplified noise in the blue part (the low S/N part) of the spectrum. For moderate S/N spectra, typical of FAST, the difference between the two techniques is small, but subtraction shows smaller redshift residuals by about 25%; while for very high S/N spectra on FAST, typical of our calibration spectra, the division method gives slightly smaller residuals.

In normal use of xcsao, only continuum subtraction is allowed, but if the parameter DIVCONT is set true (T) in the template spectrum header, the continua of both the object and template spectra will be divided out.

Many galaxy spectra show both emission and absorption lines; in general, the redshifts derived from the emission lines will be different from the absorption‐line velocities, as they originate in physically different places. To obtain a correlation velocity from an absorption‐line template for a spectrum with strong emission lines, it is necessary to suppress the emission lines; Figure 1 shows a typical spectrum and its correlation function with one of our standard absorption templates suppressing the emission lines and without suppressing them; the bottom panel shows the spectrum, smoothed and with the lines marked. The reduction with emission‐line suppression produces a believable redshift, with an r‐value (TD79) of 4.50, the reduction without it, where the r‐value is 1.88, does not.

Fig. 1.— The effect of emission‐line suppression on absorption‐line correlation redshifts. The upper panel shows the observed spectrum; the second shows the correlation function with an absorption‐line template after emission‐line clipping; the third panel shows the correlation function when the emission lines are not clipped; the bottom panel shows the observed spectrum after smoothing, and with the main lines marked.

Open New Window

Removing emission lines before correlating with an absorption‐line template was a routine feature of the TD79 software; in xcsao we have extended and generalized the procedure. Both emission and absorption lines may now be removed, and the process may be controlled by keywords in the template header. Because we now obtain emission‐line velocities using correlation methods (§ 4.1), the same emission‐line suppression cannot be used for every template, lest we remove the emission lines prior to trying to measure them.

For routine redshift reductions at the CfA, we use a subset of the xcsao capabilities, with emission lines replaced by the continuum when we correlate against an absorption‐line template, and the absorption lines replaced by the continuum when we correlate against an emission‐line template. The exact parameters used are stored in each template; normally 2 σ variations above (below) the continuum are sufficient to remove emission (absorption) lines, with the number of iterations and growing parameter as set in contpars (§ A5).

The emission and absorption‐line suppression is controlled by the parameters s_emchop (t_emchop) for the object (template) spectra. Other parameters are used as well, and when control is given to the template or object spectrum header, the interactions can be complex. They are more fully described in § A1.

In addition, xcsao permits the user to replace specific regions of the spectrum with a simple linear approximation to the continuum. This feature is typically used to prevent poor subtraction of the bright night sky lines from compromising the results.

2.3. Apodization, Zero Padding, and Fourier Filtering

Once the continuum has been removed, the spectra are apodized, zero padded, and bandpass filtered. Each of these operations has a substantial effect on the final results.

Apodization is the simplest; essentially the goal is to remove any ringing in the Fourier transform, by forcing the ends of the spectrum smoothly to zero, while suppressing as little of the actual data as possible. Our apodization is performed by a cosine taper function, which begins a symmetric set percentage of the spectrum from the ends. For both FAST spectra and the earlier Z‐Machine (Latham 1982) spectra, 0.05 is a reasonable value; i.e., the taper begins 5% from the ends (§ A1).

Zero padding the spectrum is intended to remove any artifacts caused by computation of the correlations in Fourier space, thus using a circular convolution. The primary artifact that the zero padding removes is the confusion of Hα with [O ii] λ3727 due to the wraparound of the convolution. The zero padding has two side effects that must be considered. First, the relation of the TD79 r‐statistic with error is changed, although the actual calculated error remains correct (§ 3). Second, the envelope of noise fluctuations of the correlation function, which is flat in the nonzero‐padded case, is, in the zero‐padded case, a symmetric linear function of the number of overlapping nonzero pixels and is maximum at the redshift of the template spectrum. For the case of low S/N spectra where the redshift is substantially different from the template's redshift, this structure in the noise can result in the wrong correlation peak being chosen. Zero padding can be controlled via the parameter file or the template header. We only zero‐pad spectra when correlating against the emission‐line template.

The design of the Fourier bandpass filter is critical to the optimal measurement of redshifts. As with the apodization, we use a cosine taper to suppress the ends of the (in this case) Fourier spectrum. Several other filter techniques were tried (e.g., Oppenheim & Schafer 1975), but no difference was seen for any reasonable choice of taper function (a sharp cutoff is not reasonable because of Gibbs 1899 ringing). In addition, we tested a spectral weighting function shown by Hassab & Boucher (1979) to produce the maximum likelihood estimator for the lag (radial velocity) in the limit of infinitely wide spectra; this weighting function had no positive effect on our results, and we have not implemented it.

Removing high spatial frequency information, via a high‐stop Fourier filter, is intended to increase the S/N by removing information that contains more noise than signal. The design question is where to set the high‐frequency turnoff; the method we use is to examine sets of high S/N calibration spectra, e.g., our nightly exposures of NGC 4486b. We correlate each of these against the best‐match template (in this case, NGC 7331), using a high‐pass Fourier filter that filters out all the low‐frequency information, leaving only the high‐frequency noise. If the turn‐on frequency is set too high, the correlations all give an incorrect redshift; if the turn‐on frequency is set too low, all correlations give the correct redshift. We choose the turn‐on frequency where half the redshifts are correct and half incorrect, and set this turn‐on frequency to the turnoff frequency for our high‐stop filter. This procedure gives a turnoff frequency approximately equal to that obtained by the “optimal filter” method of Brault & White (1971), which chooses that point where the power of the signal is twice that of the photon noise. This frequency is less than half that which corresponds to the projected slit width of the FAST.

The high‐stop filter is only used for absorption‐line spectra. For emission‐line spectra, the redshifts are seriously degraded if the high spatial frequencies are removed from the data. The key word FI‐FLAG in the template header controls the implementation of the high‐stop filter. In addition, as it is possible to prefilter the templates, to save unnecessary computing, this flag also controls whether and how to filter the templates. Table 2 shows all the possibilities.

Removing low spatial frequency information, via a low‐stop Fourier filter, is intended to remove any residual large‐scale systematics that remain following the continuum suppression. Essentially this can be viewed as a second continuum removal, equivalent to the continuum removal technique of La Sala & Kurtz (1985), thus giving what Kurtz & La Sala (1991) call a “reflattened” spectrum.

The design question is where to put the low‐frequency turn‐on point. The difficulty with making this decision is that there is no point where excluding all information with higher spatial frequencies does not result in a redshift (i.e., even the lowest spatial frequencies still contain accurate redshift information), and there is no reasonable point where suppressing more low‐frequency information does not result in lower residuals for high S/N sets of spectra, such as our set of NGC 4486b spectra.

We therefore set the low‐frequency turn‐on point by a simple heuristic. We estimate the scale in wavelength of the broadest spectral feature useful in estimating a redshift, in the case of FAST galaxy spectra, this is the change in the slope of the continuum around the Ca ii H + K lines, and we suppress those spatial frequencies that correspond to twice this scale, or greater.

The remaining design decision for the Fourier filter is the width of the turn‐on and turnoff ramps. It may be expected that this is only important at the low‐frequency turn‐on point, as the power there is typically 2 orders of magnitude above the high‐frequency turnoff point. The problem is that if the turn‐on is too sharp, Gibbs ringing will be introduced into the data. A full turn‐on width of 1.5% of the width of the power spectrum is sufficient to ameliorate this effect.

The exact filter implemented, especially the exact implementation of the low‐stop filter, affects the resulting redshifts, their errors, and the relation of the TD79 r‐statistic with their errors. For example, for the set of NGC 4486b spectra, the mean redshift obtained using only the lowest spatial frequencies that we include in our standard filter differs from the mean redshift obtained by using only the highest included spatial frequencies by 61 km s−1, and different reasonable choices for the Fourier filter can give redshifts that differ in the mean by 10 km s−1. These differences may be compared with a typical variation about the mean of 15 km s−1 (1 σ). The sign and amplitude of this effect changes with each object‐template pair; NGC 4486b versus NGC 7331 is a typical result. In addition, the relation of the TD79 r‐statistic with error depends on the filter (§ 3).

2.4. Cross‐Correlation, Rebinning, and Redshift Evaluation

The cross‐correlation is the normal product of the Fourier transform of the object spectrum with the conjugate of the transform of the template spectrum, as described in TD79.

The object spectrum and the template spectra need to be pairwise rebinned to have a common dispersion. The number of bins nbins is set by the user and must be a power of 2; we recommend that nbins always be larger than the number of observed pixels. The spectral region rebinned is set to obtain the maximum overlap between the template spectrum and the portion of the object spectrum between st_lambda and end_lambda in the rest frame. On the first pass, the rest frame is determined by a user guess to the redshift (czguess) or from a previous reduction (§ A1). On subsequent passes (if nzpass > 0), the rest frame is determined from the redshift obtained in the previous pass. We recommend that czguess be set to the approximate redshift expected (normally a better guess than zero), and nzpass = 2.

Next, the correlation peak is determined and fit (§ A1). The type of fit (pkmode) has little effect on the result, but the amount of the peak that is fit, pkfrac, is critical. The fit is performed from the top of the peak down to where the peak is pkfrac of the maximum; thus more of the peak is fit if pkfrac = 0.5 than if pkfrac = 0.7. Because of sidelobes in the correlation function due to the proximity of N ii to Hα, emission‐line correlation peaks cannot be fit as far down the peak as absorption‐line correlations. The template header parameter PEAKFRAC overrides pkfrac on a template‐by‐template basis; we use this to set the fit parameters for the emission‐line templates.

3. ERROR ANALYSIS

 

There are three main questions concerning the output of xcsao: Are the results reliable, can the process be automated? What is the size and nature (random or systematic) of the error? and Does what is measured correspond to the physical property the investigator wants to measure?

In this section, we examine methods for assessing the reliability and error associated with RVSAO redshifts. While the detailed results apply only to the instruments and templates we use in the analysis, the general nature of the results, and the calibration methods we employ, should have wide applicability.

3.1. Reliability

To answer the first two of these questions, we use a data set designed for this purpose; it contains 626 pairs of spectra observed with the FAST spectrograph between 1994 and 1996. Each pair consists of two independent observations of the same object; about half of these were observed to calibrate the velocity errors for the 15R survey (Geller et al. 1998), and about half are spectra that were below the quality standard and required a second integration (these would be summed in our normal reductions, but not here). All these spectra have been subjected to our normal processing (Tokarz & Roll 1997) and thus have had most of the cosmic rays removed by a labor‐intensive process.

As discussed by TD79, the r‐statistic can be calibrated as a confidence measure. We prefer to calibrate it empirically, rather than use the prescription in TD79. Figure 2 shows the results of correlating each of the 1252 spectra against each of two templates. Plotted are the absolute value of the velocity difference between the two observations versus the minimum r‐value for each pair. Each pair appears on the graph twice. The open circles are measurements using the NGC 7331 template, and the filled triangles use an emission‐line template, emtemp.

Fig. 2.— The blunder diagram; see text for discussion. The dotted lines are at and velocity differences of 300 km s−1; the dashed line is at ; both axes are on a logarithmic scale.

Open New Window

Figure 2 shows two groupings: those with the absolute velocity difference km s−1, which we will assume are reliable observations (differences ≲300 km s−1 are consistent with the expected random errors); and those where the velocity difference is greater, which we will assume are unreliable. For the moment, we will adopt an r‐value of 3.0 (rmin), above which we expect the velocity determination to be reliable. This is lower than is typically used for FAST reductions.

We would then expect no points in the upper right quadrant of the plot, where min and the km s−1; however, there are 15 points in what we will call the “blunder” region.

We will examine each point to determine which rules would catch the blunders in a fully automated reduction. The two nearly coincident circles (circles are reductions with the NGC 7331 template) at and velocity difference about 50,000 km s−1 are both objects with strong emission lines. It is known that the NGC 7331 template systematically returns velocities of 49,000 km s−1 for some emission‐line objects, so along with the high r‐value emission‐line velocity, these measures could be discarded automatically. The circle at , km s−1 is also an object with strong emission lines. A simple rule that requires that spectra with discordant, but otherwise valid, velocities be checked manually would catch this, but the reduction could not be fully automatic.

The four triangles (triangles are the emission‐line template) between and having a km s−1 all have good absorption‐line velocities and would be caught by the rule of discordance. They would be caught by another rule, however, one that does not require that the absorption‐line velocity be “good”; they are all spectra where N ii λ6583 is stronger than Hα, and the difference with the absorption‐line templates is about 948 km s−1. We adopt the rule that all emission‐line velocities that differ from an absorption‐line velocity (even if one with a low r‐value) by about 948 km s−1 must be checked manually.

The circle near and km s−1 is probably not a real blunder. Examination of the POSS prints shows that the object has two nuclei with 4 separation. We assume that the velocity difference is real and that the two observations of this object each correspond to a different nucleus.

All of the seven remaining objects in the “blunder” region of the plot are emission‐line velocities. One has the night sky line at 5577 Å mistaken for [O ii] λ5007; we can eliminate this error either by turning the bad‐lines removal feature on to replace the region around 5577 Å with the continuum or by adopting the rule that all emission‐line redshifts near 34,152 km s−1 be examined manually.

The remaining six are all spectra contaminated by cosmic rays. Five of these would be tagged by the rule of discordance, and the sixth has an absorption‐line r‐value of 2.77, so it would be tagged by only a slightly more stringent rule of discordance; it cannot be assumed, however, that the cosmic‐ray problem can be solved by looking for absorption lines. We therefore require that emission‐line redshifts must all be checked using emsao (§ A2), that at least four lines must be found that correspond to the correlation velocity, and that at least two must be fit; “blunders” are made when only three lines are found or only one line is fit. Using that criterion, all six spectra would be tagged for visual inspection, as well as 30 of the 178 emission‐line velocities ( ) that do not have a confirming absorption‐line velocity with .

For the 610 objects (of 626 total objects observed) where at least one of the two different template reductions gave a result with , xcsao yielded the correct result with no further problem for 595. Of the remaining 15 spectra, 12 are easily discovered because two valid redshift measures disagree; one is probably caused by source confusion on the sky, and two are found by emsao.

Looking at Figure 2, it is clear that many spectra where do indeed give the correct redshift. Of the 16 objects that have neither emission nor absorption reduction with , four have both emission and absorption redshifts equal (within normal errors) and could be accepted (we do not currently do this). Also 66 spectra where the emission‐line r‐value is greater than 3 have confirming absorption‐line velocities with ; if these are assumed correct (which we also do not currently do), then the number of emission‐line spectra that must be visually inspected after emsao would drop from 30 to 19. This would bring the total number that require visual inspection to 31, or 5%.

With the aid of emsao for quality control and partial manual reduction of 31 spectra, xcsao obtained the correct redshift for all 614 objects that yielded redshifts, save for the one object that was probably an observational error.

A second experiment was made using 8606 emission‐line spectra from the Z‐Machine archive, which had r‐values with the emission‐line template above 3. The results of the correlation with the emission‐line template were compared with the stored redshift in the archive, which was obtained by manually fitting the emission lines with a precursor program to emsao. After sifting the results using rules like those described above, 15 spectra (0.2%) had the wrong redshift; essentially all these spectra were the victims of very poor sky subtraction. This may be compared with the 24 spectra where the redshifts were incorrectly listed in the archive. The Z‐Machine spectra were substantially noisier than the FAST spectra; nearly 15% failed the sifting and would have had to be manually reduced.

Used carefully, the RVSAO suite provides redshifts with a very low blunder rate. The automation rate obtainable with RVSAO is strongly affected by the S/N of the observations. Absorption‐line objects must be observed long enough to have a fully reliable absorption correlation velocity (we currently use , which is conservative) or a confirming weak emission velocity. Emission‐line spectra must have a confirming weak absorption redshift, or be based on at least four lines, and at least two of the four must be fit by emsao (§ 4.1).

3.2. Error Estimation

Besides estimating the redshift of a spectrum, xcsao also estimates the error in the redshift. The error estimator can be derived analytically following the discussion in § III.c.i of TD79 with the additional assumption of sinusoidal noise, with the half‐width of the sinusoid equal to the half‐width of the correlation peak. The derived error estimator is where error is the error in a single velocity measurement by xcsao, w is the FWHM of the correlation peak, and r is as defined in TD79.

Although the assumption of sinusoidal noise with half‐width equal to the correlation peak's half‐width is reasonable, there is no compelling argument for this assumption. Therefore it is necessary to demonstrate the effectiveness of the approximation by experiment.

We will use four data sets to examine the behavior of the error estimator: (1) the 610 duplicates described above; (2) the 8606 Z‐Machine emission‐line spectra described above; (3) 7810 synthetic spectra, each identically Poisson‐sampled from a 45 Å section of a model atmosphere for a 5500 K dwarf star (Kurucz 1992), taken from the set of synthetic stellar templates used by the CfA digital speedometry program (Morse, Mathieu, & Levine 1991; Nordström et al. 1994); and (4) 50,000 synthetic spectra, using the same 5500 K dwarf star template, each with a different number of simulated photons (we confine ourselves to using the 49,880 spectra that, when correlated against the synthetic template, achieved ).

First we will look at the set of duplicate spectra. We will limit ourselves to cases where both reductions have r‐values greater than 3.5. This is spectra for the two absorption‐line combinations and 297 for the emission‐line comparison.

Figure 3a1 shows a typical result. The solid line shows a histogram of the absolute values of the differences between two observations of the same object, both reduced in the same way using the NGC 7331 template and divided by the sum in quadrature of the errors calculated by xcsao for the reductions. The smooth, thin curve is the expected Gaussian distribution; it is clearly broader than the data. xcsao overestimated the error by ∼20%.

Fig. 3.— The thick, solid histograms are the distributions of velocity differences for duplicate observations of the same object, divided by the error estimate. The thin smooth curves are the expected Gaussian distributions. In the lower panels (labeled 1), the error estimate is the 3w/8( ) xcsao internal error estimator. In the top panels (labeled 2), the error estimate is the calibrated error estimator. The left panels (a) are for reductions with the fn7331temp template, the middle panels (b) are reductions with the ztemp template, and the right panels (c) are reductions with the emtemp template.

Open New Window

Figure 3b1 is similar to Figure 3a1. Here the ztemp template was used on the FAST data; ztemp is a combination of bright galaxy spectra taken with the Z‐Machine and in use at the CfA since the days of TD79. ztemp has a restricted wavelength coverage (λλ = 4500–6200 Å) compared with the FAST spectra, has a different resolution, and has different residual systematics. In Figure 3b1 the (smooth, thin) Gaussian is narrower than the (thick, histogrammed) data. xcsao underestimated the error by ∼20%.

Figure 3c1 shows a similar histogram for the emission‐line template, emtemp. emtemp is a synthetic spectrum made before the creation of the linespec task (§ A3) to match FAST emission‐line galaxy spectra. Here the thick histogram that represents the data cannot be transformed to match the thin, smooth expected Gaussian by any multiplicative process (1.2 would be the best multiplicative factor); it would still have more power in the tail. Adding 15 km s−1 in quadrature helps remove power from the tail.

TD79, while giving a procedure to calculate the error, suggest that in practice the error be calculated by calibrating k in the equation using external comparisons. Paper 1 reiterates this suggestion, noting that the measurement of w has error, but when all reduction parameters remain fixed, w is essentially constant. Figure II of Paper I shows the effect of changing one of the reduction parameters (the low‐frequency roll‐off of the Fourier filter), which substantially changes the relation of ( ) to error, while w scales correctly so that 3w/8( ) still tracks the error.

Given a set of duplicate observations, the constant k can be determined by internal comparisons. The procedure is simply to vary k until the expected differences histogram matches the measured one. For the present case, we obtain km s−1, km s−1, and km s−1. Figures 3a2, b2, and c2 show the distributions compared with the expected Gaussians; as expected, in all cases the fit is better than with the unmodified xcsao errors. For large observing programs with stable reduction procedures, we recommend using calibrated relations to estimate the error.

The duplicate spectra, as is typical of single‐slit redshift survey data, do not show a very large range of S/N, or r‐value. This is due to the fact that one normally observes long enough to get a desired S/N and no longer. If we continue to restrict the duplicate pairs to those where both spectra achieved (still lower than is normally required for a FAST redshift), then it is not possible here to use the goodness of fit to a Gaussian to prove that the error calculated by using is any better than a constant error, for each template. If the error were dominated by systematics, one would expect the error to be approximately constant.

We define a measure of error Here are the r‐values of the two spectra (of the same object) being compared. We then vary the values of ksystematic and kstatistical and calculate the residuals of the fits to a Gaussian. Note that ksystematic is times the error in a single measurement. Figure 4 shows the results for the duplicate pairs reduced with emtemp. The contours represent lines of equal residuals in the fit to a Gaussian, in the (kstatistical, ksystematic) space. The outer contour is 20% larger than the inner contour. Clearly either a fully systematic or a fully statistical error is consistent with the data.

Fig. 4.— Residuals in fit to Gaussian error distribution for pairs reduced with the emtemp template. The contours represent lines of equal residual, and position in the (x, y) space represents admixtures of the statistical and systematic error models; see text.

Open New Window

Figure 5 shows the same diagram for the N7331 template. The inner contour level here represents an absolute error half that of emtemp, and the outer level is 33% larger than the inner level. Also here one cannot rule out either a fully systematic or a fully statistical error.

Fig. 5.— Residuals in fit to Gaussian error distribution for pairs reduced with the fn7331temp template. The contours represent lines of equal residual, and position in the (x, y) space represents admixtures of the statistical and systematic error models; see text.

Open New Window

In many cases, it is not possible to reduce hundreds of duplicate measurements using exactly the same reduction procedures to obtain an improved error estimator; in these cases, the xcsao error estimator is a reasonable choice. The 20% systematic deviations for the two absorption‐line galaxy templates are the largest we have seen, although Quintana, Ramirez, & Way (1996) suggest that for their data the error is underestimated by ∼30%, by comparison with external measurements. Although the xcsao error estimator differs systematically from the true error for a particular combination of instrumental setup, reduction procedure, and template, we have not seen any trend for this to be a systematic over‐ or underestimate.

While the 8606 Z‐Machine emission‐line spectra cannot be used to calibrate the error estimator, as to first order we are just comparing the differences in two different methods of fitting Hα in the same spectrum, we can use them to look at any differences in zero point, as a function of the fitting method. The mean difference between the correlation velocity and the velocity obtained by the semiautomated line fits is km s−1, which is small compared with the other sources of error, and may be due to small systematic differences in the wavelength calibration in the red and green parts of the spectrum (§§ 4.3 and 5).

The 7810 synthetic spectra are each correlated against the template from which they were identically randomly Poisson sampled; thus there is no spectral type difference adding to the errors. Figure 6 shows the template, a typical sampled spectrum, and a smoothed version of the typical spectrum. By calculating the rms velocity about the expected velocity (0 km s−1), we have a very accurate measure of the error in a single measurement; the mean error calculated by xcsao is 3% greater than this. The distribution of velocities is essentially Gaussian, and the deviation of the mean velocity is within 1 σ of zero.

Fig. 6.— The synthetic 5500 K spectra. The top panel is the template spectrum (Kurucz 1992), the middle panel is a typical Poisson‐sampled spectrum, and the bottom panel is the same spectrum smoothed.

Open New Window

Using the 7810 spectra, we can ask the question, “How many independent measurements of a spectrum are required to obtain a better measure of the error in a single measurement than xcsao provides?” Nordström et al. (1994), on the basis of a study of echelle spectra of rotating F stars, give this answer as ∼7. Here we take as many independent sets of N spectra as exist in 7810 spectra, for each N we calculate the sample standard deviation about the sample mean, and we compare it with the known error gotten by using all 7810 velocities. The rms of this difference is plotted as a function of N in Figure 7. Also plotted (as a straight line) is the rms of the xcsao error estimate about the true error. In this ideal case, more than 30 independent measures are required before a better error estimate is reached than the xcsao error.

Fig. 7.— Variance in error estimation as a function of the number of measurements for the set of 7810 synthetic spectra. The horizontal line is the variance (about the true value) of the xcsao error estimate. See text for details.

Open New Window

The 50,000 synthetic spectra, with 1–50,000 counts, give similar results. The mean error is underestimated here by 8%, and the mean velocity is 1.1 σ different from zero, using the 49,880 spectra where . The efficiency of ( ) in estimating S/N is demonstrated in Figure 8; here we show the ratio of the square root of the number of counts to the best‐fit linear relation with ( ). The 1 σ scatter is ∼12%, independent of N; this puts a limit on the inherent ability to estimate errors using ( ).

Fig. 8.— The relation between r and Each point represents the ratio of with the best‐fit linear relation between and .

Open New Window

3.3. Systematics

There are many factors that can cause the redshifts and radial velocities measured by xcsao to be other than those desired. Here we list several:

1. 

Errors in the wavelength calibration. For absorption‐line spectra, where the signal is averaged over all the lines in a complex spectrum, this effect should be tiny. For emission‐line spectra, where the signal comes from a couple of lines, this effect could be as large as the error in the pixel‐to‐wavelength calibration function. For typical FAST galaxy spectra, this error is ∼5 km s−1 averaged over the entire wavelength range; it could be larger in small regions (§ 4.3).

2. 

Offsets in the calibration lamp illumination. The light from the calibration lamps does not follow the same exact optical path as the light from the sky. This can cause systematic errors in the wavelength scale (§ 5).

3. 

Zero‐point errors in the template. Correlation velocities are actually velocity differences between the object spectrum and the template. Any error in determining the redshift of the template spectrum will propagate into the object being measured (§ 4.3).

4. 

Variations with the Fourier filter. As noted above (§ 2.3), different reasonable choices of the Fourier filter can change the measured velocity of a typical FAST galaxy spectrum by 10 km s−1, comparable to the error in our highest S/N observations. Unreasonable choices for the filter parameters can make more of a difference.

5. 

Spectral type mismatch. It has long been known that there are systematic velocity effects when the template spectrum does not match the observed spectrum. Nordström et al. (1994), for example, show the effect of rotational velocity mismatch for their echelle spectra. For typical FAST galaxy spectra, this effect is ∼20 km s−1 and is discussed in § 4.3.

6. 

H ii regions in spiral galaxies. The rotation velocities of disks and the finite number of H ii regions on the slit can combine to yield an Hα velocity that is different from the mean velocity of the stars in the bulge; this is especially true for galaxies that are distorted. Thus, while the measurement error in an emission‐line velocity for a particular galaxy may be substantially smaller than for the absorption‐line velocity, the systematic deviation from the desired quantity, the cosmological redshift, may be substantially larger.

7. 

Two‐lined systems. xcsao assumes that the template spectrum is a reasonable spectral match to the observed spectrum. For the case of two (or more) lined spectroscopic binaries, this condition is clearly violated. As demonstrated by Latham et al. (1996), substantially improved results may be obtained by using methods that explicitly model the two‐lined case, such as TODCOR (Zucker & Mazeh 1994). xcsao may be used to obtain the input data for TODCOR.

4. TEMPLATES

 

Accurate redshifts require the existence of very high signal‐to‐noise templates, which have well‐determined velocities and are good matches to the scientific program objects being measured. Here we describe methods of creating and maintaining systems of templates.

Templates must be very high signal‐to‐noise spectra; there are two basic ways to create them: (1) sum a number of observations, or (2) build a computer model of a spectrum. The CfA Digital Speedometry group has, over the past decade, switched entirely from using observed spectra to models (see Latham et al. 1996 and references therein for details). For galaxy redshift studies, we use both techniques.

The vast majority of nearby galaxies can be well matched by a typical absorption‐line spectrum (like NGC 7331), a typical emission‐line spectrum, or both. Unusual spectra, as one obtains for QSOs, Hδ strong galaxies, galaxies with extremely high or low internal velocity dispersions, etc., require special templates, although in most cases special templates only lower the error in the redshift; the “correct” redshift (i.e., the redshift that is related to the true redshift by the measurement error, not a “blunder” that is caused by assigning the redshift using a noise peak in the correlation function) is normally obtained using standard templates.

Creating templates is an iterative process: good templates are required as a prerequisite for making better templates. New templates were made for FAST when it saw first light, and again when it received a new thinned CCD in 1994 September. In 1997, sufficient new observations having been made, we were able to make a set of substantially improved templates.

4.1. Emission‐Line Template

Until 1995, redshifts for emission‐line galaxies were obtained with emsao and its precursor programs. Then a template was made by placing Gaussians with approximately correct line widths and line ratios at the emission‐line rest wavelengths. This template, emtemp, was tested against existing FAST observations and against the Z‐Machine archive. Typical results are for the Z‐Machine comparison with 8606 emission‐line spectra reduced by hand: the 1 σ difference is 13 km s−1, about half of the calculated error for a typical Z‐Machine emission‐line velocity. As noted above, the zero‐point offset is km s−1, about a two‐hundredths of a pixel. The median difference of 3929 FAST spectra with and the Hα velocity obtained by emsao is 0.24 km s−1.

Essentially for all spectra where emsao can obtain a redshift, xcsao obtains a redshift using emtemp. For low S/N spectra, xcsao plus emtemp is much more sensitive than emsao. For a sample of 2088 emission‐line galaxy spectra taken with FAST, only 42% of spectra where xcsao plus emtemp obtained a redshift with could be reduced automatically using emsao, 64% for , and 93% for . As noted above, the systematic effect of cosmic rays places severe constraints on the unsupervised use of xcsao with an emission‐line template.

To make a better emission‐line template, we have taken 6498 FAST spectra where , put them through emsao, and obtained 434 spectra where emsao found nine or more lines. For these we measured the ratio of line heights with the S ii λ6731 line and the line widths for each line. Then for each line we took the median of these quantities and, along with the laboratory rest wavelengths for the lines, put the data into the program linespec to produce a synthetic template that we call femtemp97 (Fig. 9).

Fig. 9.— The emission‐line template femtemp97

Open New Window

femtemp97 is indeed a better template than emtemp. Using the set of 626 duplicate spectra described in § 3.1, we can compare the velocity differences between pairs directly; the median difference using femtemp97 is ∼9% smaller than in the reduction using emtemp. The errors for femtemp97 are more Gaussian than for emtemp; Figure 10 shows the fit to the two‐parameter error model (§ 3.2) and may be compared with Figure 4. The interior contour level here is half that in Figure 4. Also note that the y‐intercept is about 32 km s−1, substantially less than the 45 km s−1 in Figure 4, implying a greater reduction in error than the direct comparison of velocity differences would indicate. The median ratio ; femtemp97 yields good velocities for a substantial number of spectra where emtemp fails. A comparison of with the 3w/8( ) error estimator shows that 3w/8( ) underestimates the error by 9%.

Fig. 10.— Residuals in fit to Gaussian error distribution for pairs reduced with the femtemp97 template. The contours represent lines of equal residual, and position in the (x, y) space represents admixtures of the statistical and systematic error models; see text.

Open New Window

4.2. Absorption‐Line Template

The first absorption‐line templates used on the FAST data were the TD79 vintage ztemp and the NGC 4486b template from the MMT spectrograph, as well as some secondary MMT templates. Over the next year, a program of template observations provided several high S/N observations of candidate templates. These were summed on an object‐by‐object basis, to form extremely high S/N templates. The best of these is the NGC 7331 template (fn7331temp) (§ 3.1).

We are now in a position to create a better template. Using fn7331temp and femtemp97, we selected galaxy spectra where the r‐value for the reduction with fn7331temp was greater than 8 and the r‐value for the reduction with femtemp97 was less than 3, which excludes spectra with emission lines strong enough to give a reliable redshift. These are 1959 moderate‐to‐high S/N absorption‐line spectra. This set of spectra still contains spectra of objects with substantial emission. NGC 7331 shows clear emission in the N ii λ6583 line, for example. We examined the differences between the velocities derived using fn7331temp and femtemp97, where (1) the difference was near zero, meaning that there was enough emission present to get a correct redshift (about 100 spectra), and where (2) the difference was near 948 km s−1, meaning that N ii λ6583 was confused with Hα (about 300 spectra); we removed those spectra from the sample.

The remaining 1489 spectra were shifted to a common rest velocity, using the fn7331temp velocity; they were normalized to the same number of counts; their continua were subtracted, using a moderate‐order spline; and finally they were summed. The sumspec task (§ A4) performed these tasks. The resulting spectrum had its residual continuum subtracted using a high‐order spline (possible because of the substantially higher S/N of the sum) and was normalized to represent the average spectrum. This is the final new absorption‐line FAST template, fabtemp97, shown in Figure 11.

Fig. 11.— The absorption‐line template fabtemp97

Open New Window

fabtemp97 is a better template than fn7331temp. The median difference between pairs of spectra from the 628 duplicates of § 3.1 is actually ∼2% larger (insignificant) using fabtemp97, because fn7331temp is better able to match emission‐line objects (it has N ii visibly in emission and other lines buried in the noise); if one restricts the comparison to pairs of spectra where for each spectrum, fabtemp97 shows a median difference ∼12% smaller than fn7331temp. As with the comparison of femtemp97 and emtemp, the differences between the new and the old templates are not large.

Figure 12 shows the fit to the two‐parameter error model in § 3.2 (compare with Fig. 5 for fn7331temp). The contours in Figure 12 are the same as in Figure 5; note that even the outer contour does not reach the y‐axis. This suggests that the error for these spectra cannot be modeled by a single number, ksystematic, independent of r; this is not true for fn7331temp, or for the emission‐line templates. A comparison of with the 3w/8( ) error estimator shows that 3w/8( ) overestimates the error by 4%.

Fig. 12.— Residuals in fit to Gaussian error distribution for pairs reduced with the fabtemp97 template. The contours represent lines of equal residual, and position in the (x, y) space represents admixtures of the statistical and systematic error models; see text.

Open New Window

The analyses of the four templates, shown in Figures 4, 5, 10, and 12, all are consistent with a model where ∼20 km s−1 constant systematic error is combined with a statistical error determined by the value of the r‐statistic.

We also attempted to make a template for narrow‐lined absorption‐line objects. We built a template by summing, after shifting to a common rest frame, several high S/N spectra from a set of M31 globular clusters; we call this template fglotemp. Next we extracted from the FAST database all spectra where and and . More than 90% of these spectra were calibration stars; we excluded them, and the globular clusters themselves, and were left with ∼300 galaxy spectra. These we shifted and summed with sumspec in the same manner as for fabtemp97, giving us a new narrow‐lined template. This template and fabtemp97 were correlated against the 300 spectra. There was no significant difference in the results. From this we conclude that for typical redshift survey spectra, from a spectrograph with resolution , there is no need to have a narrow‐lined template.

4.3. Velocity Zero Point

Correlation techniques actually measure the velocity difference between the object spectra and the template. To obtain the absolute velocity of the object, the velocity of the template must be accurately known; this provides the zero point for the radial velocity system.

We have developed a new methodology for defining the velocity zero point. Previous methods have depended on external calibrators, such as 21 cm measurements; our new methods are fully internal and should minimize velocity offsets due to spectral type differences.

Figure 13 illustrates the problem with spectral type difference. With the FAST spectrograph, we observed M31 on 116 occasions and NGC 4486b on 75 occasions. Each of these spectra was reduced using ztemp and fn7331temp. The resulting velocities were shifted so that, for each galaxy, the median redshift obtained by using the ztemp template was zero. The bottom panel of Figure 13 shows the results for M31; the solid histogram shows the result of correlating the 116 spectra with fn7331temp, the dotted histogram shows the result of correlating these spectra with ztemp. Similarly, the top panel of Figure 13 shows the results for NGC 4486b, with the solid histogram representing 75 fn7331temp correlations and the dotted histogram 75 ztemp correlations.

Fig. 13.— Spectral type differences in determining velocity zero points. The top panel shows velocities for 75 observations of NGC 4486b, the solid histogram using the fn7331temp template, and the dotted histogram using the ztemp template. The x‐axis has been shifted so that the median redshift obtained using ztemp is zero. The bottom panel shows the same information for 116 different spectra of M31. No shifting of zero points can make both sets of histograms agree.

Open New Window

No change in the zero point of either template can make both sets of histograms agree; any setting of the zero point by matching velocities for one object would make the systematic difference between different template reductions for the other object worse.

We choose to define the zero point of our velocity system to minimize the systematic difference between the two main types of galaxy spectra: emission‐line spectra and absorption‐line spectra. To insure that our internal procedure matches the “true” system, we need only make the reasonable assumption that we accurately know the rest velocities of the main spectral lines in galaxies, such as Hα and [O iii]. We must also be certain that the internal wavelength system of the spectrograph matches the external wavelength system of the sky (§ 5).

Our basic procedure is to force the median difference between the emission‐line velocity and the absorption‐line velocity, for those galaxies that strongly show both sets of features, to zero. Because the absorption‐line velocity comes primarily from the K giant stars in the bulge, while the emission‐line velocity comes mainly from H ii regions in the disk, which are moving at ∼± 200 km s−1 with respect to the central bulge, there is no reason to expect that the absorption‐line velocity and the emission‐line velocity should be the same for any particular object. These differences should be randomly distributed about zero, however.

We began by setting the velocity for the fn7331temp template so that the median velocity difference between the fn7331temp velocity and the femtemp97 velocity, for spectra where and , was identically zero. This yields a velocity of 797 km s−1 for NGC 7331, which may be compared with km s−1 from the 21 cm observations (Bottinelli et al. 1990).

As described in § 4.2, fabtemp97 was created by shifting 1489 spectra to the rest frame defined by their individual correlations with fn7331temp, then summing them. A comparison of fabtemp97 velocities with emission‐line velocities should give an indication of the stability of the zero‐point technique.

Figure 14 shows the difference in velocities between fabtemp97 and femtemp97 as a function of observation date, for 1787 FAST spectra where and . The median difference is 7.1 km s−1, with an interquartile range of 62 km s−1. Selecting subsets of the spectra with higher S/N ratios has no significant effect on this result. This difference is ∼0.1 pixel and is probably due to systematic, nonlinear errors in our wavelength scale. Figure 15 shows the difference between the fabtemp97 velocity and the velocity found by fitting Hα in emsao for a representative subset of the data. The median difference is −2.1 km s−1, consistent with zero. Figure 16 shows the difference between the femtemp97 velocity and the [O iii] λ5007 velocity, −0.16 km s−1, and Figure 17 shows the difference between the femtemp97 velocity and the Hα velocity, −9.51 km s−1. These differences are very similar to the median differences between individual lines measured in the same spectrum with emsao. Figure 18 shows the difference between fits to Hα and fits to [O iii]; the median difference is 6.4 km s−1.

Fig. 14.— The difference between velocities obtained with the femtemp97 template and the fabtemp97 template, for 1787 FAST spectra where both reductions had , as a function of observation date.

Open New Window

Fig. 15.— The difference between velocities obtained by fitting Hα and correlating with the fabtemp97 template, for 1514 spectra where and Hα was fitted, as a function of observation date.

Open New Window

Fig. 16.— The difference between velocities obtained by fitting [O iii] λ5007 and correlating with the femtemp97 template, for 3527 spectra where and [O iii] was fitted, as a function of observation date.

Open New Window

Fig. 17.— The difference between velocities obtained by fitting Hα and correlating with the femtemp97 template, for 4833 spectra where and Hα was fitted, as a function of observation date.

Open New Window

Fig. 18.— The difference between velocities obtained by fitting [O iii] λ5007 and fitting Hα, for 2008 spectra where Hα and [O iii] were both fitted, as a function of observation date.

Open New Window

Figures 1418 contain a number of similarities and differences. The scatter in Figure 14 is very similar to the scatter in Figure 15, and the scatter in Figure 16 is very similar to the scatter in Figure 18, as one might expect from the small scatter in Figure 17. In addition, the small‐scale fluctuations seen in Figure 16 are more similar to those in Figure 17 than those in Figure 18, indicating that they are due to variations in the emission‐line and correlation systems, rather than variations in the nonlinear wavelength calibration errors.

We conclude that our construction of fabtemp97 gives a velocity zero point equal to the velocity zero point of the emission‐line system; this equivalence is as accurate as our ability to define the wavelength system using standard lamps and polynomial fits.

5. TESTING SPECTROGRAPH ZERO POINT AND STABILITY USING emsao AND NIGHT‐SKY SPECTRA

 

Before 1995, we measured redshifts for emission‐line spectra using emsao. Now we cross‐correlate against an emission‐line template with xcsao (§ 4.1). emsao is used to check for errors in emission‐line correlations, interactively to reduce emission‐line spectra, automatically to measure equivalent widths, line heights, and line widths for sets of emission‐line spectra, and to do various custom projects using alternate line lists (e.g., to study QSOs).

We used emsao to calibrate the Z‐Machine and FAST spectrographs for stability and zero point by measuring the apparent velocity of the night sky lines. These are from forbidden oxygen airglow lines from the upper atmosphere and mercury and sodium emission from artificial sources such as street lights. Sodium is a blended doublet, and we do not know the effective wavelength a priori, so we cannot use it to establish a zero point. Using the zero point defined by [O i] λ5577, we measure an effective wavelength of 5891.2 Å using FAST, with the Z‐Machine data in agreement. We also obtain a systematic difference between the oxygen airglow lines and the mercury streetlamp lines (5461, 4358, and 4047 Å) of about 20 km s−1 with both spectrographs. We assume that this effect is due to differences between the effective wavelengths of Hg in streetlamps and calibration lamps. We therefore choose to calibrate our zero points with the oxygen airglow lines.

Figure 19 shows the results for the Z‐Machine. The top panel is a typical sky spectrum, the second panel shows the apparent velocity of the Hg i λ5461 line, the third panel shows the apparent velocity of the [O i] λ5577 line, and the bottom panel shows the apparent velocity of [O i] λ6300. The velocities are shown as a function of observation date, over 15 yr. Each point is the median of a single nights observations; nights with fewer than 10 observations were excluded.

Fig. 19.— emsao reduction of night sky lines from the Z‐Machine. The top panel shows a typical night sky spectrum, with the lines marked. The next panel shows the apparent velocity of [Hg i] λ5461 as a function of observation date. Each point is the nightly median, where nights with fewer than 10 observations are excluded. The next panel is similar, but shows the apparent velocity of [O i] λ5577, and the bottom panel shows the apparent velocity of the [O i] λ6300 line.

Open New Window

The three lines show similarities and differences with each other. Hg λ5461, apart from having a zero‐point offset, follows the pattern shown by [O i] λ6300 more closely than the pattern shown by [O i] λ5577. In particular, the slow increase in velocity from 1984 September to 1990 June is seen clearly in the Hg λ5461 and [O i] λ6300 lines (and also in the NaD and [O i] λ6363 lines), but it is not seen in the [O i] λ5577 line. We believe that [O i] λ5577 differs from the other lines because of the manner in which the wavelength calibrations were performed and that the behavior of Hg λ5461 and [O i] λ6300 better represents the physical changes in the spectrograph than [O i] λ5577.

All the night sky emission lines (except, after 1992, NaD) were routinely used to supplement the HeNeAr calibration lamp lines in the Z‐Machine reductions, and a seventh‐order least squares polynomial fit was used to obtain wavelength from pixel position. The night sky lines were included exactly to ameliorate the problems of matching the instrumental system to the sky. Because of the detailed distribution of HeNeAr lines, the fit in the region of 5577 Å is less dominated by the HeNeAr lines (there are relatively fewer in this spectral region) than in the vicinity of the other night sky lines. This allows the region near 5577 Å to be better calibrated to the “true” system of the sky.

The 5.75 yr, 40 km s−1 drift in the zero point shown best by [O i] λ6300 would represent a correctable error, although only now, after 25,000 spectra have been taken, were the system of the spectrograph and the system of the sky not coupled in a nonlinear way by the use of the night sky lines in creating the wavelength calibration. Because of this coupling, Falco et al. (1998), in the final reduction of the Z‐Machine data, have not attempted to perform a correction but have added an amount in quadrature to the error (35 km s−1) to account for this and other systematic errors.

Figure 20 shows the results for FAST. Again the top is a typical sky spectrum, the second panel shows the apparent velocity of the Hg λ5461 line, followed by the apparent velocity of the [O i] λ5577 line, and at the bottom the apparent velocity of the [O i] λ6300 line. The abscissa covers 4 yr, and the scale of the ordinate is half that of Figure 19, in order to show the correlated nature of the variations.

Fig. 20.— emsao reduction of night sky lines from FAST. The top panel shows a typical night sky spectrum, with the lines marked. The next panel shows the apparent velocity of [Hg i] λ5461 as a function of observation date. Each point is the nightly median, where nights with fewer than 10 observations are excluded. The next panel is similar, but shows the apparent velocity of [O i] λ5577, and the bottom panel shows the apparent velocity of the [O i] λ6300 line.

Open New Window

The night sky lines are not used in the wavelength calibration of FAST, and the HeNeArFe calibration lamp lines are fit with a third‐order polynomial. The Hg λ5661, [O i] λ5577, and [O i] λ6300 lines should be good measures of the instrumental zero point and the stability of the instrument and reduction techniques, save for the systematic offset in Hg λ5461.

All three lines show essentially the same behavior, especially since the CCD was changed in 1994 September. The scatter in these line positions is less than 20% of the scatter for the same lines in the Z‐Machine. Over the 4 yr of FAST operation, an additive offset of km s−1 brings the instrumental system into agreement with the true system of the sky.

FAST is a remarkably stable instrument. The scatter in the [O i] λ5577 line can provide a good measure of the errors due to a combination of instrumental instability, wavelength calibration, and line fitting. We measure the random error in fitting a single line by dividing the scatter in the apparent velocity difference between [O i] λ6300 and [O i] λ6363 by . Subtracting it in quadrature from the scatter in [O i] λ5577 yields the error due to the interaction of instrumental instability with the wavelength calibration. For the entire 4 yr period, this is 2 km s−1, where most of the error comes from systematic changes associated with changing the CCD, changing the line list for the HeNeArFe calibration, and changing the dewar. For the most recent 18 months, since the dewar change, the error from instrumental instability is unmeasurably small, is consistent with zero, and has a 2 σ upper limit, by bootstrap resampling, of 0.7 km s−1.

6. OTHER METHODS AND COMMENTS

 

Press (1995) has suggested a new methodology for determining redshifts. One first reduces a set of galaxy spectra (shifted to a common rest frame) to a set of orthogonal basis vectors, using singular value decomposition, SVD. Next, using the fast numerical methods of Rybicki & Press (1995), the most significant of these vectors are repeatedly fit to a spectrum with unknown redshift, with each fit being at a different redshift. For each fit, χ2 is calculated; the redshift corresponding to a minimum χ2 is the redshift of the galaxy.

Recently Glazebrook, Offer, & Deelet (1998) have developed a nearly identical scheme; the difference is that rather than calculate χ2 exactly, they use the simplifying assumption that the correlation function may substitute for χ2. Thus they can treat the (SVD derived) eigenvectors as templates in a cross‐correlation program (such as xcsao) and obtain a final redshift by summing the correlation functions in quadrature, weighted by the eigenvalues of each eigenvector‐template.

Both Press (1995) and Glazebrook et al. (1998) have suggested that the coefficients of the fits can be used to classify the spectra and that the position of a spectrum in coefficient space may be used to develop confidence measures in the derived redshift. Recently, Bromley et al. (1998) used this technique to classify spectra from the LCRS (Shectman et al. 1996).

We have not adopted the SVD method for creating templates, nor the χ2 minimization technique for determining redshifts. We have, following Press's (1995) suggestion, developed methods to use best‐fit parameters for rough classification and blunder discovery.

We do not implement the Glazebrook et al. (1998) method for a number of reasons. Their approximation assumes the variance in a spectrum is not a function of wavelength. This assumption is clearly incorrect and is one of the reasons why continuum subtraction is superior to continuum division in the low S/N regime (§ 2.2). The use of the eigenvalue as a weight in the combination of correlation functions overweights emission‐line correlations when reducing absorption‐line spectra, thus making them even more susceptible to shot noise. But finally, as the exact methods are available (Press 1995), we see no advantage in the approximation.

Press's (1995) χ2 minimization technique, enabled by the Rybicki & Press (1995) algorithm, shows substantial promise. Especially for the case where the rest frame of the template is much different from that of the unknown galaxy, we expect that a rigorous accounting of observational errors as a function of wavelength will be important; the new, deep redshift surveys of the next decade, e.g., with Hectospec (Fabricant, Hertz, & Szentgyorgyi 1994), will test whether the χ2 techniques will obtain better results for high‐redshift objects. We do not now implement the χ2 minimization because we believe that current redshift survey data, with , would not benefit; some confirmation of this view comes from the preliminary results of Press (1995; W. H. Press 1997, private communication), who finds no improved ability to obtain redshifts from low S/N spectra in the LCRS (Shectman et al. 1996) compared with the original reduction, done by H. Lin using RVSAO 1.0.

We believe our template creation techniques are superior to a SVD decomposition of a group of spectra for underlying physical reasons. Emission‐line spectra and absorption‐line spectra arise from independent physical causes in different physical locations within a galaxy; basis vectors that are admixtures of absorption‐ and emission‐line spectra make no physical sense. Additionally, emission‐line velocities and absorption‐line velocities are not identical for any particular galaxy, as a perusal of optical rotation curves (e.g., Barton et al. 1998) clearly shows.

The simultaneous use of more than two templates or basis vectors to obtain redshifts is unnecessary. We performed a SVD decomposition on the 1489 pure absorption‐line spectra used to create fabtemp97; the first eigenvector was essentially identical to fabtemp97; the next three represented small differences in the continuum subtractions; the fifth eigenvector shows the H + K lines, which are systematically weaker in the early data, before a blue sensitive chip was installed; the next couple of eigenvectors also represent continuum differences. The systematic residuals from the reduction dominate the higher dimensions of the SVD decomposition.

As a practical matter, the limiting factor in obtaining redshifts for faint objects is the ability to get redshifts for weak‐lined absorption‐line systems. Essentially, if one observes long enough with a fiber instrument to get redshifts for the absorption‐line spectra, the spectra with emission lines will all yield redshifts with almost any technique.

We therefore will use fabtemp97 and femtemp97 as our fitting functions; we fit each spectrum as a linear combination of these two templates. We take as a representative sample the last 2000 spectra observed as part of the 15R survey (Geller et al. 1998). Figure 21 shows the absorption‐line versus emission‐line fit coefficients. The range of the coefficients is a function of the normalization of each template and of each object spectrum. Because the emission‐line template has no continuum, the normalization is arbitrary. The solid dots are for and the open circles are for . There is a clear locus; outliers could be easily discovered and removed from any fully automatic data reduction process. Combined with more traditional measures, e.g., equivalent widths and line ratios, Figure 21 can form the basis for a spectroscopic classification. The x‐ and y‐axes roughly measure absorption‐line strength (metallicity perhaps) and emission‐line strength.

Fig. 21.— The coefficient of the fit to fabtemp97 vs. the coefficient of the fit to femtemp97, for 2000 spectra from the 15R survey (Geller et al. 1998). About a dozen objects have coefficients not in this range and would be looked at manually. The filled dots are spectra where , and the open circles are for spectra where . The diagram can be viewed as a crude absorption‐line strength vs. emission‐line strength classification diagram.

Open New Window

Figure 22 shows the relation between the fit residuals and the inverse of the fab statistic for the fabtemp97 reduction. The symbols have the same meaning as in Figure 21. There is a very good correlation of fit residuals with for spectra where ; for objects with emission lines, the correlation is much weaker.

Fig. 22.— The relation of the residuals to the fit to fabtemp97 to the inverse of fabtemp97. The symbols have the same meaning as in Fig. 21.

Open New Window

Figure 23 shows the relation between the fit residuals and the inverse of the fem statistic for the femtemp97 reduction; the symbols here are filled dots if and open circles if . The correlation of fit residuals with for spectra where is weak and for objects with absorption‐line velocities is nonexistent.

Fig. 23.— The relation of the residuals to the fit to femtemp97 to the inverse of femtemp97. The filled dots are for spectra where , and the open circles are for spectra where .

Open New Window

These two plots, along with the Figures 10 and 12 in § 4, which show the statistical and systematic error models, demonstrate that for typical redshift survey spectra, residuals for absorption‐line objects are mainly determined by S/N; for emission‐line objects, the residuals are mainly systematic.

7. CONCLUSION

 

We have demonstrated the techniques that we have developed in the RVSAO suite to create a system for the accurate, automated reduction of spectra for galaxy redshifts and stellar radial velocities. More than half of all published redshifts have been measured using these techniques, as well as a large number of stellar radial velocities.

The correlation method for obtaining redshifts can be successfully extended from absorption‐line spectra to emission‐line spectra, with a substantial improvement in effectiveness over the previous method for obtaining emission‐line redshifts, automated line fitting. The reduction of emission‐line spectra requires different reduction steps than absorption‐line correlations. Emission‐line correlation redshifts are susceptible to blunders due to the presence of cosmic rays. However, with the use of automated line fitting (emsao) and absorption‐line correlation velocities, the blunder rate can be kept near zero, with the degree of automation kept high.

We have developed new techniques for calibrating and characterizing the blunder rate and the individual errors in redshift measurements. The blunder rate for RVSAO reductions can be kept near zero by the use of some simple heuristics to identify possible mistakes. For typical redshift survey data from the FAST spectrograph, the automation rate is 95%. Our self‐calibrating internal error estimator is accurate to ∼20%. Large, stable surveys enable development of more accurate and stable error estimators.

We have developed new methods for creating, calibrating, and using galaxy redshift templates. We have created an emission‐line template, femtemp97, having the median properties of a large set of strong emission‐line spectra. We have created an absorption‐line template, fabtemp97, having the mean properties of a large set of absorption‐line spectra showing no sign of emission. These spectra arise from physically distinct processes and can be used to form a pair of basis vectors to perform a two‐dimensional spectral classification. We have developed a new method for establishing the zero point for redshift observations, a method that minimizes the systematic differences between emission‐ and absorption‐line redshifts. This zero point is determined as accurately as we can establish the wavelength calibration using standard HeNeArFe lamps. We have shown a technique for measuring and eliminating differences between the instrumental zero point and the true zero point.

We have shown improved techniques for a number of the substeps necessary to obtain accurate redshifts, including (1) removal of emission (absorption) lines when correlating against an absorption‐ (emission‐) line template, (2) suppression of the night sky lines, (3) suppression of the continuum, (4) design of the Fourier filter, (5) iteration of the rest frame to maximize object‐template overlap, and (5) template‐dependent zero padding of spectra.

The rapid development of large‐aperture telescopes with multiobject spectrographs presents substantial challenges for redshift and radial velocity reductions. Reducing 1 or 2 orders of magnitude more spectra of objects 1 or 2 orders of magnitude fainter while maintaining high quality control standards and minimal personnel costs is clearly a difficult problem. RVSAO provides a solid methodological and software basis to meet these new challenges.

A large number of astronomers have sent us suggestions, complaints, bug reports, and kudos over the years. We should like to thank all of them and ask that this continue. We thank W. Press and B. Bromley for permitting us to quote from their work in advance of publication. We have benefited greatly from a number of detailed scientific and technical discussions with D. Fabricant, J. Huchra, D. Latham, and G. Torres. E. Falco carefully read the manuscript and made several suggestions to improve the clarity of the text, as did an anonymous referee. Susan Tokarz has reduced tens of thousands of spectra using RVSAO; her experience, patience, and friendly collaboration have been crucial to this project. Margaret Geller's support has enabled RVSAO to achieve its high degree of robust effectiveness.

APPENDIX
THE ELEMENTS OF RVSAO

 

The RVSAO package consists of six IRAF tasks: xcsao, emsao, linespec, sumspec, contpars, and bcvcorr. Each of these tasks is controlled by a set of user‐settable parameters. In this Appendix, we list and describe all the relevant parameters for these tasks and demonstrate their use in a series of examples.

A1. CROSS‐CORRELATING A SPECTRUM IN xcsao

 

Digital cross‐correlations in the RVSAO system are performed by xcsao, basically following the prescription of TD79, with a large number of refinements and additions. xcsao is capable of reducing a wide variety of input spectra, but it must have its many parameters properly set. Parameters for xcsao are set in its parameter list (Table 1), in the contpars program's parameter list (Table 15), and by special instructions to the program in the headers of the template spectra (Table 2). In this section, we examine in detail the elements of xcsao, and we describe how the parameters are used in the code, along with the algorithmic details of their use. We follow the cross‐correlation of a single spectrum against two very different templates, one showing primarily absorption features and one showing only emission features, to demonstrate the ways in which processing varies in response to the template spectra being used.

For each object spectrum file from the input list spectra and/or each aperture specified in the aperture list specnum, the fitting subroutine is called. Spectrum files are all read from the directory specified by specdir, but full or relative pathnames may be used in spectra, and specdir may be null. Spectra that are in flux units instead of counts should be renormalized by setting renormalize to “yes.” If obj_plot is “yes,” the object spectrum is plotted as in Figure 24, and the plot is kept on the screen available for zooming and editing until a “q” is typed. If fixbad is “yes,” regions specified in the file named by badlines are replaced by straight lines connecting the adjacent pixels, and the image is plotted again if obj_plot is “yes.”

Fig. 24.— Object spectrum showing both emission and absorption lines as plotted by xcsao if obj_plot is “yes.”

Open New Window

Band tempband of each template spectrum in the list templates and list of multispec apertures, tempnum, is loaded. Template files are all read from the directory specified by tempdir, but full or relative pathnames may be used in templates, and tempdir may be null. If echelle is “yes,” tempnum is ignored, and the multispec lines used for templates track those used for object spectra. If temp_plot is “yes,” the template spectrum is plotted.

A zero‐point redshift is computed by adding the solar system barycentric velocity correction, from a source specified by svel_corr, the redshift of the template from the VELOCITY parameter in the template header, an optional template‐dependent velocity shift from the TSHIFT parameter in the header, and an optional constant velocity shift from the tshift parameter. The template spectrum's barycentric velocity correction, from a source specified by tvel_corr, is subtracted because the template spectrum's observed velocity, not its corrected one, gives that spectrum the redshift which we are comparing. An initial redshift source may be specified by vel_init: if “guess,” this is from czguess; if “zero,” the initial velocity is 0; otherwise, it can be read from an object spectrum header parameter VELOCITY (“combination”), CZXC (“correlation”), or CZEM (“emission”). If such an initial redshift has been called for, or if this is the second pass or greater (nzpass > 1), the template log wavelength limits are shifted by that initial redshift (on the first pass) or the current correlation redshift. The wavelength region over which the template and object spectra overlap is computed. If the wavelength in angstroms specified by st_lambda is greater than the blue limit of the overlap region, it becomes the new limit. If the wavelength in angstroms specified by end_lambda is less than the red limit of the overlap region, it becomes the new limit. The two spectra are rebinned into log wavelength with a number of pixels, set by ncols using an interpolation mode specified by interp_mode.

First, the continuum, and, optionally, emission and/or emission lines, are removed from the rebinned object and template spectra. Parameters for fitting the continuum are in the IRAF pset task named contpars (see § A5). Emission and/or absorption lines may be removed from either or both of each object‐template pair. s_emchop controls whether lines will be removed from the object spectrum. If the template spectrum header parameter SUBCONT is present, its value overrides that of s_emchop. s_absrej and s_emrej set the lower and upper acceptable limits for object spectrum pixels to be used in the continuum fit. If pixels outside that range are rejected, as they are when the object spectrum we are following is to be correlated against an absorption‐line template, the rejected data points are plotted as in Figure 25 if contsub_plot is “yes.” A graph of the continuum‐removed data, as shown in Figure 26, is displayed if the contsub_plot parameter is “yes.”

Fig. 25.— Possible emission lines removed from the object spectrum before cross‐correlation against an absorption‐line template spectrum in xcsao.

Open New Window

Fig. 26.— The object spectrum with (a) its continuum subtracted and emission lines removed for cross‐correlation against an absorption‐line template and (b) its continuum subtracted but its emission lines kept for cross‐correlation against an emission‐line template spectrum in xcsao.

Open New Window

Template pixels with values outside of the lower t_absrej and upper t_emrej acceptable limits, in standard deviations from the continuum fit, are replaced by continuum values in the template spectrum if t_emchop is set to “yes.” The continuum is then subtracted from the template spectra.

Template and object spectra are then apodized, tapered at each end for the fraction bell_window of the entire spectrum. If zeropad is “yes,” both spectra are padded with an equal length of null (zero) spectrum. A graph of the continuum‐removed, apodized object spectrum is displayed, as shown in Figure 27, if the apodize_plot parameter is “yes.”

Fig. 27.— The object spectrum apodized at each end, ready to be Fourier transformed for cross‐correlation against (a) absorption and (b) emission templates in xcsao.

Open New Window

The spectra are then Fourier transformed. The Fourier power spectra of the object and template spectra are displayed as in Figures 28 and 29 if fft_plot is set to “yes.”

Fig. 28.— The object spectrum Fourier transformed over the wavelength ranges of the (a) absorption and (b) emission template spectra in xcsao.

Open New Window

Fig. 29.— Fourier transforms of (a) absorption and (b) emission template spectra. in xcsao.

Open New Window

The transformed spectra are then filtered with a cosine‐bell filter. The low frequencies are filtered from low_bin to top_low and the high frequencies are filtered from top_nrun to nrun. The template header parameter FI‐FLAG controls whether the template transform is filtered and whether the high‐frequency filter is turned off for both template and object transforms to leave in emission lines. If tfft_plot is set to “yes,” the filtered Fourier power spectra of the object and template spectra are displayed as in Figures 30 and 31.

Fig. 30.— The object spectrum transform filtered for correlation against (a) absorption and (b) emission templates in xcsao. Removing high frequencies adversely affects the shape of narrow emission lines, so no high‐frequency filtering is done for emission‐line correlations.

Open New Window

Fig. 31.— Filtered Fourier transforms of (a) absorption and (b) emission templates in xcsao. No high frequencies are filtered out of the emission template transform.

Open New Window

The filtered transforms are then cross‐correlated and normalized. If uxcor_plot is “yes,” the unfiltered correlation is displayed. If xcor_plot is “yes,” this result is displayed as in Figure 32, and a specific peak may be selected using the cursor if curmode is “yes.” In that case, the maximum value within pksrch pixels of the cursor position is used. Otherwise, the highest correlation peak between the velocities minvel and maxvel is used. The redshift is calculated by fitting a parabola or similar function specified by pkmode to the portion of the peak above pkfrac of the maximum value of that peak. The R‐value and error are computed, and control returns to XCFIT to set up the template for the next pass.

Fig. 32.— Normalized, filtered cross correlations of the Fourier transforms of the object spectrum against the Fourier transforms of (a) absorption and (b) emission templates in xcsao.

Open New Window

After all of the template spectra have been correlated against an object spectrum, the template with the highest R‐value is selected. The results are displayed as text to the devices specified by logfiles in the format specified by report_mode. Figure 33 shows the default report.

Fig. 33.— xcsao report of results of cross‐correlation of the object spectrum transforms against transforms of absorption and emission template spectra

Open New Window

If displot is “yes,” the object spectrum, and, optionally, the selected correlation peak, are plotted to device in the format specified by dispmode. Figure 34 shows the summary graph if dispmode is 1. Figure 35 shows the dispmode 2 summary graph, with emission and absorption lines labeled. If hardcopy is “yes,” the same graph is sent a printer.

Fig. 34.— xcsao summary display for emission‐line cross‐correlation

Open New Window

Fig. 35.— xcsao summary display with labeled lines for emission‐line cross‐correlation

Open New Window

If nsmooth > 0, the object spectrum is smoothed by a 1‐2‐1 sliding filter nsmooth times for display purposes only. This smoothing may be changed interactively using the “g” command in cursor mode. The filtered cross‐correlation with the best R‐value is displayed centered on the redshift cvel (in km s−1) with a width in km s−1 of dvel. If cvel is INDEF, the fit redshift is used; if dvel is INDEF, the width is set to 20 times the peak width. If the correlation is not displayed, absorption lines (ablines = yes) and/or emission lines (emlines = yes) may be labeled from line lists in the directory linedir, as shown in Figure 35.

If curmode is “yes,” the user can interact with the display using the terminal cursor to zoom in on portions of the spectrum, rerun the cross‐correlation, change the display format, edit the spectrum, or several other functions. For example, Figure 36 shows the correlation result for the second‐best template, selected using the “T” cursor command.

Fig. 36.— xcsao summary display for absorption cross‐correlation

Open New Window

If save_vel is “yes,” cross‐correlation redshift results are written into the object spectrum image header in a form appropriate to the spectrum format: two entries plus one per template if multispec; otherwise, one value per keyword, as in Table 3.

A2. FITTING REDSHIFTED LINES IN A SPECTRUM IN emsao

 

While xcsao can now use emission‐line templates, it is still useful to measure redshifts of emission‐line spectra directly. For large surveys, the interactive determination of line centers and calculation of redshifts by a program like IRAF's splot is simply too slow. emsao, a companion to the cross‐correlation task xcsao, was written to find emission lines automatically, compute redshifts for each identified line, and combine them into a single radial velocity. The results may be graphically displayed or printed, saved to a file, and/or stored in the spectrum file header. emsao is designed to run with minimal human intervention, but options may be set to allow manual improvement of the line identifications and resulting redshift. The graphic cursor may be used to change fit and display parameters. Table 4 shows the full parameter list for emsao.

For each object spectrum file from the input list of spectra, and/or each aperture specified in the aperture list, the fitting subroutine, EMFIT, is called for the specified spectrum image band. Spectrum files are all read from the specified directory unless a full pathname is given in the spectrum list. Relative pathnames may be used for spectra. If a directory is not set, the spectra are expected to reside in the current working directory.

After the spectrum is loaded, it is renormalized, if the renormalize flag is “yes.” This should usually be done if the spectrum is in flux units. If the fixbad flag is “yes,” regions specified in the file named by badlines are replaced by straight lines. The spectrum is then smoothed nsmooth times. If obj_plot is “yes,” the spectrum is plotted and the plot, as shown in Figure 37, is kept on the screen available for zooming and editing until a “q” is typed.

Fig. 37.— emsao spectrum display if obj_plot is “yes”

Open New Window

To use a sky spectrum to compute noise statistics and to improve the error computations, an aperture or band must be specified. If skynum is not zero, a sky spectrum is read from that multispec aperture in the same file as specnum. If skyband is not zero, a sky spectrum is read from that multispec band in the same file. The sky spectrum, used to get the noise for error computations, is plotted if obj_plot is “yes.” Previous results that have been saved in the spectrum image header may be used by setting linefit to “no”; in that case, all of the fitting below is skipped and the results are displayed.

The continuum, computed by the IRAF curve‐fitting subroutine that is driven by the parameters set in the contpars pset task, is subtracted from the spectrum. If contsub_plot is “yes,” the spectrum is plotted with the continuum removed, as in Figure 38.

Fig. 38.— The smoothed spectrum with the continuum subtracted is displayed by emsao if contsub_plot is “yes.”

Open New Window

The source of the initial redshift is specified by vel_init. If it is “guess.” the starting redshift is read from the parameter czguess. If it is “search,” one line in the spectrum is identified by the program using the table specified by the emsearch parameter, as shown in Table 5. This table lists line centers in angstroms and the wavelength range over which each one should be the strongest line. It can be modified by the user to match the data. The brightest line in each region is assumed to be the one in the table, and its observed wavelength is saved. The redshift of the brightest of those lines is returned as an initial value to be refined by looking at more lines.

If vel_init is “combination,” the initial redshift velocity is read from the spectrum header parameter VELOCITY; if “correlation,” it is read from the spectrum header parameter CZXC; and if “emission,” from CZEM.

Search regions are read for each line in the file specified by the emlines parameter. Table 6 shows the contents of such a file. Each region is then shifted by the guessed redshift and expanded in each direction by wspan angstroms. All emission lines within each specified wavelength region are found. A spectrum pixel is assumed to be a line center if the pixel value is the maximum of the npfit neighbors on either side and greater than linesig times the square root of the average counts in those pixels. A second‐order fit is then made to the (2 ∗ npfit) + 1 points centered on the peak to refine the center and peak height. The brightest line in each region is kept unless it has already been identified. Order matters: the brightest line in a region should be listed first, so that if it is the only one present in overlapping regions, it is correctly named.

Before line profiles are fit, a copy of the spectrum is smoothed esmooth times using the same smoothing algorithm as is used before the line search is conducted. The parameter esmooth should be left at zero unless the data are especially noisy. It is best never to go above 2. Because local continuum values are important to the line fit, the continuum is removed from this copy of the spectrum using the same parameters as were used before conducting the line search. If contsub_plot is “yes,” the spectrum with the continuum removed is plotted as shown in Figure 39.

Fig. 39.— Unsmoothed spectrum with continuum removed, displayed by emsao if it contsub_plot is “yes.”

Open New Window

Each identified line is checked to see if it is part of one of the groups of close emission lines listed in the file specified by emcombine, such as that shown in Table 7. If it is, all lines in the combination will be simultaneously fit. Those members of a line combination that are not found are initialized at the redshift of the most recently found line of the group. The missing line heights are assumed proportional to that line according to the relative heights in the emcombine file. The lines are fit by one to three Gaussians, along with an optional local continuum level using zero to three additional coefficients as set by nlcont. Redshift (computed from the wavelength of the center pixel coordinate), width, height, and errors are returned for each line.

Each emission line is checked to see whether (1) it has a fit center greater than zero, (2) it is wider than lwmin times the mean width of all of the identified lines that are found and (3) narrower than lwmax times that mean width, (4) it has an equivalent width more than lsmin times the error in the equivalent width, (5) it is too close to the blue edge of the spectrum, (6) it is too close to the red edge of the spectrum, or (7) the error in the center of the Gaussian is zero. If dispmode is 1 or report_mode is 1, wavelengths and redshifts are printed for each line, with an “X” followed by the code for the test that failed at the end of entries that were omitted from the fit. If a line has been successfully fit, the rejection can be overridden interactively (if curmode is “yes”) using the “+” and “−” commands in cursor mode from the final display. If a line has been added or subtracted in cursor mode, a “+” or “−” at the end of the entry indicates that fact. A mean velocity is computed, weighted by the square of the error in the line centers, and returned. If only a single line is found, the error is set to sigline, which should be set to the uncertainty in the dispersion function in angstroms. This is usually significantly greater that the error in the fit to the center of the Gaussian.

After all of the lines have been fit and a combined velocity has been computed, a zero‐point redshift is computed by adding the solar system barycentric velocity correction, from a source specified by svel_corr.

The results are displayed as text to the devices specified by logfiles in the format specified by report_mode. Options include (1) one line per found emission line under a self‐documenting summary, (2) a single line report listing the names of the lines that are found, and (3) a single line report listing a velocity for each searched‐for line. Table 8 shows the mode 1 report for the spectrum we are following.

If displot is “yes,” the spectrum is plotted to device in the format specified by dispmode. Figure 40 shows the mode 1 report; Figure 41 shows the mode 2 report. If hardcopy is “yes,” the same graph is automatically sent to a printer as well.

Fig. 40.— The emsao summary display shows the spectrum and line fit results if dispmode is 1

Open New Window

Fig. 41.— The emsao summary display shows only the spectrum if dispmode is 2.

Open New Window

If nsmooth is greater than zero, the displayed spectrum is smoothed by a 1‐2‐1 sliding filter that many times. Absorption lines listed in the file ablines are labeled if dispabs is “yes,” and emission lines listed in the file emlines are labeled if dispem is “yes.” Both files are found in the directory linedir.

If curmode is “yes,” the user can interact with the display using the terminal cursor to zoom in on portions of the spectrum, identify lines and refit the emission lines, change the display format, edit the spectrum, or perform several other functions.

If save_vel is “yes,” emission‐line redshift results are written into the spectrum image header in a form appropriate to the spectrum format: two entries plus one per line if multispec; otherwise, one value per keyword, as shown in Table 9.

A3. CREATING A SPECTRUM WITH linespec

 

linespec is an IRAF task for making a spectrum from a list of emission and/or absorption lines that is driven by the parameters shown in Table 10 and uses the same line models, and hence the same line‐defining parameters, as emsao. The IRAF task artdata.mk1dspec allows more complicated line and continuum models, but linespec eases tuning of a line spectrum to a particular spectrograph. It was written to create templates for use by xcsao. An example of such a line list, with center wavelengths, widths, and heights for each line that may appear in the spectrum, is shown in Table 11. The line list is read from the file specified by the linefile parameter in the directory designated by the parameter, linedir. If linedir is null, the file is assumed to be in the current working directory.

A blank (all zero) spectrum, with the object name specobj and file name specfile in the directory specdir is created. It is linear in wavelength, with a resolution in angstroms given by pix_lambda and a range from st_lambda to end_lambda. Spectral world coordinate system information is written to the header, with the standard FITS parameters CRPIX1, CRVAL1, and CDELT1. If verbose is “yes,” files specified by logfiles are opened and a header is written.

The center of each line in the table is redshifted according to the parameters zspec and velspec. velspec (cz = apparent Doppler‐shifting velocity) is used unless zspec (Δλ/λ) is not zero, in which case zspec is used for the redshift. The linewidth, if it is tabulated in km s−1, is converted to angstroms at the shifted line center. The line width is also broadened appropriately if the line is redshifted. The redshift velocity is put in the spectrum header using the VELOCITY keyword. If maxwidth is “yes,” the linewidth parameter is used as the width of the line if it is greater than the tabulated width. If maxwidth is no, the width from the table is used, and the linewidth smoothing is done later. For each line, a Gaussian at the shifted center wavelength, half‐width, and tabulated height is added to the spectrum. As each line is added to the spectrum, the line definition is written to the spectrum's header, producing a table such as that shown in Table 12. After all of the lines are computed, a continuum level specified by the parameter continuum is added to the spectrum.

The spectrum is plotted to a graphics device, as shown in Figure 42, if spec_plot is “yes.” After the graph is displayed to the terminal, if spec_int is “yes,” the user must interact with it using single‐character cursor commands, such as the “z” for zoom command, the result of which is shown in Figure 43. If the “@” command is given in this mode, the graph is immediately sent to the device specified by plotter. The “q” command must be typed to leave the graph and proceed.

Fig. 42.— Graph of an emission‐line spectrum produced by linespec

Open New Window

Fig. 43.— Graph of Hα region of the same spectrum displayed using the “z” cursor command in linespec.

Open New Window

If maxwidth is no, the resolution of the spectrograph is simulated by convolving the entire spectrum with a Gaussian of height 1.0 and sigma (half‐width at half of maximum) of linewidth. If spec_plot is “yes,” the spectrum is displayed again.

The linespec version and the date the program is being run are written to the spectrum image header, the spectrum is written, and the file is closed.

A4. CREATING A COMPOSITE SPECTRUM WITH sumspec

 

sumspec is an IRAF task for making a composite spectrum that is adjusted for velocity. It rebins input spectra to a common wavelength or log wavelength range and resolution, with optional continuum removal and normalization. While similar in some ways to the IRAF tasks onedspec.sarith and onedspec.scombine, it includes use of known velocity components in aligning spectra at specified redshifts and other features, such as line rejection, which are useful in making templates. In addition to making templates, sumspec has turned out to be very useful for rebinning individual spectra as well as combining multiple spectra. The parameter list for sumspec is shown in Table 13.

The spectra to be added are specified by the list of files in the parameter spectra and/or the list of spectrum numbers in specnum. Unless a pathname is specified as part of the file name in spectra, each spectrum is expected to be found in the directory specdir. All spectra are taken from the band specband in multispec spectra. Files specified by logfiles are opened to receive logging information. If debug is “yes,” additional information about the input files and the progress of the program is written to STDERR.

If either st_lambda or end_lambda is set to INDEF, all of the spectrum headers are read to find the limits of the overlap of their wavelength coverage, and the missing limit(s) of the output spectrum are set accordingly.

Before the first spectrum in the list is added, a blank (all zero) spectrum, with the object name compobj and file name compfile in the directory compdir is created. If complog is “yes,” it is linear in log wavelength; otherwise, it is linear in wavelength. Either way, the output spectrum contains npts data points between st_lambda, or its log, and end_lambda, or its log. World coordinate system information is then written to the header.

Spectra are read from the input list one at a time and renormalized if renormalize is set to “yes.” This should be done if the spectra are fluxed or whenever pixel values are much bigger or smaller than their variation. If spec_plot is “yes,” each spectrum is plotted to the device specified by the parameter device, as in Figure 44. If spec_int is “yes,” the user can interact with that graph using single character cursor commands and must type a “q” to proceed with the task. If the “@” command is given in this mode, the graph is sent to the device specified by plotter.

Fig. 44.— Input spectra for sumspec

Open New Window

Each spectrum is rebinned using an interpolation mode specified by interp_mode into an npts pixel wavelength‐linear spectrum covering the range computed above. If complog is “yes,” the rebinned spectrum is linear in log wavelength; if complog is “no,” the rebinned spectrum is linear in wavelength. If specplot is “yes,” the rebinned spectrum also is plotted.

If s_contin is set to “subtract” or “divide” instead of “no,” the IRAF curve‐fitting subroutines are used to fit a continuum to each input spectrum. The CONTSUM task, described in the following section, sets all of the appropriate parameters.

As it is rebinned, each input spectrum is redshifted according to the parameters ztemp ( ), if it is not zero, or veltemp (cz = apparent Doppler‐shifting velocity). A velocity correction to the solar system barycenter is removed according to the svel_corr parameter. Set it to none if the input spectra have not been shifted. If “file,” BCV is used if present in the file header, or else HCV. If “hfile,” the header parameter HCV is always used. If neither is found, no correction is made. If “heliocentric” or “barycentric” corrections are chosen, position and time parameters are read from the spectrum data file header, and the correction is computed as described in § A6 below.

Emission and/or absorption lines may be removed from each input spectrum if reject is “yes.” If the spectrum header parameter SUBCONT is present, its value overrides that of reject. The parameters absrej and emrej set the lower and upper acceptable limits for input spectrum pixels in standard deviations of the continuum fit to the spectrum.

Graphs of the continuum‐removed, apodized data are displayed, as in Figure 45, if the cont_plot parameter is “yes.”

Fig. 45.— Rebinned sumspec input spectra with the continuum removed

Open New Window

The composite spectrum is plotted to the device specified by the parameter device, as in Figure 46, if comp_plot is “yes.” If comp_int is “yes,” the display is held to allow the user to interact with it using single character cursor commands. If the “@” command is given in this mode, the graph is sent to the device specified by plotter. A “q” command allows the task to proceed.

Fig. 46.— Output composite of two input spectra plotted by sumspec

Open New Window

After each spectrum is added, the sumspec version and the current date are written to the output spectrum image header and that spectrum image file is written. After the last spectrum is added, the file is closed. Table 14 shows the parameters added to the spectrum header. EXPTIME is the total exposure time from all of the input images. A one‐dimensional image is produced, so DISPAXIS is always 1. DC‐FLAG is 1 if the output spectrum is log wavelength, in which case the log wavelength of the first pixel is given in both CRVAL1 and W0, and the log wavelength per pixel is given in CDELT1 and WPC. The VELOCITY is set by the velcomp or zcomp parameter. If no velocity is specified, it is left at 0. If savenames is “yes,” the instrument, filename, input redshift velocity, and barycentric velocity correction are written to the header.

A5. REMOVING A SPECTRUM'S CONTINUUM USING contpars OR contsum

 

contpars and CONTSUM are IRAF parameter files that are used by a single subroutine to fit and remove a continuum from object or template spectra. The subroutine uses the IRAF interactive curve‐fitting subroutine package and is based on the technique used in the IRAF RV package (Fitzpatrick 1993). The parameters that are set are shown in Table 15. The same subroutine is called prior to cross‐correlation by xcsao, an emission‐line search by emsao, and before or after summation of spectra by sumspec. If the SUBCONT header keyword in a template spectrum is set to “F,” no continuum is fit or removed from that template spectrum. If the keyword is “T” or not present, the continuum is fit and removed.

For a spectrum, such as that in Figure 47, the continuum is fit using the IRAF curve‐fitting subroutine. It may be used interactively by setting the c_interactive parameter to “yes.” In that case, a graph of the fit and the spectrum is displayed as shown in Figures 48 and 49. The result may be examined or refit after individual pixels are removed. A “q” command is needed to exit from the interactive graph.

Fig. 47.— Input spectrum before the continuum is removed according to contpars

Open New Window

Fig. 48.— On this display shown during interactive continuum fitting, spectrum pixels rejected from the fit are marked with diamonds.

Open New Window

Fig. 49.— Here is a closeup of the region around Hα, showing the continuum that has been fitted and giving a better view of the rejected pixels.

Open New Window

The c_function parameter specifies the function to be used to fit the continuum. If the continuum is irregular, it is best to use the spline3, but spline1, Legendre, or Chebyshev polynomials may be used as well. The parameter order sets the order of the fit. The portion of the spectrum to be fit is specified, in image section format, by the parameter sample. Before the fit, the spectrum is averaged in groups of naverage pixels.

Pixels with values outside of the limits s_low_reject and s_high_reject for spectra and t_low_reject and t_high_reject for templates are rejected. Separate limits are used because composite template spectra will have much higher S/N levels than object spectra and limits may need to be set differently to fit a good continuum. These limits are specified in standard deviations of the fit. To make sure that the wings of lines are thoroughly eliminated, this rejection process is repeated niterate times. Each time, a number of pixels, specified by the grow parameter, are rejected on each side of the rejected pixel.

Normally, the continuum which is fit is subtracted from the spectrum, preserving the signal‐to‐noise ratio of a spectrum when it is dominated by photon noise. To remove the continuum of both template and object spectra by division, the template spectrum header keyword DIVCONT is set to “T.” In that case, noise in portions of the spectrum with low S/N levels will be amplified. If emsao.contsub_plot, xcsao.contsub_plot, or sumspec.cont_plot is “yes,” the spectrum is plotted with the continuum removed, as shown in Figure 50.

Fig. 50.— Fig. 47 spectrum with its continuum removed

Open New Window

A6. CORRECTING RADIAL VELOCITIES TO THE SOLAR SYSTEM BARYCENTER

 

To compare redshift velocities observed when the Earth is at different positions in its orbit, the velocity of the Sun relative to the Earth is added to the redshift velocity. This heliocentric velocity correction is more accurately done to the true stationary point of the solar system barycenter, its center of mass. This barycentric velocity correction is shown in Figure 51. A second, smaller correction is added for the motion of the observer relative to the center of the Earth as it rotates. This vector velocity is projected in the direction of the observed object, and that component is saved.

Fig. 51.— The correction of a radial velocity observed from the Earth to that which would be observed from the center of mass of the solar system is shown at two different times in the Earth's orbit.

Open New Window

The xcsao, emsao, and sumspec tasks of the RVSAO package compute this correction using subroutines that read the time of observation, object position, and observatory position from the spectrum image header. Several common alternative keywords are built into those subroutines. RA, DEC, and EPOCH give the right ascension, declination, and equinox of the observed object. SITELONG, SITELAT, and SITEELEV for longitude, latitude, and altitude, respectively, give the observer's location on the earth. OBS‐DATE yields the observation date and UTMID the midtime of the observation. If UTMID is not found in the header, UT, assumed to be the end time of the observation, or UTOPEN, the start time of the observation, are used in combination with EXPOSURE or EXPTIME, which give the duration of the observation in seconds to compute a midtime. The requested heliocentric or barycentric correction is then calculated at that midtime.

Since there is really no standard for the meaning of these keywords, a separate task, bcvcorr, has been added to RVSAO to allow several alternate ways of specifying these three major pieces of information. bcvcorr can write its result to the header of the image that it is processing; the other RVSAO tasks will use this value when their svel_corr parameter is set to “file.” First, the sky direction of the object is set. The right ascension, declination, and equinox of the object's sky position are read from the keywords specified by the keyra, keydec, and keyeqnx parameters. If those are not all found, the position is read directly from the parameters ra, dec, and equinox.

The time for the velocity correction is read from a Julian date keyword specified by keyjd. If none is found, the date is taken from the header keyword specified by keydate. The UT midtime is taken from the keyword in keymid. If that is not found, the observation midtime is computed by finding the exposure duration from the keyexp keyword value and adding half of it to the start time from the keystart keyword value or subtracting half of it from the end time in the keyend keyword value. If no time can be found in the header, the Julian date is read directly from the gjd parameter, if it is greater than zero, or the midtime is taken from the year, month, day, and ut parameters.

The heliocentric Julian Date (HJD) is the time at which the light from the object for this observation reaches the Sun. It is needed if multiple radial velocity observations of an object are to be compared accurately. The HJD is not used in computation of the velocity correction, which is dependent on the time of the observation at the Earth. It may be read from the header keyword keyhjd or set by the parameter hjd (only if greater than 0). If none is set, the HJD is computed from the observation time.

If the obsname parameter is “file,” the observatory name and position is read from the image header using the keywords in keyobs, for the name; keylat, for the latitude; keylong, for the longitude; and keyalt, for the altitude. Otherwise, the value of obsname is used to get a position from IRAF's observatory database. If the string is not found there, the longitude, latitude, and altitude are given by the parameters obslong, obslat, and obsalt.

If verbose is “yes,” the information in Table 16 is printed to standard output. If savebcv is “yes,” several spectrum image header keywords are set. BCV is the barycentric velocity correction in km s−1 (see Table 17). The midtime of the observation is stored three ways: UTMID, the universal time; GJDN, the geocentric Julian Date, and HJDN, the heliocentric Julian Date. These times are displayed by xcsao and emsao as the time of observation.

TABLE 1
TABLE 1 Parameter List for xcsao

Open New Window

TABLE 1
TABLE 1 (Continued)

Open New Window

TABLE 2
TABLE 2 Template Spectrum Header Parameters

Open New Window

TABLE 3
TABLE 3 Object Spectrum Header Parameters Added by xcsao

Open New Window

TABLE 4
TABLE 4 Parameter List for emsao

Open New Window

TABLE 5
TABLE 5 Emission Lines for emsao Initial Redshift Estimate

Open New Window

TABLE 6
TABLE 6 Emission Lines in Galaxy Spectra for emsao

Open New Window

TABLE 7
TABLE 7 Emission‐Line Combinations for Galaxy Spectra for emsao

Open New Window

TABLE 8
TABLE 8 emsao Report for Emission‐Line Galaxy

Open New Window

TABLE 9
TABLE 9 Keywords to Save emsao Results

Open New Window

TABLE 10
TABLE 10 Parameter List for linespec

Open New Window

TABLE 11
TABLE 11 List of Emission Lines from linespec Input File

Open New Window

TABLE 12
TABLE 12 Multiple‐valued Keywords

Open New Window

TABLE 13
TABLE 13 Parameter List for sumspec

Open New Window

TABLE 14
TABLE 14 Parameters Set in sumspec Output Spectrum

Open New Window

TABLE 15
TABLE 15 contpars Parameters as Shown by IRAF

Open New Window

TABLE 16
TABLE 16 bcvcorr Parameters as Shown by IRAF

Open New Window

TABLE 17
TABLE 17 Output from bcvcorr

Open New Window

REFERENCES

 
© 1998. The Astronomical Society of the Pacific. All rights reserved. Printed in U.S.A.