Modeling the Echelle Spectra Continuum with Alpha Shapes and Local Regression Fitting

TL;DR

This paper introduces two novel, data-driven algorithms for continuum normalization of echelle spectra, improving accuracy over polynomial fitting and addressing detector artifacts in high-precision spectrometry.

Contribution

The paper presents two new algorithms, AFS and ALSFS, for spectrum continuum normalization that outperform polynomial methods and incorporate lab reference data.

Findings

Improved continuum normalization accuracy over polynomial fitting.

Effective mitigation of CCD detector artifacts.

Validated on simulated spectra and real spectrometer data.

Abstract

Continuum normalization of echelle spectra is an important data analysis step that is difficult to automate. Polynomial fitting requires a reasonably high order model to follow the steep slope of the blaze function. However, in the presence of deep spectral lines, a high order polynomial fit can result in ripples in the normalized continuum that increase errors in spectral analysis. Here, we present two algorithms for flattening the spectrum continuum. The Alpha-shape Fitting to Spectrum algorithm (AFS) is completely data-driven, using an alpha shape to obtain an initial estimate of the blaze function. The Alpha-shape and Lab Source Fitting to Spectrum algorithm (ALSFS) incorporates a continuum constraint from a lab source reference spectrum for the blaze function estimation. These algorithms are tested on a simulated spectrum, where we demonstrate improved normalization compared to polynomial regression for continuum fitting. We show an additional application, using the algorithms for mitigation of spatially correlated quantum efficiency variations and fringing in the CCD detector of the EXtreme PREcision Spectrometer (EXPRES).

Tables

Table 1. Table 1: Unit: 1 × 10 − 3 1 superscript 10 3 1\times 10^{-3} . Medians of RMSE over 1000 samples for orders B and R with varying S/N for AFS, ALSFS, and the iterative method. Medians listed in this table are corresponding to distributions of RMSE in Figure 12 .

	Order B			Order R
S/N	AFS	ALSFS	Iterative	AFS	ALSFS	Iterative
300	7.58	4.31	12.44	4.06	2.70	8.10
150	7.71	5.67	15.79	4.80	5.31	8.25
50	8.94	9.93	22.78	8.07	13.82	11.49

Equations

$$\sum _{\begin{array}{c}{\lambda }_j\in N_{m_0}({\lambda }_i)\end{array}}{\omega }_i({\lambda }_j)(y_j-f_{{\beta }^{({\lambda }_i)}}({\lambda }_j){)}^2,$$

$$y_{corrected,l}=w_l{y}_{1,l}+(1-w_l)y_{2,l},l=1,\dots ,m,$$

$$threshold({\lambda }_i,t)=fi{t}_{polynomial}^{(t)}({\lambda }_i)-\frac{0.5}{t+1},$$

Full text

instrumentation: spectrographs; techniques: spectroscopic, radial velocities; methods: statistical, data analysis

††journal: AJ

1 Introduction

Spectroscopy is a powerful observational technique for understanding fundamental astrophysics. A high dispersion stellar spectrum contains detailed information about individual atomic transitions that enables the derivation of parameters such as effective temperature, surface gravity, elemental abundances, and the measurement of Doppler shifts in stellar spectra reveals the presence of stellar and planetary companions (Fischer et al., 2016, and references within). Most spectroscopic analysis techniques require a flat, continuum normalized spectrum (Blanco-Cuaresma et al., 2014). For example, the precision of equivalent width measurements for abundance analysis or cross-correlation for exoplanet detection is very sensitive to even small errors in continuum normalization (Torres et al., 2012, and references within).

An echelle spectrograph disperses light so that high spectral orders can be recorded. The higher orders have greater dispersion and therefore provide higher resolution spectra over a broad range of wavelengths. However, most of the brightness of a spectrum is concentrated in the zeroth and lower spectral orders; higher orders are intrinsically fainter. The optical grating in an echelle spectrograph has angled facets designed to shift the intensity envelope of dispersed light to high spectral orders. When this phase shift is introduced, the grating is said to be blazed and with cross dispersion, several dozen spectral orders can be stacked onto the detector. Each order has an intensity distribution characterized by the blaze function so that the continuum intensity is strongest in the center of the order and drops off steeply toward the edges of the order. The term “blaze function” is often used to describe the shape of the continuum across an echelle order and indeed, the blaze distribution is the dominant effect.

One approach for flattening the spectrum is to divide by the theoretical blaze function (Barker, 1984), which depends on grating parameters (incident angle, the angle of the blazed facets, and the grating facet width and spacing) and can be calculated for each order of the spectrum. However, this method for normalization or flattening of the spectrum will leave residual variations in the continuum because of departures from the theoretical blaze function: manufacturing defects in the grating, chromatic aberrations in the optics of a spectrograph, or quantum efficiency variations in the electronic detector. In addition to the instrumental blaze, additional wavelength-dependent intensity variations will occur because of the black-body temperature of stars or calibration lamps.

Another approach for normalizing the continuum is to divide by an extracted flat-field calibration source (e.g., Skoda et al., 2008), such as a quartz lamp. However, the quartz lamp will also have a black-body curve superimposed on the blaze function. One of the most common methods for normalizing the continuum of a spectrum is to fit a polynomial to high points along the order. This is a completely agnostic approach that does not require prior knowledge about the blaze function and while it works fairly well, polynomial fitting can fail near broad and deep lines, especially if they are close to the edges of the orders. Here, we describe a new approach for continuum fitting and compare our method with polynomial fitting in Section 3.

