A NEW METHOD FOR DETECTING VELOCITY SHIFTS AND DISTORTIONS BETWEEN OPTICAL SPECTRA

The DC method directly compares two spectra of the same object at the same wavelengths using a χ² minimization to find the best-fit velocity shift between them. It can be used to compare different exposures from the same telescope or spectra of the same object observed with different telescopes. The ultimate goal of this method is to find any velocity distortions or velocity shifts between spectra and remove these velocity distortions over the whole spectrum. This method only requires three user-defined parameters to be set: the scale of a simple Gaussian smoothing applied to the spectra, the size of the spectral regions to be directly compared (and for which an individual velocity shift measurement will be made), and the significance threshold above which regions containing reliable velocity shift information will be selected. The distribution of these individual velocity shifts can be used to decide the best method to correct for velocity shifts and distortions across the whole spectrum. Since there is no a priori knowledge of what the functional form of the distortions will be, the functional form of the final correction will depend on the results measured by the DC method. In addition, the DC method produces a reliable velocity shift and a reliable uncertainty estimate. It also has the advantage that one does not first need to identify features in the spectra, or which features are sharp and narrow enough to provide the best results.

Consider two spectra to be directly compared, f₁(i) and f₂(j), which are dispersed onto different arrays of pixels, i and j, and do not necessarily have the same velocity dispersion or resolution. To compare the two spectra and allow for a sub-pixel velocity shift between them, we used the barak package¹ to convolve the spectra with a Gaussian of constant FWHM in velocity space. Smoothing with a velocity kernel roughly equal to the instrument resolution allows us to increase the signal-to-noise ratio (S/N) as high as possible without significantly degrading the resolution of the spectra or losing spectral information. The size of this smoothing kernel is the first of the tunable parameters for the DC method. A natural choice for the FWHM is roughly equal to the resolution element of the telescopes, i.e., 2–3 pixels.

After smoothing, we interpolate the flux and error array of one spectrum with a cubic spline (we use "spline" colloquially as a verb). Therefore, we convert the flux, f₂(j), and error, σ₂(j), from a discrete spectrum to continuous functions of velocity, v, and user-specified parameters p—f_m(v, p) and σ_m(v, p)—where the extra parameters p are the parameters of the χ² minimization procedure discussed below, including the main parameter of interest here, the velocity shift between the spectra. It is ultimately better to spline the spectrum with more "information" in it so as to maximize the accuracy of the spline. For pairs of spectra with similar resolutions, it is preferable to spline the higher S/N spectrum, but if one spectrum has a smaller velocity dispersion, it may be better to spline it instead.

The spectra are broken into velocity "chunks"—sub-sections of the spectra of user-specified velocity length—and the chunk size is the second tunable parameter in the DC method. Chunks may overlap with each other to any user-specified extent. One might consider overlapping the chunks if looking for distortions over short wavelength scales at the expense of having independent measurements. However, in this work we only consider independent, non-overlapping chunks. We discuss the velocity length in the context of real examples in Sections 3.2.1 and 3.2.2. Each of the chunks in one spectrum is compared to the corresponding chunk in the other. We then minimize a modified χ² statistic to determine the best-fit velocity shift between corresponding chunks:

$$\begin{equation} \chi ^2= \sum \limits _{i}\frac{\left[f_{1}(i)-f_{m}(i)\right]^2}{\sigma _{1}^2+\sigma _{m}^2}. \end{equation} \tag{ 3 }$$

In our implementation we minimize χ² via a Levenberg–Marquardt method where the free parameters are a velocity shift, a relative flux scaling, and a tilt between the pair of chunks. In this equation, f₁ refers to the flux of the non-splined spectrum, f_m is the flux of the continuous, splined spectrum, while σ₁ and σ_m are their corresponding error arrays. After χ² is minimized we get a value and an uncertainty for the velocity shift, amplitude, and tilt, calculated from the relevant diagonal terms of the covariance matrix. The parameter in which we are most interested is the velocity shift; however, the other free parameters allow for possible differences in how the data were normalized and reduced.

The decision about how to set the velocity length of the "chunks" is based on several considerations. One determining factor is to minimize any differences in the shape of the continuum over the chunk length so that only first order corrections (scaling and tilt) are necessary to accurately compare the chunks. Smaller chunks help in this regard. However, the minimum chunk size should be larger than the extent of typical absorption features so that the velocity shift and slope parameters fitted in the DC method are not degenerate. Finally, in trying to decide a maximum size for velocity chunks, we aim to have as many independent chunks as possible. That is, while, in principle, the user may choose a very large chunk size, that may limit the understanding gained about any short-range or even some longer-range wavelength distortions present in the spectra. A secondary consideration is that, at lower S/Ns, a larger chunk contains more potential for the noise fluctuations in the flux to influence the velocity shift estimate and its uncertainty—see Section 2.6 for an exploration of low-S/N extremes. Another major benefit to breaking the spectra into chunks is that it allows the DC method to be applied to the whole spectrum without first identifying regions of absorption features. Once the process is complete, regions dominated by well-defined features have small error bars relative to those with fewer or broader features (Figure 7).

To test the reliability of the DC method, we produced synthetic spectra and ran them in pairs through a series of Monte Carlo simulations. Consider two such spectra with the same absorption line profile but shifted with respect to each other. We represent these spectra as single velocity chunks of unit continuum with a single Gaussian absorption feature. Each simulation involves two spectra with the same feature but with different random Gaussian noise. We then apply a sub-pixel shift between the two spectra followed by the DC method to see how well we can recover that shift. Figure 2 shows an example of the simulated spectra and the resulting distribution of velocity shifts from a Monte Carlo test with 5000 realizations. These spectra have an S/N of 100 per pixel in the continuum, a velocity dispersion of 1.5 km s⁻¹ pixel⁻¹, are Gaussians with a σ of 2.0 km s⁻¹ and are offset by 0.3 km s⁻¹. The reason that we choose a σ for the Gaussian that is so small is that many metal lines in high-resolution spectra are unresolved. Therefore, it is natural to choose a width of simulated spectra to be similar to what we would see as an unresolved metal line. The right-hand panel of Figure 2 shows that there is a roughly Gaussian distribution of recovered velocity shifts around the correct shift of 0.3 km s⁻¹. However, we find unreliable uncertainties on the measured velocity shifts. This problem occurs because of correlations between the fluxes of neighboring pixels, but it is simple to address and we detail the small corrections required in Section 2.3. Even with this problem with the uncertainty, this simple test demonstrates that the DC method recovers (1) the correct sub-pixel velocity shift and (2) a Gaussian distribution of velocity shift measurements, at least at relatively high S/N.

