A Robust Method to Measure Centroids of Spectral Lines

When analyzing data cubes produced by radio telescopes or integral field units, astronomers frequently aim to produce a map of velocity as a function of pixel coordinate or time. In some cases, the spectral dimension is dominated by a single line and the operation of inferring a velocity is simplified to a question of determining the centroid of the line. For this problem, there exist multiple popular methods (see, for example de Blok et al. 2008), however the most common are the "intensity weighted average velocity" (also known as the "first moment" method) and the velocity of the maximum intensity, called the "ninth moment,"³ as they are implemented by CASA (McMullin et al. 2007). Even when the emission is dominated by a single line component, these methods do not offer the accuracy or precision required to extend studies to small scale structures.

The "first moment" approach is strongly affected by the noise in the spectrum and any asymmetry in the line profile. Sigma clipping is typically used to circumvent these issues, however this can result in the addition of spurious features or the masking of astrophysically interesting features in the resulting velocity map. "Ninth moment" maps are more robust to these issues, however the precision of this approach is limited by the velocity resolution of the data.

Other approaches involve fitting a model line profile (e.g. a Gaussian) to the data, but these methods are both more computationally expensive and more sensitive to the details of the line shape. In the case of a double-Gaussian profile, due to radial infall, for example, the line center for a single Gaussian component (or a higher order expansion) will therefore be offset to account for this double component.

In this note we present a new method for inferring line centroids that is simple, computationally efficient, and robust to noise and errors in line shape.

It has been demonstrated that the centroid of photometric sources can be determined with near-optimal precision by fitting a quadratic surface to the brightest central pixels of the point-spread function (Vakili & Hogg 2016). Here we apply this method in one dimension to estimate the centroid of spectral lines. The fundamental idea is that we can fit a quadratic model to the brightest pixel and one pixel on either side to estimate the curvature of the line profile near maximum.

The procedure is as follows:

• 1.  
  We find the pixel of maximum intensity in the spectrum and extract that pixel value I₀ and the one to the left I₋ and right I₊.
• 2.  
  We model these three intensities as a quadratic function

  $$\begin{eqnarray}&&I(x)={a}_{0}+{a}_{1}\,(x-{x}_{0})+{a}_{2}\,{(x-{x}_{0})}^{2}\end{eqnarray} \tag{ 1 }$$

  where x₀ is the pixel coordinate of maximum intensity and x is in units of pixels. Given the intensities in the three target pixels, we find

  $$\begin{eqnarray}&&{a}_{0}={I}_{0}\end{eqnarray} \tag{ 2 }$$

  $$\begin{eqnarray}&&{a}_{1}=\displaystyle \frac{1}{2}({I}_{+}-{I}_{-})\end{eqnarray} \tag{ 3 }$$

  $$\begin{eqnarray}&&{a}_{2}=\displaystyle \frac{1}{2}({I}_{+}+{I}_{-}-2\,{I}_{0}).\end{eqnarray} \tag{ 4 }$$
• 3.  
  The maximum of Equation (1) is at

  $$\begin{eqnarray}&&x={x}_{\max }={x}_{0}-\displaystyle \frac{{a}_{1}}{2\,{a}_{2}}\end{eqnarray} \tag{ 5 }$$

  $$\begin{eqnarray}&&=\ {x}_{0}-\displaystyle \frac{{f}_{+}-{f}_{-}}{2\,({f}_{+}+{f}_{-}-2\,{f}_{0})},\end{eqnarray} \tag{ 6 }$$

  and this is our estimate for the centroid of the line in pixel coordinates.

We can also estimate the statistical uncertainty on x_max by linearizing and propagating the uncertainty from the fluxes to the centroid estimate. If we assume that the intensity uncertainties are normal, independent, and homoskedastic with standard deviation σ, we find the following approximation for the statistical uncertainty on x_max

$$\begin{eqnarray}&&{{\sigma }_{{x}_{\max }}}^{2}=\displaystyle \frac{{\sigma }^{2}}{8}\,\left(\displaystyle \frac{3}{{{a}_{2}}^{2}}+\displaystyle \frac{{{a}_{1}}^{2}}{{{a}_{2}}^{4}}\right).\end{eqnarray} \tag{ 7 }$$

We emphasize that this only includes the statistical uncertainty and any application of this method to real data will also require a treatment of systematic uncertainties.⁴

The top row of Figure 1 compares different methods for estimating a line centroid the case of an asymmetric line profile, modeled as the sum of two Gaussian components. Both the first moment and Gaussian fit results are biased by the asymmetry of the line. Although the "ninth moment" approach avoids this bias, its precision is limited by the velocity resolution of the data.

Zoom In Zoom Out Reset image size

The bottom row shows velocity maps of ¹³CO from the protoplanetary disk, HD 135344B (van der Marel et al. 2016, ALMA Project 2012.1.00158.S), imaged at a velocity resolution of 330 m s⁻¹. With a 2σ clip applied to the data, a "first moment" approach results in a noisy map, making it hard to distinguish real structures from noise. As with the asymmetric spectrum, the "ninth moment" is less affected by the noise but the resulting map loses fidelity due to the coarse velocity resolution of the data. The quadratic approach does not suffer from noise features while maintaining a high precision on the inferred velocity and achieving a statistical uncertainty of 47 m s⁻¹, or 15% of the velocity resolution.

The code used to generate this figure and an open source implementation of this method is available at https://github.com/richteague/bettermoments. Version v1.0 used in this note is available at https://doi.org/10.5281/zenodo.1419754.

Software: bettermoments (Teague & Foreman-Mackey 2018).