2 Description of methods

In this work, the goal is to estimate the blaze function of a target spectrum so that it may be removed, leaving a flattened spectrum. Two algorithms are proposed to accommodate different scenarios: (i) the Alpha-shape Fitting to Spectrum algorithm (AFS) (the baseline approach), (ii) the Alpha-shape and Lab Source Fitting to Spectrum algorithm (ALSFS) (when a lab continuum source is available)111The implementation code of AFS and ALSFS, as well as example data can be found and downloaded from: https://github.com/xinxuyale/AFS or https://zenodo.org/badge/latestdoi/173169370..

In the proposed algorithms, we first fit an alpha shape (Edelsbrunner et al., 1983), which is a polygon enclosing a dataset, to capture the general shape of the blaze function of a target spectrum. Then we use this preliminary estimation to select a set of pixels that are ultimately used to fit the final blaze function model. The pixels are selected so that they are generally on or near the continuum and not in the absorption lines. With the selected set of pixels, we use local polynomial regression (Cleveland, 1979) to estimate the blaze function. The two algorithms are introduced next, followed by a discussion on how to select the tuning parameters.

2.1 AFS Algorithm

The proposed AFS algorithm is a versatile algorithm that can be used to remove the blaze function whether or not a corresponding lab source spectrum is available. In this algorithm, an alpha shape is used to obtain a preliminary estimation of the blaze function’s high-level shape. An alpha shape is a generalization of a convex hull, but is not required to be a convex set; it is a region bounded by a set of segments generated from a set of points. To generate an alpha shape, imagine a piece of paper with a plot of spectrum printed on it. Consider a special paper cutter that can only cut out circles with radius $\alpha$, known as $\alpha$-balls. An incomplete circle is allowed, but the cutter cannot cut anything from the spectrum. In Figure1a, the blue circle is an $\alpha$-ball that can be cut out from the paper, but we cannot move the $\alpha$-ball any lower vertically since it would cut the spectrum. Continue to cut as many $\alpha$-balls as possible, and the remaining paper is called an alpha hull. Figure1a shows the resulting alpha hull with $\alpha=5$. By connecting the points where the alpha hull touches the spectrum with straight segments, it becomes an alpha shape, as displayed in Figure1b. An alpha shape can capture the general shape of a spectrum, and its upper boundary is used as a starting model for the blaze function estimation.

Another technique used in the AFS algorithm is local polynomial regression (Cleveland, 1979), which is a non-parametric method for estimating functions given a set of points. At each point, $(\lambda_{i},y_{i})$, a p-degree polynomial model is fitted to a subset of neighboring points of $\lambda_{i}$. The p-degree polynomial regression is fitted by weighted least squares, giving more weight to points closer to $\lambda_{i}$ and less weight to points further away. Consider a dataset $(\lambda_{1},y_{1}),(\lambda_{2},y_{2}),\dots,(\lambda_{n},y_{n})$, and a p-degree polynomial function in the neighborhood of $\lambda_{i}$ is $f_{\beta^{(\lambda_{i})}}(\lambda_{j})=\beta_{0}^{(\lambda_{i})}+\beta_{1}^{(\lambda_{i})}(\lambda_{j}-\lambda_{i})+\beta_{2}^{(\lambda_{i})}(\lambda_{j}-\lambda_{i})^{2}+\dots+\beta_{p}^{(\lambda_{i})}(\lambda_{j}-\lambda_{i})^{p}$, where $\lambda_{j}$ is in the neighborhood of $\lambda_{i}$. Then the estimation of $\beta^{(\lambda_{i})}$, $\hat{\beta}^{(\lambda_{i})}=(\hat{\beta}_{0}^{(\lambda_{i})},\hat{\beta}_{1}^{(\lambda_{i})},\dots,\hat{\beta}_{p}^{(\lambda_{i})})$, is obtained by minimizing

[TABLE]

where $N_{m_{0}}(\lambda_{i})$ is the neighboring set of $\lambda_{i}$ containing the nearest $\lfloor m_{0}n\rfloor$ pixels to $\lambda_{i}$, $m_{0}$ is a smoothing parameter, and \omega_{i}(\lambda_{j})=K\Big{(}\frac{|\lambda_{j}-\lambda_{i}|}{\max\limits_{\begin{subarray}{c}\lambda_{j^{*}}\in N_{m_{0}}(\lambda_{i})\end{subarray}}|\lambda_{j^{*}}-\lambda_{i}|}\Big{)}, with $K(x)=(1-x^{3})^{3}$ as a common weighting function. The local polynomial estimate at $\lambda_{i}$ is $\hat{\beta}_{0}^{(\lambda_{i})}$. An advantage of local polynomial regression over ordinary polynomial regression is its ability to adapt to local characteristics of a dataset rather than fitting all the data points using one single model. This flexibility makes it useful for modeling complex data sets where regular polynomial regression fails.

Let $\{(\lambda_{i},y_{i})\}_{i=1}^{n}$ be an observed spectrum, where $\lambda_{i}$ is the wavelength of pixel $i$, and $y_{i}$ is the intensity of pixel $i$. Our method is summarized in Algorithm 1 and illustrated in Figure 2 with the details of the AFS algorithm presented next.

