CORRECTING FOR TELLURIC ABSORPTION: METHODS, CASE STUDIES, AND RELEASE OF THE TelFit CODE

All ground-based astronomical spectra suffer from contamination by the Earth's atmosphere, which introduces so-called telluric lines into the spectrum. The wavelengths of the telluric lines change very slightly with wind along the line of sight and pressure shifts in water lines, and the relative line strength can vary a great deal with the observatory location, weather, and the airmass of the observation. The telluric lines must be removed from the spectrum to retrieve many stellar features redward of about 500 nm, which is a nontrivial task.

Often, telluric lines are removed by observing a rapidly rotating hot star (spectral type A or late B) near the same time and airmass as the target star and using it as a template for the telluric spectrum. This approach has several disadvantages: (1) it can be difficult to find a suitable standard star near the same time and airmass as the observation, so approximate corrections are made to the line strengths using Beer's law (Beer 1852); (2) the water vapor content changes much more rapidly than the other absorbing species in the atmosphere, so scaling the whole empirical telluric spectrum is incorrect; (3) the standard star has strong hydrogen and helium lines and weak metal lines that can further contaminate the science spectrum; (4) observing a standard star can take a significant amount of precious telescope time, especially for high signal-to-noise work.

A better solution is often to generate a theoretical telluric absorption spectrum from a line list and the observing conditions. Recently, several groups (Seifahrt et al. 2010; Bertaux et al. 2014; Husser & Ulbrich 2013; Cotton et al. 2014; Gullikson & Dodson-Robinson 2013) have used the LBLRTM code¹ (Clough et al. 2005; Line By Line Radiative Transfer Model) for this purpose. However, the interface for the LBLRTM code can be difficult to learn, and can not directly fit an observed spectrum.

In this paper we describe Telfit,² a code that acts as a wrapper to LBLRTM and allows for easy fitting of the telluric spectrum in astronomical data. We compare the results to those obtained from empirical telluric correction in optical echelle spectra, and find that model-fitting produces similar telluric line residuals and is in fact better for dim targets or high signal-to-noise ratio work. We only demonstrate the use of TelFit in optical spectra in this work, but note that an earlier version of this code corrected telluric methane absorption near 2.3 μm to ∼5% of the continuum level (Gullikson & Dodson-Robinson 2013). Seifahrt et al. (2010) provide a case study of telluric modeling in several bands in the near-infrared using a different code, and find similar results.

We describe the observations and data reduction used in this paper in Section 2. We then describe the TelFit code and give a brief description of the telluric fitting procedure in Section 3. Finally, we demonstrate the use of the TelFit code to correct several optical telluric bands in Section 4.

We observed representative spectra for the early and late-type stars listed in Table 1. The A-stars HIP 20264 and HIP 25608 were observed with the CHIRON spectrograph on the 1.5 m telescope at Cerro Tololo Interamerican Observatory. This spectrograph has a spectral resolution of R = 80,000 from λ = 460 to 860 nm, and uses 3x1 binning along the spatial direction. The data were bias-subtracted, flat-fielded, and extracted using standard IRAF tasks. The reduced spectra were wavelength-calibrated using a ThAr lamp exposure from immediately before the observation of the star. In order to reach a high signal-to-noise ratio and avoid detector saturation, we co-added seven spectra after reduction and telluric correction (see Section 3).

In addition, several M-type stars were observed with the Magellan Inamori Kyocera Echelle (MIKE) optical echelle spectrograph on the Clay telescope at the Magellan Observatory as spectral type standards for a study of young stars (Kraus et al. 2014). We used the 07 slit, which yields a spectral resolution of R = 35,000 from λ = 335 to 950 nm. The pixel scale oversamples the resolution with the 07 slit, so we used 2x binning in the spatial and spectral directions. The spectra were reduced using the CarPy pipeline (Kelson 2003).³ The reduced spectra were wavelength calibrated using a combination of a ThAr lamp exposure and the 760 nm telluric A band (see Kraus et al. 2014 for details).

