Unbiased flux calibration methods for spectral-line radio observations
Key Words.:Methods: observational – Techniques: spectroscopic
iint \savesymboliiint \savesymboliiiint \savesymbolidotsint \restoresymbolAMSiint \restoresymbolAMSiiint \restoresymbolAMSiiiint \restoresymbolAMSidotsint
Position and frequency switching techniques used for the removal of the bandpass dependence of radio astronomical spectra are presented and discussed in detail. Both methods are widely used, although the frequency dependence of the system temperature and/or noise diode is often neglected. This leads to systematic errors in the calibration that potentially have a significant impact on scientific results, especially when using large-bandwidth receivers or performing statistical analyses. We present methods to derive an unbiased calibration using a noise diode, which is part of many heterodyne receivers. We compare the proposed methods and describe the advantages and bottlenecks of the various approaches. Monte Carlo simulations are used to qualitatively investigate both systematics and the error distribution of the reconstructed flux estimates about the correct flux values for the new methods but also the ‘classical’ case. Finally, the determination of the frequency-dependent noise temperature of the calibration diode using hot–cold measurements or observations of well-known continuum sources is also briefly discussed.
Spectroscopic data obtained from radio astronomical observations typically need to be post-processed in several steps before one can proceed with the scientific analysis. The uncalibrated flux values not only have to be converted from device units (counts) to physical quantities (e.g. Kelvin or Jy/beam), but one must also deal with the frequency-dependent gain (bandpass curve), radio frequency interference, removal of any residual baselines, and so on. While several methods are discussed in the literature to solve the bandpass issue (Heiles, 2007, and references therein), the flux calibration is often handled using only a very simplified approach. In particular, the widely used position and frequency switching techniques (see Heiles, 2007, for a review) suffer from an improper treatment of calibration. In this paper we describe the underlying problem, present possible solutions, and assess the quality of the various approaches in terms of flux calibration accuracy. We compare the proposed methods with the ‘standard’ approach, and show that the latter is subject to bias effects that can have a large impact on scientific results.
The paper is organised as follows. Section 2 introduces the basic equations and definitions necessary to describe the problem. Sections 3 and 4 present the position and frequency switching procedures, along with several methods to correctly calibrate the measured spectra. An alternative approach is shown in Section 5 which may be of interest in cases where neither of the switching schemes are applicable. To illustrate the calibration methods, artificial spectra were generated with well-defined input quantities. Section 6 discusses the various approaches and their ability to reconstruct the simulated flux values. The flux calibration presented here relies on the signal of a noise diode fed into the receiver. Two possibilities for measuring the spectrum of this calibration signal are discussed in Section 7. We conclude with a summary in Section 8.
2 Basic equations and definitions
The output signal, , of a spectroscopic (heterodyne) receiving system can be described with the simple formula
where the telescope itself and the receiving system exhibit a frequency-dependent gain, , on the incident radiation. Often referred to as , , is the antenna temperature caused by the flux density of the astronomical source of interest, which potentially incorporates line emission and a continuum contribution, i.e. . Apart from the flux of the astronomical source, , that we are interested in, there are several noise contributions, which we subsume in the so-called system temperature, .111Note that sometimes the system temperature is defined differently to include . This consists of
Standing waves (SW) are produced when the incident radio signal is reflected at the aperture of the antenna leading to interference patterns that often exhibit a sinusoidal shape. For a more detailed discussion, we refer the reader to Rohlfs & Wilson (2004).
In this paper we make extensive use of the calibration signal , which is a part of many heterodyne radio receivers. Its purpose is to provide a well-defined intensity standard that can be used to calibrate the measured spectra in terms of the antenna temperature. In practice, there are two common techniques to make use of . One is to rapidly switch the signal on and off on time-scales of the order of tens of milliseconds up to a few seconds throughout the observation, leading to two so-called switching phases cal and non-cal. This method is used at the 100-m Effelsberg telescope, for example. In the following, we adopt the notation and when referring to non-cal and cal phases respectively, whilst indicates either of the two. The alternative method in use slowly cycles between the two phases during a dedicated measurement of a calibration source (or blank sky) every few hundred seconds, e.g. at the Arecibo telescope. In this paper, we assume a rapid switching for the sake of a simpler presentation of the calibration methods. However, our work may easily be extended to the slow cycling scheme. We would also like to point out that many (sub-)mm-wave telescopes make use of chopper wheels or hot-cold loads to improve the calibration. It should again be straightforward to adapt our methods to these techniques, though this is beyond the scope of this paper.
A spectroscopic back-end usually records the spectral density function, , in arbitrary units (e.g. counts) such that the raw spectra need to be calibrated in terms of both flux density values and also with respect to the frequency-dependent bandpass shape . We note that for simplification we incorporate the coefficient that converts the antenna temperature to counts into the gain factor . Using a heterodyne receiver, one has to distinguish between the gain in the receiver/radio frequency (RF) and in the intermediate frequency (IF) part of the receiving system.
One major problem in the reduction of spectroscopic data is that the bandpass curve(s), , must be disentangled from the input signal, . From Eq. (1), it is clear that one has to divide the measured signal, , by the bandpass, , to obtain the input temperature spectra. Unfortunately, the bandpass curve is unknown222Even if it had been measured, it is still subject to instabilities, e.g. caused by temperature drifts, such that frequent re-evaluation would be needed.. To cope with this issue, one often applies a switching technique to obtain a reference spectrum, e.g. position or frequency switching as explained below. Dividing by this reference removes the frequency-dependent gain. Residual baselines caused by SW, for example, may then be subtracted. We note that one can neither bandpass-calibrate a measurement by simply subtracting a reference spectrum nor remove residual baselines by division of a baseline model.
Even after the (successful) removal of the gain curve, the extracted spectrum is not equal to the true brightness temperature, , of the astronomical object. Both the atmospheric dampening and aperture efficiency, , which act on prior to amplification, must also be taken into account
where is the zenith opacity and the elevation-dependent airmass. 333Note that does not accurately describe the airmass at low elevation angles and must be replaced by a site-dependent model, including vertical values of density and refraction index (e.g., Maddalena & Johnson, 2005). The aperture efficiency, , is the product of several other efficiencies. The antenna pattern and blockage of the aperture play a role, as well as the surface accuracy and to some extent ohmic losses. This cannot be calculated theoretically, especially for the large single-dish telescopes. It is possible to model opacity (see Kraus, 2009). In the best case, an atmospheric model makes use of simultaneous measurements of water vapour to achieve greater accuracy (e.g. using a water vapour radiometer, Rottmann & Roy, 2007). For the remainder of this paper, we neglect these corrections for simplification as they are completely independent of the calibration methods presented here.
3 Position switching
In position switching, one points the telescope first to a reference position (Off) before measuring toward the true target (On). Applying Eq. (1), we obtain
The aim of position switching is to remove the bandpass curve. In the ideal case, is equal for both positions and dividing both spectra leads to
This can also be written as
where we introduce as
the difference between the system temperatures in the On and Off positions. The receiving system is usually stable for the duration of a single position switch and only a few contributors must be taken into account. While and are elevation dependent, may be a function of the incident continuum flux. In some cases, SW would only be noticeable if one pointed toward strong continuum sources, e.g. calibrators, and in other cases during daytime because of the sun, regardless of the observed source. The Galactic background continuum radiation depends slightly on the pointing position of the telescope. Using an appropriate Off position, for instance at constant elevation angle, one should usually be able to minimise . The Off position should obviously not contain a continuum source itself. The equations could otherwise still be interpreted in a way where is not the true continuum flux of the On position. Re-arranging Eq. (7), we obtain
It is important to realise that is a frequency-dependent function. One can ensure a proper flux calibration over the complete spectrum only if this is determined correctly. In the remaining part of this section, we develop two methods to determine . The first is easier to apply and more robust. The second may be of interest under particular circumstances when using frequency switching and is also discussed in the case of position switching for completeness. Both approaches make use of the signal of the calibration diode (or noise tube). Furthermore, we discuss the ‘classical’ method, which assumes that and , as a simplification of the first method.
3.1 Setting up simulations
In the following subsections, we present the three alternative calibration methods. Artificial spectra were generated to enable a direct comparison of input and output spectral line intensities. In Section 6, we statistically assess the quality of these methods to evaluate their general robustness and error distribution.
The test spectrum uses a relatively simple set-up. The continuum flux of the Off position is described as a single-slope power law , where . This is a strong simplification of reality where several different power laws (e.g. from the galactic background, sky, and ground radiation) add to a more or less ‘arbitrary’ temperature contribution from the receiver. For the astronomical source, , a set of three Gaussians superposed on the power law (where ) was used. Each Gaussian has equivalent line parameters (3 K amplitude, 1.4 MHz line width, FWHM) but are added at different frequencies. These are used as ‘probes’ to measure the flux calibration quality at different frequencies after data reduction. We chose to model a possible frequency dependence of the noise diode using a shallow power law (where ). In reality is usually a much more complicated function especially at higher frequencies, but this setup is suitable for the main purpose of our simulations, which is to illustrate the various steps in the data processing as simply as possible. Section 3.6 illustrates a more realistic test case with respect to the functional form of the system temperature and calibration diode. The resulting temperature spectra are plotted in Fig. 1 (upper panel).
Gaussian noise was added to each component according to the root mean square (RMS) noise calculated using the radiometer equation (see Rohlfs & Wilson, 2004)
where is a constant factor depending on, e.g. the input quantisation of the signals, is the bandwidth of a spectral channel, and is the integration time. In the simulations, we assume the total number of spectral channels to be 16k over a bandwidth of and an integration time (per spectral dump or measurement phase) of . For example, a system temperature of 25 K would result in (per spectral channel). Eq. (11) implies that noise is a function of frequency because the temperature values, , for each measurement phase vary over the observed band. Likewise, the additional temperature components of the continuum emission of the observed source and the calibration diode increase the noise. It is also possible to calculate the theoretical RMS level expected in the final (reduced) spectrum using the noise values of the individual input spectra (see Appendix A).
Finally, for a more complicated test case, we simulated (mono-modal) SW by adding a sine wave component to (see Fig. 1, lower panel).444Note that this set-up does not account for all possible situations. For example, if strong continuum sources alone produce SW, then they contribute to only. However, this is only a minor problem as the computation of is not affected. Nevertheless, one has to account for such a case if one uses continuum sources for calibration. For the position switching test, .
The input temperature spectra were multiplied with a simple bandpass curve shown in Fig. 2. The resulting power spectral densities (in arbitrary units) are plotted in Fig. 3. The superposition of the SW and the ripples in the bandpass shown in the lower panel of Fig. 3 creates a rather complex pattern.
3.2 Using the Off position to obtain
leads to an equation to infer the system temperature of the Off position, which depends only on the temperature of the noise diode, . If the latter is known, then
However, the noise diode temperature is a frequency-dependent quantity. It can be measured with good precision using the hot–cold method (compare also Section 7.1 and Appendix B.1), which, unfortunately, is a time-consuming procedure. One solution is to establish a catalogue of astronomical calibrators, i.e. bright continuum sources, which serve as reference to (re-)calibrate the spectrum on appropriate timescales (see Section 7.2 and Appendix B.2). This approach was also proposed by Maddalena & Johnson (2005). Since is a time-dependent quantity (e.g. owing to a change in the environmental conditions), it is always a potential source of (systematic) error. Nevertheless, a temporal stability on the order of 1% can typically be expected on the timescale of one hour.
Unfortunately, even if one has a good model of , the measured quantity is still subject to noise and will substantially increase the noise in the final reduced spectra. One solution is to suppress the noise in the obtained spectrum before substituting it into Eq. (9).
In simpler cases, namely in the absence of a standing wave contribution and a relatively flat , the quantity
is approximately proportional to , which can be described by a power law. Hence, a low-order polynomial might already suffice to describe the quantity .
Unfortunately, at higher frequencies in particular, it is almost impossible to engineer a sufficiently flat . In such cases, a filtering approach or the use of high-order polynomials is required to suppress noise (see Section 3.6 for a more realistic example). In most cases, it is easier to model for numerical stability, as usually . The noise-free (or low-noise) model can then safely be inverted to obtain .
Fig. 4 illustrates for the two test cases with and without SW. In the former case (upper panel), a third-order polynomial was used to describe . For the latter, a more complicated model is necessary
where are polynomial functions of degree . This result was obtained only after providing suitable initial fit parameter values and is shown in the lower panel. An entirely automated procedure would most likely have difficulties in handling SW — even for the simplest scenario of monochromatic SW.
We note that the continuum flux of the source is also reconstructed. To obtain an optimal signal-to-noise ratio (S/N), the reduced cal and non-cal spectra should be averaged.
After applying Eq. (16) to the simulated data (using a polynomial model for ; see Fig. 4), the spectra shown in Fig. 5 are obtained. To measure the spectral line intensities, we fit a Gaussian superposed on a third-order polynomial to each spectral line. We note that Fig. 5 shows the reduced spectra for the SW case. The non-SW scenario produces very similar results. The upper panel shows the direct result and the lower panel contains the spectrum after baseline removal. The noise values and the recovered line intensities match the expected results. The residual noise is higher than one might naively expect (the RMS in the final spectrum is around the same level as for each of the four contributing phases). This is due to the division by the reference spectrum that does not contain signal but adds noise (see also Appendix A.1).
It is possible to smooth the reference spectrum prior to division greatly reducing the RMS in the final spectrum by almost a factor of . However, as Braatz (2009) points out, this procedure can have unwanted side-effects, for example, degradation of the baseline or emphasis of narrow features present in the reference, e.g. radio frequency interference (RFI), and results in (locally) correlated noise.
One drawback of the proposed method is that the Off position alone is utilised to infer . As a consequence, a factor of in sensitivity is lost in calculating . However, the final result is unaffected so long as is described by a (smooth) model. If sensitivity is not a major concern, even the original (noisy) spectrum, i.e., without modelling, may be used. This is advantageous since the solution is then unaffected by the choice of a given model.
3.3 Using On and Off positions to obtain
The second method utilises Eq. (7) and uses continuum contributions only, to define the quantities
Since the On spectra are involved, it is important that the fitting algorithm excludes spectral lines. In simple cases, the use of relatively low-order polynomial functions should suffice to provide a good model. If SW occur and only few modes (i.e. SW frequencies) are present, one should be able to model these with superposed sine-waves. In principle, one could even try to use a filtering approach to obtain a noise-free model .
The fits calculated for the example spectrum are shown in Fig. 6. For the non-SW case, a third-order polynomial model is used. The SW case uses
where are polynomial functions of degree .
such that is a function of . Eq. (9) now becomes
Using the calculated model fits and applying Eq. (22), one obtains the reduced spectrum. Fig. 7 shows the standing wave example where the upper panel shows the reconstructed spectrum, , and the lower panel contains the result after baseline subtraction. The continuum contribution, , of the source is correctly reproduced. In some cases however, the baseline may contain a residual imprint caused by standing waves, which affects the quality of the model fitting to some extent, especially in the presence of noise. In the worst case, this might also affect the flux calibration. The more complex the required model, e.g. in the presence of SW, the larger the potential error in flux calibration.
We note that this method may not be applicable in all cases. Inspecting Eq. (22) in detail reveals that the denominator of the term
can become zero. This occurs if , i.e. for sources with low continuum flux. In these cases, a nearby calibrator may be used for an independent reference position in order to obtain a well-behaved correction term.
3.4 Determination of neglecting any dependence on frequency — the ‘classical’ approach
The ‘classical’ case where and are both treated as constants is clearly a simplified case of the first method (presented in Section 3.2) and of course causes errors in the flux calibration. This simplification was justifiable to some extent until about two decades ago when broadband receiving systems started to become more commonly used. Today, this approach is still widely used, e.g. within the online reduction pipeline at the 100-m telescope in Effelsberg (see Kraus, 2009), in order to provide a fast and robust online display of measured spectra.
Eq. (12) may be rewritten as
In principle, this provides a direct means of calculating from . Unfortunately, a simple evaluation of the above equation is numerically unstable because the denominator can have values close to zero or even become negative owing to noise. If we treat and as a constant with respect to time and frequency, we find that
where is obtained by calculating the average of the inner 50% of the difference spectrum in the online reduction pipeline at the 100-m telescope at Effelsberg. The current data reduction pipeline at the GBT (GBTIDL; Braatz, 2009)555http://gbtidl.nrao.edu/ uses the equivalent method except using the inner 80% of the spectra. For future versions of GBTIDL, a vectorised approach is planned whereby averaging is performed on certain bins.
The parameter is calculated by the same averaging procedure, such that
Using a scalar value for , roughly known for each receiver from Eq. (25), we obtain the ‘calibrated’ spectrum
Both, cal and non-cal spectra, should again be added to reduce noise in the final spectrum.
Although this approach is much easier to implement in software, the reader should be warned that neglecting the frequency dependence of introduces systematic errors, which can be a serious problem. This is especially the case when computing line ratios or dealing with statistical analyses of large samples of sources at various radial velocities/frequencies (see Section 6).
Figure 8 (upper panels) shows the resulting spectra after applying Eq. (27). We note that the baseline does not match the input baseline at all. For reference, we subtracted a baseline (which proved to be rather complicated in the SW case) from the resulting spectra. The result clearly shows a systematic flux calibration error and is displayed in Fig. 8 (lower panels).
3.5 Computing the bandpass curve
From Eq. (5), it directly follows that
where can be calculated, for example, via Eq. (14), , or using Eq. (20) and (21). This can be of great practical use. While the inferred from a single measurement (see Fig. 9) is rather noisy, one may use several observations, even of different sources, to average the gain curve and suppress noise. The calculated could then in turn be applied to all involved data sets. In practice, the gain is usually imperfectly stable in time. Eq. (28) enables us to monitor the time evolution of the receiving system. It is even likely that temporal drifts can be approximated to some degree with a low-order polynomial (spectral channel-wise) and can then be used to improve the accuracy of the data reduction.
3.6 A realistic case
We present the results of a more complex simulation that resembles realistic observations. In comparison to the previous example, both and now have more structure. Both quantities were modelled to imitate the outcome of a measurement with the 1.9-cm primary-focus receiver at the 100-m telescope. The resulting spectrum is illustrated in Fig. 10. shows a steep increase towards higher frequency, attributed to the water vapour line at (during the observations, weather conditions were rather bad).
In this instance, the first method as described in Section 3.2 is applied using spectral filtering instead of polynomial models to describe . The result is shown in Fig. 11. The second approach requires spectral windows to be set around the (expected) emission lines such that a simple filtering approach would not work. It is therefore neglected here.
Fig. 12 displays the resulting spectrum. Despite the very complex set-up, the calibration worked well and the underlying continuum flux of the observed source could be reconstructed. Fig. 13 shows the result of the classic approach. For purposes of comparison, a constant continuum contribution was chosen. The resulting baseline is very complex and is in fact the incorrectly calibrated continuum flux of the source. In a worst case scenario, this can lead to a misinterpretation of the results depending on the line widths of the expected astronomical features. For example, Fig. 13 suggests that the middle spectral line has a blue-shifted wing.
4 Frequency switching
In frequency switching, one uses two different local oscillator (LO) frequencies to provide a reference spectrum. The shifted spectrum serves as a reference to remove the bandpass dependence in the non-shifted spectrum. In practice, a symmetric LO shifting pattern is often used, such that Eq. (1) may be written as
Together with the cal and non-cal phases, the LO switching leads to a total of four different phases of sig+cal, ref, sig, and ref+cal, which are usually rapidly stepped through in order to avoid time-dependent instabilities.
Using a shorter notation and assuming , which should be fulfilled since switching is performed on short timescales
We first assume that , which is fulfilled only in the very rare cases that . Then
where . We note that is equal in the cal and non-cal phases, while is different ( in the non-cal phase and in the cal phase). Both quantities depend only on the slope of and because all contributors to remain approximately equal during two adjacent switching phases. Using Eq. (33), it follows that
To obtain an optimal S/N, one should not only average cal and non-cal phases but also shift (to match RF frequencies) and average and
This is clearly a much more complicated procedure than for the position switching case.
4.1 Bandpass ghosts
In Eq. (36), apart from the desired signal, , so-called ‘spectral-line ghosts’ (sometimes also referred to as bandpass ghosts), , appear in the resulting equation. They have half of the amplitude of the original signal. Furthermore, there is a residual additive666This is an important point, as it allows for a simple subtraction of a baseline fit, without destroying the calibration. baseline depending on the relative changes in and with frequency. These can make the final baseline fitting a challenging task, especially in the presence of standing waves (SW).
One thing to note is that usually we are unable to determine and can only compute which then is substituted into the right-hand side of Eq. (34) and (35). This is an essential drawback. Consequently, the estimated are wrong for all frequencies where is not fulfilled (see Fig. 14).
Luckily, this usually only happens for frequencies at which the spectral-line ghosts appear (except for complicated cases; see Section 4.6). If one uses both Eq. (34) and (35) and applies shift-and-averaging as described above, this is not a problem. However, if just one of these two equations is used and instead a so-called ‘fold’ procedure (e.g., the Class task fold from the Gildas package777http://www.iram.fr/IRAMFR/GILDAS) is applied, one will achieve incorrect results if (Fig. 14, bottom panel). This is not a rare case, especially for observations in the decimetre-wavelength regime. The term folding denotes the use of spectral-line ghosts themselves to add to the positive emission line, by changing the sign of and both shifting and averaging the signals
4.2 Setting up simulations
The input spectra for frequency switching are similar to those shown in Fig. 1, except that each measurement phase contains the contribution. They are also shifted according to an LO frequency of before multiplying by the bandpass (). A different standing wave frequency, , is also used for the method that we describe in Section 4.4 (since must be a multiple of for this second method). Fig. 15 shows the resulting input spectra as they would have been measured by the backend. We note that for frequency switching the measured sig and ref phases belong to different RF frequencies. To clarify whether an operation is performed in either the IF or RF domain, we assume a hypothetical IF centre frequency of 150 MHz and plot figures accordingly.
4.3 Using the two switching phases individually to obtain
As in the previous sections, the signal of the calibration diode is used. Computing
One needs to neglect the line emission contribution, , on the right hand side of the equations, but this has no serious impact, except that spectral-line ghosts have unexpected amplitudes.
Fig. 16 displays the spectra, along with appropriate fitting models. Using these models leads to the calibrated spectra shown in Fig. 17 (top panels). Whilst the more simple case without SW shows no residual baseline owing to the specific form of Eq. (36) in the resulting spectrum, the SW case displays a strong residual pattern. For the lower right panel in Fig. 17, we subtracted a sine-wave signal modulated with a polynomial.
4.4 Using both switching phases together to obtain