In step 1, the intensity vector, $y=(y_{1},y_{2},\dots,y_{n})$, is rescaled by multiplying by a value $u=\frac{\max(\lambda)-\min(\lambda)}{10\times\max(y)}$. This $u$ corresponds to the $\alpha$ value for the alpha shape. Since the construction of an alpha shape depends on coverage areas of small circles, the relative range of $\lambda$ and $y$ affects results: if the range of $y$ is too large or too small compared to the range of $\lambda$, larger alpha values should be used for the alpha shape. For a common blaze function, $u=\frac{\max(\lambda)-\min(\lambda)}{10\times\max(y)}$ works well by scaling the range of $y$ and $\lambda$ to be $1:10$, in coordination with a recommended value for $\alpha$ later. In step 2, the alpha shape, $AS_{\alpha}$, is constructed with radius $\alpha$, which is an infinite point set containing all the points within the boundary of the alpha shape (see Figure2a). In $AS_{\alpha}$, for each $\lambda_{i}$, there are infinite $y^{*}_{i}$ such that $(\lambda_{i},y^{*}_{i})\in AS_{\alpha}$. The upper boundary of $AS_{\alpha}$ is defined as $\widetilde{AS}_{\alpha}$, which is a finite point set only including the largest $y^{*}_{i}$ such that $(\lambda_{i},y^{*}_{i})\in AS_{\alpha}$ for each $\lambda_{i}$, as displayed in Figure2a. In step 3, we fit a local polynomial regression using all the points in $\widetilde{AS}_{\alpha}$, denoted as $\hat{B}_{1}$, such that $\hat{B}_{1}$ is a smoothed version of $\widetilde{AS}_{\alpha}$ and is the initial model of the blaze function (see Figure2b). $\hat{B}_{1}$ is not an accurate estimate of the blaze function, but rather an approximation of its shape. Next, $y$ is divided by $\hat{B}_{1}$ to get the first estimate of the flattened spectrum, denoted as $\hat{y}^{(1)}$, displayed in Figure2c.

In step 4, we denote the intersection of $\widetilde{AS}_{\alpha}$ and the spectrum $\{(\lambda_{i},y_{i})\}_{i=1}^{n}$ as $W_{\alpha}$, shown in Figure2b and2c. Notice that $W_{\alpha}$ contains points mostly near the continuum. The blue line in Figure2b, which connects the points in $\widetilde{AS}_{\alpha}$, can be thought of as a collection of segments, and points in $W_{\alpha}$ are vertices of those segments. Each point in $W_{\alpha}$ is a vertex of a window where splits are defined - the spectrum is cut into small windows by the points in $W_{\alpha}$ in order to get the local quantiles. The purpose of the next steps is to select a subset of points that do not fall into an absorption feature. This is accomplished by selecting $\hat{y}_{i}^{(1)}$ in the upper quantiles of these windows. In particular, in the j-th window, we select the points whose $\hat{y}^{(1)}$ values are larger than the $q$ quantile of the $\hat{y}^{(1)}$ in the window, and the selected points make up the set $S_{j,\alpha,q}$. For $j=1,2,\dots,|W_{\alpha}|-1$, where $|W_{\alpha}|$ is the number of points in $W_{\alpha}$, the combined set of $S_{j,\alpha,q}$ for different j’s is defined as $S_{\alpha,q}$, displayed in Figure2c. $S_{\alpha,q}$ contains points that are locally in the upper $1-q$ quantile, which will be used to estimate the blaze function since these points generally do not fall into an absorption line. In step 5 we run a local polynomial regression on $S_{\alpha,q}$ and fit it to the whole spectrum. This regression is our final estimate of the blaze function, denoted as $\hat{B}_{2}$. The red line in Figure2d shows the final estimate of the blaze function. Then, in step 6, $y$ is divided by $\hat{B}_{2}$ to get the blaze-removed spectrum, denoted as $\hat{y}^{(2)}$ and shown in Figure2e.

2.2 ALSFS Method

The ALSFS algorithm is a method for removing the blaze function when a lab source spectrum, such as an LED or quartz lamp spectrum, is available as a reference. Generally, when a reliable reference spectrum is available, it can be used as the preliminary estimate of the blaze function shape. The use of a reference spectrum can be particularly advantageous in situations when the science spectrum contains wide absorption lines, which often make estimation of blaze function shape especially challenging.

The AFS algorithm in section 2.1 can be adapted to take advantage of this additional information in the reference spectrum. To a first approximation, a lab source continuum spectrum should trace the instrumental blaze function for each order. However, the calibration source will typically have some effective blackbody temperature – that is, its intrinsic intensity will peak at a specific wavelength. Over the limited wavelength range covered by a single spectral order, this effective blackbody function is approximately linear. Therefore, an observation of a lab source spectrum can be modeled as a linear transformation of the instrumental blaze function, and therefore the difference between the blaze function and the corresponding lab source spectrum is only a location-scale transformation. Thus, the blaze estimation problem can be translated into an optimization problem to find the best intercept, scale, and slope of the lab source.

2.2.1 ALSFS Algorithm

The ALSFS algorithm is initialized in the same manner as the AFS algorithm in steps 1 and 2, but it differs slightly in steps 3 and 4 and greatly in step 5. Let the reference lab source spectrum be $\{(\lambda_{i},l_{i})\}_{i=1}^{n}$. The ALSFS algorithm is summarized in Algorithm 2.

