Abstract
Comparison of quasar (QSO) absorption spectra with laboratory spectra allows us to probe possible variations in the fundamental constants over cosmological timescales. In a companion paper we present an analysis of Keck/HIRES spectra and report possible evidence suggesting that the finestructure constant, α, may have been smaller in the past: over the redshift range . In this paper we describe a comprehensive investigation into possible systematic effects. Most of these do not significantly influence our results. When we correct for those which do produce a significant systematic effect in the data, the deviation of from zero becomes more significant. We are led increasingly to the interpretation that α was slightly smaller in the past.
Introduction
The idea of varying fundamental constants is not a new one: Milne (1935, 1937) and Dirac (1937) independently suggested a varying Newtonian gravitational constant. More modern theories, such as superstring theory, naturally predict variations in fundamental coupling constants (e.g. Forgács & Horváth 1979; Marciano 1984; Barrow 1987; Li & Gott 1998), strongly motivating an experimental search. Several authors have demonstrated the idea (Varshalovich & Potekhin 1995 and references therein), first used by Bahcall, Sargent & Schmidt (1967), of using spectra of highredshift gas clouds seen in absorption against background quasars (QSOs) to constrain possible variation of the finestructure constant, , over cosmological timescales. This method is based on a comparison of the observed transition wavelengths of alkali doublets and recent applications of it have yielded constraints on , where α_{z} and α_{0} are the values of α at the absorption cloud redshift z and in the laboratory respectively (e.g. Cowie & Songaila 1995; Ivanchik, Potekhin & Varshalovich 1999; Varshalovich, Potekhin & Ivanchik 2000). Recently, we have increased the level of precision obtained with the alkalidoublet method to using higherquality data (Murphy et al. 2001b).
However, a new method offering an order of magnitude improvement in precision was suggested by Dzuba, Flambaum & Webb (1999a,b) and was demonstrated by Webb et al. (1999, hereafter W99). The increase in precision derives from the fact that the relativistic corrections to the energy levels of different ions vary from species to species by up to an order of magnitude. A comparison of the observed Mgi, Mgii and Feii transition wavelengths in 30 absorption systems allowed W99 to suggest tentatively that the mean value of over the redshift range : the finestructure constant may have been smaller in the past.
The present work is the companion paper of Murphy et al. (2001a, hereafter M01a) in which we discuss our techniques in detail and present the results of recent work. We summarize the main points of that paper here. We obtained Keck/HIRES spectra of 28 QSOs and analysed 49 absorption systems lying along their lines of sight. These systems included those systems considered in W99 and also a higherredshift sample of damped Lymanα systems (DLAs) containing absorption lines of ions such as Mgi, Mgii, Alii, Aliii, Siii, Crii, Feii, Niii and Znii, so that our redshift range now covers , corresponding to lookback times of billion years . Our absorption systems fell conveniently into two subsamples: the systems comprised the lowredshift sample and the DLAs comprised the highredshift sample . There were, however, two absorbers within the highz sample that lay at low redshift but we have included them in the highz sample since the observations and data reduction for the two subsamples were carried out independently by two different groups.
The analysis is based on the following equation for the wavenumber, ω_{z}, of a given transition at a redshift z:
where ω_{0} is the laboratory rest wavenumber of the transition, q_{1} and q_{2} are relativistic coefficients, and . Note that the second and third terms on the righthand side only contribute if . The values of q_{1} and q_{2} have been calculated using manybody techniques in Dzuba et al. (1999a,b, 2001) and we tabulate those values relevant to our analysis in table 1 of M01a.The values of ω_{0} must be known to a high precision in order to place stringent constraints on . A precision of corresponds to an uncertainty on (∼0.3kms^{−1}).^{1} This corresponds fairly well to the data quality since we can centroid an absorption feature to approximately onetenth of the HIRES resolution, ∼7kms^{−1}. Such precision has not been previously achieved for most resonance transitions of ionized species, and so measurements of ω_{0} for all relevant species were carried out by several groups (Nave et al. 1991; Pickering, Thorne & Webb 1998; Pickering et al. 2000; Griesmann & Kling 2000) using laboratory Fourier transform spectrometers. We have summarized and tabulated these measurements in M01a.
Knowing all values of ω_{0}, q_{1} and q_{2}, we then deconvolve each absorption system into individual Voigt components and perform a nonlinear leastsquares χ^{2} minimization to find the bestfitting value of at each absorption redshift. For the lowz sample we confirm the W99 result, our new results giving a weighted mean of . For the new, highz sample, we find a similar result: . Both results independently suggest a smaller value of α in the past. Taking the sample as a whole, we find , a 4σ result. These new results are also summarized in Webb et al. (2001, hereafter W01).
Given the potential importance of such a result for fundamental physics, a thorough investigation of possible systematic problems in the data or in the analysis is essential, and this is the purpose of the present paper. In Section 2 we list all effects that one might consider as being able to mimic a systematically nonzero . We are able to exclude many of them with general arguments. Sections are then given over to a more detailed analysis of five specific effects: wavelength miscalibration, ionic line blending, atmospheric dispersion effects, differential isotopic saturation and possible isotopic variation. Section 8 describes another test that is sensitive to a simple monotonic wavelength distortion, specifically for the highz sample. We make our conclusions in Section 9.
Summary of Potential Systematic Errors
In M01a and W01 we found that values of were systematically negative and that this trend was evident in two different samples of data occupying different redshift regimes. We note that the lowz sample is more susceptible to systematic errors than the highz sample. The Mgii lines act as anchors against which the larger Feii shifts can be measured. All of the Feii lines have similar magnitude shifts in the same direction (i.e. q_{1} is large and positive in all cases). Furthermore, all the Feii lines lie to the blue of the Mgii doublet. If some systematic error in the wavelength scale is present, which effectively compresses the spectrum, this would mimic a negative shift in the mean value of (if it were not degenerate with redshift). The highz sample is far more complex in this regard since it contains ions not only with positive coefficients, but also with negative q_{1} coefficients (i.e. all Crii lines and Niiiλλ1741 and 1751). Also, the shifts for the various ions are of various magnitudes, making it more complicated to predict the effect of a systematic error on our measured values of .
Laboratory wavelength errors
Errors in the values of ω_{0} will lead directly to errors in any single determination of . If, for example, the laboratory wavelengths of the Mgii lines were all shifted to the blue or red, then this would systematically bias the values of for the lowz sample. In M01a we summarize the measurements of all values of ω_{0}. Most measurements were repeated independently and good agreement was obtained in those cases. A typical measurement accuracy from these laboratory measurements ω_{0} of is ∼, which could only introduce a maximum error in around an order of magnitude below the observed deviation from zero. In practice, since we use several transitions (particularly for the highz sample), any associated errors are likely to be smaller. We therefore consider it unlikely that the systematic shift in is due to errors in the values of ω_{0}.
As also noted in M01a, some laboratory spectra were calibrated with the Arii lines of Norlén (1973) while others were calibrated with those of Whaling et al. (1995). The Whaling wavenumbers are systematically larger than those of Norlén such that To deal with this discrepancy we have rescaled those wavenumbers which had been calibrated using the Whaling scale (Alii, Aliii, Siii) to the Norlén scale. We stress that the difference between the Norlén and Whaling wavenumbers is a scaling only, and not a distortion. Therefore, is it not possible to introduce an apparent deviation in from zero, since small zeropoint offsets in the wavelength scale are degenerate with redshift.
Heliocentric velocity variation
During a 1h exposure, representative of the quasar integrations for the sample we have used, the heliocentric velocity may change by as much as ∼0.1kms^{−1}. This will act to smear any spectral feature: the instantaneous spectrum is convolved with a tophat function in velocity space. However, since we fit a redshift parameter to the spectrum when determining , then a velocity space smearing is completely absorbed into the redshift parameter. Heliocentric smearing will have no effect on measurements of .
Differential isotopic saturation
The ions used in our analysis, with the exception of those of Al, have naturally occurring isotopes. Thus, each absorption line will be a composite of absorption lines from all the isotopes of a given species. In all cases other than Mg and Si, the laboratory wavenumbers we have used correspond to composite values. The composite wavenumbers will only strictly be applicable in the optically thin regime (linear part of the curve of growth). As the column density increases, the strongest isotopic component begins to saturate and the line centroid will shift according to the natural abundances (cf. Section 2.4) of the other isotopes.
The mass isotopic shift, being the dominant splitting mechanism for light atoms, is proportional to for m the atomic mass. The wavelength separation between the isotopes is therefore largest for the lighter (low q_{1}) ions. Isotopic separations have been measured for the Mgiλ2853 transition by Hallstadius (1979) and for Mgiiλ2796 by Drullinger, Wineland & Bergquist (1980), and the isotopic structure of all Mg lines is presented in Pickering et al. (1998). However, similar measurements do not exist for the other transitions of interest. For Si we used estimates of the isotopic wavenumbers based on a scaling from the basic Mg structure by the mass shift. For all other species, measured (composite) wavenumbers were used. See M01a for full details.
We discuss and quantify this effect in Section 6.
Isotopic abundance variation
In the above we have assumed that the isotopic abundances in the absorption systems are equal to terrestrial values (Table 3). However, if the isotopic abundances at high redshift are significantly different, then this may lead to apparent shifts in the absorptionline centroids and this would cause a systematic error in . Again, the most susceptible transitions are those of Mg and Si since (i) these ions will have the largest isotopic separations (i.e. m is low) and (ii) transitions in these atoms are prominent in our data. If, for the lowz sample, the ^{25}Mg and ^{26}Mg abundances were zero, we would incorrectly infer a more positive by assuming terrestrial isotopic abundances. In the simple case of finding from Mgiiλ2796 and Feiiλ2344, we would expect to find in the absence of any real variation in α.
We consider this possibility in detail in Section 7.
Hyperfinestructure effects
Hyperfine splitting of energy levels occurs in species with odd proton or neutron numbers. Different hyperfine components will have different transition probabilities but the composite wavelength of a line will be unchanged by the splitting (i.e. the centre of gravity of the hyperfine components is constant with increased splitting). However, a similar differential saturation effect will occur for the hyperfine components as discussed in Section 2.3 for the isotopic components.
The worst case here is that of ^{27}Al – the only stable isotope of Al – since it has an odd proton number and a very large magnetic moment. In the laboratory experiments of Griesmann & Kling (2000), the hyperfine structure of the Aliii doublet was clearly resolved and can be seen to be quite prominent. However, as noted in M01a, we only include Aliii lines in the analysis of six absorption systems of our highz sample, and in these cases we see very little or no saturation. The effect of removing the individual Aliii transitions from those systems on can be seen in Fig. 4. In the case of Aliiλ1670, both the ground and excited states have zero spin and the pwave hyperfine splitting in the excited state is very small (cf. ^{25}Mg in Pickering et al. 1998). Thus, we do not expect any Al hyperfine saturation effects on .
All other dominant isotopes in our analysis have even proton numbers, and where odd isotopes occur these have very low abundances (see Table 3). The ^{25}Mg hyperfine structure has already been taken into account in our analysis but the hyperfine separations are so small that neglecting it would make no difference to even upon saturation. The odd isotopes of Si, Cr, Fe, Zn and Ni have very low abundances. Also, Si, Cr, Fe and Ni will have very small hyperfine splittings due either to low magnetic moments or to the nature of the ground and excitedstate wavefunctions.
We therefore expect that differential hyperfinestructure saturation has not affected in our analysis. An unlikely interpretation, given the above, is cosmological evolution of the hyperfinestructure constants, but very large relative variations would be required.
Another possibility is that the populations of the two hyperfine levels of Aliii are not equal in the absorption clouds. Since the gas clouds have very low density, the equilibrium between the two hyperfine levels will be maintained predominantly by interaction with cosmic microwave background (CMB) photons. For the Aliiiλ1862 transition, this implies a lower limit on the relative populations of the two levels of , leading to a shift in the line centroid of ∼0.01cm^{−1} (∼. This is a potential problem for the six systems in which we fit Aliii lines. We do not consider this effect in detail, but we note that Fig. 4 in Section 4.2 shows that potential lack of thermal equilibrium has no significant effect on .
Magnetic fields
Shifting of ionic energy levels due to large magnetic fields may also be a possibility. If very largescale magnetic fields exist in the intergalactic medium and quasar light is polarized, then this could result in correlated values of in neighbouring regions of the Universe. However, the QSOs studied in M01a and W01 cover a large region of the sky. Furthermore, even in Abell clusters, intracluster magnetic field strengths are typically ∼μG (e.g. Feretti et al. 1999), roughly nine orders of magnitude below the strength required to cause substantial effects. We consider this possibility to be very unlikely.
Kinematic effects
An assumption in the analysis of M01a and W01 is that the absorption velocity structure is the same in different species. Of course, if there were some nonzero velocity between, say, the Mgii and Feii components of a given cloud, then we would incorrectly determine . However, it is difficult to imagine a mechanism by which the Mgii is systematically blue or redshifted with respect to the Feii when averaged over a large sample of absorption systems. Suppose, for example, we have an asymmetric, expanding (or contracting), rotating cloud in which all the Feii and Mgii lie at different ‘radii’. Any large sample of sightlines through such a cloud would yield an average value of . Furthermore, if this was a major effect, then the observed scatter in the points would be greater than expected on the basis of our error estimates, and in M01a we showed this not to be the case. No matter how contrived the velocity structure and ionic segregation, when averaging over a large number of values of , it is not feasible to generate a systematic offset from zero.
Wavelength miscalibration
There may be inaccuracies in the wavelength calibration software used to reduce the spectra. Also, human error cannot be ruled out. W99 stated that a drift in the wavelength calibration by an amount corresponding to roughly ∼2.5 times the mean wavelength–pixel residuals over the range of the chargecoupled device (CCD) could result in a significant error in the lowz results. A considerably more complicated form of the miscalibration would be required to produce a systematic error in the highz results.
The wavelength calibration of the QSO CCD images is done by comparison with thorium–argon (ThAr) lamp exposures taken before and after the QSO frames. We consider errors in the laboratory wavelengths of the selected ThAr lines to be negligible, as their relative values are known to much greater precision than our values of ω_{0} (Palmer & Engleman 1983). However, ThAr line misidentifications may lead to serious miscalibration of the wavelength scale over large wavelength regions. Such misidentifications can also be applied to the rest of the spectra taken over the same series of observations (i.e. applied to the spectra of many QSOs). Thus, a priori, we cannot rule out the possibility that the wavelength scale has been set improperly in this process and that a systematic shift in α has not been mimicked. Such a potential effect needs careful investigation.
In Section 3 we discuss and quantify this effect for our data set.
Airvacuum wavelength conversion
Most ThAr line lists are presented as air wavelengths. The usual data reduction procedure is to carry out the wavelength calibration fits to the data using air wavelengths and then convert the calibrated spectrum to vacuum wavelengths. For example, iraf uses air wavelengths for a selection of lines from the thorium spectral atlas of Palmer & Engleman (1983) and from the argon lines of Norlén (1973).^{2} The same list is used in makee– the HIRES data reduction package written by T. Barlow. makee was used to calibrate 12 of our highz sample and the remainder were calibrated within iraf. The latter were actually calibrated with a set of vacuum wavelengths, λ_{vac}, that iraf states are derived from the air wavelengths, λ_{air}, by using the Edlén (1966) formula for the refractive index, n,
at 15°C and atmospheric pressure. Instead of using , iraf makes the approximation . makee converts the final air wavelength−calibrated QSO frames to vacuum using the Cauchy dispersion formula (Weast 1979), at 15°C and atmospheric pressure. λ_{vac} and λ_{air} are both measured in angstroms in the above equations. This difference begs the question: ‘What is the absolute and relative accuracy of these conversion formula?’First let us refine this question. The measurements reported by Palmer & Engleman (1983) and Norlén (1973) were all carried out in vacuum. The air wavelengths they quote are calculated using the Edlén (1966) formula. Thus, iraf and makee will only produce ‘correct’ vacuum wavelength scales if they invert the Edlén formula to calculate the refractive index. As stated above, this is not the case. Therefore, to restate the above question: ‘What error does the iraf and makee air–vacuum wavelength conversion introduce into the resulting QSO wavelength scale?’ To illustrate the answer to this we plot the various conversion formulae commonly used in the literature, including those discussed above, in Fig. 1. We have also included the original Edlén (1953) formula for comparative purposes.^{3}Peck & Reeder (1972) published a fit to revised refractive index data, proposing a fourparameter formula:
From Fig. 1 we see that the Cauchy formula seriously deviates from the other formulae. Furthermore, the difference between the Edlén and Cauchy dispersions is strongly wavelengthdependent. Thus, spectra calibrated with makee may have systematically nonzero values of as a result of this distortion. As we stated in M01a, we have corrected the 12 affected spectra in our (highz) sample by converting the wavelength scale back to air wavelengths using the (inverted) Cauchy formula and then back to vacuum using the (inverted) Edlén formula. From Fig. 1 we can also see that there is only a small distortion to the wavelength scale introduced due to the approximation made within iraf. This difference is below our level of precision.
Temperature changes during observations
The refractive index of air within the spectrograph depends on temperature (and also on pressure, but this is a smaller effect). If a QSO spectrum is calibrated with only a single ThAr exposure taken at a different spectrograph temperature, we expect a systematic miscalibration of the QSO frame's wavelength scale. Consider a QSO frame taken at 15°C that is calibrated with a ThAr exposure taken at 0°C and similar pressure. Any wavelength separation in the QSO exposure will be overestimated. For example, the difference between the wavelength separations of vacuum λλ4000 and 7000 at 15°C and 0°C is (using the tables in Weast 1979). A comparison with values of Δλ in table 1 of M01a shows that this would result in a significantly positive for the lowz sample, the situation not being so clear for the highz sample.
Note that the above effect can only mimic a systematically nonzero if the spectrograph temperature is systematically higher or lower for the ThAr frames compared with the QSO frames. This is possible if the ThAr exposures were always taken before or after the QSO frames and the temperature evolves monotonically throughout the night. However, the effect is greatly reduced if ThAr exposures are taken near the time of the QSO observations. This was the case in our observations. We have used image header information to calculate the QSO–ThAr temperature difference, , inside HIRES for both the low and highz samples. Here the average of T_{QSO} and T_{ThAr} is taken over all exposures for each object. We find mean values of and for the low and highz samples respectively. Taking into account that, on average, restframe separations between transitions are ≲300Å then we can see that temperature variations could not have mimicked any significant shift in in either sample.
Line blending
The errors in presented in M01a and W01 take into account errors from signaltonoise ratio (S/N) and spectral resolution considerations and the velocity structure of the profile fits. The errors are also reduced when more lines are fitted simultaneously. However, we have assumed that we have deconvolved each absorption system into the correct number of velocity components. There may have been weak, interloping, unresolved lines, which, if the interloping species were in the same absorption cloud, could have produced a shift in the fitted line wavelengths of all velocity components of one or more transitions.
We distinguish between random blends and systematic blends. Random blends may occur if many absorption clouds at different redshifts intersect the line of sight to a single QSO. A systematic blend will occur when two species are in the same cloud and have absorption lines with similar rest wavelengths. Such an effect could mimic a systematic shift in α.
The importance of such an effect diminishes as the number of transitions used in each fit is increased. Therefore, the effect, if present, would be expected to be smaller in our highz sample.
In Section 4 we describe the results of a detailed search for atomic (and not molecular) interlopers for all transitions of interest, and place constraints on their strengths and positions.
Atmospheric dispersion effects
As noted in M01a, all of the lowz sample and many of the highredshift sample of QSOs were observed before 1996 August, at which time an image rotator was fitted to HIRES. Objects observed prior to 1996 August were therefore not observed with the slit length (spatial direction) perpendicular to the horizon. Atmospheric dispersion leads to a stretching of the target spectrum relative to the calibration spectrum.
An additional consequence is the potential wavelengthdependent truncation of the QSO seeing profile on the slit jaw edges. This can introduce a wavelengthdependent asymmetry in the point spread function (PSF).
Together, these effects can conspire to mimic a nonzero . We provide detailed calculations for these effects in Section 5 and show that they will tend to produce an apparent positive, and so cannot mimic our results.
Instrumental profile variations
If the instrumental profile (IP) of HIRES shows significant intrinsic asymmetries then we expect absorptionline centroids to be incorrectly estimated. If the asymmetry varies with wavelength, this could mimic a nonzero . Valenti, Butler & Marcy (1995) have determined the HIRES IP for several positions along a single order and they find that the IP is indeed asymmetric and that the asymmetry varies slightly along the echelle order. Valenti et al. (1995) do not quantify the asymmetry variation across orders, although it is likely to be comparable to the variation along the orders.
As a byproduct of the detailed investigation into possible wavelength calibration errors, we find that this effect is negligible (see Section 3 for details).
We have therefore eliminated, with some reliability, potential systematic errors due to: laboratory wavelength errors, heliocentric velocity variation, hyperfinestructure effects, magnetic fields, kinematic effects, air–vacuum wavelength conversion errors and temperature changes during observations.
In the following sections we investigate in considerable detail the remaining effects: wavelength miscalibration (and instrumental profile variations), line blending, atmospheric dispersion effects, differential isotopic saturation and isotopic abundance variation. We also consider a simple test for any simple but unidentified systematic errors for the highz sample in Section 8.
Detailed analysis: wavelength calibration errors
In this section, we investigate the possibility that wavelength calibration errors could have led to an apparent nonzero . To quantify this directly, we analyse sets of ThAr emission lines in the calibration lamp spectra in the same way as each set of QSO absorption lines has been analysed. If no calibration error is present, then we expect that for all clouds (i.e. for all z).
Method
For each QSO spectrum, we selected ∼15Å sections of ThAr spectra at wavelengths corresponding to the observed absorption lines in the QSO spectra. Each ThAr ∼15Å section contains several (typically up to ∼10) lines. Error arrays were generated assuming Poisson counting statistics. We fitted continua to these ThAr sections, dividing the raw spectra by the continua to obtain normalized spectra.
The main steps in the procedure are then as follows. We select several independent sets of ThAr lines.^{4} Each set of ThAr lines therefore corresponds to, and has been selected in an analogous way to, one QSO absorption system.
In order to obtain the best estimate of any wavelength calibration errors, the average is determined for each ThAr set (typically sets for each QSO absorption system). The process of defining sets in this way also provides a direct additional estimate of the errors on .
To fit the ThAr spectra, a modified version of vpfit was used (to fit Gaussian profiles, see M01a), adopting ThAr laboratory wavelengths from Palmer & Engleman (1983) for values of ω_{0} in equation (1). The q_{1} and q_{2} coefficients for the QSO absorption lines are applied to the corresponding ThAr lines. That is, we treat the set of ThAr lines as if they were the QSO lines. The free parameters involved in each fit were linewidth (assumed the same for all ThAr lines in a given set, although the results were insensitive to this assumption), ‘redshift’ (numerically close to zero) and peak intensity.
Results and discussion
We applied the above method to each available ThAr spectrum corresponding to each QSO absorption cloud considered and the results are illustrated in Fig. 2. Our first observation is that the scatter in these results is an order of magnitude less than that of the QSO absorption results themselves. This clearly demonstrates that wavelength miscalibration of the CCD has not mimicked a shift in α to any significant degree. The results above also show that the data quality permit a precision of .
Secondly, many points in Fig. 2 deviate significantly from zero. There are three possible reasons for this: weak blended emission lines, ThAr line misidentifications and errors in the ThAr laboratory wavelengths. Inspecting the ThAr spectra does indeed reveal typically ∼20 easily identifiable, but very weak, lines per 15Å, and we would expect this to add to the scatter in the final results for . Furthermore, this effect would not produce correlated deviations in , as observed. Consistent with this interpretation, we found that the value of χ^{2} per degree of freedom for a fit to a given set of ThAr lines was ≫1. Note that, if ThAr line misidentifications had occurred, one may expect systematically correlated deviations. Finally, we do not expect any significant errors in the ThAr laboratory wavelengths (as discussed in Section 2.8).
On inspection of some of the ThAr spectra, we noted that some ThAr lines towards the blue end of the spectrum were slightly asymmetric. This may indicate a degradation in the polynomial fits near the edges of the fitting regions. There may also be intrinsic asymmetries in the IP (Section 2.13). If such asymmetry variations exist, then they should also apply to the QSO spectra. However, the results in Fig. 2 show that such variations do not produce a systematic effect in at any detectable level.
We therefore conclude this section by stating that the results in Fig. 2 show unambiguously that wavelength calibration errors and variations in IP asymmetries are not responsible for the observed shifts in α in the QSO absorptionline results.
Detailed analysis: systematic line blending with unknown species
In this section we suppose that there is an unidentified transition, arising from some species in the same gas cloud as is being studied. We further assume that this interloping transition is not strong enough and/or sufficiently displaced from one of the lines in the analysis to be detected directly. We describe two approaches to explore this: in Section 4.1 we attempt to identify candidate blending species and in Section 4.2 we describe the results of a test where we remove one transition or one species at a time to investigate the impact on .
Search for potential interlopers
There are two main parameters that will determine whether or not any particular interloper can give rise to a significant effect in terms of column density and the position in wavelength with respect to the ‘host’ line of interest. In order to search for a potential interloper, we explore the range in parameter space through numerical simulation, generating trial interlopers, blended with a single stronger ‘host’ line, and examining the resulting shift in the combined centroid in terms of . Having established the allowed characteristics of a potential interloper, we use photoionization equilibrium models to search for candidate atomic species. We ignore molecular species.
Characteristics of a candidate interloper
We restrict the following discussion to blending with Mgiiλ2796, since is particularly sensitive to blending with this anchor line in the lowz sample. We later generalize our discussion to include blending with all other relevant transitions.
Consider an interloper – a spectral line of a blending species X – which is separated from the fitted line centroid wavelength, λ_{0}, by −Δλ in the frame of the cloud. Let the separation between the actual Mgiiλ2796 wavelength and the fitted line position be dλ (i.e. in the frame of the cloud. That is, the line is fitted as a single component, despite the presence of the interloper. The interloper has a column density N(X). However, as the species X is not known a priori, we may only determine the quantity N(X)f_{X} where f_{X} is the oscillator strength of the interloping transition.
We generated a synthetic optically thin Mgiiλ2796 line and matched this against a composite Mgiiλ2796 profile, which included an additional weak blended component, representing our candidate interloper. The model was fitted to the ‘data’ using vpfit. The synthetic spectrum was generated using a column density of and a b parameter of 5.0kms^{−1}, as these were representative values for the Mgii lines in M01a and W01 (Churchill 1997).
We restricted the range of b parameters for the interloper using the following physical considerations. If we consider all lines to be thermally broadened, then the b parameter varies inversely as the square root of the atomic mass. Any species lighter than Mg will have a larger b parameter. Since we are considering X to have transitions in the ultraviolet (UV), we may assume that it probably has a mass greater than or comparable to that of Mg. We therefore place an upper limit of . Also, we did not consider any species with atomic number greater than ∼100 and so we adopted a lower limit of .
By comparing sets of model and synthetic blended lines in this way, we can constrain the range of possible values of dλ, given some apparent .
We present the results of the simulation in Fig. 3, where we define and . The ripples on the three surface reflect errors (≲0.1dex in Γ_{X,Mgii}) introduced due to the restrictive upper limit placed on b(X).
Using Fig. 3, we can now specify a certain value of (i.e. we specify dλ) and an interloping transition of a species X (i.e. we specify Δλ) and find the lower bound on the column density which that species must have in order to have caused the said shift in . For example, if we are concerned with then , and the section of Fig. 3 with contains information that restricts the candidate interlopers.
Photoionization modelling
Given a species X, what is an upper limit on the column density in any given cloud? To answer this, we modelled the absorption clouds with a grid of cloudy models (Ferland 1993). cloudy estimates the photoionization conditions in absorption clouds by assuming that the cloud can be approximated as a series of zones, in a plane parallel geometry, illuminated by UV radiation from one direction. The fact that we are trying to model a (presumably) spherical cloud with UV flux incident from all sides is expected to introduce abundance errors of the order of a factor of 2 (Ferland 1993).
Photoionization equilibrium is characterized by (a) the shape and intensity of the UV spectrum, (b) the number density of hydrogen atoms in the cloud, n_{H}, (c) the size of the cloud, which we assume is the same for all species [the size can be found by specifying n_{H} and N(Hi), the neutral hydrogen column density] and (d) the metallicity, Z, of the cloud. It is also useful to define the ionization parameter as , where n is the number density of ionizing photons incident on the cloud. Here we follow the work in Bergeron & Stasińska (1986), Bechtold et al. (1987), Steidel & Sargent (1989) and Bergeron et al. (1994) and choose a powerlaw spectrum, , for f the incident flux at frequency ν. The low and highenergy extremes are given by (for and (for .
All the above parameters [U, n_{H}, N(Hi) and Z] vary from cloud to cloud and are poorly known quantities. We therefore use a grid of cloudy models to find an upper limit on the column density of a given ionic species using ranges in these parameters as follows. For U we adopt the range in line with the measured range of parameters for those clouds studied in M01a and W01 (Churchill 1997). Similarly, for N(Hi) we use . For Z we adopt a conservative range, (following Morris et al. 1986; Bergeron & Stasińska 1986; Bechtold et al. 1987; Steidel & Sargent 1989; Bergeron et al. 1994; Churchill 1997). The range of parameters we explore here is characteristic of the lowredshift sample. For a related discussion of DLAs, see Prochaska (1999) and Prochaska & Wolfe (2000).
The parameter n_{H} is the most difficult to constrain. To obtain an estimate of the lower bound, one may use the method, first suggested by Bahcall (1967), whereby one places an upper limit on n_{p}, the number density of protons in the cloud, using the equation
Here, N(Cii*) is the column density of Cii that has been collisionally excited to the state Cii*; n_{e} is the electron number density. We assume that and that , and so, using and from Songaila et al. (1994), we estimate that . Under these simple assumptions, this value is consistent with Morris et al.’s (1986) estimate, , and we adopt this as a conservative upper limit. A lower limit may be obtained by using upper limits on the cloud sizes. We used results from Bergeron & Stasińska (1986) to place a lower limit of where we have used cloudy to estimate that . We compare our parameter range estimates above with those assumed in Lopez et al. (1999) and find that they are consistent after considering the differences between the assumed radiation field.
We present the results of the grid of cloudy models in Table 1, where we only give the value of , the maximum value of Γ_{X,Mgii} that occurred in the grid. We also present the model responsible for this value of Γ_{X,Mgii}. All values of should be treated as having an error of at least 0.3dex. We find that these results are consistent with the predictions of model 1 in Lopez et al. (1999) and are also consistent with their measured values.
A search for interlopers and results
Having obtained estimates of the column density upper limits (Table 1), we can use the results of Fig. 3 to restrict the possible values for Δλ for any of the species above. Fig. 3 shows the value of this test, since, given a value for Γ_{X,Mgii}, we can place very stringent constraints on the allowed range in Δλ. Then, we search tabulated transitions for each species to see if any candidate interlopers exist.
Table 1 and the results expressed in Fig. 3 show that any interloper must have a wavelength differing by no more than ∼0.1Å from that of Mgiiλ2796. We searched published energy level tables –Moore (1971) and the Vienna Atomic Line Database (VALD, Piskunov et al. 1995; Kupka et al. 1999)– to find potential lines of those species shown in Table 1 which satisfy the above criteria. We searched only for E1 transitions, as all other transitions are associated with a large factor of suppression. For example, E2 and M1 transition probabilities are suppressed by factors ∼α^{−2} and ∼ respectively (Z_{n} is the nuclear charge). Firstforbidden E1 transitions can also be considered since the oscillator strengths are typically ∼ times that of a normal E1 transition in a related or similar ion.
Oscillator strengths, f_{X}, are not generally available for the potentially interloping transitions, so we must adopt an upper limit by using other (stronger) transitions for the same multiplet (for which oscillator strengths are known).
Given upper limits on f_{X} and Γ_{X,Mgii} and a value for Δλ, we can estimate dλ and hence derive a conservative upper limit on the interloper's effect on .
Although Fig. 3 relates specifically to Mgiiλ2796, similar simulations were carried out to search for interlopers near all other ‘host’ absorption lines (i.e. those in table 1 in M01a). The model parameters were also extended to DLAs (e.g. We have not been able to identify any candidate interlopers satisfying all of the criteria above.
We note that cloudy only computes column densities for the most abundant species. However, it is possible that other lowabundance species could mimic a shift in α. Therefore, we searched for interloping lines – in Moore (1971) and VALD (Piskunov et al. 1995; Kupka et al. 1999) – of ionic species of elements such as Sc, Ti, V, Cr, Mn, Co and Zn and found only one possibility: a Crii transition with that lies 0.13Å to the red of the Zniiλ2026 ‘host’ line. This line has an oscillator strength of 0.047 according to VALD. We simulated its effect on for the highz sample.
We obtained the largest value of directly from the observed QSO absorptionline data. We then generated a synthetic, singlevelocity component spectrum containing Zniiλ2026 and the Crii interloper. We also generated several other lines seen in the highz sample. We produce synthetic spectra with several different values of and then used vpfit to measure . The result was that the measured value of was always consistent with the nominal value – that is, we could detect no effect due to the interloping Crii transition.
In summary, we have found no interloping species that could have caused a significant shift in . We do note, however, that the ionic energy level data in Moore (1971) and VALD are incomplete for many highly ionized species, and so we cannot rule out the possibility of blending with undiscovered transitions. Also, as noted previously, this analysis only takes into account blends with atomic, not molecular, species. First, typical oscillator strengths for molecular transitions are much smaller than for resonance atomic transitions. Also, the typical column densities of the lowabundance elements in Table 1 suggest that very few molecules will exist in the QSO absorption clouds. Therefore, we consider the possible effect of any molecular blending to be negligible.
Line removal as a test for systematic blends
If a particular host transition (e.g. Mgiiλ2796) were systematically blended, then removing that transition from the analysis in M01a and W01 should result in a different value of . Line removal is also a rough test for isotopic ratio or hyperfinestructure effects as introduced in Sections . We therefore remove a particular transition in an absorption complex only if doing so still allows us to obtain meaningful constraints on . For example, we do not remove the Mgiiλ2796 line from a fit to an system if the Mgiiλ2803 line is not present. We also removed entire species where possible. A small difficulty occurs with Crii and Zniiλ2062, whose profiles are sometimes blended together (i.e. when the velocity structure is broad). If the profiles were blended then we simply removed both transitions simultaneously, only removing them separately when the lines were completely separated. We obtain with the same procedure described in M01a.
Removing a transition from the fit may result in a slightly modified estimate for the velocity structure, revising the estimate for for that system. The question is, how robust is the estimate for , averaged over the entire sample, to this line removal process?
We present our results in Fig. 4. For each transition removed, we compare the weighted mean of (and the associated 1σ error) before and after the line removal. It must be stressed that both the value and the errors are not independent, and so interpreting Fig. 4 is difficult. However, it seems that there are no extreme deviations after line removal. Thus, we might conclude from such a diagram that systematic blending has not occurred and that we have no reason to suspect that one or more of the laboratory wavelengths is incorrect.
Concentrating on the heavier ions, we might also conclude that there is no effect on due to possible isotopic ratio or hyperfinestructure effects. We cannot draw similar conclusions for Mg since we cannot remove all Mg lines simultaneously (this would prevent us from obtaining meaningful constraints on for the lowz sample). We also see that any lack of thermal equilibrium in the absorption clouds, leading to unequal population of the hyperfine levels of Aliii (see Section 2.5), has had no significant effect on . The large error bars on the Aliii points reflect the fact that those lines are only present in six of our highz systems.
Qualitatively, this line removal test suggests that there is no one single ‘host’ transition or species that is systematically blended and has therefore affected our results in any significant way. The results also suggest that isotopic ratio or hyperfinestructure effects have not significantly affected the heavier ions in our analysis. This point is more fully explored in Sections 6 and 7.
Detailed Analysis: Atmospheric Dispersion Effects
The problems that arise due to atmospheric dispersion effects can be understood with reference to Fig. 5. If we observe a point source with the spectrograph slit at some angle θ to the vertical (as projected on the sky), then the object is dispersed with some component along the spectral direction of the slit. The optical design of the Keck/HIRES combination at the time at which a large fraction of our data was taken was such that for ξ the zenith angle of the object (Tom Bida and Steve Vogt, personal communication). As we noted in Section 2.12, this leads to two effects on the wavelength scale:

(i)
The first effect is a stretching of the spectrum (relative to the ThAr calibration frames) as a result of the angular separation of different wavelengths entering the slit. Consider two wavelengths, λ_{1} and λ_{2}, falling across the slit as in Fig. 5. If we were to measure the spectral separation between these two wavelengths on the CCD, Δλ′, we would find that
where a is the CCD pixel size in angstroms, δ is the projected slit width in arcseconds per pixel (for HIRES, per pixel) and Δψ is the angular separation (in arcseconds) of λ_{1} and λ_{2} at the slit. Δψ is a function of the atmospheric conditions along the line of sight to the object and can be approximated using the refractive index of air at the observer and the zenith distance of the object (e.g. Filippenko 1982). Note that equation (6) assumes that the seeing profile is not truncated at the slit edges. If tracking errors or seeing effects cause truncation, equation (6) provides an upper limit to the measured Δλ′. 
(ii)
If truncation of the seeing profile does occur, it may be asymmetric (as the diagram illustrates). In the case of HIRES, the optical design is such that a blue spectral feature will have its red wing truncated and a red spectral feature will be truncated towards the blue. This is the case shown in Fig. 5. Thus, an asymmetry is introduced into the instrumental profile (IP) and this will be wavelengthdependent. The extent of the asymmetry will depend on Δψ (which depends most strongly on ξ).
Note that when we centroid spectral features to find their separation, both effect (i) (shifting) and effect (ii) (IP distortion) effectively ‘stretch’ the spectrum. For the lowz sample, this should lead to systematically positive values of . That is, atmospheric dispersion effects cannot explain the average negative value of seen in the QSO data. Removing these effects from our data would make the results more significant. Although there is no obvious correlation between and ξ as would be expected in this scenario (see Fig. 6), the potential effect is quantified in our analysis below.
For all objects that were observed without the use of an image rotator, we have used the observational parameters reported in the image header files to find average values of ξ (see Fig. 6), temperature, pressure and relative humidity for each object. This allows us to correct for the stretching of the spectrum implied by equation (6). We have recalculated in the same manner as described in M01a and we plot the results in Fig. 7. For several QSOs, only a single image header (i.e. for one of the QSO exposures out of a series) was available. In these cases we have added (in quadrature) an error of . This is an average value determined by calculating for several values of the zenith distance for those objects concerned and corresponds to an uncertainty in the average ξ of about ±20°.
We have also tested the effect of the IP distortions [effect (ii) above] on . We generated high S/N, singlecomponent absorption spectra (with for several redshifts. We included different combinations of transitions in the different spectra so as to make a representative estimate for both samples of QSO data (particularly the highz sample, which contains many combinations of different transitions). A wavelengthdependent PSF was constructed using the following procedure. We assume that the seeing (taken as 08) and tracking error (taken as 025) generate Gaussian profiles across the slit and that a wavelength λ_{0} falls at its centre. Other wavelengths are refracted to different parts of the slit depending on the zenith distance. To compute this, typical Keck atmospheric conditions (temperature, pressure, relative humidity) were used.
We truncate the profile by multiplying by a tophat function (whose width is equal to that of the slit) and then convolve the result with a Gaussian IP of width (Valenti et al. 1995). This final PSF is then convolved with the synthetic absorption spectra. The value of λ_{0} is ∼5500Å (i.e. the TV acquisition system is centred at about this wavelength). Using this value for λ_{0} and we find to be quite insensitive to the IP distortions, particularly in the highz regime. The effect on may be as large as for the lowz regime but roundoff errors in vpfit become important at this level of precision. To exaggerate the effect, we also chose extreme values of λ_{0} to be 4000 and 7000Å. The effect on in these cases was still negligible in the highz regime but increased to the level of for the lowz spectra. This latter value corresponds to ∼ of the correction due to the shifting effect. The sign of the above corrections to are such that applying them to the QSO data make the observed average more negative.
In summary, we find to be most sensitive to the shifting effect described by equation (6) and reasonably insensitive to the IP distortion effect. Thus, we do not apply a correction for the latter in the results shown in Fig. 7. The results show the values of once the shifting effect has been removed from the data. As expected, we see a significant decrease in in the lowz sample. The highz sample is insensitive to this systematic effect, as can be expected due to the many different transitions used with q_{1} coefficients of different signs. We summarize the results in Table 2 and compare them with those found in M01a and W01. We stress that the corrections above (equation 6) are upper limits. Seeing and tracking errors will diminish the effects substantially.
Detailed Analysis: Differential Isotopic Saturation
In this section we quantify the effect of differential isotopic saturation on our measured values of . As discussed in Section 2.3, this comes about when the strongest isotopic component of an absorption line saturates. The weaker isotopic components may still be on the linear part of their curves of growth. Using a single weighted wavelength to represent the composite profile position (i.e. the profile centroid) can therefore, in principle, result in systematic errors in .
The various isotopic abundances of concern are shown in Table 3. If only composite wavelengths were used then, if there exists a saturated dominant isotope in the QSO data, this may result in systematic errors in . However, composite wavelengths were only used for the four heaviest atoms (i.e. Cr, Fe, Ni and Zn) and, of these, only some Fe transitions appear to be saturated in some of our QSO spectra. The Cr, Ni and Zn lines were always very weak (typically absorbing only ∼20 per cent of the continuum) and so using composite transition wavelengths is justified.
It is therefore left to quantify the effect of differential isotopic saturation for the Fe lines in our spectra. Ideally, if we knew the isotopic wavelength separations for the various Feii transitions, then we could do this by simply fitting our spectra with this structure to find . The difference between the weighted mean and that found in M01a and W01 would indicate the size of the effect. In reality, we do not know the isotopic separations. We can, however, place an upper limit on this effect since we know (or can estimate) the isotopic separations of the Mg and Si anchor lines. If we fit our spectra with only the composite Mg and Si wavelengths then we will obtain maximal changes to since (i) these atoms have strong lines that saturate more often than the Feii lines and (ii) they will have the largest isotopic separations (i.e. largest mass shift).
We have carried out the above recalculation of using the composite wavelengths for the Mgi, Mgii and Siii lines. In the lowz sample, the fact that the Mg lines have their weaker isotopes at shorter wavelengths means that should be systematically more negative if the composite wavelengths are used. The situation is again more complicated in the highz sample. We find only very small differences between the results and those of M01a and W01, the weighted means for the low and highz sample being and respectively [compared to and from M01a and W01].
Can we understand the magnitude of these corrections? Imagine a single velocity component and consider finding from, say, the Mgiiλ2796 and Feiiλ2344 lines only. One might estimate the maximum effect of differential isotopic saturation by assuming that, when all isotopes of Mgii are saturated, the line centroid would lie at the unweighted mean isotopic wavelength. This would result in a correction to of ∼ (i.e. the observed value would become more negative when the correction is applied).
However, several factors will reduce this correction in reality. First, once a line begins to saturate, then because the intrinsic linewidth (i.e. the b parameter) is far greater than the isotopic separations, the weaker isotopes are swamped by the dominant one. This reduces the effect on . Secondly, if the velocity structure is more complex (i.e. if it contains more than one velocity component), the unsaturated components, which do not suffer from the differential isotopic saturation problem, will provide much stronger constraints on . Finally, as more unsaturated transitions (e.g. Feii) are incorporated into the fit of a particular QSO absorption system, the effect is further reduced since some constraints on come from shifts between the Feii lines.
It is therefore clear that differential saturation of the isotopic components of the various Feii lines cannot have led to any significant effect on . Our results also indicate that the weighted mean for the highz sample will not be very sensitive to errors in our estimate of the Siii isotopic separations.
Detailed Analysis: Isotopic Abundance Variation
The analysis we have carried out so far assumes terrestrial isotopic abundances. However, it is quite possible, even likely, that these highredshift gas clouds have rather different abundances. As we noted in Section 2.4, a substantial systematic shift in α might be mimicked and will be particularly pronounced in the lowz sample since the Mg isotopic separations are the largest. No observations of highredshift isotopic abundances have been made, so we have no a priori information on the QSO absorption systems themselves. However, observations of Mg (Gay & Lambert 2000) and theoretical estimates for Si (Timmes & Clayton 1996) in stars clearly show a decrease in isotopic abundances with decreasing metallicity. For example, at relative metal abundances , i.e. about 20 per cent below terrestrial values (see Table 3). Theoretical calculations (Timmes, Woosley & Weaver 1995) suggest even larger decreases.
Estimates have also been made of the metal abundances in a subset of the QSO absorption systems used in our analysis. The DLAs have a mean metallicity of (Prochaska & Wolfe 2000) and the Mg/Fe systems span the range to 0.0 (Churchill et al. 1999).
Therefore, we explored the effect of isotopic abundance evolution on for the entire range in isotopic abundance from zero to terrestrial values (i.e. in the case of Mg, the limiting case is . This allows us to derive a secure upper limit on the magnitude of the effect on .
We have recomputed for the entire sample after removing the following isotopes: ^{25}Mg, ^{26}Mg, ^{29}Si and ^{30}Si. The wavelengths for the dominant isotopic components, ^{24}Mg and ^{28}Si, for the various transitions of interest, are given in table 3 of M01a and are summarized here in Table 4.
The results are presented in Table 5 and Fig. 8. These show, as expected, that the lowz sample is more sensitive to this systematic error than the highz sample. We stress that the signs of the above corrections are such that applying them to the QSO data makes the observed average more negative. Since we have no information on the isotopic separations for the other atoms of interest, we cannot perform a similar analysis for all species. However, the corrections due to variation in the isotopic abundances of Mg and Si should be dominant (as we discuss in Section 2.4) so the corrections shown in Table 5 and Fig. 8 are conservative upper limits.
A General Test for Systematic Errors
Although we have already convincingly ruled out any significant wavelengthscale distortion using ThAr spectra (Section 3), the QSO and ThAr light does not follow precisely the same light path. Indeed, for these reasons we had to consider some of the previous effects such as atmospheric dispersion. Thus, it is conceivable that some residual longscale wavelength distortion or other simple systematic error remains in the QSO spectra.
The dependence of on the q_{1} coefficients is simple for the lowz sample but considerably more complicated for the highz sample. This can be seen by inspecting table 1 in M01a. If systematic errors are responsible for the nonzero that we measure, it is therefore surprising that is so consistent between the two samples. A priori, it would be even more surprising if we were to take subsets of the QSO data, grouped according to the sign and magnitude of q_{1}, and found a consistent value for .
We apply this idea to the highz sample only, since it contains transitions with a wide range of q_{1} values. Of the 21 absorbers in this sample, we can find a subset of 13 systems that include at least one anchor line (transitions with small q_{1}), at least one positiveshifter and at least one negativeshifter . If some kind of general longscale distortion of the wavelength scale was present, we might expect to find a change of sign in the average when comparing, for example, the negativeshifters or positiveshifters against the anchors (see equation 1).
From these systems we remove all transitions with q_{1} coefficients of ‘mediocre’ magnitude so as to delineate clearly the three types of transition (anchor, positive and negativeshifter). We then calculate for this sample. We then remove each type of transition separately and calculate again. There is a small complication due to the presence of the occasional blend between Crii and Zniiλ2062 since these lines shift in opposite directions. In cases where these two lines were blended together, we removed both of them from all fits since they are of different ‘type’ and cannot be removed individually.
We present our results in Table 6. Again, precise, quantitative comparison of the values of obtained with the various types of transitions removed is difficult since neither the values nor the 1σ error are independent. Also, when we remove the positiveshifters, the uncertainty in is quite large. However, we do not find any evidence to suggest that is somehow dependent on the sign of q_{1} from this test. That is, the highz results seem to be robust against the presence of simple systematic errors. Indeed, this is corroborated by the results of removing atmospheric dispersion effects (Table 2 and Fig. 7).
Conclusions
The most important conclusion of the work described here is that, after a detailed search for instrumental and astrophysical effects, we have found no systematic effect that can easily mimic a negative .
Of those effects we have identified and explored, two systematic effects have been found to be potentially significant: atmospheric dispersion and isotopic evolution. If we were to remove these effects from the data, the resulting would actually be more negative than the uncorrected effect (see Figs 7 and 8). These corrections are not applied to the data in M01a and W01 because (i) our estimates provide upper limits on the magnitude of the effects, and (ii) we prefer to be as conservative as possible as to the significance of the deviation of from zero.
Tighter constraints on α are possible for realistic S/N from current telescopes/spectrographs. The analyses above indicate that another order of magnitude precision can be gained using our methods before wavelength calibration and laboratory wavelength errors become significant. We hope that M01a, W01 and the present work encourage further wellcalibrated, higher S/N observations of QSOs and that an independent check on our results can be made.
Acknowledgments
Data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.
We are very grateful to Tom Bida and Steven Vogt who provided much detailed information about Keck/HIRES. We thank Ulf Griesmann, Sveneric Johansson, Rainer Kling, Richard Learner, Ulf Litzén, Juliet Pickering and Anne Thorne for much detailed information about their laboratory wavelength measurements. We also acknowledge helpful discussions with Michael Bessel, Bob Carswell, Roberto De Propris, Alberto FernándezSoto, John Hearnshaw, Alexander Ivanchik, Jochen Liske and Geoff Marcy.
We are grateful to the John Templeton Foundation for supporting this work. MTM and JKW are grateful for hospitality at the IoA, Cambridge, where some of this work was carried out. CWC received partial support from NASA NAGS6399 and NSF AST9617185.