TelFit can handle rapidly rotating early-type stars without preprocessing the spectra by fitting a Savitzky–Golay smoothing filter (Savitzky & Golay 1964) as a pseudo-continuum for the telluric lines. For more feature-rich, late-type stars, we remove an approximation to the stellar spectrum using the procedure below before fitting the telluric lines.

• 1.  
  Make an initial guess telluric spectrum with TelFit. The guess spectrum intentionally has slightly higher molecular mixing ratios than the correct values so that it over-corrects the telluric lines.
• 2.  
  Divide the observed spectrum by the guess telluric spectrum. This leaves a spectrum with stellar lines in absorption and telluric residuals appearing like emission lines.
• 3.  
  Find the best-fit stellar model to the absorption lines using a grid of PHOENIX model spectra (Hauschildt et al. 1999) with effective temperature and log(g) near the expected values for the star.
• 4.  
  Divide the original data (before the telluric over-correction in step 1) by the best-fit normalized stellar spectrum. This leaves a spectrum that is mostly telluric lines, but which may have strong stellar line residuals.
• 5.  
  Enumerate the wavelength ranges that are strongly affected by stellar line residuals, and tell TelFit to ignore them. We ignored all wavelengths from 817.5 to 820 nm in our M-star spectra, which have strong residuals from the sodium lines.

The fit to the primary star does not need to be perfect. Indeed, if we could perfectly model stellar spectra there would be little reason to observe them! The purpose of the preprocessing described above is to ensure that telluric lines dominate the spectrum before fitting. Any strong residuals from that process can and should be ignored in the telluric fit. While a physically better solution would be to fit the stellar and telluric spectra simultaneously, the method described above will work well enough for most purposes and is much faster than a simultaneous fit.

The TelFit code performs a least-squares fit using a constrained Levenberg–Marquardt algorithm (Marquardt 1963). It is extremely flexible; the user can decide which atmospheric parameters to fit, which spectral regions to use in the fit, how the detector point-spread function (PSF) is modeled, the order of the continuum fit used in each iteration, and whether the data wavelengths are calibrated to fit the model or vice versa. The algorithm will optimize the fitting variables, as well as the detector resolution and the wavelength solution. LBLRTM requires an atmosphere profile giving the temperature, pressure, and the abundance of several molecules as a function of height. We provide a default atmosphere profile⁴ with the TelFit code that is acceptable for mid-latitude observatories, but provide the ability to easily give an alternate atmosphere profile for more accurate results.

While TelFit explicitly fits the temperature, pressure, and abundances at the observatory altitude, it scales the quantities at all atmosphere layers. For the pressure and temperature, it finds the difference between the requested pressure (temperature) at the requested altitude, and the pressure (temperature) in the atmosphere profile. If P_r is the requested pressure at altitude z, and the atmosphere profile pressure is P₀ at that altitude, then the pressure P_i at each atmosphere layer z_i is scaled as

$$\begin{equation*} P_i \rightarrow P_i - (P_r - P_0) e^{\frac{-(z_i - z)^2}{2\cdot (10 \rm km)^2}}. \end{equation*}$$

The temperature is scaled in the same way as the pressure above, and is done in this way so that the quantities more than about 10 km above the observatory are effectively unchanged. The mixing ratio of all molecules is scaled in a much simpler way: if A_r is the requested mixing ratio at the telescope altitude, and A₀ is the mixing ratio in the atmosphere profile, the mixing ratio at each layer A_i is given by

$$\begin{equation*} A_i \rightarrow A_i \frac{A_r}{A_0}. \end{equation*}$$

Because of the way TelFit scales the temperature, pressure, and molecular mixing ratios in every atmosphere layer, the optimized values should not be taken as a measurement of the actual surface mixing ratios.