Zoom In Zoom Out Reset image size

The Gaussian smoothing of the spectra causes the uncertainty in individual velocity shift measurements (derived from the χ² minimization) to be smaller than the Monte Carlo distribution of velocity shifts in Figure 2 (right panel). There are two factors contributing to this problem.

First, by smoothing the flux arrays of the spectra, the statistical noise per pixel is reduced, and this must be reflected in a correction to the error arrays. Without that correction, the smoothing leads directly to a reduced χ², $\chi ^2_{\nu }$, from the DC method which is too small. By repeating the simulations represented in Figure 2 with different S/N and FWHM of the smoothing kernel, we find that, to first order, $\chi ^2_{\nu }$ does not depend on the S/N and only depends on the FWHM. As such, we can model $\chi ^2_{\nu }$ as a function of smoothing size, as shown in Figure 3), to derive a correction factor, ξ_χ, for the error arrays of the spectra. From Figure 3 we find best the approximation for the underestimation of $\chi ^2_{\nu }$ to be

$$\begin{equation} \xi _\chi = \exp (-0.16x^2-0.62x-0.76), \end{equation} \tag{ 4 }$$

where x is the FWHM of the smoothing kernel in number of pixels. After finding the correction value, ξ_χ, we then multiply the corresponding error array by the square root of ξ_χ prior to performing the χ² minimization in the DC method. This modifies Equation (3) to

$$\begin{equation} \chi ^2= \sum \limits _{i}\frac{\left[f_{1}(i)-f_{m}(i)\right]^2}{\xi _{1\chi }\left(\sigma _{1}\right)^2+\xi _{m\chi }\left(\sigma _{m}\right)^2}\,, \end{equation} \tag{ 5 }$$

where ξ_1χ and ξ_mχ are the corresponding correction factors to the un-splined error array and the splined error array, respectively. With these corrections made, we obtain a reduced χ² value of approximately unity when analyzing smoothed, simulated spectra.

Zoom In Zoom Out Reset image size

The second correction we must make is to the uncertainty measured in the velocity shift between chunks, σ_Δv, derived from the χ² minimization. The correction in the previous step (ξ_χ) only compensates for the overestimation of the error arrays because of the smoothing applied to the flux arrays. However, smoothing the fluxes also introduces extra covariance between pixels that is not accounted for by the previous correction. These correlations between pixels cause the root-mean-square (rms) flux variations to be slightly smaller than implied by the error array, even after correction in the previous step. To determine a correction for this, we model the ratio of the uncertainty calculated from the χ² minimization to the 1σ width of the Monte Carlo velocity shift distribution, σ_Δv/Δv_MC, for a given S/N, as a function of x, the FWHM (in pixels) of the smoothing kernel. Unlike the previous correction, we find that this ratio depends on both x and the instrumental resolution, ${\rm FWHM}/2\sqrt{2 ln(2)}$, in units of the spectral dispersion, which we define as y, i.e., $y\equiv \left({\rm FWHM}/2\sqrt{2 ln(2)}\right)/{\rm dispersion}$. We find that the dependence of the ratio σ_Δv/Δv_MC on x and y is separable, that is, it can be expressed as ξ_σ × ξ_FWHM where ξ_σ is a function of x only and ξ_FWHM is a function of y only. We find that the best-fitting expressions for ξ_σ and ξ_FWHM are

$$\begin{equation} \xi _{\sigma } = \exp (-0.05x^2-0.17x-0.28), \end{equation} \tag{ 6 }$$

and

$$\begin{equation} \xi _{\rm FWHM} = -0.15y+1.22. \end{equation} \tag{ 7 }$$

Figure 4 shows values of σ_Δv/Δv_MC for a range of x which is appropriate for applications to real spectra while holding the resolution of the simulated spectra constant. Over-plotted is the total correction, ξ_σ × ξ_FWHM, as a function of x, for the three different dispersions. The three different curves represent the three values of y. For the typical values of y explored here (0.8 ≲ y ≲ 2.0) we recover a DC method uncertainty that is correct to within a few percent; this is apparent in the relatively small scatter around the three curves in Figure 4.

Zoom In Zoom Out Reset image size

It is worth pointing out that the corrections from Equations (4) and (7) (seen in Figures 3 and 4) do not appear to approach unity exactly as the smoothing kernal size approaches 0 km s⁻¹. If there were no smoothing at all, we would expect there to be no need for the correction functions, though we do not force our fits to obey this boundary condition. However, we have chosen only to model the corrections needed over a range of practical smoothing values. Smoothing kernels smaller than 1 pixel are not recommended as there is no additional information present on a sub-pixel level. Likewise, smoothing kernels larger than a few pixels are also not recommended because there will be a greater loss of spectral information.

With this correction, we can convert the uncertainty measured from the χ² minimization, σ_Δv, into a more accurate and statistically meaningful uncertainty on the velocity shift between two spectral chunks:

$$\begin{equation} \sigma ^\prime _{\Delta v} = \frac{\sigma _{\Delta v} \times \sqrt{\chi ^2_{\nu }}}{\xi _{\sigma } \times \xi _{\rm FWHM}}. \end{equation} \tag{ 8 }$$