Steps 1 and 2 are the same as the AFS algorithm. Since prior knowledge of the shape of blaze function is available, fewer points are needed for the second local polynomial fitting. Step 3 is similar to step 3 of the AFS algorithm, but we make set $S_{\alpha,q}$ more accurate by also calculating the $2q-1$ quantile of $\hat{y}^{(1)}$, denoted as $Q_{2q-1}$. $Q_{2q-1}=Q_{1-2(1-q)}$, which means the upper $2(1-q)$ quantile of $\hat{y}^{(1)}$, displayed in Figure3a. Compared to the upper $1-q$ quantile of $\hat{y}^{(1)}$ within each window, a smaller quantile is used for the global quantile; otherwise, too many points will be excluded if the same quantile is used. Step 4 also departs from the AFS algorithm: for the j-th window, we select the points $\lambda_{i}$ where both $\hat{y}^{(1)}_{i}\geq$ the $q$ quantile of the window and $\hat{y}^{(1)}_{i}\geq Q_{2q-1}$ into set $S_{\alpha,q}$, shown in Figure3a. In step 5, we use the lab source spectrum’s intensity curve as a reference model, displayed in Figure3b, to find the best linear coefficients for blaze function estimation using the points selected in the previous step. We apply a linear transformation on the lab source spectrum: $\hat{l}_{i}(a,b,c)=a+bl_{i}+c\lambda_{i}$, where $a$, $b$, and $c$ are intercept, scale, and slope parameters, respectively. Since our ultimate goal is to have a flat spectrum without the blaze function, we use an objective function \sum\limits_{\begin{subarray}{c}i\end{subarray}\in S_{\alpha,q}}\Big{(}\frac{l_{i}}{\hat{l}_{i}(a,b,c)}-1\Big{)}^{2}, which measures the total distances from the removed spectrum to the constant 1. We minimize this objective function on the set $S_{\alpha,q}$ to get estimates for $a$, $b$, and $c$. Then the modified lab source spectrum $\hat{a}+\hat{b}l_{i}+\hat{c}\lambda_{i}$ is our final estimate for blaze function, displayed in Figure3b. Step 6 is the same as the AFS algorithm, $y$ is divided by $\hat{B}_{2}$ to get the blaze-removed spectrum.

2.2.2 Lab Source Smoothing by AFS Algorithm

The process of flat-fielding spectra is necessary for removing pixel-to-pixel quantum efficiency (QE) variations in charge coupled devices (CCDs) that are used as detectors in astronomical spectrographs. In the case of fiber-fed echelle spectrographs, flat fielding can be carried out by extracting a featureless calibration spectrum and dividing the extracted science spectrum order-by-order. However, the flat-field source is generally not perfectly uniform in intensity over a large wavelength range, and therefore it will typically have an effective black-body temperature that does not match the stellar effective temperature. As a result, this division leaves behind residual trends. By using the AFS algorithm a model can be fitted to each order of the flat field calibration spectrum. Division of the flat field echelle orders by this fitted model will then produce a normalized spectrum that can be used to divide out the QE variation. Figure 4 shows an example of how this process was used to create a normalized flat in red orders that exhibit fringing from interference of red wavelengths in thinned silicon detectors. This fringing can be removed by dividing stellar spectra with this normalized flat.

In some orders a cosmic ray or a pixel with very low QE will create an upward or downward spike that is confined to one or a few pixels. In this case, the AFS algorithm is slightly modified to iteratively reject these pixels using outlier rejection. Let the original lab source spectrum be $\{(\lambda_{i},L_{i})\}_{i=1}^{n}$, and $\Delta L=\{|L_{i}-L_{i-1}|$, $i=2,\dots,n\}$. Let $Q_{q_{s}}$ be the $q_{s}$ quantile of $\Delta L$, where $q_{s}$ could be a number between $0.95$ and $0.99$. In the beginning, the $0.99$ quantile of $\Delta L$ is denoted as $Q_{0.99}^{(0)}$. In the $j$-th iteration, remove pixels where $|L_{i}-L_{i-1}|>Q_{0.99}^{(j-1)}$ and calculate the $0.99$ quantile of the new $\Delta L$, denoted as $Q_{0.99}^{(j)}$. Continue the iteration until $Q_{0.99}^{(j)}<Q_{q_{s}}$. The remaining pixels are used for smoothing by the AFS algorithm, displayed in Figure 5.

2.3 Parameter Selection

In the proposed algorithms, there are several parameters that need to be selected by users. In the AFS algorithm, there are three parameters: $\alpha$ for the alpha shape, quantile $q$ for the point selection, and $m_{0}$ in two local polynomial regressions. The $\alpha$ determines how many windows a spectrum is cut into, because the number of windows is determined by smoothness of the alpha shape, which is controlled by $\alpha$. Using these windows, $q$ determines how many points are selected in each window. After selecting the points, $m_{0}$ determines the smoothness of local polynomial fitting on these selected points. Practically, the three parameters are robust within appropriate ranges. For the example spectrum in Figure6a, we recommend a set of standard parameters $(\alpha=\frac{1}{6}\times$wavelength range, $q=0.95,m_{0}=0.25)$ as a default. In Figures 6-8, we show the blaze estimates using extreme parameter choices (top panels) and its resulting normalized spectra compared with the standard parameter choice (bottom panels).

