Wide Bandwidth Desired
Wide bandwidth high-resolution spectroscopy is one of the most valuable diagnostic tools for astronomy, from Doppler radial velocimetry in the search for exoplanets to studying the characteristics of stars through the size, breadth, and position of spectral features. A wide bandwidth is needed in radial velocimetry to collect as many photons as possible to reduce Doppler velocity noise. Wide bandwidth is desired for analysis of species population because the characteristic wavelengths are often distributed over a wide band.
However, it is difficult to achieve both wide bandwidth and high spectral resolution for a dispersive spectrograph, since this requires not only a large and expensive instrument but also a detector with many pixels. This is especially so for infrared wavelengths since the spectrograph size scales with resolving power times wavelength.
Stability Against Point Spread Function Drifts Desired
Second, a serious challenge for high-resolution spectrographs is the unwanted shifting of the focal spot (point spread function, PSF) along the dispersion direction, such as due to air convention or thermomechanical drifts. While use of a calibrant spectrum measured in parallel or comingled with the stellar spectrum can remove some lower order components of the drift, the higher order components corresponding to gaps between strong features in the calibrant spectrum can remain, producing a wavenumber error. These environmental insults are frequently the limitation, rather than photon noise, to achieving the highest desired Doppler precision, or sufficiently accurate spectra to use as a template for subsequent Doppler calculations.
Externally Dispersed Interferometry
A technique called externally dispersed interferometry1 (EDI) can dramatically reduce the necessary size of spectrographs for Doppler radial velocimetry1–9 and high-resolution spectroscopy.10–15 Other workers have adopted the EDI method from our laboratories and demonstrated a Doppler planet detection.5 We have demonstrated EDI in both applications in the near-infrared at the Hale 5-m telescope at Mt. Palomar Observatory (Fig. 1) with a version called TEDI, (the T denotes the TripleSpec16 spectrograph to which it was coupled). An earlier report2 describes the Doppler velocimetry, and this report describes the spectroscopy.
This is the first time multiple delays have been used on starlight to recreate a high-resolution spectrum. Previous EDI use on stellar spectroscopy used a single delay.11 Multiple delays had been used previously12 only on laboratory sources.
The TEDI project was designed primarily to test Doppler measurements, and the selection of interferometer glass etalons and observing schedule were prioritized for this purpose. These were not optimal for testing high-resolution spectroscopy (having a gap in coverage of delay space when contiguous coverage is needed to avoid ringing in the lineshape). Nevertheless, the preliminary results for high-resolution spectroscopy using multiple delays are a resounding success in two important areas.
Resolution, Robustness Demonstrated
We have demonstrated on starlight that the EDI technique can produce both high resolution, a factor of to beyond the native spectrograph resolution of , and have extremely wide simultaneous bandwidth, limited only by the bandwidth of the native spectrograph, in this case (0.95 to ).
Second, the native spectrograph in our study suffered from extremely large and irregular PSF drifts that would normally have precluded the high resolution we obtained, even if by some other means the slit width and focal spot had been reduced by several times and the detector pixel density increased. We realize now that a great advantage of the EDI technique, especially when used with multiple delays rather than a single delay as used previously, is an order of magnitude or more improvement in the stability against PSF drifts.
We observe an approximate 20-fold decrease in the translational reaction of a ThAr lamp line to a translational insult along the dispersion direction. This robustness to PSF errors is explained theoretically and demonstrated with measurements. A recently realized method of further reducing the PSF shift sensitivity to zero is described theoretically and demonstrated in a simple simulation which produces a 350× times reduction.
Third, we demonstrate EDI’s dramatic robustness to fixed pattern (FP) noise such as from bad detector pixels. These pollute the native spectrum but not the EDI-derived spectrum measured simultaneously.
The TEDI data set is instrument limited rather than photon limited. A theoretical description of EDI photon-limited sensitivity compared to other spectroscopy techniques is presented in a companion paper.
How Externally Dispersed Interferometry Works
The EDI technique uses a field-widened Michelson interferometer (Fig. 1) crossed with (in series with) a dispersive spectrograph to heterodyne unresolvable narrow features to a detectable low frequency moiré pattern. Detailed instrument description is in Ref. 2 and instrument theory for a single delay spectroscopy is in Ref. 11. The heterodyning effectively shifts the native spectrograph sensitivity peak to higher frequencies (Fig. 2)
Early EDIs used linear spectrographs in which the crossing dimension was orthogonal to the dispersion, where interferometer phase varied along the slit length. In more recent EDIs using an echelle spectrograph (having the important advantage of extremely wide bandwidth), the spacing between orders is limited and the crossing dimension is time. A series of phase-stepped exposures are taken with a sequence of one to a few fixed delay values, which map out delay-space contiguously up to a maximum value , currently 3 cm, but which could easily be higher.
After Fourier processing to shift the frequencies to their original value, a spectrum with a resolution 4 to 10 times higher than the existing grating spectrograph has been recovered. Higher boost factors are possible with purchase of more glass etalons to allow a larger contiguous delay range. The EDI resolving power is not determined by the focal spot size, but the maximum delay relative to the wavelength, , and thus can greatly exceed the ceiling of the conventional spectrograph. In contrast with a dispersive spectrograph, it is now independent of the slit width, focal blur, pixel Nyquist limitations, and grating dispersion.10,12,13
Externally Dispersed Interferometry Has Corrugated Point Spread Function
The EDI PSF10 (instrument response) is a corrugated peak or wavelet (bold curve of Fig. 3), having the envelope set by the native spectrograph (dashed curve), multiplied by a sinusoid set by the interferometer delay chosen by the user. The key point is that now the PSF slope can be much higher than the low slope of envelope (native spectrograph), being proportional to the interferometer delay. It is easy to implement the 1- to 5-cm delay needed to match the slope of a high-resolution spectrograph.
The PSF slope is now much more mathematically controlled, largely independent of the envelope shape. The interferometer has only three degrees of freedom (phase, period, and intensity offset), rather than the many degrees of freedom typical of a grating PSF (contributed by the aggregate of tens or hundreds of grating grooves). Second, since the interferometer component is small (), it is much easier to put in a stable environment than a stand-alone disperser, which for a 50,000 resolution can be several meters in size.
Combining Several Delays Makes a Peaked Point Spread Function
For use in high-resolution spectroscopy, the subject of this report, the many redundant peaks in the corrugated PSF for a single delay would be a problem. The solution is to use several delays. The central area (cross-hatched) is made unambiguous by combining several wavelets of different periodicities (delays), as demonstrated in simulation in Fig. 4, and in practice in Fig. 5. The latter shows measurements of two ThAr lamp lines resolved by the EDI technique, but which otherwise are unresolved by the native spectrograph.
Data Analysis Overview
The data analysis pipeline is described in several sections:
• (Stage 1) Order extraction [converts the two-dimensional (2-D) echelle data into a manifold of one-dimensional (1-D) spectra],
• (Stage 2) Phase-stepping arithmetic (separates fringing and nonfringing components),
• (Stage 3) Rebinning and regularization (linearizes horizontal axis to wavenumber and compensates for etalon glass dispersion),
• (Stage 4) Spectral reconstruction (shifts data in Fourier space, sums over etalons, and equalizes to final resolution).
Stage 1: Order Extraction
Two Stripes (A, B) for Each Echelle Order
After optional flat fielding to remove pixel-to-pixel gain variations, the 2-D data from echelle spectrograph detector, which are an FITS file, are converted to 16 1-D spectra () by collapsing vertically across each of four stripes for each order. There are four useful orders, enumerated A to D starting from the top (high pixel), and the fifth (E) is weak and mostly redundant in range with the D. These have wavenumber (wn or ) ranges of A: 4100 to , B: 5400 to 7100, C: 6700 to 8800, D: 8100 to 10,600, and E: 9400 to 10,700. These are wavelength ranges of A: 2.44 to , B: 1.85 to 1.41, C: 1.49 to 1.14, D: 1.23 to 0.94, and E: 1.06 to 0.93.
In the first version of the TEDI apparatus,6 both complementary interferometer outputs were sent through the spectrograph and detected. However, that interferometer design suffered from noncommon path errors. We rebuilt it in 2010 (version 2) using collimated beams, which eliminated the path errors. We also decided to use one of the two interferometer outputs to monitor the flux for maximizing starlight alignment at the input side (beam path described in Ref. 2). Although this prevents summing the two complementary outputs to create a nonfringing ordinary spectrum for each exposure, the latter is created by summing the 10 phase stepped exposures in the phase-stepping math process (Sec. 5) to create one nonfringing spectrum per set of exposures. The increase in flux due to better alignment of the starlight into the input fiber made up for the loss of half the outputs.
Hence we recorded four stripes per order, which we denoted sky A, fiber A, fiber B, and sky B. But since the companion sky stripes (which sensed blackbody background) were soon optionally subtracted from the data-carrying fiber stripes, we effectively had two stripes per order. (For the bluer orders, we could omit the subtraction.) The fibers A and B could contain either ThAr spectral lamp or starlight combined with ThAr light. We explored several data taking configurations. Typically, one fiber contained ThAr-only and the other ThAr + star. Often these roles were alternated during the night.
Fringe Comb Visible for Small Delay
Figure 6 shows a portion of an FITS echelle spectrograph output around the C-order, and Fig. 7 shows three closeups showing the effect of interferometer delay stepping. In Figs. 6 and 7, note the periodic fringes imprinted on the stellar spectrum by interferometer acting on the continuum portion of the input spectrum. A small interferometer delay (, “E1”) was used here, which produces fringes wide enough to be resolvable by the native spectrograph (res ). The sinusoidal fringes are not required to be resolvable, and for most of the higher delays (), are not. The heterodyning occurring here creates moiré patterns that manifest high frequency information beaten down to lower, resolvable frequencies.
Figure 7 is a closeup in C-order of three exposures of the same set as Fig. 6 but having different phase. The top fiber (A) has only ThAr spectral lamp, and the bottom fiber (B) is a sum of stellar and ThAr. The change in fringe phase between the two exposures can be seen in the different intensities of the ThAr lines, comparing top to bottom panel, and in the apparent sideways shift of the periodic comb in fiber B stripe.
Summing Across Stripes
For each of the four stripes of an order, the counts are summed across 15 vertical pixels to form a set of 1-D spectra outputted from stage 1. The pairs of white and red curves in the figures demark the 15-pixel vertical width of the stripes.
The stripe shape does not change significantly from night to night, and it merely needs to be offset vertically. The vertical summing does not correct for the effect of the slight slant in the spectrograph slit across the order, but the blurring this creates is insignificant since most of the counts are concentrated in 3 or 4 vertical pixels. However, the effect of the slant between the A and B stripes is more significant. Our routines compensate for this slight wavenumber shift between the two fibers in the pixel to wavenumber mapping done in stage 3, by having individual calibrations for each fiber. (However, our sky subtraction is currently unshifted.)
Detector Artifact Cleaning and Flatfielding
An optional first step, before order extraction, is to run a routine that removes electronic ghosting artifacts specific to this detector. This routine was supplied by TripleSpec spectrograph group at Cornell.16 Then an optional flatfielding exposure divides out some FP noise due to pixel gain variation. Flatfielding is not critical for the fringing component since a built-in flatfielding-like process occurs already during phase-stepping mathematics that rejects stationary patterns. However, flatfielding benefits the ordinary spectrum.
Stage 1: Output Is Two-Dimensional Array Stack Versus Delay and
For each of the stripes, the 10 sequentially recorded phase-stepped 1-D spectra are stacked into an array versus exposure number along -axis. Since exposure number is essentially delay, which is related to phase, this forms a pseudomultiphase format (Fig. 8). These arrays of 1-D uniphase data are thus in the same format as 2-D multiphase data of earlier EDIs (such as Fig. 5 of Ref. 1), where the phase spatially varies along slit length in each exposure. This is a visually intuitive and algorithmically convenient way of presenting the phase set data, and routines written earlier for multiphase data can be applied.
Echelle spectrographs require uniphase mode because they usually do not have enough vertical pixels across an order to support varying phase spatially. Instead, the phase is varied temporally via interferometer mirror movement in 10 exposures. But there are some advantages to the uniphase mode that makes it attractive in any case. In uniphase mode, the same pixel is used repetitively, so FP noises are easily suppressed (because each pixel has a constant gain). Second, it allows using the spatial dimension to splay across the target, for imaging spectroscopy.
Comparing the two different delays shown in Figs. 8(a) and 8(b) note that the pattern of phase shifts for the ThAr lines is different. These constitute a “fingerprint” that reveals during analysis the precise delay that was used for each exposure, so we do not need to rely on knowledge of the commanded interferometer mirror position.
Cosmic Ray Removal
For EDI, it is more reliable to do the cosmic ray (CR) removal when the data are in pseudomultiphase format, rather than earlier when individual exposures are in 2-D format. (In the latter, a bright ThAr line can be mistaken as a CR by a median filter.) We fit a sinusoid along each wavenumber column versus the phase. Unusually large residuals are deemed CRs and their values replaced by the sinusoidal fit value.
Stage 2: Phase-Stepping Math to Separate Fringing from Nonfringing
Reference 11 describes single delay EDI theory regarding spectroscopy. Fringing and nonfringing components can be computed from a series of phase stepped exposures of regular phase (evenly distributed around the phase circle) by
We now describe the phase-stepping algorithm we use to separate fringing () and nonfringing () spectra from the spectral array data () of Fig. 8. This algorithm can handle unknown and irregular steps and intensities, and therefore it is useful to other interferometry applications generally (the first author uses it fruitfully for analyzing Doppler velocimetry interferograms having irregular visibilities and intensities in shock physics17). This distinguishes our algorithm from most others in the literature,18 which have trouble handling unknown and irregular steps.
Alternate solving X and Y dependencies
The problem we face is to solve both and behaviors. The axis is wavenumber , and the axis is delay step position . We alternate between analyzing data along a column using provisional knowledge of delay steps to produce “template” spectra ( and ), and then analyzing along rows using provisional knowledge of the template spectra to produce detailed delay step values (flowchart Fig. 9).
The observed TEDI delay steps [Fig. 10(b)] are reasonably uniform so the initial calculation of the provisional templates by assuming constant steps produces such an accurate result for W and B that the whole process converges rapidly in only 1 to 3 iterations, using only a single order. (A much tighter constraint on the delay step would result with the increased bandwidth of multiorder processing in parallel, but this has yet to be implemented.)
Figures 10 and 11 show example results for and for star *11UMi plus weak ThAr lamp, for delay E4 (). The magnitude of the whirl is typically just a few percent of the ordinary nonfringing spectrum for the higher delays where the sinusoidal comb is unresolved. For the lowest delay (E1, ) where the sinusoidal comb is fully resolved the magnitude can approach 100%.
Pathologies in the iterative process can be detected by monitoring the plot of in Fig. 10(b). Since the typical TEDI delay step is fairly uniform, if the observed is erratic this indicates a problem in the data. Often this is due to an uncorrected CR.
Phase step irregularity from wide bandwidth
The phase step can be irregular because even with a regular delay step , the phase varies across the wide band asFig. 10(c)], that also changes across the band. For example for the A-order, varies . Our nominal delay step produces a phase change on the red side of 0.102 cycle and blue side of 0.132 cycle. After 10 steps they will be mismatched by 0.3 cycle, which is significant.
Lastly, mechanical vibrations, air convection, and atmospheric seeing can potentially alter step delay and exposure flux, and smear time averaged fringe visibility, so for most generality we assume that the phases and visibilities are unknown and irregular.
Fitting Along a Column Assuming Y-Knowledge
Sine fitting versus vector fitting
For each (wavenumber) column of Fig. 8, the intensity versus phase is essentially fitted to a sinusoid. The sine and cosine amplitudes represent the imaginary (blue) and real (red) parts of the complex spectrum [Fig. 11(b)], called “whirl,” which represents the moiré pattern. The vertical offset to the sine fit represents the contribution of the nonfringing ordinary spectrum , also called the native spectrum [Fig. 11(a)].
In the case of uniform visibility fringes and knowledge of phase step positions, a standard sine fitting routine can be used. But if the visibility is nonuniform [as in Fig. 12(a)], such as under atmospheric seeing induced intensity variations, then a standard sine fitting routine can give an erroneous result. (Magnifying a weak signal also magnifies the photon noise, so fixing severe drops in fringe visibility only through normalization is not recommended.)
Rather than sine fitting, we use a vector combination method instead [Fig. 12(b)], which has greater flexibility in handling irregular fringe visibility and irregular phases, and has the ability for massively parallel computations over many columns simultaneously—useful for rapidly processing the wide bandwidths of our data. To separate the fringing and nonfringing component, the vector sum must be zero, and weightings are applied to achieve this.
The method of choosing the optimal weightings is discussed in Sec. 5.3.3 on “healing.” Notably, a “clumping” technique is used so that if some of the vectors have near zero magnitude, they do not spoil the value of the clump. So dropouts in signal due to a passing cloud, e.g., are tolerated well in our scheme.
It is useful to introduce a set of designating vectors (Fig. 13), whose phase and magnitude represent the interferometer comb phase and visibility for a given pixel or group of pixels. Whatever operation is performed on the designating vectors is also performed on the underlying data it represent.
We visualize Eq. (4) for the ordinary spectrum as a sum of vectors [Fig. 13(d)] such that the fringing component is cancelled because the vector sum is zero. We are allowed to weight the vectors. We find a combination of weights (there are many) that produces vector cancellation, and then apply the same weights and arithmetic to the data to produce the ordinary nonfringing spectrum.
Analogously, we visualize Eq. (3) for the fringing component as a series of rotations , tailored to each exposure, to align all the vectors to point in the same direction. These software rotations are ideally exactly counter-rotating to the delay-induced rotations. Now, it is the nonfringing terms which appear to rotate and vector cancel, leaving the purely fringing component.
Conversion of irregular to regular phase steps
Since there are more phase steps than degrees of freedom, there are many allowable weighting schemes that solve the phase-stepping vector problem. The scheme described below works well and is programmatically straightforward, but it is not a unique solution, nor necessarily universally optimal.
We handle irregular phase steps by first regularizing the data by an arithmetic recombination we call “healing,” which is a scheme that converts irregular to regular data using only adjustment of weightings and no rotations. Then to this ersatz regular data, we can apply any of a wide variety of methods designed for regular phase steps, such as the direct equations Eqs. (3), (9), and (10), or by applying a standard sine fitting algorithm column by column through the array.
We use a rule of thumb that we never rotate real data (e.g., a single exposure), but only adjust its weight. This is because real data contains at least two phasors of opposite polarity, one clockwise (negative frequencies) and the other counterclockwise (positive frequencies), as in Eq. (5) of Ref. 11, which is
Figure 13 shows the healing process that adds or subtracts other uniphase data from the same multiphase array, to each uniphase row. Fundamentally, this is allowed mathematically only if the nonfringing component has already been removed from each uniphase data, so that the so-modified data have a sinusoidal behavior exactly proportional to its designating vector, i.e., no intensity offset. So this healing process takes place after an ordinary spectrum is calculated (described in next section), then subtracted from each of the uniphase rows to produce a purely fringing array of
For each column (-pixel), the phase (and underlying data) is sorted so it is monotonically increasing. Then we create a regular set of goal vectors , evenly distributed about the circle having constant magnitude [Fig. 13(d)]. We compute the difference between the goal and nearest actual vectors and call these healing vectors [red vectors in Fig. 13(c)]. These represent how much fraction of other channel data we need to add to a given channel data to bring it to regularity. Note, we are not allowed to rotate a given channels data directly, and we can only modify it by adding or subtracting various amounts of other data channels.
It is best to select data vectors that lie in similar direction to the needed healing vector. We create three vector primitives, labeled , , [Fig. 13(b)], by averaging together those vectors pointing in similar directions, grouped in thirds of a circle [Fig. 13(a)]. The advantage of using clumps is that it avoids using a single data that by chance is weak and noisy.
Then we form each healing vector from a linear combination of two of the three clump vector primitives. The weightings are then calculated easily since all the angles and magnitudes are known.
Once all the arithmetic operations have been identified that convert the irregular to regular vectors, the same operations are applied to the data. This forms a corrected multiphase array having regular phase steps . Then one applies a sine fitting or Eqs. (3), (9), and (10) to form the whirl and ordinary spectrum .
Calculation of from three clumps
The used in Eq. (9) was calculated by weighting and then summing the data associated with the three clump primitive vectors [generalization of Eq. (4)]. With one vector held constant, the two weightings of the other vectors needed to produce vector cancellation can be found directly through algebra, since there are two equations ( and ) and two unknowns.
Delays relative to first exposure
The first exposure (bottom row of the multiphase array) is the designated representative exposure for the set, when finding the absolute delay , which would be , or outputting the SR spectrum. We may use and to denote the change offset from the first exposure, e.g., and .
(Out of habit, we may use term “phase step” generically where “delay step” is more proper. Fundamentally, it is the delay which is commanded to change and which is constant across the band. Also, “steps” may refer to delay position, but in other places refers to change in delay position.)
ThAr provides step and calibration info
We can determine the individual exposure’s relative from either stellar or ThAr spectra, but it is more accurate to get it from ThAr since it has narrow isolated lines that have very high visibility, near unity, rather than the few percent typical of the stellar whirl. So we process the ThAr-only fiber data first, and use that ThAr step information to process the other fiber containing StarThar. The presence of some starlight with the ThAr does not hinder finding accurate phase step.
Fitting Along a Row Assuming X-Knowledge
Fitting data to twisted whirl
The purpose of this section is to fit each row of 1-D uniphase data to a three component equation, where the unknowns are the relative delay change and magnitude of the whirl, plus the amplitude of the ordinary nonfringing spectrum.
Anticipating that the whirls for higher delays than the 0th are twisted, we form a twisted template whirl
An initial value for is used in a direct evaluation of Eq. (12) by the projection method described next (Sec. 5.4.2) to yield . We decrease our guess for by about and repeat. This converges rapidly to zero phase remainder and our output value for . (This iteration loop is inside the box labeled “Y-dependence discovery” of Fig. 9 and separate from the overall iteration loop.) This produces the two template spectra and B for the next section used to determine for every channel .
Projection method for determining individual step phase
Now that we have found templates and that characterize the set of channels, the next step is to find the phase , fringe visibility , and amplitude of ordinary spectrum of each individual channel (each exposure). This is done using a method of projections.
This method only works when the templates are approximately orthogonal to each other, meaning their overlap integral is small relative to the self overlap integral. This occurs when the bandwidth is sufficiently large so that and have dissimilar shapes.
The guiding equation Eq. (14) is re-expressed into a linear model having three components: two from the fringing and one from the nonfringing
We use overlap integrals to generate from Eq. (17) three equations for the three unknowns. Let the integral between two functions and be represented by the more compact expression
Note that because the equations include a third component for the nonfringing , these can be used either before or after the subtraction of from each , in Eq. (9). Also, it is okay if the templates and contain multiple components, such as (star + ThAr), since all the components are locked in phase with respect to each other, regarding how they respond to an interferometer delay change.
As seen in the flowchart Fig. 9, a channel-dependent normalization () divides the data to become more uniform against fluctuations such as passing clouds. The challenge is to learn the source intensity variation that is apart from the sinusoidal variation seen naturally from fringes. There are two methods of finding .
The first method uses the amplitude from the three component decomposition (Sec. 5.4.2). This method is the best choice when the spectral features are sparse, such as for the ThAr lamp. For ThAr, the second method of simple summation is much worse because there are too few features to have the statistical benefit of averaging. Then summation does not sufficiently diminish the net fringe visibility, so the sinusoidal nature “prints through” to the and produces a wobble-like effect on the whirl, which prevents perfect separation of ordinary and whirl components in the two outputs.
Initially, we also applied this first method to the Star + ThAr data, but discovered that it produced a slight residual of the ordinary spectrum in the whirl, which produced a ringing artifact. The remedy was to use a second normalization method, which is simple summation over all X-pixels. This method works well for the Star component because the number of features is high so the net fringe is washed away statistically. In the case of StarThar data, we temporarily mask away the bright ThAr lines prior to the summing.
Normalization problems lead to pollution of whirl with ordinary
We found that if the normalization was not working properly, this harms the phase-stepping math, and pollutes the whirl with a small fraction of ordinary component, causing the whirl to have a pedestal appearance [Fig. 14(c)], seen between the real and imaginary parts and having the same shape as the ordinary spectrum. This is not critical for the ThAr because it is large. But for the Star whirl which is only a few percent, a small fractional pollution of the larger ordinary component produces a significant absolute amount of pedestal.
This pedestal is a DC-like (zero frequency) component that gets upshifted during heterodyning reversal to contribute a ringing artifact [Fig. 14(a)]. This was most easily detected in places of the telluric spectrum which are heavily absorbing and thus should be relatively flat (zero intensity). When we eliminated the pollution of the whirl by changing normalization schemes, we got a beautiful ring-free result [Fig. 14(b)]. (We could also solve this by manually subtracting an adjustable amount of ordinary spectrum from the whirl, but that did not address the root problem.)
The pollution and ringing problem disappeared when we changed the method of computing normalization for the StarThar to the simple sum method (after masking away ThAr lines).
Normalizing challenge: Star twinkles but ThAr does not
The ThAr lamp intensity is steady with time, but the stellar intensity is not, due to atmospheric seeing fluctuations. For our present algorithm and StarThar data, we do not attempt the mathematically challenging task of normalizing both simultaneously. Instead, we first use the ThAr-only signal to provide the calibrations for phase step values and for equalization used for all signals. Then we process the Star-ThAr signal, normalizing the stellar component, since it is the weaker and more valuable signal. This is adequate.
Stage 3: Rebinning and Regularization
Stage 3, called rebinning and regularization, linearizes the horizontal axis, allowing comparison of the data to theory when plotting in mathematically more convenient wavenumbers (). The horizontal variable is rebinned more finely, to a uniform grid of so that it can hold the high frequency information of the eventual high-resolution output spectrum. We also remove the phase twist versus due to the changing refractive index of glass etalons comprising each delay. This regularization is applied to each of the eight etalon channels of the data, to form two output arrays sent to the next stage 4 (spectral reconstruction), which are a stack of whirls and a stack of ordinaries , eight high versus etalon #. Two arrays are thus produced for each of source: ThAr, StarThar, calculated ThAr, and calculated Telluric + ThAr.
The latter two sources are simulations of data which can be injected into the data analysis pipeline here at the beginning of stage 3, on simulated pixels. This allows comparison of data to theory both in detector pixel space and in uniform wavenumber space. The comparisons are used on ThAr data to compute the glass dispersion detwisting, and tweaks to the pixel to wavenumber calibrations.
The horizontal pixels are rebinned more finely , to a uniform grid of . The criteria are that half the reciprocal of the final grid spacing (called the Nyquist frequency) is 10 cm, and this needs to be larger than the largest delay we intend to use (4.6 cm of E8).
We found that a standard cubic spline does a better job than our custom-written FFT method, regarding avoiding any overshoot artifacts produced by the narrowest features.
Pixel to Wavenumber Map Calibration
During rebinning, we map the data from pixels to wave numbers using an improved calibration curve (map) which we fine-tune here. We use individual calibrations for fibers A and B that accounts for the slit slant between those fibers. We add to the original polynomial fit a line segmented “tweak map” which forces the map to exactly pass through each important ThAr line or other feature. We primarily use the magnitude of the whirl for the comparison because it is much cleaner than the ordinary spectrum (since the former is insensitive to FP noise). In sparse ThAr regions we use telluric features.
Bulk Shift of Map versus (Time or Etalon)
We find that the map has an etalon-dependent offset, which we call shear. This is due to a time-dependent drift, probably thermomechanical from the changed gravity vector after repointing the telescope. The magnitude of the shift is over 1 or 2 h.
After removing the monolithic shift in each ordinary spectrum, we notice that there remains a significant misalignment which changes its amount rapidly along the band, and even changes polarity. Examples are shown in the instrument noise section.
Lowpass Filtering to Prevent Nyquist Artifacts
Figures 15(a) and 15(b) show aliasing artifacts in data from etalon E5 (1.27 cm), that occurs when the detector pixel density is too sparse. Energy from the continuum fringe comb spike leaks beyond the Nyquist frequency (, or ), wrapping to the other frequency branch on the negative frequency side. This causes an erroneous signal.
Our remedy [Fig. 15(c)] is to restrict frequencies to , then apply a standard cubic spline interpolation routine to rebin. The latter has good behavior similar to a Hanning-shaped lowpass filter.
Detwisting the Whirl, Removing Glass Dispersion
The interferometer delay is a function of wavenumber due to the refractive index dispersion of the glass etalon. This creates a phase twist versus wavenumber, , associated with each etalon. Data whirls are untwisted via . Otherwise, the FFT will not produce a narrow spike [Fig. 16(b)] for the fringe comb.
For each etalon, we find [Fig. 17(a)] by comparing measured to calculated phases of ThAr lines. Fitting a polynomial to produces a baseline fit good to cycle in chunks of 20 to 200 pixels. The precision is further improved later to better than 0.01 cycle using the ThAr or telluric features of each target’s data (Fig. 18).
For the spectral reconstruction, we only need this phase curve. However, it is useful to produce a delay versus curve [Fig. 17(b)] by taking the -derivative. This is useful for detecting skipped integer fringes of phase, lost spanning band regions that lack features.
Fine Delay Determination
Due to our wide bandwidth, by observing the twist in phase along the band, we could determine the delay unique to the nearest fringe (nearest ), i.e., with no integer ambiguity. Then within that fringe we further determined the delay to a fraction of a fringe and removed the last fractional phase shift.
Due to mirror repositioning errors, the delays would change a few microns from night to night; we could also detect a delay difference of a fraction of a micron between A and B fibers. Since these were so small, we kept the original delay value and instead slightly adjusted the whirl phase.
The values of the eight delays in our quiver for 2010 September were for E1-E8 0.083, 0.34, 0.66, 0.96, 1.27, 1.75, 2.92, 4.63 cm, and 2.4 cm for E6.5, added in June 2011. (The latter was used to fill the gap between 1.7 and 3 cm, and swapped for the E1 delay on the eight-position filter wheel.)
Calculating Phase Difference Between Theory and Data
We have two methods for calculating the phase difference between theory and data. (1) The dot product method is automated and thus used initially. However, it has the drawback that it blindly uses the entire theory spectra without regard to any incorrect components to it. (2) Then we fine-tune the result using a manual Lissajous method which has the ability for the user to ignore components suspected to be incorrect.
Dot product method of phase evaluation
The method of dot products uses an overlap integral
We calculated simulated whirl data from theoretically supplied high-resolution spectra of our sources, such as ThAr from the Nat. Inst. Standards and Techn. (NIST) group20 or a telluric model from Roe.19
Fine tuning phase and , alignment via Lissajous
A fine tuning alignment of data to theory uses a Lissajous presentation to compare two wavelets for a small region of the bandwidth, data versus theory. This style of presentation has the advantage of showing small changes in both phase and somewhat independently. The user observes the Lissajous shape (Fig. 18) while adjusting phase and translations. We use a rainbow color change of the curve to indicate the approximate , so that features can be identified. Potentially this method could be automated but is currently manual.
An important advantage of this style of presentation is the ability for the user to disregard features in theory or data which are suspected to be unreliable. Examples include blended lines, containing multiple species such as both Ar and Th within the resolution element. We have discovered (Sec. 8.4.1) that our ThAr lamp differs from the NIST ThAr lamp—our Th species appears several times lower in magnitude, whereas our Ar species has the expected magnitude.
Figure 18 shows that a small misalignment in phase manifests a different appearance than a misalignment in . Thus both can be corrected independently. (a) In perfect-alignment each loop becomes a line. The phase misalignment (c) makes the individual loops more open and elliptical. The misalignment (b) tends to spread the loops out like a fan.
Stage 4: Spectral Reconstruction
Summary of Spectral Reconstruction
In stage 4, “spectral reconstruction,” the whirls are numerically shifted in frequency space (heterodyning reversal, Fig. 19) to undo the analog frequency shift imparted by the interferometer. This produces a wavelet signal for each etalon, covering a different spatial frequency range centered on its delay value. These wavelets are summed to produce a provisional output spectrum . This is equalized to remove the bumps in the net frequency response, so that the latter is Gaussian with a user-chosen final resolution, to form a final spectral output having minimal ringing. Comparison of ThAr data to theory for an isolated peak provides the equalization curve used for both ThAr and Star. The choice of final resolution is limited by the maximum delay which produces a contiguous range of delays, as .
Bell-Like Weighting to Reduce Net Photon Noise
Prior to heterodyning reversal, frequencies outside a set value (0.52 cm) were deleted, since they contain little signal. In addition, a bell-shaped weighting was sometimes applied .
For photon-limited noise, the ideal shape is that of the native sensitivity peak, but for instrument limited noise (the TEDI situation) there is little benefit. If we used a weighting, it was a Hanning function since it goes to zero absolutely at its edges (set to 0.52 cm).
Reversal of Heterodyning
Figure 19(b) shows the effect of reversal of the heterodyning, where the signal has been shifted in the direction that moves the continuum comb to zero. The effect is illustrated in frequency space, but is more accurately computed in the software in space as
Our convention has the continuum fringe comb starting at a positive frequency . Hence during numerical heterodyning reversal, Figs. 19(a) and 19(b), the signals shift to left toward negative . This convention is arbitrary since after the real part is taken (d) both frequency branches (c) have energy.
Wavelets Produced from Real Part of Shifted Whirl
The output spectrum we are reconstructing must be purely real. We take the real part of the heterodyne reversed signal (in space) to form what wavelet signalFigure 19 shows that for a small delay of E2 has similar frequency content as the native spectrum, for the whole A-order.
In practice, the above equation might actually require using the imaginary part of , depending on how the phase is defined. The test is that the continuum comb portion of (produced by a small delay) should be, after heterodyning reversal, mostly real and positive.
Summation of Wavelets to Form Provisional Spectral Reconstruction
The wavelets are summed to form a provisional SR spectral output called7.5.2).
Ordinary spectrum useful for photon limited performance
Each of the etalons produces an ordinary spectrum . We average them together to form a cleaner version . This is equivalent to an ordinary spectrum exposed over the total exposure time of all etalons. Hence we expect this to have times less fractional photon noise than the signal from the E1 delay which covers a similar range of low frequencies. Hence for a photon-limited situation, it is optimal to weight it stronger in the sum, by so we have
Ordinaries not included in TEDI output spectra
Initially, we included the ordinary native spectrum in the output spectrum, to provide the lowest frequencies, but we discovered for TEDI data we get cleaner results NOT including the ordinary spectrum. Instead, we rely on fringing from the low delay etalons E1 and E2 for the lowest frequencies information. (However, we still include it in plots of the SR output, for comparison.)
Figure 20 demonstrates in stellar data that the fringing-derived spectrum (red curve) agrees extremely well with the telluric model (black curve). In some cases, it is more accurate than the native spectrum (green curve), such as when FP noise is present. The EDI is extremely robust to FP noise because only changes between exposures contribute to the whirl.
Each native resolution element is independent
The heterodyning and SR process occurs independently for each native resolution element, massively parallel across the whole bandwidth of the instrument. Thus noise or a data defect such as a CR in a given resolution element only effects that resolution element.
In Fig. 6 of Ref. 15, similar to Fig. 20 here, we noted that the bumps of the conventional signal looked different than the model. Now we realize that this was due to our use of a lowpass filter on the conventional signal having a sharp cutoff frequency (i.e., rectangular) which introduced ringing. The lowpass filtering prevented aliasing artifacts when shifting the native signal to compensate for spectrograph drifts versus time. Since we now usually omit the ordinary signal in the SR sum we also omit lowpassing it.
Modulation Transfer Function Displays Net Frequency Response
The net frequency response of the EDI system is probed by a modulation transfer function (MTF), which is the ratio out versus in of complex amplitude, for each frequency
Using a narrow bandwidth for modulation transfer function
Evaluating the MTF over the full bandwidth of the order has the disadvantage of producing an MTF that has too many points (thousands) in its curve. We found it advantageous to evaluate the MTF for an extremely small bandwidth, centered over a single prominent feature such as an isolated ThAr line. We use a subset spectrum 100 final-pixels spanning . (By final-pixel we mean the rebinned points). We then scale this 100 point MTF over any size bandwidth up to the bandwidth of the entire order, which has 240 times as many points, by linear interpolation in frequency space.
Measured MTF of externally dispersed interferometry data
Red curve in Fig. 21 shows an MTF measured for a star. Note that the positions of the peaks (associated with each delay) are not ideal. We desire shoulder-to-shoulder coverage to produce a ring-free net lineshape. Unfortunately, the set of TEDI delays had a gap between E7 and E8 (3 and 4.6 cm), since TEDI delays were selected for Doppler radial velocimetry purposes, not spectroscopy, and there were only eight positions in the rotary etalon holder.
The presence of the 3 to 4.6 cm gap limits the effective resolution, because attempting to use a wide Gaussian MTF that extends over this region (such as the curve in Fig. 21) would produce ringing in the final lineshape in -space. If additional etalons filled the gap between 3 and 4.6 cm, a resolution of 45,000 at is possible.
For data taken prior to June 11, 2011, there was a gap at 2.4 cm between E6 and E7. This was remedied by a purchase of a new etalon, called E6.5, of 2.4 cm delay which occupied the position in the eight-position “filter wheel” previously held by the smallest etalon E1. (E2 can cover the range of E1.) Only a few exposures have been taken with the E6.5, during what turned out to be the last run of TEDI, June 11 to 19, 2011.
Equalization Applies Taper to Modulation Transfer Function
The purple curve in Fig. 21(b) is the equalization curve used to reshape the MTF into a smoother Gaussian. It is the ratio between the desired MTF (black curve) and the actual MTF (red curve)
The EQ is applied to the FFT of the data, via
The choice of goal resolution for the Gaussian in Eq. (33) is a tradeoff between resolution and ringing. In the presence of gaps in the delay coverage, increasing resolution increases the ringing.
Selecting etalon values in a future instrument
As the average increases, the width of each etalon peak in the MTF narrows, so a designer needs to select etalon values to produce contiguous coverage at the blue end of the native spectrograph bandpass. It is also beneficial to have enough overlap of etalons to produce two independent sets of SR results, one for the odd numbered etalons and the other for even numbered (using different equalizations). These can be used to confirm fidelity of the SR process, since it is unlikely an artifact would have the same appearance for different delays.
Example Reconstructed Spectra
Observations of many stars were performed with the TEDI interferometer producing reconstructed resolution as high as ( boost), using 3 cm maximum delay with TripleSpec having a native resolution of . Figure 22 is an example of resolution boosting. Figure 23 shows the extremely wide bandwidth we can achieve, covering all four orders of the native spectrograph from 4100 to . Figures 24 and 25 show closeups for the reddest order.
Figure 22 shows 10-fold boost in resolution (2700 becoming 27,000) observing mostly telluric features in the spectrum of star CrB along with ThAr calibration lamp emission lines, observed June 2011. Seven contiguous delays were used [as seen in MTF Fig. 21(a)], including a new etalon 2.4 cm E6.5 that filled a gap between 1.7 and 3 cm delays allowing a higher goal resolution without excessive ringing. The E1 (0.1 cm) position on the eight-position rotary holder was used for E6.5. [The E1 is not critical since low frequencies can be supplied by E2 (0.3 cm) or ordinary spectrum.] If additional delays were available to contiguously fill the gap between E7 and E8 (3 to 4.6 cm), then the equalized output resolution could be as high as 36,000 [dotted curve of Fig. 21(a)] or more without significant ringing.
Figure 23 shows the wide bandwidth capabilities: star HD219134 measured simultaneously over the four orders, 4100 to . The EDI resolution scales with wavenumber (or as ), so the resolution is not constant across this wide bandwidth.
Figure 24(a) shows the whole A-order at resolution boost and Fig. 25 shows closeups. The black curve is the theoretical telluric + ThAr model. The data (red and green curves) are attenuated significantly on the red end of the spectrum, probably due to transmission loss of the native spectrograph.
Note how the broad features of the red and green curves follow each other accurately, even though the red curve is obtained entirely from the fringing signal component and the green curve from the nonfringing portion. These are independent measurements extracted from the same data set. The nonfringing data were not included in the sum which forms the red curve.
Comparing Two Fibers on Same Star
Figures 26 and 27 show EDI measurements comparing the two fibers A (blue curve) and B (red curve) measuring the same source, their difference curve (gold), and compared to the telluric model (gray). The two fibers used different detector pixels in the direction and were exposed at times different because the role of the fibers was being alternated between Star and Star + lamp. Hence the fibers experienced different PSF distortions to the native spectrograph (see Sec. 9.4). Yet the resulting absorption lines yield reasonably accurate lines because it sums over several wavelets.
Measurement Comparing Star to Literature Spectrum
Figure 28 compares our EDI measurement of HD219134 with a spectrum obtained by another group with a conventional spectrograph, the NIRSPEC at Keck Observatory.21 All the stellar NIRSPEC features are seen in the EDI output, but many are unresolved in the TripleSpec native spectrum (green curve).
Measurement Comparing Two Stars: Unexpectedly Deep Lines and Doppler Shift Value
Figures 29 and 30 compare two stars measured on same night, star HD219134 (blue curve) and HD962 (red curve). We observed the latter star because the SIMBAD database described it as a featureless A-type, and thus could indicate the telluric component of the spectrum. However, our observations disagree with SIMBAD in two ways for HD962, while agreeing for HD219134.
First, the depth of narrow stellar absorption features of HD962 is similar to that of HD219134, which is a K-type. So the designation of A-type for HD962 is suspect. This is seen clearly in Fig. 31 which compares the native spectra of HD962 to several A- and B-type stars (HR1178, Gam Oph, HD7891) and the K-type star HD219134.
All were measured on the same night September 19, 2010. Strikingly, HD962 has more similarity to the K-type star than the A- or B-type stars, in that it shares the narrow stellar absorption lines at , which are absent for the three telluric calibrator stars. Interestingly, the three telluric calibrator stars all show (Fig. 31) a broader absorption feature at not seen in the K-type, and the dispute star HD962 also shows this, although weakly.
Second, we observe a relative Doppler shift between the stellar lines of the two stars, but it is of the wrong polarity. According to SIMBAD, we should be seeing HD962 blue shifted relative to HD219134. Instead we see HD962 redder. The Doppler velocity of HD219134 that we measure agrees perfectly with the SIMBAD database when the barycentric motion is accounted for. The two stars were measured within an hour of each other, using the same data analysis procedure, and the telluric and ThAr features match up perfectly (see dips at 9188 and ). We checked the RA and Dec values of the FITS file header to confirm the telescope was pointed at the right star. Thus we believe the SIMBAD database for HD962 is in error.
Some details of the Doppler calculation–Solar spectrum from FTS measurements at NIST are shown shifted by velocities 5.7 and to match the position of absorption lines. We observe geocentric blueshifts of 5.7 and for HD962 and HD219134, respectively. Given Earth’s motion of 11.679 and toward HD962 and HD219134 at UT times 10:49 and 9:50, respectively, this implies heliocentric blueshifts of 5.7 to for HD962 and 26.4 to for HD219134. SIMBAD blueshift values are 23.0 and for HD962 and HD219134. Thus they agree nicely for HD219134, but are off by for HD962.
High Dynamic Range and No Fixed Pattern Noise
Figure 32 shows the extremely high dynamic range and robustness to FP noise of which EDI is capable. These are of our ThAr lamp to a boosted resolution of 19,000 in the D-order. We resolve weak peaks of heights of the brightest line in that order. (Note the weak but distinct peaks at 9171 and , for example.) The data are normalized so that the brightest line of the order, at , defines unity height. The noise at the base of the 9188 line (used to calibrate the equalization) is a dynamic range of , which is excellent. There is very little ringing ().
This data also illustrates EDI’s robustness to FP noise. The inset is a small section of the 2-D raw detector image at the D-order associated with this data, showing a couple of bad pixels. These produce erroneous bumps in the native spectrum (green curve) at 9174 and . The red EDI output spectrum does not react to these bad pixels, because they are constant and thus do not respond to the sinusoidal interferometer phase dithering.
Individuality of each ThAr hollow lamp
Interestingly, we see a discrepancy between the measurement of our lamp, and the NIST measurement20 of their lamp, in the very weak lines (heights ). Using the line identifiers (purple letters) provided by NIST, this is explained by the Th of our lamp being weaker to the Ar lines by , compared to NIST. This is consistent with the reports20 that each ThAr hollow cathode lamp can be very individual in character.
This stresses the importance of measuring each ThAr lamp and demonstrates that EDI is scientifically useful for this.
We now understand earlier noticed discrepancies in fringe phase between data and “theory” (NIST measurement) for some weak features, such as during calibration of our glass etalon. For blended lines, changes in relative amplitude of components (Th relative to Ar) can create phase errors.
Instrument Noise Behaviors
The TEDI data set is dominated by instrument noises at . Photon noise, at , is about 30 times smaller for our typical exposures. The detector also suffers from occasional bad pixels which give false peaks on the native spectra (as in Figs. 32 and 33). Figures 34–38 describe the extreme PSF drifts that afflict the native spectrograph. The EDI is robust to all these types of insults to the native spectrograph. This is evident by the high-resolution spectra it produces, which would not be possible with a conventional spectrograph under the same conditions. We will also explain the theoretical basis of the robustness.
Fixed Pattern Noise Rejection Advantage
Additive fixed pattern noise
Figures 32 and 33 are examples of EDI’s extraordinary rejection of FP noises—only intensity changes synchronous to the phase dithering contributes to the whirl. Any constant intensity pixel will cancel in Eqs. (3) or (5) for . Even if the user decides not to increase the resolution, i.e., a boost of , the EDI result can be cleaner than the native spectrum. (If the photon noise dominates over FP noise, the native spectrum information should be included with the EDI fringing information.)
Multiplicative fixed pattern noise
If the data are not properly flatfielded, there can be a pixel-to-pixel gain variation artifact that changes the magnitude of the whirl. [However, the phase is not affected because of the constant error term cancels in Eqs. (3) or (5) for discussed above]. This is the same magnitude error imposed on the ordinary spectrum. For different delays, this magnitude error pattern will get shifted to different final frequencies. For Doppler velocimetry, a purely magnitude error is less harmful than a purely phase error, since the location of an isolated spectral feature is mostly determined by the phase of the wavelets.
Severe Distortions of Native Spectrograph
Figure 34(a) shows that high-resolution spectra can still be produced in spite of native PSF drifts that are extremely large and irregular—locally as much as as shown in Fig. 35, and averaged over in Fig. 37. Figure 35 is a stack of the eight native spectra showing drifts along the dispersion direction. There is a bulk drift versus time (Fig. 36) which we call shear, which is removed prior to heterodyning reversal. Note the significant residual irregular component, which even causes the shear polarity to flip over small wavenumber distances. Figure 38 shows how this horizontal noise is equivalent to an intensity noise of .
Robustness to Point Spread Function Widening
The EDI is very robust to both changes in the PSF position and width, such as due to atmospheric seeing disturbances or thermomechanical drifts. Let us first discuss the effect of an unexpected change in the width of the native lineshape. For concreteness, let us assume the net focal spot size is a sum in quadrature of the native slit width and some fixed additional blur, fixed in detector space, and that we have a boost so that the slit width of the EDI native spectrograph is wider than the conventional spectrograph it is trying to match.
This is numerically simulated in Fig. 39. The result is that EDI lineshape (red) is dramatically more robust than the conventional spectrograph (black) because EDI data are sensed at spatial frequencies defined by a physical delay, not controlled by focal spot size or slit width. The conventional spectrograph lineshape (black solid to dash) doubles and its MTF is halved in width. Whereas the effect on the EDI MTF is just to have slightly narrower bumps [solid to dashed red curves of (b)]. Since the original bumps are equalized away, the EDI final resolution is mostly unchanged by the blurring.
Robustness to Point Spread Function Drift of Native Spectrograph
Translational reaction coefficient for single delay externally dispersed interferometry
Now we will discuss, at length, the EDI robustness to an unwanted translation or “insult” (in wavenumber units) of the native spectrograph PSF along the dispersion direction on the detector, producing a perceived shift in the final output spectra. We define a translational reaction coefficient (TRC) . For a conventional dispersive spectrograph, .
For a single delay EDI, the TRC is smaller than unity by a factor of several, and can be zero, depending on the shape of the spectrum where TRC is being evaluated. For example, at an isolated narrow feature such as a ThAr line, it is zero. This is because the moiré pattern has a phase which has no slope along the dispersion direction, so translating the moiré pattern does not change the perceived phase, and hence perceived .
However, for a pair or more of lines that are partially blended together by the native blur, there can be a slope to the moiré pattern. This slope, acting against a translation , can produce a phase change which leads to a nonzero TRC (but still much smaller than unity by a factor of several). How TRC behaves on average for an arbitrary spectrum is best seen in Fourier space, by multiplying the peak shape at a delay of by a line running through the center of the peak [Fig. 40(b)]. A similar line runs through the native spectrograph in (a) to calculate its reaction to . We have colored green similar regions of in the panels to suggest how the net reaction for the native spectrograph (a) is several times that of the EDI peak in (b), because the location of multiplier line is at much higher for the EDI, because the heterodyning effect moves it, along with the peak, to be located at .
Figure 40 is notional and does not yet multiply the peak shape against the line. That is done in Fig. 41(a). This latter graph actually depicts the net reaction (bold black curve) for multiple delays, but if you look at the reaction curve under one of the two highest delays, E7 or E8, they are isolated, and thus depict what happens for a single delay. This net reaction curve (bold black) in Fig. 41 is a phase change in space, which manifests translations in space. The black and gray dashed slopes depict how that phase would change for and , respectively. Note that the reaction curve for E7 at 3 cm is roughly below the slope, and would be under for the highest delay E8 at 4.5 cm. The TRC scales as for a single delay, and would be of order 0.15 for the smaller delay E4 at 1 cm.
Translational reaction coefficient for multiple delays
We find that with the use of multiple delays, and without any special lineshape reshaping, we can reduce the per by roughly a factor of 20. The idea is that having multiple overlapping peaks will cause the negative reaction on the left side of one peak to cancel with the positive reaction on the right side of a neighboring peak. This is notionally shown in Fig. 40(c), and calculated for actual TEDI peak positions and original shapes in Fig. 41(a). The bold black curve is consistent with in the region of the overlapped peaks E1 to E6, 0 to 2 cm. This agrees with the simulation of Fig. 42, even though this was done in a different manner—it calculates movement of a specific feature rather than estimate for an arbitrary spectrum, and does it in wavenumber space rather than Fourier space. It follows a ThAr line in TEDI data and computes when all etalon rows of the moiré data are artificially shifted by to the left. The result is a shift, for a . This means the Doppler noise would be smaller than if using the native spectrograph alone (if is the only noise insult type).
New Crossfading Method to Cancel Point Spread Function Drift
The factor of improvement in the robustness to a PSF shift is remarkable, since it can potentially reduce radial velocity Doppler noise by the same factor (if is the only noise insult mode). What is even more remarkable is that with a simple reshaping of the lineshapes, we can theoretically achieve . We have realized this only recently. We present the following theoretical explanation and simple simulation below.
General Lineshape for Two Crossfading Peaks
To engineer perfect cancellation of the net reaction of overlapped peaks, the shape and width of the peaks are modified into a so-called “crossfading” shape. The original psf is reshaped to a new psf by multiplying the data in Fourier space with a multiplier , as in Fig. 43(a) from a typical peak with gradual wings into a triangle crossfading shape having a well-defined width. For the dual peak crossfade discussed here (mixing more peaks is possible), the most important consideration is that the edge of the (called a “foot”) goes to zero magnitude at the center position of the neighboring , and stays zero beyond that.
From the isolated peaks E7 and E8 of Fig. 41(a), one can see that each peak produces an odd symmetry sine-like contribution to the reaction curve, with the right side being positive and the left is negative. The new idea is to reshape the peaks to overlap so that the positive contribution of one peak perfectly cancels the negative contribution of the neighboring peak. This is shown in Figs. 43(b) and 43(c) for the simple case of a triangle lineshape. Between the two peaks, i.e., between and , a crossfading is said to occur.
There are a family of curves that satisfy this condition of producing in this crossfade region between the peaks. Let us derive the equation that describes them. The net reaction to a horizontal PSF drift is the product of the line through the center of the peak, times each peak. For peaks this is
Case of Triangle Lineshape
The simplest case is a triangle lineshape, separated by distance 1 and height 1, which thus has linear crossfading of unity slope between the peaks. Then the right side of the left peak is , so the first term of Eq. (36) is , and the left side of the right peak is so the second term is . Hence Eq. (36) is satisfied since , and in the region between the two peaks .
We apply this crossfade reshaping to every pair of overlapping peaks to fill a contiguously overlapped region entirely with . For the case of irregular , the peaks will be asymmetrical, treated as half-peaks, in order to always place the foot on the center of the neighboring peak and vice versa. That is, the left half of peak at is the reflection of right half of peak at , and these can be a different shape than neighboring pairs of half-peaks. So one should view a manifold of peaks as instead a manifold of crossfading pairs of half-peaks.
Other Crossfading Lineshapes
Colleague Eric Linder pointed out22 that a sinc function, , for just the region between its first zeros, also satisfies Eq. (36). This produces , which has the required property that shifting by 1 flips its polarity. This has the advantage of a rounded peak that avoids some possible ringing produced by the sharp tip of a triangle. Interestingly, the sinc function (but with its many satellite peaks) is the Fourier transform of a rectangle and hence is naturally the psf of a slit opened much wider than its spectrograph focal blur. We are studying the utility of other crossfading shapes that satisfy Eq. (36).
Crossfaded Multiple Delays Produce TRC0
Figure 41(b) is the same calculation of net reaction as (a), but using sinc function peaks instead of original peak shapes, and having asymmetrical shapes to satisfying the foot to center condition for the irregular delay spacing. Remarkably, this produces over the contiguously overlapped region from E1 to E6 (up to 1.8 cm). This means that for resolutions that only need 0 to 2 cm of frequency ( boost, or in the A-order), the interferometer is completely decoupled from the native spectrograph and not disturbed by changes in the native PSF location. (In principle, three new etalons at 2.3, 3.4, and 4.1 cm to fill gaps between E6, E7, and E8 could be added to extend the contiguous delay region to 4.6 cm, greatly increasing resolution for operation to ).
Crossfading can include the native peak, but that is not done here since the smallest delay E1 provides a cleaner signal of the lowest than native spectrograph peak (robust to the FP noise seen in the native).
Crossfading has not yet been applied to bulk TEDI data. How well one can model the actual peak shape, and how much it varies across the band is under study.
Simple Crossfading Simulation Using Two Delays
Figure 44 shows a simple demonstration using a synthetic source spectrum consisting a group of evenly spaced lines (a), so that a reasonably narrow frequency range around 1.4 cm is produced, so that only two peaks that straddle this frequency, at 1.2 and 1.6 cm, are needed. The source spectrum and moiré pattern are calculated twice, one in the original position, and again after artificially shifted by to the left.
A key point is that the moiré pattern changes its slope (phase versus wavenumber) when the delay passes through 1.4 cm, from 1.2 cm (b) compared to 1.6 cm (c). Since the moiré flips its slope polarity, the idea is to mix the two results to produce a net cancellation to the insult of a translation. The mixing happens not to the data in moiré form, but after it has been processed into an output spectrum (d). Both a triangle and sinc function crossfading shape were tested (sinc results are shown), and the triangle produces slightly more ringing.
The average of the 1.2 and 1.6 outputs form the final outputs shown in (e), for original (blue dash) and forcibly shifted (black) cases. These two cases overlap in position extremely well, producing a perceived shift of only for a , which is remarkably small.
Correlated Point Spread Function Drifts Needed
The above crossfading method presumes that the same PSF drift applies to the two overlapped delay peaks. Otherwise the cancellation will be incomplete. This is naturally satisfied by long-term drifts in a TEDI style apparatus which measures delays sequentially on the same pixels using different glass etalons that are rotated into the interferometer beam, where “long” means longer than the time to measure each delay (typically several minutes).
We propose that shorter term drifts would be more correlated between peaks using a stepped mirror (Fig. 45) which can be designed to manifest several delays simultaneously from the same source, although on different detector pixel rows. Then comparison of original (having shift) to crossfaded data (having no shift) yields information on the PSF shift , which thus could be at short times scales as well as long. However, this would be susceptible to any row dependence to the drifts (such as thermally induced on the detector, and thus likely to be long time scale).
By combining in the same apparatus a stepped mirror with a sequentially changed etalon similar to TEDI, information on short timescale variations from the stepped mirrors but on different pixel rows could be joined with that sequentially measured on the same rows. We speculate that this could produce a near zero TRC apparatus for Doppler velocimetry or high-resolution spectroscopy that was robust both on long and short timescales. The details of how to combine the data to best tease out the desired information are still under investigation.
Implications for Precision Radial Velocimetry
In light of the great advantage of having at least two overlapping delays to cancel the effect of PSF drift, we encourage users of single delay EDI for velocimetry to consider modifying their apparatus to measure two or more delays quasi-simultaneously, such as in Fig. 45. Or if their existing single delay overlaps their native peak, they can apply crossfading to mix native and single delay signals, for frequencies up to but not beyond the delay. The decrease in TRC should directly reduce a horizontal (wavelength) noise, translating into reduced Doppler velocity noise if PSF translation is the dominant type of insult, and if not dominated by photon or read noise. (Example noise sources not addressed by the TRC could be a telluric contribution that is time varying and blended into the spectrum, and spatial drifts in the reference ThAr beam relative to the stellar beam, causing slight differential phase shifts in the interferometer.)
Discussion and Conclusions
The TEDI is dominated by instrumental noise of about 3% due to severe drifts of the native PSF, which would otherwise prevent high-resolution spectroscopy without the interferometer. The photon shot noise is about 30 times smaller at about 0.1%. A companion paper23 describes an EDI theory of photon noise for various configurations of multiple delays and final resolutions.
An important ability of EDI is its great robustness to distortions to the native spectrograph instrument lineshape PSF. This robustness is valuable because in many high-resolution spectroscopy applications, particularly radial Doppler velocimetry, the limiting noise is not photon noise but fluctuations produced by instrumental irregularities. This is why many radial velocity spectrographs are put in vacuum tanks, built with thermomechanically stable materials, and fiber optically scrambled (in spite of the flux loss at fiber injection). The deleterious effects of instrument noises on the conventional spectrograph are so significant that it is worth the costs of these remedies.
The instrument error that results is a product of two factors: the magnitude of applied insult (e.g., ) times the sensitivity (e.g., TRC) of the instrumental technique to such insult. Optimally, both of these factors should be reduced for a precision radial velocimeter or high-resolution spectrograph, but we have focused our attention in this project on the latter (sensitivity) factor. Previously, it was not possible to reduce the TRC of a conventional dispersive spectrograph—it is fixed at unity, so naturally the effort by the community was put into reducing the former factor, such as by environmental stabilization and fiber scrambling. Remarkably, we have shown using EDI with multiple delays we can reduce TRC by at least a factor of 20 using our standard processing, with the promise of an order of magnitude further reduction of TRC possible with strategic reshaping of the data in a so-called crossfading method.
Disadvantages of Externally Dispersed Interferometry Implementation
Disadvantages of implementing EDI include increase number of optics that the beam must pass through or reflect, and thus increases the parasitic loss, which decreases the flux. This is only relevant for photon-limited applications. We note that researchers have passed starlight through fiber scramblers24 to reduce focal spot mode variation, while accepting some flux loss such scrambling entails. The increase in PSF stability is worth some flux loss, because it is critical to attaining precision high-resolution spectroscopy or Doppler velocimetry.
Another EDI disadvantage is that there are two interferometer outputs, which (unless polarization is used) necessarily travel on two separate beam directions. So additional optical engineering is required to detect both outputs, otherwise a 50% flux hit is suffered. This could be with optics that redirects the two outputs to different portions of the same spectrograph slit, or to two different spectrographs. For TEDI version 1 both outputs were detected,6 but for version 2 the complementary output was found2 more useful as a monitoring port to maximize the amount of flux entering the input fiber.
Third, additional read noise is generated because multiple exposures (at least three) are needed. However, this disadvantage may often be immediately outweighed by the advantage of FP noise rejection that such phase stepping creates (where two exposures recently taken are compared for differences), since FP noises can often be more significant than photon noise. (Read noise behaves like detector noise, which is described in the companion paper.23) Also, a recent paper9 proposes using a readout noise-free detector (EMCCD) with fast scanning of interferometer phase to improve Doppler precision.
Advantages of Externally Dispersed Interferometry
Summarizing, the EDI method offers three significant advantages: (1) resolution increases beyond classical limits (Fig. 46), and (2) dramatic robustness to PSF irregularities (Figs. 41 and 42), and (3) superb rejection of FP noises such as from bad pixels (Fig. 33). While the resolution increase tends to be discussed most often, what may prove to be the greatest advantage is the robustness to PSF irregularities, since this is the most serious practical limitation to many precision Doppler velocimeters or high-resolution spectrographs in astronomy today.
In summary, the EDI allows one to make high-resolution measurements not possible or practical with your current apparatus due to physical laws of optics, such as a nonzero focal spot size, sparse pixel spacing, bad pixels, and severe drifts and irregularities of the focal position and diameter. The bandwidth is set by the native spectrograph, which tends to be wide because it has lower resolution. Hence the bandwidth resolution product can be several times larger than in a conventional spectrograph, which is limited by the number of pixels. Spacecraft and aircraft have severe weight restrictions that make large spectrographs impractical. Since the cost and weight of a dispersive spectrograph grows in a power law with resolution, any technology that reduces the native spectrograph size by a factor of several to achieve a given resolution goal, offers a significant opportunity.
This material is based upon work supported by the National Science Foundation under Grant Nos. AST-0505366 and AST-096064, NASA Grant No. NNX09AB38G, and by Lawrence Livermore Nat. Lab. under Contract DE-AC52-07NA27344. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Thanks to Ed Moses for his valuable support during the genesis years. Thanks for high-resolution spectra from Henry Roe (UC Berkeley), and Florian Kerber (NIST) and Kevin Covey (Lowell Obs). Thanks to Jason Wright, Daniel Harbeck, Alex Kim, Ron Bissinger for astrophysics guidance and applications. Thanks to students Sam Halverson, Tony Mercer, Danny Mondo, and Agnieszka Czeszumska. Thanks to TripleSpec PI Terry Herter and Cornell staff Charles Henderson and Stephen Parshley. Thanks to Palomar Observatory and UC Berkeley Space Sciences staff including Mario Marckwordt, Michael Feuerstein, and Gregory Dalton. Thanks to Eric Linder for recent analysis and discussions on optimal EDI lineshape for enhanced stability.
David J. Erskine received his PhD from Cornell in 1984. He is an experimental physicist at Lawrence Livermore National Laboratory with experience in femtosecond lasers, semiconductor physics, superconductivity, diamond anvil cell high-pressure physics, shock physics, high-speed recording techniques, Doppler interferometry, white-light interferometry, digital holography, Fourier signal processing, image reconstruction, and phase-stepping algorithms for interferogram analysis. He collaborates with astronomers on interferometric techniques for Doppler planet search and high-resolution spectroscopy. He is a member of SPIE.
Edward H. Wishnow received a PhD in physics from the University of British Columbia. He is now a research physicist at the Space Sciences Lab at UC Berkeley. He is working on stellar interferometry and spectroscopy in the mid-infrared and visible.
Martin Sirk has over 30 years experience in the design, construction, calibration and science analysis of astronomical instrumentation. This has included working with Digicon detectors on the Hubble Space Telescope, CCD detectors on ground based telescopes, micro-channel plate detectors and optics on 6 NASA missions (EUVE, FUSE, ORFEUS, CHIPS, SPEAR, and ICON), and photographic plates at Lick Observatory.