If there were no correlations between pixels (such as in an un-smoothed simulation) we would expect $\sigma ^\prime _{\Delta v}$ to equal σ_Δv. However, since there is always some correlation between pixels in real spectra, even before we smooth the flux arrays, we must also multiply our final uncertainty estimate by $\sqrt{\chi ^2_{\nu }}$. Thus, Equation (8) provides the final velocity uncertainty estimate from the DC method.

Figure 5 shows another Monte Carlo simulation with all of the same input parameters as those used to create Figure 2, with the exception that the corrections from Equations (5) and (8) have been applied. This gives a reduced χ² of $\chi ^2_\nu \approx 1$ as well as an uncertainty roughly equal to the distribution of the measured velocity shifts.

Zoom In Zoom Out Reset image size

The third and final parameter that must be tuned in the DC method is a means of selecting which chunks are dominated by noise and which contain reliable measurements. If the S/N of the spectra, which clearly plays a major role in determining the uncertainty on the Δv measurements, were constant across the entire wavelength range of interest, we could select such regions simply by choosing a maximum velocity shift uncertainty and rejecting all Δv measurements with uncertainties larger than this threshold. However, in reality, the S/N of echelle spectra changes significantly over the optical wavelength range, showing a particularly strong fall-off at progressively bluer wavelengths than ∼4000 Å. A chunk must contain enough sharp features to return a velocity shift uncertainty somewhat smaller than the "saturation value" for that chunk size, the point where the noise fluctuations are contributing a comparable amount to χ² as the features. In our implementation of the DC method, we calculate a "tracker point" using simulations for each chunk to objectively determine whether that chunk contains enough spectral information—enough sharp features with a high enough S/N—to provide a reliable velocity shift measurement with a reliable uncertainty estimate. We describe this technique below.

We first determine how much information is in each chunk via a method outlined by Liske et al. (2008). They provide a means of calculating the uncertainty on Δv between two (chunks of) spectra, determined by how much spectral information is present in the chunk. This velocity uncertainty, $\sigma _{v_i}$, is calculated by the following equation for each pixel, i:

$$\begin{equation} \sigma ^2_{v_{ i}} = \left(\frac{1}{\frac{dS_i}{dv}} \right)^2 \left[ \sigma _{1i}^2 + \sigma _{2i}^2 + \frac{\left(S_{2i} - S_{1i} \right)^2 }{ {\left(\frac{c}{v}\right)^2} \left(\frac{dS_i}{dv}\right)^2} \sigma ^2_{S^{\prime }_i} \right]. \end{equation} \tag{ 9 }$$

This is dependent on the uncertainty in flux, the gradient in flux and in the uncertainty of the gradient of the flux at each pixel in the pair of spectra. S is the flux, dS_i/dv is the average gradient in the flux at the ith pixel, and $\sigma ^2_{S^{\prime }_i}$ is the average uncertainty on dS_i/dv between the two spectra. The Liske et al. uncertainty on Δv for an entire chunk is then the sum over all of pixels in that chunk:

$$\begin{equation} \sigma ^2_{v} = \frac{1}{\sum \limits _{i} \sigma ^{-2}_{v_{ i}}}. \end{equation} \tag{ 10 }$$

Calculating this Liske et al. uncertainty requires that both spectra be on the same pixel grid, so it is therefore always calculated after smoothing the spectra and the spline interpolation has been performed. We express dS_i/dv as an averaged, one-sided finite-difference derivative:

$$\begin{equation} dS_i/dv = \frac{\left(\frac{\left(S_{i+1}-S_{i} \right)_1}{\sigma ^2_1} + \frac{\left(S_{i+1}-S_{i} \right)_2}{\sigma ^2_2} \right)}{\left(\frac{1}{\sigma _1^2} + \frac{1}{\sigma _2^2}\right)}/dv. \end{equation} \tag{ 11 }$$

Propagating the errors in the above equation gives an expression for $\sigma ^2_{S^\prime _i}$,

$$\begin{equation} \sigma ^2_{S^{\prime }_i} = \left(\frac{1}{\sigma _{1i}^2} + \frac{1}{\sigma _{2i}^2} \right)^{-2} \left(\frac{\sigma _{1i+1}^2 + \sigma _{1i}^2}{\sigma _{1i}^2} + \frac{\sigma _{2i+1}^2 + \sigma _{2i}^2}{\sigma _{2i}^2} \right), \end{equation} \tag{ 12 }$$

and the following expression for the Liske et al. uncertainty on any pixel:

$$\begin{equation} \sigma ^2_{v_{i}} = \left(\frac{v_{i+1}-v_{i}}{dS} \right)^2 \left(\sigma _{1i}^2 + \sigma _{2i}^2 + z\right), \end{equation} \tag{ 13 }$$

where

$$\begin{equation} z=\frac{\left(S_{2i} - S_{1i} \right)^2 }{\left(\frac{S_{1{i+1}} - S_{1{i}}}{\sigma ^2_{1i}} + \frac{S_{2{i+1}} - S_{2{i}}}{\sigma ^2_{2i}} \right)^2} \left(\frac{\sigma _{1i+1}^2 }{\sigma _{1i}^2} + \frac{\sigma _{2i+1}^2}{\sigma _{2i}^2 } + 2 \right). \end{equation} \tag{ 14 }$$

The sum over all the pixels in the chunk, as per Equation (10), then yields the Liske et al. error for the chunk.

Our aim is to use the above Liske et al. formalism, applied to synthetic realizations of a chunk, to determine whether that chunk contains enough spectral information to yield a reliable Δv measurement, i.e., whether or not that chunk's spectral features, if any, are sharp and/or numerous enough to dominate contributions to χ² over those from noise fluctuations alone. From Figure 6 it is clear that the Liske et al. velocity shift uncertainty shows the same behavior as that determined using the DC method: it "saturates" at S/N ≲ 7 pixel⁻¹, indicating that, for a feature of relative depth 0.7 (see Figure 2 left-hand panel), the noise, not the feature, dominates χ² at lower S/N. We simulate what an average Liske et al. uncertainty would be for a given pair of chunks with no feature present and compare this value to the Liske et al. uncertainty calculated from the real pair of chunks. This allows us to select any chunks with Liske et al. uncertainty significantly lower than measurements from corresponding simulated chunks of continuum. These selected points will have reliable measurements of velocity shifts and accurate uncertainties.