First, the $\alpha$ can be selected based on the shape of a blaze function. For example, in Figure2a, we find that the blaze function increases first and then decreases. Also, it is convex first, then becomes concave, and turns to be convex in the end. We want to select an $\alpha$ that can capture the convex parts (not too large), but will not go too deep into absorption lines (not too small). The choice of $\alpha$ is mainly determined by the shape, e.g., curvature and concavity, of a blaze function. However, when there is a wide absorption line, a larger $\alpha$ may be needed than the shape would suggest. For echelle spectra orders, since each convex portion generally takes up about $\frac{1}{6}$ of an order, we recommend selecting an $\alpha$ that is $\frac{1}{6}$ of the wavelength range of the order, but have found empirically that $\alpha$ is rather robust from $\frac{1}{12}$ of the wavelength range to $\frac{1}{3}$ of the wavelength range. Figure6a shows an example of a very large $\alpha$ equal to the order’s entire wavelength range. This $\alpha$ value is too large for the $\alpha$-balls to capture the shape of the spectrum near the boundaries; in the blaze-removed spectrum, the left part drops downward like a wide absorption feature because points in that portion of the spectrum were not selected in $S_{\alpha,q}$. Figure6b shows an example of a small $\alpha$ of $\frac{1}{50}$ of the wavelength range. With such a small $\alpha$, the $\alpha$-balls fit into absorption lines so that points inside an absorption are selected into set $S_{\alpha,q}$. In the removed spectrum, there are regions that are above the reference line $y=1$, compared with the cyan spectrum where $\alpha=\frac{1}{6}\times$wavelength range.

The parameter $q$ depends on the S/N of a spectrum and the amount of absorption. The goal when selecting $q$ is to find points on the spectrum that do not drop into absorption lines, but instead are on the continuum. After flattening the spectrum using a preliminary estimate of its shape, we select points in the local upper $1-q$ quantile for the set $S_{\alpha,q}$ to be used in the final estimation. If we happened to know the true blaze, then after dividing the spectrum by the blaze we would expect to see points randomly scattered around $1$. In the proposed algorithms, these points are approximated by set $S_{\alpha,q}$. If the S/N is high or there is a large amount of absorption, a larger $q$ is needed to select points in $S_{\alpha,q}$ so that these points do not fall in absorption lines. Conversely, if S/N is low or there is minimal absorption, a smaller $q$ is needed to get enough points into set $S_{\alpha,q}$222Alternatively, an adaptive $q$ can be used for each small window to capture the noise more accurately. An adaptive $q$ can incorporate the overall S/N and the average intensity in each window to characterize the noise variance better.. For an order with a similar amount of absorption as the one displayed in Figure6a, a $q$ from $0.95$ to $0.99$ works for S/N 300, a $q$ from $0.85$ to $0.95$ works for S/N 150, and a $q$ from $0.5$ to $0.85$ works for S/N lower than 150. Figure7a shows the effects of a small $q$ such as $0.7$: too many points are selected into $S_{\alpha,q}$ so that the blaze function estimate is dragged downward by the points in absorption lines, and in the removed spectrum, the left and right boundary regions are above the reference line $y=1$. In contrast, Figure7b shows a large $q$ of $0.999$. In this example $S_{\alpha,q}$ contains too few points, which makes the estimate sit almost completely above the spectrum. In the removed spectrum, almost all the pixels are under the reference line.

The smoothing parameter $m_{0}$ depends on the distribution of points in $S_{\alpha,q}$ along the spectrum, which is determined by the amount of absorption. If there are many absorption lines or any absorption lines that are wide, the set $S_{\alpha,q}$ has large gaps between pixels and thus a large $m_{0}$ is needed to get a good estimate. If there are few absorption lines or absorption lines that are narrow, a small $m_{0}$ is needed so that the estimation better adapts to local regions. For an order with a similar amount of absorption as the one displayed in Figure6a, an $m_{0}$ value from $0.15$ to $0.3$ has worked well empirically. In Figure8a, $m_{0}$ is set to be $0.5$ (too large) and the estimate is off on the left part of the spectrum, where the blaze-removed spectrum rises to above $1.5$. Figure8b shows the results of a too small $m_{0}$ value of $0.1$: the blaze estimate has some small bumps, but the blaze-removed spectrum looks reasonable to the eye. However, since the true blaze function does not have small bumps, this blaze estimate is not as good a fit as it might appear at first glance.

In general, the blaze function estimate is more sensitive to small changes in $q$ than to $\alpha$ and $m_{0}$. For an echelle spectrum order with a similar amount of absorption as the one displayed in Figure6a, we can start with an $\alpha$ equal to $\frac{1}{6}$ of the wavelength range of the order, an $m_{0}$ equal to $0.25$, and tune the parameter $q$ within the range according to its S/N and amount of absorption. A more detailed set of recommendations for parameters in different situations is provided in the Appendix.

The ALSFS algorithm has the same parameters: $\alpha$, $q$ and $m_{0}$. Because the final estimate is the reference lab source spectrum with linear modification, $m_{0}$ has less influence on the results than for the AFS algorithm. In the lab source spectrum smoothing process, $\alpha$, $q$, and $m_{0}$ operate the same as in the AFS algorithm. The quantile parameter $q_{s}$ in $Q_{q_{s}}$ depends on the particular appearances of spikes. If spikes are long, use a smaller $q_{s}$ such as $0.95$; if spikes are very small, use a larger $q_{s}$ such as $0.98$ or $0.99$.

2.4 Complications and Corrections