In each iteration for the main variables (temperature, pressure, telescope zenith angle, and molecular abundances), TelFit refines the wavelength solution and detector resolution. The wavelength fitter performs a third-order polynomial fit to adjust the wavelengths of the model such that the telluric lines match the data. Alternatively, the data wavelengths can be adjusted to fit the model. The wavelength solution of the data should be very close, as TelFit will not shift any wavelengths by more than 0.1 nm. By default, TelFit fits the spectrograph PSF as a Gaussian, and fits the detector resolution by finding the best Gaussian width to convolve the model against. Alternatively, TelFit can do a nonparametric fit to the spectrograph PSF using the singular value decomposition described in detail in Rucinski (1999).

We fit the spectra for the rapidly-rotating A0V star HIP 20264 with a two-step approach. In the first step, we fit the observatory temperature and humidity in the echelle orders dominated by the water bands at 590 nm, 650 nm, 700 nm, and 730 nm. For the rest of the fit, the humidity and temperature were fixed at the χ²-weighted average of the fits for each band. Next, we fit the O₂ abundance with the γ and B bands near 630 nm and 690 nm, respectively. Because the blue end of the B band is extremely strong and absorbs nearly all the light, we only use the red half in the fit. The O₂ abundance was fixed at the χ²-weighted average of the individual fitted values for the two echelle orders. We then applied the best-fit parameters to every echelle order, allowing the fitter to adjust the model wavelengths and detector resolution separately for each order. For all telluric fits in this paper, we adjusted the temperature, pressure, and water vapor atmosphere profile with sounding data from the Global Data Assimilation System (GDAS) meteorological archive.⁵ As stated in Section 2, we fit each frame of HIP 20264 separately, and co-added the telluric-corrected spectra. The fitted value of the relative humidity varied from 18.7% to 22.1%, and the airmass of the star increased from 1.15 to 1.52 (a change of 19° in zenith angle) over the course of the 7 frames, making the individual fits crucial.

We show the telluric correction for the observation of the A0V star HIP 20264 in four water bands and two O₂ bands in Figures 1 and 2, respectively. The telluric model shown is the average telluric model of the 7 individual exposures, since we corrected each one separately to better account for the changing water vapor content and telescope zenith angle. We apply a Savitzky–Golay smoothing filter to the corrected data to remove any broad features in the stellar spectrum. The water features are corrected to near the noise level of the spectrum, with the exception of the very strong telluric water band near 730 nm. The poor correction may be due to a slightly incorrect temperature and water-vapor atmospheric profile, which causes the line profile to change appreciably for strong lines. In addition, small errors in the line strength parameters are more noticeable for strong telluric lines. The correction in the O₂, B, and γ bands (Figure 2) is somewhat worse. The same O₂ mixing ratio undercorrects the γ band and slightly overcorrects the B band. This systematic error may be from incorrect line strengths in the HITRAN database for the two bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed the A1V star HIP 25608 30 minutes after HIP 20264 and at similar airmass, and so use it as a telluric standard star to directly compare our method to the empirical telluric correction method. To have comparable S/N ratios as in Figures 1 and 2, we co-added all seven frames of both the target star (HIP 20264) and the standard star (HIP 25608). We account for the different column density of absorbers in the target and standard star observations using Beer's law (Beer 1852) as follows. For each order, we scaled the normalized flux of the standard star (f_s) with

$$\begin{equation} f_s \rightarrow f_s^{\gamma } \end{equation} \tag{ 1 }$$

where γ is determined using the normalized flux of the standard star and target star (f_t) and is defined as

$$\begin{equation} \gamma = {\rm median} \left(\frac{\log {f_t}}{\log {f_s}} \right). \end{equation} \tag{ 2 }$$