Zoom In Zoom Out Reset image size

In practice, our approach is as follows. For each chunk of the real spectra we create 50 synthetic realizations of that chunk. The flux in each of these synthetic spectra is generated as random Gaussian noise with an average value of unity. The standard deviation of the Gaussian noise, as well as the error array, is set to a constant value which is the median value of the error array from the corresponding chunk in the real spectra. The mean Liske et al. uncertainty and rms deviation then provide a "tracker point" for that chunk: if the Liske et al. uncertainty in the real chunks is ⩽N times the rms below the tracker point, we accept that chunk as having sufficient spectral information to provide a reliable Δv measurement. The number, N, times the rms of the tracker points is the final parameter set by the user. Therefore, choosing a higher N will be more selective in determining which features to use, while a lower N-value will accept chunks with less information.

We test this method with simulations that can be seen in Figure 7. We simulate spectra in the same fashion as before but this time apply a 1.0 km s⁻¹ velocity shift between the two spectra and allow for chunks of different depths. The middle panel of Figure 7 shows the "tracker point" calculated via the simulations for each chunk compared with the Liske et al. uncertainty from the real spectra. As the feature depth increases, the spectral information increases and the Liske et al. uncertainty decreases. Once it decreases below N = 4 times below the Liske et al. uncertainty range from the simulations, the lower panel demonstrates that the velocity shift measurement becomes more reliable (and precise).

Zoom In Zoom Out Reset image size

In previous sections we address single, unresolved features in simulated spectra. However, real spectra can contain much more complicated structures. In this section we show that the uncertainty estimate on the velocity shift found by the DC method is robust to more complex absorption structures. From Equation (9), in the case of a single feature per chunk, we can see that as the feature becomes broader the Liske et al. uncertainty becomes larger, in linear proportion, because the gradient in flux decreases. Conversely, as the feature increases in depth (and the S/N in the continuum remains the same) the Liske et al. uncertainty becomes smaller, again in linear proportion until the features saturate, after which no additional information is added. Similarly, in the case of multiple, unblended absorption features, the uncertainty decreases with the square root of the number of features. These scalings stem from the fact that the Liske et al. uncertainty accurately reflects the spectral information available. Also, as demonstrated by Liske et al. (2008), this is true even for much more complex spectra with many blended features of varying depth and width, e.g., real Lyα forest spectra.

Figure 8 demonstrates that the uncertainty measured from the χ² minimization in the DC method traces the Liske et al. uncertainty for varying feature width, depth and degree of blending of two spectral features. In the case of two independent spectral features (far right side of Figure 8) we see that the uncertainties are relatively small, reflecting the increased information present, and very similar to the Liske et al. uncertainties. The chunk on the far left shows the two features entirely blended together, again resulting in small DC method and Liske et al. uncertainties. Because the features are on top of each other, the depth of the feature increases (lowering the uncertainty) but also broadens slightly (increasing the uncertainty). We can see that, combined, this leads to an uncertainty slightly larger than two independent features; again, this closely reflects the behavior of the Liske et al. uncertainties. In the central chunks of Figure 8 the two measures of uncertainty increase together to a maximum when the two features are blended such that they produce a broad feature without gaining much depth.

Zoom In Zoom Out Reset image size

These simulations demonstrate that the DC method uncertainty is accurate because it traces well the Liske et al. uncertainty when features broaden, deepen, and blend. The former should therefore prove reliable once absorption structure is generalized even further to complex structures like the Lyα forest spectra upon which the Liske et al. approach has previously been demonstrated. It also supports our use of the Liske et al. uncertainty in simulations for selecting chunks with sufficient spectral information for reliable velocity shift measurements.

There are limitations to the DC method at low S/N. To demonstrate this, 10,000 realizations of pairs of chunks were created with S/N ranging from 3 ⩽ S/N ⩽1000. Figure 6 shows that the uncertainty on Δv is underestimated with the DC method for a chunk with S/N ≲ 7 per pixel in the continuum. Therefore, the DC method cannot deliver reliable information about the velocity shifts between two spectra at these low-S/N values. However, the tracker points described above are designed exactly for this reason—to only select chunks with accurate uncertainties. It is worth recalling that the other methods, like line-fitting and cross-correlation, will also fail at low S/N, so the DC method is not alone in this shortcoming.

While it is not recommended to trust any chunk with an underestimated uncertainty (i.e., falling within N times the simulation rms of the tracker points), there is an even more extreme breakdown of the DC method at lower S/Ns. The left-hand panel of Figure 9 shows a Monte Carlo simulation of 10,000 pairs of spectra at an S/N of 10 per pixel in the continuum. While a chunk with S/N of 10 will likely be rejected, it is interesting to note that we still recover a roughly Gaussian distribution around the correct (0.3 km s⁻¹) velocity shift. However, in the case of very low S/N (right-hand panel of Figure 9) we see that even the distribution is non-Gaussian. Rather, there seems to be an artificial preference for half-pixel shifts. This is most likely due to the need for interpolating the flux and error arrays in one spectrum; in our implementation of the DC method, we used a spline. As the splined error array is an overestimate of the actual uncertainty away from pixel centers, there will be a small dip in χ² when half-pixel shifts are applied between the pair of (chunks of) spectra. Regardless of the cause, it is important to note that chunks with such low S/N will never be accepted by even the most lenient selection of an N-value for the tracker points.

Zoom In Zoom Out Reset image size