We have found that the proposed algorithms work well in most cases. However, there are several special cases that can result in poorer estimates of the blaze function, for which we have developed corrections to mitigate these issues. Since the AFS algorithm relies only on the science spectrum itself, it is more susceptible to complications than the ALSFS algorithm.

2.4.1 Boundary Correction

An order is normalized by dividing by its estimated blaze function. Since the blaze function approaches zero near the edges of the order, small errors in the blaze shape will be magnified in the divided spectrum. The ALSFS algorithm is less susceptible to this problem because of the strong constraints provided by the lab source, but other blaze estimation methods, including the AFS algorithm, can be strongly affected.

The AFS algorithm also has difficulty in this region in cases where the edge of an order splits an absorption line. Fortunately, neighboring orders often have some region of overlap that can be used to correct the boundaries. A weighted average of the blaze-removed spectrum of the two orders can be used as an estimate of the blaze function in the overlapping region. For example, Figure 9 shows two neighboring orders that share an overlapping region. Figure9a shows the right boundary of the left order, which looks good as a blaze-removed spectrum using AFS algorithm. Figure9b shows the left boundary of the right order, which spuriously rises above $1$. We correct the overlapping region using the following:

[TABLE]

where $y_{1}$ and $y_{2}$ are the intensities for overlapping regions from the left and right orders, respectively, $y_{corrected}$ is the array of intensities for the corrected overlapping region, $m$ is the number of pixels in the overlapping region, and $w_{l}=1-\frac{l}{m}$, for $l=1,\dots,m$. The result of the correction is shown in Figure9c. The boundary-corrected spectrum is much better than the original estimate, and we can achieve further improvement by changing the definition of $w_{l}$. For example, if it is known that one of the two orders has a better estimate on its boundary, we can assign more weight toward the better order. This might be the case for a pair of orders with a broad spectral feature that is cut-off on only one of the orders.

2.4.2 Wide Absorption Lines

Sometimes an order has wide absorption lines that influence the performance of the blaze function estimation. Wide absorption lines can significantly influence the AFS algorithm performance, but only slightly influence the ALSFS algorithm. Figure10b shows the order containing the two deep Na D lines, which are each so broad that the spectrum does not fully return to the nominal continuum level between them. When attempting to fit the blaze function, the alpha shape will dip into wide features like these and pull the final blaze function estimate downward. In Figure10a, the spectrum (same as in Figure10b) before the blaze function removal is displayed. Despite failing to return to the proper continuum level, this attempt at normalization looks reasonable to the eye. Without prior information, the proposed data-driven AFS algorithm cannot consistently determine the continuum level over regions of greatly extended absorption like this one. If we know there is a wide absorption feature before applying the algorithm, the regions can be masked in step 4 of the AFS algorithm, excluding those points from the estimate.

2.4.3 Continuous Opacity

If an absorption region is wider than half of the wavelength range of the target order, we refer to it as continuous opacity. This problem impacts both the AFS algorithm and the ALSFS algorithm. Since the region is so wide, masking is not an option for the AFS algorithm. The only way to deal with it is to use prior information about continuous opacity: location and intensity. Then the spectrum can be adjusted by accounting for the opacity to get a good estimation. The ALSFS algorithm can address this by connecting and adjusting neighboring orders. Since there are overlaps between neighboring orders, one can search neighboring orders until finding an order that does not contain continuous opacity as a reference. Then the orders with continuous opacity can be adjusted to the level of the reference order. An example is shown in section 3.4 to illustrate the continuous opacity correction.

3 Simulations

In our simulation study, we use an integrated-disk solar flux atlas spectrum (Wallace et al., 2011) produced by the National Solar Observatory (NSO), obtained with the McMath-Pierce Solar Telescope’s Fourier transform spectrometer. Since the spectra in the atlas were obtained with a Fourier transform spectrometer rather than an echelle spectrograph, they have no intrinsic blaze function. The atlas has been approximately continuum normalized. The spectral resolution of the atlas ranges from 350000 to 700000, and the spectra are essentially noiseless.

To mimic the data characteristic of EXPRES, we use the same wavelength ranges as EXPRES to divide the NSO spectrum into artificial orders. We impose a shape based on a blaze function estimated from a B-star spectrum onto each order and use our algorithms to remove the blaze function. Then simulated photon noise (Gaussian white noise, which well approximates Poisson noise for high S/N) is added corresponding to a S/N of 300. The noisy blaze-imposed spectrum is divided by the true blaze function of produce a benchmark flattened spectrum. In this simulation study, we use two orders: one from the bluer end of the spectrum and one from the redder end. The AFS algorithm and the ALSFS algorithm are tested on the two orders, respectively. Additionally, we compare our algorithms with the commonly used iterative method. The iterative method is introduced next.

3.1 Iterative Method

The iterative method is commonly used to remove blaze function from a spectrum. A polynomial model is fit to a spectrum order and the fit is considered as the starting estimation for the blaze function. Next, the method iterates. In each iteration, there is a threshold curve obtained by:

[TABLE]