In Equation (2), we only use the points with f_t < 0.95 in the median computation, so that the noise does not dominate in the orders with few telluric lines. These empirical telluric corrections are shown in the top row of Figure 3. The corrections are quite poor in this case, mostly because the airmass of the target (standard) star changes from 1.15 to 1.52 (1.31 to 1.86) from the first to last frame, and co-adding the spectra amounts to a flux-weighted average over airmass that the correction in Equation (1) cannot capture.

Zoom In Zoom Out Reset image size

We also test an empirical telluric correction using only the last frame of the target star and the first frame of the standard star, and show the result in the bottom row of Figure 3. In this case, the airmass can be treated as approximately constant and Equation (1) does a much better job in accounting for the small airmass difference between the target and standard star spectra, with the exception of the nearly saturated oxygen lines (bottom left panel). The water line correction results in similar line residual amplitudes to the telluric modeling method (see Figure 1 for comparison).

We now turn to a more typical-use case for the TelFit code: correcting the telluric lines in late-type stars near a feature of interest. For this, we use a series of M-type stars (see Table 1) with observations near the 819 nm sodium doublet, which is used as a gravity-sensitive age indicator for late-type stars (Slesnick et al. 2006). As the stellar age and therefore surface gravity increases, the line strength increases. Amongst main sequence M-stars, later spectral types have higher surface gravity and so we expect the equivalent width of the sodium lines to increase as we go toward later spectral types.

Figure 4 shows the telluric correction for the M1V star GJ 908. For this and all M-star spectra in this work, we fit the water vapor, temperature, and O₂ mixing ratio at the same time for the spectral order covering 819 nm. The bottom panel of Figure 4 shows that each of the spectral lines after telluric correction come from the star itself, and that the telluric contamination is reduced to near the noise level of the spectrum. Importantly, the telluric lines that fall within the sodium doublet line profile no longer affect the profile shape. We can now directly compare the line profile of the sodium doublet lines as a function of spectral type, without the influence of telluric lines. Figure 5 shows the evolution of the doublet lines. The central depth stays about the same throughout the sequence, but the later spectral type stars show significantly broader line wings, leading to the known sequence in equivalent width (Slesnick et al. 2006). The robust recovery of the sodium line strengths and profiles demonstrate that TelFit can accurately remove telluric lines, even in feature-rich spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We presented the TelFit code, an object-oriented Python code capable of accurately fitting the telluric spectrum in ground-based spectra. We use a high signal-to-noise ratio echelle spectrum of the A0V star HIP 20264 to demonstrate the fit quality in Figures 1 and 2. Contrary to expectation, the water lines are typically fit to much higher precision than the O₂ telluric lines. This is likely coming from a systematic error in the HITRAN line strength database, since the same O₂ mixing ratio overfits the B band and underfits the γ band. We compare our code to an empirical telluric correction of HIP 20264, and find that the model is at least as accurate. In fact, TelFit is significantly more accurate than the empirical method when several frames of both the target and standard star are co-added.

We also demonstrated the use of TelFit for in-depth study of spectral features in late-type stars with a series of M-star observations near the 819 nm sodium doublet. The telluric lines were removed to near the noise level of the observations, allowing for analysis of the sodium line profiles and recovery of the known sequence of increasing equivalent width with later spectral types. Regions contaminated by telluric lines are often ignored in optical spectral analysis; accurate correction of telluric features could help open these regions up for further analysis.

This code was mostly developed and tested for correction of optical spectra. However, an early version of this code was used in Gullikson & Dodson-Robinson (2013) to correct for telluric methane absorption in B-star spectra, with similar line residuals to those that Seifahrt et al. (2010) found in the same spectral region. We encourage the use of TelFit for correcting telluric absorption in near-infrared as well as optical data.

We thank Andreas Seifahrt for his help in the early stages of code development. This project was funded by a UT Austin Hutchinson fellowship to Kevin Gullikson and start-up funding to Sarah Dodson-Robinson from the University of Texas. This work makes use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), as well as SciPy, NumPy (Oliphant 2007), and of course LBLRTM. We would like to thank the developers of all of those packages.