Summarizing the results of these tests, the DC method appears to be a robust tool that can measure sub-pixel velocity shifts between spectra of the same object. It is model-independent, meaning that there is no bias from user input, and does not require any prior demarcation of spectral features. Since it can be applied quickly and automatically across a pair of spectra in chunks, it can be used to find velocity shifts between each pair of chunks, as well as changes in these shifts as a function of wavelength, i.e., spectral distortions. Moreover, the DC method provides a robust estimate of the uncertainty in its velocity shift measurements. However, the DC method does have its limitations at S/N lower than ∼7 per pixel in the continuum. As long as this is kept in mind, the DC method can be applied to real spectra to better understand systematic errors between telescopes or epochs of observation. The next section explores the DC method as applied to the QSO J2123−0050.

J2123−0050 is a bright (r = 16.45 mag) quasar at redshift z_em = 2.261 which was identified in the Sloan Digital Sky Survey data release 3 (Abazajian et al. 2005). The presence of so many spectral lines of varying degrees of depth and width make J2123−0050 an ideal target on which to test the reliability of the DC method. The HIRES observations were performed by Milutinovic et al. (2010) in 2007 with a seeing of 0.3–0.5 arcsec, resulting in a spectrum that covers 3071–5869 Å and has an S/N in the continuum of ∼25 per 1.3 km s⁻¹ pixel at ∼3420 Å and roughly 40 at ∼5000 Å. The UVES spectrum was observed by van Weerdenburg et al. (2011) in 2008 with wavelength coverage of 3047–9466 Å and an S/N in the continuum of roughly 65 per 1.5 km s⁻¹ pixel at ∼3420 Å and 60 at ∼5000 Å. There are 86 H₂ lines and seven HD lines over the wavelength range 3071–3421 Å.

In the application of the DC method to J2123−0050, we only consider the final one-dimensional spectra, i.e., those resulting from the combination of many individual exposures. For each of the final HIRES and UVES spectra, the contributing exposures were combined without concern for overall velocity shifts or velocity distortions between them. The position of the quasar in the slit during a given exposure imparts an approximately constant velocity offset, and this offset will vary between exposures. These velocity offsets are only of secondary importance for varying fundamental constant measurements, as can be seen from Equations (1) and (2), and generally only lead to a very slight broadening of narrow spectral features. However, velocity distortions between exposures would be a much more important systematic effect for such measurements. The DC method could be used to identify and quantify such velocity shifts and/or distortions between individual exposures from a single spectrograph. However, for the purposes of demonstrating the DC method in this paper, and for assessing the systematic effects on the previous measurements of Δμ/μ from the HIRES and UVES spectra of J2123−0050, instead we only apply it to the final, combined one-dimensional spectra here.

If there is a variation in μ, the redder Lyman H₂ transitions should shift in the opposite direction to (and by similar amounts as) the bluer Lyman transitions, i.e., the redder transitions' K-coefficients have opposite sign to those of the bluer transitions.² A long-range velocity distortion in one spectrum would, therefore, have a very similar effect on the relative velocities of the H₂ transitions as a varying μ and spuriously give a non-zero Δμ/μ. However, the transitions of atomic hydrogen in the Lyα forest should not vary with μ, and give us a means to break the degeneracy between a variation in μ and long-range distortions. The presence of so many narrow, and likely unresolved (∼6 km s⁻¹ wide), H₂ transitions which constrain Δμ/μ, and the Lyα forest of relatively broad H i lines which do not constrain Δμ/μ, make J2123−0050 a particularly interesting target on which to utilize the DC method.

While Malec et al. (2010) and van Weerdenburg et al. (2011) used the spectra only for measuring Δμ/μ from the H₂ and HD transitions in the Lyα forest, there are also well-defined metal transitions present redward of the Lyα emission line that could, in principle, be used to measure Δα/α in the damped Lyα system at z_abs = 2.059. There are also metal-line transitions at other redshifts, which would not usually be used to measure Δα/α. The presence of both these types of lines allows the DC method to be applied redward of the Lyα emission line. Thus, the Keck/HIRES and VLT/UVES spectra of J2123−0050 allow a comprehensive test for velocity distortions using the DC method over the wavelength range ∼3000–5200 Å. Finally, one unusual feature of this QSO is that it has time-variable C iv absorption systems falling just blueward of the C iv emission line in the wavelength range 4840–4925 Å (Hamann et al. 2011). Therefore, the region containing these time-variable absorption systems must be masked out when calculating possible velocity distortions between the spectra.

3.2.1. Redward of Lyα Emission

Figure 10 presents our results for the red portion of J2123−0050. The top panel shows the Keck and VLT spectra in the wavelength region containing a series of absorption features, ∼4650–5700 Å. Note that while there are more features in J2123−0050 both to lower and higher wavelengths than those shown in Figure 10, only this wavelength region contains overlapping spectra from both telescopes. The HIRES spectrum has a resolving power of R ≈ 110, 000 while the UVES spectrum has R ≈ 60, 000 in the red, meaning that the sharper metal lines are substantially less resolved in the UVES spectrum. This can clearly be seen in Figure 10 where some features appear slightly shallower in the UVES spectrum than in the HIRES spectrum. Ideally, the DC method would be applied to pairs of spectra with the same resolution, so as not to lose much information in the convolution process. However, since our aim is to measure velocity shifts between the spectra (cf. measuring metal column densities or Doppler broadening parameters), the information lost from smoothing should have only minor effects on the results. For the purposes of our demonstration of the DC method on J2123−0050, we do not attempt to quantify this here. Also, this region contains the aforementioned C iv variable absorption features as well as a UVES chip gap, both of which have been masked in the analysis. At 5000 Å the S/N for the HIRES spectrum is ∼40 per 1.3 km s⁻¹ pixel and for UVES it is ∼60 per 1.5 km s⁻¹ pixel. These S/N values, coupled with the sharpness and multicomponent nature of the absorption features in this wavelength range, yield an average uncertainty in velocity shifts calculated by the DC method of ∼0.18 km s⁻¹. In this portion of the spectrum we find that most features are well defined and, at most, ∼150 km s⁻¹. Because of this, we chose a chunk size that is slightly larger, 200 km s⁻¹, allowing us to cover whole regions of absorption but also maximize the number of measurements we can make in a limited wavelength range.

Zoom In Zoom Out Reset image size