where $fit_{polynomial}^{(t)}(\lambda_{i})$ is the 7th-order polynomial regression estimate at $\lambda_{i}$ and $t$ is the iteration time. As t increases, the threshold curve increases. The method builds a subset $M_{t}$ containing wavelength $\lambda_{i}$’s with intensity larger than $threshold(\lambda_{i},t)$. Then in the next iteration, it only uses spectrum points whose wavelength values are in $M_{t}$ to fit a polynomial model and the fit is a new estimation for continuum. To stop the iteration, a stopping time, denoted as $T$, can be set at a particular iteration. Another option is to set a standard deviation value $sd$, such that the algorithm stops when the standard deviation of the blaze-removed spectrum is smaller than $sd$. The $sd$ cannot be too small (much smaller than the standard deviation of the true spectrum without the blaze function), otherwise the iteration will not stop. In this work, we stop the iteration if either it arrives at the $T$-th iteration or the standard deviation of the blaze-removed spectrum is smaller than $sd$. Empirically, an $sd$ value around $0.05$ works well, but the performance of the iterative method is influenced by the choice of $T$. A higher S/N requires a larger $T$ and, in general, we have found a $T$ value equal to S/N seems to well.

3.2 Simulated Spectra

The first order (blue) has wavelengths ranging from 4478 to 4528 Å, called “order B”, and the second order (red) has wavelengths ranging from 6154 to 6221 Å, called “order R”. For this simulation, we require a realistic blaze function to add to the NSO spectrum. To do this, we estimate a blaze function of a B-star, HR 5501, observed with EXPRES (Jurgenson et al., 2016). We first use ALSFS algorithm, with the corresponding LED spectrum (after using AFS algorithm on it) as the lab source, on the B-star spectrum to get the estimate of the blaze function. Then we apply this estimate as the blaze function to the two orders. The raw spectrum of the B-star spectrum is then used as the lab source reference for the ALSFS algorithm. The AFS algorithm, ALSFS algorithm, and iterative method are applied to estimate blaze functions of the two orders and residuals are calculated. The residuals of order B and R with S/N$=300$ are displayed in Figure11a and 11b, respectively. Overall, ALSFS has the smallest residuals that are consistent across the whole order. AFS has larger residuals than ASLFS, which increase in the boundary regions, but the iterative method has larger residuals than AFS in both boundary regions and middle regions.

3.3 Results

Figure 11 displays a single realization of the noisy spectrum. We repeat the procedure by adding different realizations of noise 1000 times. For each of the 1000 realizations, the AFS algorithm, ALSFS algorithm, and the iterative method are applied to the two orders. This is carried out for an S/N 300, 150, and 50. The results are displayed in Figure 12 and Table 1. The root mean squared error (RMSE) is calculated as $\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}r_{i}^{2}}$, where $r_{i}$ is the residual at pixel $i$. The ALSFS algorithm has the smallest median RMSE for three scenarios: S/N 300-order B, S/N 300-order R and S/N 150-order B. For the other three scenarios, the AFS algorithm has the smallest median RMSE. Except for S/N 50-order R, the iterative method has the largest median RMSE.

3.4 Correction for Continuous Opacity

An example of continuous opacity is illustrated using the same simulation setup: we apply a blaze function obtained from a B-star spectrum to the NSO spectrum, and then we attempt to recover the underlying spectrum. Here we examine five consecutive artificial orders. A region of continuous opacity spans about 80 Å, which is captured in parts of the third, fourth and fifth orders, displayed in Figure13a. After applying the ALSFS algorithm, the resulting blaze-removed spectra are displayed in Figure13b. Although the blaze-removed spectra are flat, the ALSFS algorithm does not capture the continuous opacity correctly on its own. We know that the first two orders are not affected by the continuous opacity, and so they are used as the reference to correct the other orders: each order is linearly adjusted in intercept and slope to be align with its left neighboring order using the points within the overlapping region whose normalized intensities are in the top $\tilde{q}$ quantile. Ideally the combined segment should recover the continuous opacity in the normalized intensities after the blaze function removal. However, error in the slope estimation of the first order is amplified: the combined spectrum in Figure13c is not perfectly horizontal and it goes upward from left to right. We fit another linear regression using the points in the combined spectrum whose normalized intensities are in the top $\tilde{q}$ quantile to remove the extra slope. The resulting spectrum recovers the continuous opacity well, as displayed in Figure13d.

3.5 Cosmic Rays

The proposed methods can also be used for spectra with cosmic rays. In particular, the modified AFS algorithm for lab source smoothing, described in section 2.2.2, can be used directly to deal with the presence of cosmic rays. To demonstrate this, an order with simulated cosmic rays is displayed in Figure14a, which contains two upward spikes. The blaze function estimate from the AFS algorithm is shown in Figure14a, and the blaze-removed spectrum is displayed in Figure14a in comparison to the true spectrum. The ALSFS algorithm can be modified similarly to deal with cosmic rays.

4 Applications and Discussions

The AFS algorithm can be useful for studying telluric and micro-telluric absorption lines. Telluric lines, originating in the Earth’s atmosphere, create time-varying and humidity-dependent perturbations to the shapes of stellar lines (Leet et al. 2019, in prep). It is a particularly acute problem in the field of high-precision exoplanet radial velocity detection where uncorrected telluric lines contribute a radial-velocity error of $\sim$ 0.2 to 1 meters per second in optical wavelengths (Cunha et al., 2014) and as much as a few meters per second in near infrared wavelengths (Bean et al., 2010).

Telluric lines can be measured using spectroscopic observations of B-stars, which are bright, rapidly rotating young stars whose spectra are devoid of all but the strongest absorption lines because of extreme rotational broadening. A B-star acts as a background against which narrow telluric lines can be observed as a calibration tool for radial velocity measurements of other stars. Fitting and removing the blaze function and continuum of B-stars allows the depth of the telluric lines to be measured, which folds into current and proposed methods to mitigate their effects (e.g., Leet et al. 2019, in prep).

