A constrained linear model for continuum normalization of stellar spectra

Inferring stellar parameters and chemical abundances by forward modeling stellar spectra usually requires a spectral synthesis code, or an emulator constructed from a curated training set. In these situations continuum normalization is often implemented as a pre-processing step that is independent of stellar parameters. This leads to results that are biased, or inconsistent across signal-to-noise ratios. A more justified approach is to forward model spectra with all nuisances simultaneously, but in practice this can be an expensive or non-convex optimization procedure. Here we describe a constrained linear model that can fit stellar absorption, telluric transmission, the joint continuum-instrument response. Stellar absorption and telluric transmission are each modeled by factorizing a grid of rectified theoretical spectra into two non-negative matrices with a chosen number of basis components. This model characterizes all possible spectra in many fewer parameters than comparable data-driven models. The non-negativity constraint ensures basis vectors are strictly additive, which limits rectified flux to less than or equal to unity, such that we can distinguish normalized spectra from the joint instrument-continuum response. The model requires no initial guess, and the linearity ensures that inference is convex, stable, and fast. This model allows us to reliably fit nuisances (e.g., tellurics, continuum), and is readily extensible to radial velocity and rotational broadening, without any prior knowledge about the fundamental stellar properties. We demonstrate our method by fitting ESO/HARPS high-resolution echelle spectra of BAFGKM-type stars. With repeat observations of $α$-Centauri A we present results that are best in class: consistent across time to 0.2% at S/N ~ 100, and to better than 0.5% at S/N ~ 30.

💡 Research Summary

The paper presents a novel framework for simultaneously normalizing stellar spectra and modeling nuisance effects such as telluric absorption and the combined continuum‑instrument response. Traditional continuum‑normalization methods treat the continuum as a preprocessing step independent of the stellar parameters, which introduces biases that depend on signal‑to‑noise ratio (S/N) and stellar type. The authors instead formulate a forward model in which the observed flux at each wavelength is the element‑wise product of three physical components: (i) a continuum‑normalized stellar line spectrum, (ii) a telluric transmission spectrum, and (iii) a smooth continuum‑instrument response. By taking the natural logarithm of the observed flux, the multiplicative model becomes an additive one, allowing the problem to be expressed as a linear combination of basis functions.

The core technical innovation is the use of non‑negative matrix factorization (NMF) to compress large libraries of theoretical, continuum‑normalized stellar spectra and telluric spectra into a small set of non‑negative basis vectors. The authors first compute a grid of high‑resolution synthetic spectra (the BOSZ grid) covering a wide range of stellar parameters, convolve them to the resolution of the HARPS spectrograph (R≈115 000), and apply a negative‑log transformation to obtain a sparse matrix. NMF then factorizes this matrix into K stellar basis vectors (F) and corresponding weights (W★) for any test spectrum, and similarly into L telluric basis vectors (G) and weights (W⊕). The non‑negativity constraint guarantees that the reconstructed rectified flux never exceeds unity, preserving the physical bound that absorption can only reduce flux.

The continuum‑instrument response is modeled with a set of Fourier basis functions (H) whose coefficients (γ) are solved simultaneously with the stellar (α) and telluric (β) weights. The full log‑space model is:

Y = log y ≈ W★ α + W⊕ β + H γ + noise,

where Y is the log‑transformed observed spectrum and the covariance matrix of Y is approximated via a second‑order Taylor expansion. Because the model is linear in the unknown coefficients and all coefficients are constrained to be non‑negative, the inference problem is convex; a global optimum can be found efficiently using non‑negative least squares (NNLS) or similar algorithms. No initial guess is required, and the computational cost scales linearly with the number of pixels.

The authors demonstrate the method on a large set of HARPS echelle spectra spanning B‑ through M‑type stars. Training the NMF model on hundreds of thousands of synthetic spectra takes about 30 minutes on a single CPU core; fitting a single HARPS spectrum (≈3 × 10⁶ pixels) takes only a few seconds. Validation on repeated observations of α Centauri A shows remarkable stability: at S/N≈100 the continuum‑normalized spectra agree to better than 0.2 % across epochs, and at S/N≈30 the agreement remains better than 0.5 %. These figures outperform traditional sigma‑clipping or hand‑crafted continuum methods, which typically show larger scatter and systematic trends with S/N.

Additional experiments explore the method’s robustness to rotational broadening by convolving the learned stellar basis vectors with a rotational kernel before fitting; the approach accurately recovers broadened line profiles without degrading continuum consistency. The paper also discusses limitations: the synthetic grid must adequately span the stellar parameter space of interest; the Fourier basis may struggle to capture sharp instrumental features; and residual systematics can appear at very high S/N where subtle mismatches between the synthetic and real spectra become significant. Potential extensions include incorporating more flexible continuum models (e.g., splines or Gaussian processes), jointly fitting radial velocities and rotational velocities, and applying the framework to other spectrographs (e.g., ESPRESSO, NIRSPEC) by swapping the instrument‑specific convolution kernel.

In summary, the authors introduce a constrained linear model that leverages NMF to compress theoretical spectra, transforms the multiplicative physical model into an additive linear one via logarithms, and solves for stellar, telluric, and continuum components in a convex, fast, and initialization‑free manner. The method delivers high‑precision, S/N‑independent continuum normalization across a broad range of stellar types, paving the way for more reliable automated pipelines in large spectroscopic surveys and high‑precision stellar abundance studies.