The middle panel of Figure 10 shows how the reliable velocity shift measurements—i.e., those associated with features containing sufficient spectral information—were discerned from those unlikely to be reliable. Liske et al. uncertainties (blue points) for the real spectra that fall below five times the rms of the tracker points (green points) correspond to reliable chunks to use for Δv measurements. It is reassuring to note that chunks that are selected via this method do correspond to regions with strong features in both spectra.

The bottom panel of Figure 10 summarizes the results from the red portion of the J2123−0050 spectra. All velocity shifts and their 1σ uncertainties determined from the DC method are plotted, while the measurements selected as reliable are represented as red stars. It is immediately clear that the points identified as reliable show a much smaller scatter in velocity than those associated with featureless continuum regions of the spectra. The reliable points show a general consistency in their velocity shift values, with a possible trend toward larger velocity shifts between the spectra at longer wavelengths. A weighted linear regression was fit to these selected points, giving an average offset of 0.31 ± 0.05 km s⁻¹ and a slope of (2.1 ± 1.0) m s⁻¹nm⁻¹ with a χ² per degree of freedom ($\chi ^2_{\nu }$) of 2.6. Because $\chi ^2_{\nu } > 1$, a linear fit may not be the best model for the measured velocity shifts. Therefore, to obtain a more representative error estimate for the slope, we added a constant value in quadrature to the uncertainties of all the points such that $\chi ^2_{\nu }$ reduced to unity. This effectively marginalizes over the extra scatter around the best-fitting line. Doing this gives us a best-fitting line with a slope of (0.4 ± 1.7) m s⁻¹nm⁻¹, corresponding to a null detection of distortions between the spectra. However, this may not be the full story. If the four features blueward of ∼5000 Å are excluded from the fit then we measure a slope of (4.9 ± 2.2) m s⁻¹nm⁻¹ over a range of 450 Å, which is marginally non-zero. It is therefore possible that a velocity distortion does exist between the Keck and VLT spectra in this case, though it may have a complex variation with wavelength. Finally, the overall offset of 0.31 km s⁻¹ is likely due to the quasar, on average, being in a slightly different part of the slit when observed with Keck/HIRES and VLT/UVES and does not directly impact measurements of varying constants (see Equations (1) and (2)). The positive sign of this average offset means that the HIRES spectrum is redshifted relative to the UVES spectrum, on average.

Although we find no significant velocity distortion in our very simple (linear) fit above, it is interesting to understand how precisely the distortion measurement could constrain systematic errors in a varying-constant analysis. As an example of a typical Δα/α analysis, we imagine an absorber with metal-line transitions spanning the ∼4650–5100 Å wavelength range used in measuring the most extreme slope between the HIRES and UVES spectra (shown in Figure 10), and consider Mg ii λ2796 and the two Fe ii transitions with rest frame wavelengths of 2586 and 2600 Å. If we assume that Mg ii λ2796 falls at the red edge of the range in Figure 10, implying an absorption redshift of z_abs = 0.815, then the Fe ii λ2586 line would fall at the blue edge. These transitions therefore form a representative combination for Δα/α studies from the wavelength range in Figure 10. These three transitions are commonly seen in quasar spectra and, because Mg ii λ2796 is insensitive to Δα/α (q = 211 cm⁻¹) and the Fe ii lines are strongly dependent on α (averaged to give q = 1410 cm⁻¹) this combination of lines normally contributes strongly to a constraint on Δα/α. The (4.9 ± 2.2) m s⁻¹nm⁻¹ slope corresponds to a total velocity distortion of 182 ± 81 m s⁻¹ between the UVES and HIRES spectra over the 370 Å wavelength range from the Fe ii λ2586 transition to the Mg ii λ2796 line. Equation (1) implies that a distortion of this magnitude and sign would lead to a measurement of Δα/α being 10.0 ± 4.5 ppm greater in HIRES than UVES.

One particularly interesting aspect of these results is the precision with which we can measure/limit the distortion and its effect on Δα/α. In this example, that precision (i.e., 4.5 ppm) is smaller than the typical statistical uncertainty on Δα/α for Mg/Fe ii absorbers in Murphy et al. (2004) and King et al. (2012) with spectra of similar S/N to those studied here. That is, the DC method is generally capable of utilizing enough spectral information in a quasar spectrum to put interesting limits on, or identify and measure, systematic velocity distortions which are important for varying α measurements. With a larger sample of pairs of quasar spectra, it would therefore be possible to identify systematic errors that are likely shared by all the pairs and understand their effect on Δα/α with higher precision.

Of course, the main caveat to this conclusion is that this Mg/Fe ii absorber is just one example of transitions used to measure Δα/α. If instead we had imagined a higher-redshift absorber, quite different transitions, perhaps with quite a different set of ΔQ values, would have fallen in the wavelength range covered by Figure 10. The systematic effect on Δα/α caused by the velocity distortion (assuming it is real; recall that it only has a 2σ statistical significance) would then be quite different, perhaps even of the opposite sign. Nevertheless, the above discussion highlights the invaluable information about systematic errors that it is possible to derive using the DC method.

3.2.2. Blueward of Lyα Emission