To illustrate the applicability of the proposed algorithms, the AFS and, for comparison, the iterative method are applied to a B-star spectrum, HR 8634, which was observed with EXPRES (Jurgenson et al., 2016) on July 7, 2018. A blue order (order B: 4376 to 4431 Å) and a red order (order R: 6085 to 6160 Å) of this spectrum were selected to show the effect of the blaze function removal. The flattened spectra are displayed in Figure 15, with a reference line at normalized intensity $=1$. Though the ground truth is unknown, it appears that the proposed AFS works well on flattening the spectra, while the spectra flattened by the iterative method have a zigzag pattern.

The ALSFS method can be useful for estimating the blaze function of the more complex spectra of late-type stars by incorporating information from a lab source spectrum to obtain an initial guess. Late-type stars are the primary targets of Extreme Precision Radial Velocity (EPRV) planet searches, which aim to suppress radial velocity measurement errors below $\sim$ 1 meter per second. Among the greatest challenges hindering EPRV is the problem of stellar activity: magnetically-driven motions within the stellar atmosphere lead to time-varying features, such as spots and faculae, that create line-profile distortions that skew the measured centroids of the lines leading to imprecise RV measurements (summarized in Fischer et al., 2016, Section 4.2).

Work to address the problem of stellar activity is ongoing, but one encouraging approach used by several teams is to investigate the sensitivity of individual spectral lines to activity, in order to obtain activity-free RVs (Davis et al., 2017; Wise et al., 2018; Dumusque, 2018). These methods all utilize some sort of continuum-fitting method, because the depths of the lines must be known with precision in order to search for correlations with activity over time. By providing flatter, more uniform blaze function estimates, the ALSFS algorithm will permit more precise measurements of the individual line depths and line profile shapes that are correlated with stellar activity.

The ALSFS and iterative method are applied to a blue order (order B: 4473 to 4529 Å) and a red order (order R: 6147 to 6223 Å) of the star 51 Pegasi, observed with EXPRES on July 8, 2018. Amplified figures of the flattened orders are displayed in Figure 16. For order B, the spectrum from ALSFS is mostly flat, while the spectrum from the iterative method has much higher intensities on boundary regions. For order R, both methods works well in flattening, but the iterative method is inaccurate in scale so that the peaks are above the reference line at Normalized Intensity $=1$.

We did not test the methods on cooler stars, where the continuum is poorly defined, such as M dwarfs. Since the sun is a relatively metal rich star, as is 51 Peg (metallicity of +0.2; Frasca et al. 2009), our methods are expected to perform at least as well on stars with lower metallicities and higher temperatures.

5 Conclusions

In this work, we presented two data-driven algorithms, AFS and ALSFS, for removing the blaze function from spectra obtained from echelle spectrographs. The key aspects of the algorithms are the use of alpha shapes to provide an initial guess of the blaze function’s shape, and the use of local polynomial regression to refine this guess. The two algorithms are designed for two scenarios: the AFS algorithm operates without a reference spectrum and may be applied directly to stellar spectra containing even a high number of absorption lines, while the ALSFS algorithm also incorporates additional information from a reference continuum spectrum to inform its initial guess. As an application of the AFS algorithm, a continuum lab source reference spectrum - such as an LED or quartz lamp spectrum - could be corrected and smoothed to be used in the ALSFS algorithm.

A simulation study was presented to illustrate the performance of the proposed algorithms compared to the commonly used iterative method for spectral normalization. In general, our algorithms have smaller RMSE than the iterative method. Overall, the ALSFS algorithm has the smallest median RMSE when S/N is high. Moreover, our algorithms are able to capture the edge effects better than the iterative approach. ALSFS is relatively robust to edge effects, and we have also developed a method of boundary correction for the AFS algorithm. Furthermore, detailed discussion regarding the applications of the algorithms was presented with examples of B-star and star 51 Pegasi spectra, which are observed with EXPRES.

This work proposes methodology to correct the continuum of an echelle spectrum by modeling the blaze function of individual orders. A flattened echelle spectrum obtained from the proposed methods works better than its original form in studying physical and astronomical properties of a star, e.g., blaze-removed B-star spectra for understanding telluric lines, more precise absorption line depths for studying stellar activity.

The selection of the parameter $\alpha$ is discussed in section 2.3 so in this section we focus on the selection of $q$ and $m_{0}$. The $m_{0}$ depends on the amount of absorption of an order, and $q$ depends on both the S/N and the amount of absorption. While S/N can be estimated, the amount of absorption is influenced by multiple factors such as wavelength, temperature, surface gravity, and stellar metallicity. Instead of recommending parameter values based on these factors individually, we provide several example orders to give an idea of the rough ranges of parameters to use. The examples below are all echelle spectra, and an $\alpha=\frac{1}{6}\times$wavelength range is used for all of them. The selected $q$ and $m_{0}$ values are listed under each figure. Figure 17 displays two orders of EXPRES spectra for the G8 V star 55 Cancri, which has a temperature of about 5165 K and a metallicity of +0.27 dex (Marcy et al., 2002). Figure 18 displays six simulated orders from the NSO solar spectrum, as described in section 3. The noise is added to each order with S/N set to be $300$, $150$, or $50$. We suggest these examples be used as a guide for parameter selection.