In Figure 11 we present our results for the blue portion of the spectra. The top panel covers all of the overlapping Lyα forest from the UVES and HIRES spectra, a wavelength range of 3071–3964 Å, a comparable range as covered in the red portion (Figure 10). The HIRES spectrum has a resolving power of R ≈ 110, 000 while the UVES spectrum has R ≈ 53, 000 in the blue, meaning that some of the sharper features, particularly metal lines, are less resolved in the UVES spectrum and consequently appear slightly broader and shallower than in the HIRES spectrum. This is less of a concern in the blue portion of the spectra than in the red (where it did not seem to have a substantial effect on the results) because the relatively broad Lyα forest lines are completely resolved in both spectra. At 3420 Å the S/N for the HIRES spectrum is ∼25 per 1.3 km s⁻¹ pixel and for UVES it is ∼65 per 1.5 km s⁻¹ pixel. These S/N values, coupled with the broader but more numerous features than in the red portion, yield an average uncertainty in shifts calculated by the DC method of ∼0.06 km s⁻¹. As with the red portion of the spectra, we want to use as large a chunk size as possible without causing the velocity shifts between the pair of spectra to be impacted by noise in the flux. Firstly, there are more spectral features in the forest, so the chunks can be larger before they begin to be dominated by noise. Secondly, the many forest lines blend with each other and often produce spectral features that are much broader than metal-line complexes seen in the red part of the spectra. However, simple visual inspection reveals that ∼500 km s⁻¹ chunks contain blends that are complex enough, with enough internal structure (e.g., sharp changes in flux), and almost always included some flatter, continuum-like regions, to reduce any risk of significant degeneracy between the parameters in the χ² minimization, particularly the velocity shift and tilt parameters. Therefore, we chose a chunk size of 500 km s⁻¹ blueward of the Lyα quasar emission line. For higher (lower) redshift Lyα spectra, where the number density of forest lines is larger (smaller) than in our spectra, one would expect to use larger (smaller) chunks.

Zoom In Zoom Out Reset image size

The center panel of Figure 11 shows the Liske et al. velocity uncertainties and their corresponding tracker points. However, the choice of N for the tracker point selection process took into account several aspects of the Lyα forest region which differ significantly from the red portion. Since the Lyα forest region contains many broad features (from atomic hydrogen) and also numerous narrow features (metals and, for this QSO, H₂ and HD) in every chunk, there is a much smaller dynamic range in the Liske et al. uncertainties. That is, the most reliable and least reliable chunks are less distinguished from each other. In principle, almost all chunks offer a reasonably reliable velocity shift estimate (and uncertainty). However, one may choose to only select the most reliable chunks, but this requires a slightly more careful tuning of the N parameter.

The bottom panel of Figure 11 shows our results from the blue portion of J2123−0050. Just as with the red portion of the spectra (Figure 10), the points identified as reliable show a much smaller scatter in velocity than those associated with regions of the spectra containing less spectral information. The reliable points show a general consistency in their velocity shift values, with a possible trend toward larger velocity shifts between the spectra at longer wavelengths. We find a slight offset between UVES and HIRES, similar to that found for the red portion, 0.44 ± 0.03 km s⁻¹. As described for the red results, we fitted a line to the velocity shifts found for the selected (i.e., reliable) chunks and increased the uncertainties (by adding a constant value in quadrature to them) to obtain $\chi ^2_{\nu } \approx 1$. This gives a slope of (1.7 ± 1.8) m s⁻¹ nm⁻¹. This result is <1σ from zero and is consistent with no velocity distortion between the UVES and HIRES spectra in the Lyα forest.

It is interesting to consider the precision with which we have constrained a linear velocity distortion above (i.e., 1.8 m s⁻¹ nm⁻¹) and its implications for varying-μ analyses. This portion of the spectra (3071–3421 Å) contains the H₂ and HD transitions at z_abs = 2.059 which Malec et al. (2010) and van Weerdenburg et al. (2011) used to measure Δμ/μ. Over this wavelength range, our precision would correspond to a velocity distortion of ∼63 m s⁻¹. For the Lyman H₂ lines over this wavelength range, the K-coefficients—their sensitivity to variations in μ—vary fairly smoothly from −0.03 for the reddest transitions to +0.02 for the bluest transitions.³ Using Equation (2), we find that this maximum distortion would correspond to a Δμ/μ measurement on UVES being more positive than its corresponding HIRES measurement by 4.2 ppm.

This maximum (1σ) systematic error of 4 ppm is similar to the statistical uncertainty on Δμ/μ derived from the lower S/N of the two spectra, the HIRES spectrum, by Malec et al. (2010). That is, the DC method, applied to the Lyα forest, utilizes enough spectral information to identify and measure systematic velocity distortions which are important for varying-μ measurements. Secondly, the size of the possible systematic error in Δμ/μ created by the marginally significant velocity distortion we find, is consistent with the very small (and insignificant) difference between the Δμ/μ values found by Malec et al. (2010) and van Weerdenburg et al. (2011): +5.6 ± 5.5_stat ± 2.9_sys and 8.5 ± 3.6_stat ± 2.2_sys ppm respectively. That is, in this case, the DC method is consistent with the DC of the Δμ/μ values derived from the same pair of spectra. However, we emphasize that, by utilizing additional spectral information, not just the transitions used to derive Δμ/μ, the DC method enables the discovery of systematic errors to which a simple comparison of the Δμ/μ values is not sensitive.

Recently, Rahmani et al. (2013) compared VLT/UVES asteroid (i.e., reflected solar) spectra with a higher-resolution Fourier transform spectrum of the Sun (Kurucz 2005, 2006) and identified velocity distortions over the wavelength range ∼3300–3900 Å with significant slopes, ranging between 2 and 6 m s⁻¹ nm. If distortions as strong as 6 m s⁻¹ nm were present between the HIRES and UVES spectra studied here, over the same wavelength range, they would have been detected by the DC method in Figure 11. Nevertheless, it is worth noting that the sign of the (albeit statistically insignificant) distortion we identified, is the same as those identified by Rahmani et al. (2013) in the their UVES asteroid exposures (assuming that any distortions in the HIRES spectra are negligible).

3.2.3. Intra-order Distortions

From the analyses of Griest et al. (2010) and Whitmore et al. (2010) we expect there to be some velocity distortions present on the scale of echelle orders (as opposed to the much longer-range distortions considered above). Griest et al. (2010) measured peak-to-peak intra-order distortions in individual HIRES quasar exposures ranging from 300 m s⁻¹ to 800 m s⁻¹ and Whitmore et al. (2010) found peak-to-peak intra-order distortions in individual UVES exposures from 100 m s⁻¹ to 200 m s⁻¹. Therefore, it is quite possible that pairs of spectra, each comprising many exposures, such as those of J2123−0050 considered here, could also contain these intra-order distortions. We search for evidence of intra-order distortions below using the DC method.

Griest et al. (2010) and Whitmore et al. (2010) used spectra of their objects observed through an iodine (I₂) cell to establish two different wavelength scales: that from the usual ThAr comparison lamp exposures and that from the dense I₂ forest of lines when compared to a laboratory Fourier transform spectrum. This allowed them to determine, with high precision, velocity shifts between the two wavelength scales as a function of position along each echelle order falling within the wavelength range covered by the I₂ transitions (∼5000–6200 Å). By comparison, the DC method only utilizes the spectral information in the object spectra themselves, not the very dense information contained in ThAr or I₂ spectra. Therefore, searching for intra-order distortions within individual orders of our J2123−0050 spectra is not really possible; there is simply not enough spectral information, even in the Lyα forest, to provide enough reliable velocity shift measurements across any given echelle order.

To remedy this lack of features we "stack" the velocity shift information for all echelle orders in each portion of the spectra (blue and red). To do this we take the residuals between the reliable Δv values and the lines of best fit in Figures 10 and 11, and calculate how far away, in velocity, each chunk is from the center of its Keck/HIRES echelle order. This information is plotted in Figure 12 for the blue and red portions of the spectra, separately. We find no obvious evidence in Figure 12 for intra-order distortions which have a similar shape and amplitude in all the Keck echelle orders, in either the blue or red portions. While Figure 12 assumes that the distortions are coherent among the echelle order structure of HIRES, we found similar results when stacking the velocity shift measurements with respect to their positions on the VLT/UVES orders.

Zoom In Zoom Out Reset image size

There are several possibilities for why Figure 12 shows no distortions but Griest et al. (2010) and Whitmore et al. (2010) find clear evidence for them. Firstly, these authors examined individual echelle orders in individual object exposures while S/N considerations force us to work with many exposures combined into a single spectrum. It is possible that when the orders are combined to make the final spectrum, some of the distortions are dampened by the overlap of neighboring orders. Another possibility arises from comparing spectra from different spectrographs. Because the order lengths are different in each spectrum, any coherent intra-order distortions across each spectrum will be differently phased and may substantially cancel each other out when analyzed with the DC method. On the other hand, they may reinforce each other, making the effect more pronounced.

To test whether the DC method is sensitive to intra-order distortions in our spectra of J2123−0050, given the echelle order structures of HIRES and UVES, we introduced a sawtooth distortion in all of the HIRES orders. We took the extracted, wavelength-calibrated echelle orders before they were combined into a single spectrum and introduced a distortion with a peak-to-peak amplitude of 800 m s⁻¹ into each echelle order. This distortion causes the edges of the orders to have a negative velocity shift and the centers of the orders to have a positive velocity shift. We then recombine the spectrum in the same manner as described by Malec et al. (2010) and reapply the DC method in the same fashion used to produce Figure 12. Figure 13 shows that we recover the additional 800 m s⁻¹ velocity distortion in the center of the HIRES orders. However, we do not recover any distortions at the edges of orders. The regions that show little sign of distortion at velocities (relative to the order centers) where we expect neighboring orders to overlap. There, the sawtooth nature of the artificial distortion causes overlapping orders to substantially cancel out the distortions. Therefore, if real intra-order distortions do exist in the HIRES spectra, but not in the UVES spectra, we would expect to see them in Figure 12 after applying the DC method to the combined, multiple-exposure Keck and VLT spectra of J2123−0050.

Zoom In Zoom Out Reset image size

Figure 13 demonstrates that the DC method can detect intra-order distortions, allowing us to set a rough upper limit to any intra-order distortions between HIRES and UVES which have a similar shape and amplitude across all orders. From Figure 13, it is clear that such intra-order distortions with a peak-to-peak amplitude ≳500 m s⁻¹ would be detectable, so we treat that as a conservative upper limit. Since we expect overlap of neighboring echelle orders to diminish the effects of intra-order distortions, it follows that we would see less evidence for intra-order distortions in the blue end of the spectrum, compared to the red, because the overlap is much greater for bluer orders (both UVES and HIRES are grating cross-dispersed spectrographs). This implies that there should be little impact on the measurements of Δμ/μ, because such measurements usually utilize the bluest set of echelle orders to access the Lyman and Werner H₂ transitions. However, it is important to note that, because the DC method can only detect distortions between pairs of spectra, it is possible that intra-order distortions with the same (or similar) shape and amplitude in both HIRES and UVES would not be detected. Having said that, if the two telescopes have similar distortions, then similar samples of absorbers on both telescopes, from which (spurious) measurements of Δμ/μ or Δα/α have been made, would be affected in the same manner, and this would not lead to different average measurements from the two different telescopes (like, for example, those leading to evidence for a variation in α across the sky).

Despite not finding any overall evidence for intra-order distortions, we do see one example of a possible coherent, short-range distortion. At the bluest edge of the red portion of spectra shown in Figure 10, from 4665 Å to 4700 Å, we noticed a negative trend in velocity shift between the HIRES and UVES spectra. We have further explored this region in Figure 14 by using a smaller chunk size (100 km s⁻¹) than in Figure 10 (200 km s⁻¹). This ∼30 Å region has an almost 500-m s⁻¹ distortion across it, including some of the only negative values of velocity shift (HIRES blueshifted relative to UVES) measured in the whole spectrum. We fitted a weighted linear regression to the velocity shifts of the selected features in the bottom panel of Figure 14. Despite being able to see clear, coherent distortions by eye, it is obvious that a linear fit does not describe these distortions well. It is unclear what causes this effect. Being in the red portion of the spectra, this would not have had any effect on the measurements of Δμ/μ performed by Malec et al. (2010) and van Weerdenburg et al. (2011). However, because only a 350 Å range of these spectra contain H₂ and HD which were used to constrain Δμ/μ, short-range distortions like the one seen at 4665–4700 Å could have a marked effect on Δμ/μ if present in the blue portion of one or both spectra.

Zoom In Zoom Out Reset image size