© The Authors 2024

1 Introduction

Cosmology measures and interprets the evolution of the whole Universe. To probe its dynamics and understand the nature of dark energy, observers need to compute distances at different epochs from the light they receive in telescopes. The evolution of cosmological distances with time indicates how dark energy, dark matter, and matter interact and how they can be modelled.

Optical surveys use magnitude and colour comparisons to build a relative distance scale. For instance, type Ia supernovae (SNe Ia) revealed the presence of a dark energy component because they appeared fainter in the early Universe than was thought (Riess et al. 1998; Perlmutter et al. 1999; Betoule et al. 2014; Scolnic et al. 2018). More precisely, because SNIa colours shift with the expansion of the Universe, high-redshift supernovae were fainter in red bands than what can be inferred from low-redshift supernovae observed in blue bands. This case underlines that colours need to be accurately calibrated in an optical survey to display the dynamics of the Universe (see e.g. Betoule et al. 2013). Every chromatic effect that alters the astral light distorts our dynamic perception of the expansion of the Universe, such as the galactic dust, the instrumental response, or the local atmospheric conditions.

In this paper, we present a forward-modelling method to analyse and extract data gathered with a dispersing element (grating or hologram) in the filter wheel of a telescope. We label our approach forward modelling because we implement a numerical simulation of the data-taking procedure that includes as much a priori knowledge as available, and then estimate model parameters from likelihood maximisation. This method is fundamentally different from the traditional flux-weighted sum orthogonal to the dispersion axis (Horne 1986; Robertson 1986) or the algebraic method that uses multiple images (Ryan et al. 2018) because it relies on physical modelling to directly describe the footprint of the spectrum on the imaging sensor. Deconvolution techniques and point spread function (PSF) modelling have been explored for optical fibre spectrographs (Bolton & Schlegel 2010; Li et al. 2019), but in our forward model, we proceed to build a physical model for the extraction of the spectra, and in particular, the atmospheric transmission from spectra. Our approach was inspired by the forward modelling developed in Outini & Copin (2020), and we applied it for punctual sources with the ultimate goal of measuring atmospheric transmission.

The scientific context of our work is the study of the atmospheric transmission variability via the repeated observation of stable (standard) stars. Burke et al. (2010) opened the path to controlling the optical survey photometry with dedicated measurements of atmospheric components by observing standard stars. Then, Burke et al. (2017) reached a 5 per mill (i.e. per thousand) relative photometric calibration between filters that covered the full optical and the near-infrared (near-IR) range and accounted for linear temporal variations of the atmosphere transmission over the nights. Because of the scale of the new SNe Ia surveys, the number of observed objects now reaches a point where even such an exquisite calibration becomes the dominant source of systematics. One of the challenges of modern SNe Ia cosmology thus is to be able to accurately measure and estimate the chromatic variations of the atmospheric transmission at the per-thousand level that allow us to probe for systematic variations, either nightly, seasonal, or directional.

While spectrophotometry at the required precision has been hinted at in The Nearby Supernova Factory (2013); Rubin et al. (2022), it relies on the dedicated use of a specifically designed integral field spectrograph. Our current approach instead focuses on exploring the spectrophotometric possibilities offered by a much simpler design: We consider a slitless spectrograph, where a disperser (either a grating or a hologram) is inserted in the converging beam of a telescope in a regular filter wheel. This implementation is used in the Rubin Auxiliary Telescope (Aux-Tel) (Ingraham et al. 2020), as well as on the Star Direct Illumination Calibration Experiment (StarDICE) telescope. StarDICE is an experiment aiming at transferring to stars the unit of optical power (watt) defined at the National Institute of Standards and Technology (NIST) with a reference cryogenic radiometer, the Primary Optical Watt Radiometer (POWR) (Houston & Rice 2006; Hazenberg 2019; Betoule et al. 2023).

This paper describes the preparatory work for these projects and presents the analysis procedure we developed and tested on a few nights of data gathered at the Cerro Tololo Inter-American Observatory (CTIO). It describes some implementation choices and demonstrates that the forward modelling approach allows us to incrementally build a detailed understanding of the data that in the end can permit the direct extraction of the parameters used to describe the atmospheric transmission variability.

The first section of this paper describes the theory of slitless spectrophotometry, the basic implementation in Spectractor, and the data and simulation sets we used to assess the quality of the algorithm. Section 3 details the different ingredients of the Spectractor software, the assumptions and the implementation choices. In particular, we detail the regularised deconvolution technique at the heart of the process to obtain a prior for the forward model and to qualify the code on simulations. The application of Spectractor to extract spectra from on-sky CTIO data is described in Sect. 4, while Sect. 5 focuses on the measurement of the atmospheric transmission. Discussions and summaries conclude the paper in Sect. 6.

2 Forward modelling of a slitless spectrogram

There are many different possible configurations for gathering spectroscopic data without the use of a slit. The slitless spectrograph configuration that we considered (a grating or hologram in a converging beam) can be implemented in different ways that could require special care in the forward-modelling analysis. For example, the field of view can be small and may contain only one star, or it can be crowded, with many source-dispersed images (so-called spectrograms hereafter) super-imposed on each other. Different detectors might span the field of view, with responses that need to be mapped and gaps that need to be accounted for. We did not try to abstractly solve all different situations. We therefore defer all the technical issues that we did not encounter to further work and concentrate on those we did and solved.

Within these restrictions, we consider in the following only the case of point sources and do not discuss, except for a passing mention in the next section, the case of extended sources such as galaxies or resolved planetary nebulae, nor the deblending of these extended sources with point sources.

2.1 Description and geometry of a slitless spectrograph

The slitless spectrograph we consider can be seen as a grating with N grooves per millimetre, placed in a telescope beam at a distance D_CCD from a charge-coupled device (CCD) sensor. In the following, the positions in the sensor plane are parametrised with the coordinates r = (x, y), and the z-axis points orthogonally toward the CCD. The disperser can be more or less complex and may be used in a convergent or in a parallel beam, and its resolution can vary, but in the end, it will spread the source light in different diffraction orders superimposed on the CCD, with a sky background that is also diffracted. Depending on the choice of N, on the size N_CCD (in pixels) of the sensor, and on D_CCD, different diffraction orders will finally be recorded by the sensor.

The special case of the zeroth order is worth mentioning. While its presence on the image is not mandatory, knowing its centroid position r₀ can significantly facilitate the setting of the zero of the wavelength calibration.

The positive and negative diffraction orders are placed on each side of the zeroth order on a line forming an angle α with the x-axis. We parametrised the position along this dispersion axis with the coordinate u and transversally with the coordinate υ. The zeroth order stands at coordinates (u₀, 0) in the (u, υ) coordinate system. If the instrument wavelength coverage spans more than one octave in wavelength, the different diffraction orders are superimposed on each other.

The source spectrogram is the total 2D CCD image formed by the cummulation of the dispersed light. All the notations are illustrated in Fig. 1.

For a periodic grating placed inside a convergent beam instead of a parallel beam, this optical system response is astigmatic, that is, the image of a point source such as a star is not point-like on the sensor. Usually, the image is elliptical, and the redder the wavelength, the wider the ellipse (see e.g. Moniez et al. 2021). However, the centroid of these ellipses is still given by the classical grating formula (see e.g. Murty 1962; Hall 1966; Schroeder & Inc 2000),

$$\sin {\theta }_p(\lambda )-\sin {\theta }_0=p{N}_{\text{eff}}\lambda ,$$(1) $$\tan \theta (\lambda )=\frac{u(\lambda )}{D_{\text{CCD}}},\tan {\theta }_0=\frac{u_0}{D_{\text{CCD}}}.$$(2)

The angles are those of the projection in the plane perpendicular to the grating lines (see Fig. 1); θ₀ is the angle of the projected telescope beam axis with respect to the normal to the grating surface, ρ is the diffraction order, θ_p(λ) is the corresponding projected diffracted angle, and N_eff is the effective spatial frequency of the grating lines at the position of the central ray¹ of the light beam (hereafter called chief ray).

The observer can choose the focus of the telescope. One common choice is to set the focus using the zeroth order, but if this is done, the spectrogram can be affected by defocusing effects at increasingly redder wavelengths (for periodic gratings). To minimise the defocusing effect and increase the spectrograph resolution, the focus for a particular wavelength can be optimised, but it is then more difficult to set the zero of the wavelength calibration with a defocused zeroth-order image. In the following, we assume that the focus has been made on the zeroth order unless otherwise specified.

2.2 Theoretical model of a spectrogram

A theoretically perfect image of a light source can be modelled by its spatio-spectral flux density C(r, λ). For a point source, we can separate the spectral and spatial distributions as

$$S_{\ast }(r,\lambda )=S_{\ast }(\lambda )\times \delta \left(r-r_0\right),$$(3)

with S_*.(λ) spectral energy density (SED) of the astrophysical object, and δ the Dirac distribution. The observed spectral and spatial distribution is then

$$C_p(r,\lambda )=T_{\text{inst},p}(\lambda )T_{\text{atm}}\left(\lambda \mid P_a\right)S_{\ast }(r,\lambda ),$$(4)

where T_atm(λ|P_a) is the atmospheric transmission, depending on a set of atmospheric parameters P_a, and T_inst,p(λ) is the instrumental transmission (including the CCD quantum efficiency) for the diffraction order p.

In the case of a monochromatic point source, we have C₀(r, λ) = Α δ(λ – λ₀) × δ (r – r₀), with A the received source flux at λ₀. The PSF describes the optical response of the telescope and of the atmospheric turbulence on a sensitive surface such as a CCD (see Fig. 2a). It depends a priori on the wavelength λ and can be modelled by a function ϕ₀(r, λ) whose spatial integral is normalised to one. Therefore, by definition of the PSF, the image of a monochromatic point source centred on r₀ on the CCD can be described as a convolution product,

$$\begin{array}{l}{I}_0(r)={{\displaystyle \int }}^{\text{}}\text{d}\lambda {\displaystyle \int {\displaystyle \int {\text{d}}^2{r}^{\prime }{C}_0\left(r^{\prime },\lambda \right){\varphi }_0\left(r-r^{\prime },\lambda \right)}}\hfill \\ =A{\varphi }_0\left(r-r_0,{\lambda }_0\right).\hfill \end{array}$$(5)

With a slitless spectrograph, the mechanism is the same, but the incoming light is dispersed in several diffraction orders p, and light from all diffraction orders of the sky background is superimposed. The position of the point-source image depends on the wavelength and on the order p. Moreover, the shape of the PSF itself can depend on the order p and on the wavelength λ. Here and everywhere else, the dispersed-imaging PSF integrates both the seeing and the instrumental PSF. We introduce the dispersion relation Δ_p(λ) as the 2D vectorial quantity that describes the position of the PSF centroids {x_c,p(λ) – x₀,y_c,p(λ) – y₀} on the CCD with respect to the zeroth-order position (x₀, y₀), for a diffraction order p and a wavelength λ. This quantity can be computed by applying the classical grating formula (1). With a point-like monochromatic source, the image recorded on the CCD is modelled as

$$\begin{array}{l}I(r)={\displaystyle \underset{p}{{{\displaystyle \sum }}^{\text{}}}}{{\displaystyle \int }}^{\text{}}\text{d}\lambda {\displaystyle \int {\displaystyle \int {\text{d}}^2{r}^{\prime }{C}_p\left(r^{\prime },\lambda \right){\varphi }_p\left(r-{\Delta }_p(\lambda )-r^{\prime },\lambda \right)}}\hfill \\ ={\displaystyle \underset{p}{{{\displaystyle \sum }}^{\text{}}}}{A}_p{\varphi }_p\left(r-{\Delta }_p\left({\lambda }_0\right),{\lambda }_0\right),\hfill \end{array}$$(6)

with Δ₀(λ₀) = r₀ and Α_p the flux density at wavelength λ₀ for the order p. On the image, we expect a series of spots, one per order p, of different intensities, with different sizes, but containing the same spectral information S_*(λ₀) about the source.

The observed sources are naturally polychromatic. For pointlike sources, the theoretical description above holds at each wavelength, and the image can be described as

$$I(r)={\displaystyle \underset{p}{{{\displaystyle \sum }}^{\text{}}}}{{\displaystyle \int }}^{\text{}}\text{d}\lambda S_p(\lambda ){\varphi }_p\left(r-{\Delta }_p(\lambda ),\lambda \right),$$(7)

with

$${\text{Δ}}_p\left(\lambda \right)=\left\{x_{c,p}(\lambda )-x_0,y_{c,p}(\lambda )-y_0\right\},$$(8) $$S_p(\lambda )=T_{\text{inst },p}(\lambda )T_{\text{atm }}\left(\lambda \mid P_a\right)S_{\ast }(\lambda ).$$(9)

Therefore, the spectrogram of a polychromatic source can be viewed as a stack of monochromatic images with different centroids, or as a dispersed image with a very chromatic PSF (see Fig. 2b). This description of the image and of the spectrogram formation is the base of our forward model for slitless spectrograph data.

To obtain the S_p(λ) spectra, a process to un-stack the monochromatic images spread across the image by the disperser is needed (Fig. 2c). This in turn requires that ingredients such as the PSF ϕ_p(r, λ) and the dispersion relation Δ_p(λ) are sufficiently well known, either a priori, through a specific data analysis, or directly fitted on data. The hardest point usually is the determination of the PSF kernel as a function of wavelength. As illustrated in Figs. 1 and 2, in the general case, the PSF is blurred and defocused. A simple Moffat profile can be used to model the PSF as long as the defocus is small; we used this approximation throughout this paper and leave the general defocus case for future work. Nevertheless, we studied and discuss two cases: The use of a standard grating, and the use of an holographic disperser that limits the defocus.

In this paper, we describe our implementation of this process in the form of the python3 Spectractor code (Neveu et al. 2021) (see Sect. 3). We also show how we tested it on simulations and data and how it could be used in order to ingest lacking modelling information from specific data analysis.

2.3 Spectractor

We call Spectractor²,³ the computer suite we wrote to analyse the future AuxTel images and the images obtained on the Cerro Tololo Inter-American Observatory (CTIO) 0.9 m telescope. It was trained on CTIO data, but with the purpose of being easily configurable for slitless spectrophotometry with other telescopes. The main steps, inputs, and outputs of the extraction part are illustrated in Fig. 3 and are described in detail in Sect. 3. To facilitate reading of the following, we summarise the main Spectractor steps below.

Zeroth-order centroid. The main inputs of Spectractor are a pre-processed image (overscan subtracted, debiased, and spatial flat fielded) obtained with a slitless spectrograph, and a configuration file setting the main geometrical and spectrographic properties of the instrument (D_CCD, N, telescope diameter, pixel size, the PSF model, etc.) At the time of writing this paper, the PSF models we implemented were either a Gaussian profile, a Moffat profile, or a Moffat minus a Gaussian profile. Additional and more detailed PSF models are planned to describe the AuxTel data more accurately as they are analysed. The centroid of the zeroth order is contained in the image. We thus implemented a searching procedure to inform us on the origin of the location of the spectrum and set the zero of the wavelength scale.

Rotation. Because of the geometry of the spectrograph and the dispersion properties of the grating, the spectrogram was cropped from the image and was de-rotated so that the dispersion axis was along the horizontal axis of the cropped image. The rotation angle is fitted later in the full model without explicitly rotating the image.

Geometric wavelength calibration. On this image, a first fit of a 1D sliced PSF model transverse to the dispersion axis is performed. The PSF shape parameters are represented by a polynomial function⁴ as a function of the distance to the zeroth order.

PSF deconvolution. The procedure is continued by a deconvolution that uses a 2D PSF model and the ID result as a prior to regularise it.

Wavelength calibration. A first wavelength calibration is performed using the detection of the principal absorption or emission lines in the extracted spectrum (astrophysical or telluric lines).

Flux calibration. The spectrum flux in ADU is converted into erg s⁻¹ cm⁻² nm⁻¹ units using the telescope collecting area, CCD gain, and exposure time.

Full forward model. Finally, given all these prior ingredients, a full forward model is initiated on the pre-processed⁵ but not rotated exposure using a 2D PSF model, and a model for the atmospheric differential refraction (ADR)⁶. The PSF shape parameters as well as the wavelength calibration are refitted in the process. The main output is a calibrated spectrum in wavelength and amplitude, but Spectractor also returns a host of useful fitted parameters such as the PSF shape chromaticity or the D_CCD distance to perform extraction quality analyses, and in the end, to improve the forward modelling or its initialisation.

All steps but the last are implemented in order to provide the required ingredients to the full forward models, and they are thus completely contingent to the availability of understanding the instrument. This is again be addressed below when we discuss how the optical transmission of the instrument can be obtained.

2.4 Data examples

While we use the natural ability of the forward-model implementation to easily provide simulations to validate the code, we also present the use of Spectractor on real data.

2.4.1 CTIO data

In order to test the approach of slitless spectroscopy, in particular, using a holographic disperser, we benefited of a run of 17 nights in May-June 2017 at the Cerro Tololo Inter-American Observatory (CTIO) 0.9 m Cassegrain telescope (f/13.7, scale at focal plane 60 μm″⁻¹). This telescope is equipped with a cooled Tek2K CCD device of 2048 × 2046 pixels, read by four amplifiers⁷. Two filter wheels are installed. The first wheel in the light path was used to insert broad-band filters, and the second wheel holds different dispersers.

While many gratings were tried, we focused on two of them here: a Thorlabs blazed grating (300 lines mm⁻¹) ref. GT50-03, and an amplitude holographic optical element (around 350 lines mm⁻¹) especially designed for this telescope. This hologram is fully described and analysed based on the CTIO data in Moniez et al. (2021). Its main advantage is that the defocusing described in Fig. 1 is very limited, which allowed us to model its chromatic PSF with simpler mathematical models.

By using these dispersers in the upstream filter wheel, we readily transformed the CTIO 0.9 telescope into a spectrophotometric instrument with a resolution of about 150-200 (Moniez et al. 2021).

Figure 4 shows an example of the data we obtained: The dispersion axes are nearly horizontal along the x-axis of the CCD, and for an optimal focusing of the amplitude hologram, the target star was placed around pixel coordinates (750, 700). The spectrum covers two amplifiers. Field stars are present, and the sky background is also dispersed (brighter in the middle).

Dome flats were taken with two different filters, blue (λ < 550nm) and red (λ > 715nm). We saw no difference in the extraction of the spectra using one or the other. In order to be more precise about the extraction of absorption bands in water in the red part of the spectrum, we chose to use only the red flat to flatten out the pixel-to-pixel fluctuations in our exposures. Combined bias was taken at the beginning of each night for the bias subtraction. We made sure that the meta data contained informations about the properties of the CCD and the on-site meteorological station.

We mainly analysed the performances of the holographic element we brought during this run. Fortunately, we had one very good night on 2017 May 30 with very stable conditions in terms of temperature and seeing, which we used to estimate the atmospheric transmission. During that night, we essentially monitored the CALSPEC⁸ (Bohlin et al. 2014, 2020) star HD111980 with an amplitude hologram. The main characteristics of these data are summarised in Table 1.

2.4.2 Simulations

To test the Spectractor pipeline, we used the full forward model for CTIO spectrograms (see Sect. 3.8) to simulate observations of CALSPEC stars (in particular, HD111980).

The simulation used in Sect. 3 shares the same known characteristics as the real-data image presented in Table 1, but with variations in the unknown parameters such as the PSF model, the amount of second-order diffraction contamination, and the atmospheric parameters.

For the simulations, a 2D Moffat circular PSF kernel φ(x, y, λ) was chosen. To model the widening of the PSF due to defocusing or chromatic seeing effects, the shape parameters γ and α evolved with wavelength as a n_PSF order polynomial function,

$$\varphi \left(x,y\mid r_c,P\right)=A{\left[1+{\left(\frac{x-x_c}{\gamma (z(x))}\right)}^2+{\left(\frac{y-y_c}{\gamma (z(x))}\right)}^2\right]}^{-\alpha (z(x))}$$(10) $$A=\frac{\alpha -1}{\pi {\gamma }^2}\text{\hspace{1em} with \hspace{1em}}\alpha >1\text{, }$$(11) $$\gamma (z)={\displaystyle \underset{i=0}{\overset{n_{\text{PSF}}}{{{\displaystyle \sum }}^{\text{}}}}}{\gamma }_i{L}_i(z),\alpha (z)={\displaystyle \underset{i=0}{\overset{n_{\text{PSF}}}{{{\displaystyle \sum }}^{\text{}}}}}{\alpha }_i{L}_i(z).$$(12)

The integral of this PSF kernel is exactly A, and its centre is at r_c = (x_c, y_c). The PSF shape parameters (γ(x), α(x)) are themselves sets of polynomial coefficients γ_i and α_i, respectively. The L_i(x) functions are the order i Legendre polynomials. We call x_min and x_max the left and right pixel positions of the spectrogram edges on the x-axis. The parameter z ∈ [−1, 1] was rescaled proportionally on the desired pixel range [x_min, x_max], set to approximately encompass the wavelength range [350 nm, 1100 nm] with the formula

$$z(x)=\frac{x-\left(x_{\max }+x_{\min }\right)/2}{\left(x_{\max }-x_{\min }\right)/2}.$$(13)

Parametrisation with Legendre polynomials has the advantage to give an equal weight to all polynomial coefficients during χ² minimisation, regardless of the degree of the polynomial functions. The chosen n_PSF, γ_i, and α_i values are quoted for each simulation. The simulation suite is fully available in the Spectractor code.

3 Forward-model extraction of a spectrum

In this section, we follow and describe in detail the steps of the Spectractor pipeline in order to obtain the first-order spectrum S₁ (λ) of a point source, calibrated in flux and wavelength. These steps cover the orange boxes of Fig. 3. The process starts from a preprocessed image containing the 2D image that is formed by a star observed through a slitless spectrograph, which we call a spectrogram.

3.1 Uncertainty evaluation

When it is not given, we must start to build the uncertainty map of the exposure when the exposure is pre-processed. The uncertainties on the pixel values are estimated using the CCD gain G_CCD(x, y) (in electrons ADU⁻¹) and its read-out noise σ_ro (in electrons). Other sources of noise, such as those coming from the flat-fielding or the dark current, are subdominant. The correlations between pixels are also negligible. The exposure unit is considered to be in ADU at that point. The uncertainty σ(x, y) on the pixel value I(x, y) is then

$$\sigma (x,y)=\frac{1}{G_{\text{CCD}}(x,y)}\sqrt{{\sigma }_{\text{ro}}^2+G_{\text{CCD}}(x,y)I(x,y)}$$(14)

because we assume that the number of photoelectrons in each pixel follows a Poisson distribution. The uncertainty map can be inverted to obtain a weight map, on which we can superimpose a mask to remove bad pixels. Assuming no correlations between pixels, we then assemble the weight matrix W as the diagonal matrix of the inverted pixel variances.

The computed noise variance uses the pixel value itself, which incorporates the noise fluctuation. Using weights incorporating the fluctuations can introduce a bias on the recovered quantities because in the case of positive fluctuation, the weight is increased, while in the case of negative fluctuation, the weight is decreased. The bias is smaller with spectrograms with a high signal-to-noise ratio (S/N). For this reason, we only considered the case with a high S/N in simulations and data, and we avoid the bias due to this uncertainty evaluation (even at the spectrum edges, where flux is low). The case with a low S/N is postponed for future developments.

3.2 Zeroth-order centroid

Because the zeroth order is included in the observed images, we used its centroid to set the zero of the wavelength scale. Therefore, an error in the determination of the zeroth-order centroid (x₀, y₀) causes a systematic shift of the wavelength calibration.

A subset of different situations encountered in the CTIO data is presented in Fig. 5. To obtain a high S/N in the spectrogram, the exposure time was set at such a value that the zeroth order is saturated, causing bleeding spikes. If the exposure has astro-metric coordinates in the World Coordinate System (WCS), then we can obtain the precise position of the star on the CCD. However, in most images we obtain from CTIO, the WCS associated with the images is incorrect because the mount calibration is not correct.

While not directly part of the extraction procedure, we present the algorithm we implemented to find the zeroth-order centroid in these difficult cases as a useful example of how to improve the starting point from a preliminary analysis of the data.

First, the image is cropped around the supposed position of the zeroth-order image, as close as needed so that the targeted star is the brightest object. Then the cropped image is projected onto the x- and (y-axis: The maximum of the two projections sets a new approximation of the zeroth-order position. From there, saturated pixels are detected, and a null weight is associated with them. A 2D second-order polynomial background with a 3σ outlier removal is fitted and subtracted from the cropped image. Finally, a 2D circular Moffat profile is fitted on the weighted pixels: Only the crown of non-saturated pixels counts and locks the fit. This last step is then repeated a second time on a new cropped image for which the width and height are divided by two, centred on the last fitted centroid. This step is illustrated in Fig. 6. We tested this process on many images, most of which were pathological, and we visually confirmed that the accuracy of this algorithm is finer than the pixel size on CTIO images.

Another issue we solved is that when a disperser is added to the telescope beam path, the WCS associated with the image can be shifted or distorted. In the case of a crowded field, for faint objects, or for a very pathological zeroth order, we developed a method for estimating the WCS using the field stars, the library astrometry.net⁹ (Lang et al. 2010), and the Gaia DR2 catalogue. The process is described in Appendix A, but it was not used by default for the CALSPEC stars we observed. However, when the centroids obtained with both methods are compared, their difference has a scatter of 0.″15 (around one-half of a pixel) on CTIO data, which confirms that both methods are accurate at the 0.″15 level at least (Fig. 7).

This accuracy is converted into a prior on the zeroth-order shift δu₀ used during the wavelength-calibration process to account for a mistake in the evaluation of the zeroth-order centroid.

3.3 Rotation

Unless special care has been taken to that end when mounting the disperser in the filter wheel, the spectrogram image can be tilted (possibly intentionally) with respect to the x-axis, with an angle, which can be poorly known depending on the mounting of the disperser into the telescope beam.

This dispersion direction can either be known a priori or fitted in the full forward-model step. We needed a supplementary step to estimate a good starting point for this angle because our case is the latter case.

In addition, we found it extremely useful that the spectrogram was roughly aligned with the x-axis of the exposure in our procedure, with the wavelength increasing with x. To this end, the exposure must be flipped and rotated accordingly before the process of extracting more information could be continue. That was useful both for diagnosing the data quality and for refining the starting point of the full forward model.

The spectrogram of sources that are sufficiently continuous in wavelength, such as the thermal emission component of stars, displays filament shapes on the 2D image that can be detected using a Hessian analysis inspired by the one developed in Planck Collaboration Int. XXXVIII. (2016). The advantage of this technique is that it comes with an analytical expression of the angle of the detected shape with respect to the horizontal or vertical axis of the CCD grid.

The Hessian matrix H(x, y) of the image is computed for each pixel value I(x, y) as

$$H(x,y)=\left(\begin{array}{ll}{H}_{xx}\hfill & H_{xy}\hfill \\ H_{xy}\hfill & H_{yy}\hfill \end{array}\right)=\left(\begin{array}{cc}\frac{{\partial }^{2}I}{\partial x^2}& \frac{{\partial }^{2}I}{\partial x\partial y}\\ \frac{{\partial }^{2}I}{\partial x\partial y}& \frac{{\partial }^{2}I}{\partial y^2}\end{array}\right).$$(15)

The two eigenvalues of the Hessian matrix Η are calculated as

$${\lambda }_{\pm }(x,y)=\frac{1}{2}\left(H_{xx}+H_{yy}\pm h\right),$$(16)

with $h=\sqrt{{\left(H_{xx}-H_{yy}\right)}^2+4H_{xy}^2}$. The eigenvalue λ_ is associated with the eigenvector directed along the spectrum dispersion axis, and λ₊ corresponds to the eigenvector with the largest change in intensity value, that is, transverse to the dispersion axis. The orientation angle of these eigenvectors with respect to the x-axis can be computed analytically. For instance, we have for λ_

$$\begin{array}{l}\alpha (x,y)=\arctan \left(\frac{H_{yy}-H_{xx}-h}{2H_{xy}}\right)\hfill \\ =\frac{1}{2}\arctan \left(\frac{2H_{xy}}{H_{xx}-H_{yy}}\right),\hfill \end{array}$$(17)

with the trigonometric formula tan 2α = 2 tan α/ (l – tan² α). After selecting the 5% pixels with the highest λ_ value above a reasonable threshold, the median α of the remaining α(x, y) values gives the orientation of the spectrum with respect to the x-axis. A linear fit can also be performed across the selected pixels, and the slope gives an angle very close to the one estimated with the median of the angle values. This process is illustrated Fig. 8.

Because of the atmospheric differential refraction (ADR), the spectrogram can be sheared transversally to the dispersion axis. The ADR shear is about 2″ across the visible spectrum, which is 5 pixels at CTIO, while the spectrogram is ≈700 pixels long, and it is thus neglected at this step of this analysis. On the other hand, it is fully accounted for in the full forward model.

3.4 Background estimation

The background of the image must be carefully subtracted to avoid bias in the estimation of the spectrum flux and its PSF. However, due to optical vignetting, it can be non-flat, and it is dispersed. It can also contain additional field stars and their corresponding spectrograms.

To estimate the spectrogram background, we first selected two lateral bands above and below the spectrogram region that had the same length and width $N_y^{(\text{bgd})}$ (see Fig. 9). First, we masked the sources detected above a 3σ threshold. Then we divided the two lateral regions in a few square boxes of size $\left(N_y^{(\text{bgd})}/2,N_y^{(\text{bgd})}/2\right)$. The box dimensions must be larger than the typical size of a field star PSF or of a spectrogram width, but small enough to account for the spatial variations in the background.

To obtain a first estimate of the background, we used a Python wrapper of the SExtractor¹⁰ algorithm for background extraction, which is roughly based on the bilinear interpolation of the median value inside the boxes, after a sigma-clipping rejection of the outlier pixels (more details in Bertin & Arnouts 1996). This process is illustrated in Fig. 9.

In a second step, we analysed the distribution of the background residuals normalised with their uncertainties in the two regions, with the sources being masked. When the histogram of the residuals normalised by their errors (the pull distribution) departs from a distribution of zero mean and standard deviation equal to 1, we refined the background estimate by dividing the box size by 2, and the process was continued iteratively until the mean was below 1 and the standard deviation was below 2, or until the box size decreases below a threshold of 5 pixels.

At the end of this process, the estimated background was interpolated between the two lateral bands to obtain the background below the spectrogram, and it was finally subtracted. We call B(r) the background map. The background RMS was also evaluated by SExtractor and was added quadratically as a background uncertainty to the error budget of the spectrogram.

At this stage of the pipeline, with α and r0 given and using a first geometrical wavelength calibration to approximately define the left and right margins of the spectrogram, we cropped the exposure to extract a background-subtracted spectrogram. This is extremely useful for diagnosis purposes and can also be directly studied with an atmospheric forward model (see Sect. 5).

3.5 First spectrum extraction

The next step of the pipeline we implemented is devoted to extracting the first-order spectrum S₁(λ) and estimating the wavelength-dependent PSF. The extraction of the spectral information S₁(λ) from the spectrogram image is a delicate process. The intuitive and traditional way to extract slitless spectra at this stage of the pipeline is to sum over the cross-dispersion direction, possibly with weights, to form what we call a cross-dispersion spectrum. Horne (1986) and Robertson (1986) presented an optimal method to achieve this. However, this method leads to distorted spectra in case of a wavelength-dependent PSF (neighbouring wavelengths contaminate each other on the sensor), which is not an issue for spectroscopy, but is problematic for spectrophotometry. In the following, we first describe a method for deconvolving the spectrum for the PSF (Sect. 3.5). The products of that process are then very useful to inform the full forward model, finalising the unbiased extraction of the spectrum (Sect. 3.8).

3.5.1 First-order model of the spectrogram

To start, we only consider the spectrogram of the first diffraction order, potentially with the superposition of a second diffraction order, as provided by the previous step of the pipeline. The cropped spectrogram has the shape (N_x, N_y) pixels.

Inspired by Eq. (7), we model the spectrogram as a discrete stack of N_x 2D PSF realisations of amplitude A_i, separated by one pixel along the x-axis,

$$I_1\left(r\mid A,r_c,P\right)={\displaystyle \underset{i=0}{\overset{N_x}{{{\displaystyle \sum }}^{\text{}}}}}{A}_i\varphi \left(r\mid r_{c,i},P_i\right),$$(18)

with r = (x, y) the vector of the pixel coordinates, A the amplitude parameter vector along the dispersion axis of the spectrogram, and φ(r|r_c,i, P_i) the 2D PSF kernel whose integral is normalised to one. This kernel depends non-linearly on the shape parameter vector P_i and on the centroid position vector r_c,i = (x_c,i·, y_c,i), where only the y_c,i coordinate is considered unknown. The x_c,i coordinate is set directly to the pixel index i. This choice of implementation can be discussed and changed, but it was found to be practical because the PSF is then well sampled by the pixel grid. In theory, we could choose another sampling for x_c,i, however, to increase the speed of the spectrum extraction or to enhance the spectral resolution.

In some way, the array of vectors r_c is a sampled precursor of the dispersion relation Δ₁(λ), and the vector A is the flux density S₁(λ) integrated within the pixels. When we index all the N_xN_y spectrogram pixels as a long vector Z = (ζ₁,…, ζN_xN_y), the Eq. (18) takes a matricial form,

$$I_1\left(Z\mid A,r_c,P\right)=M\left(Z\mid r_c,P\right)A,$$(19)

with

$$\begin{array}{l}M\left(Z\mid r_c,P\right)=\hfill \\ \left(\begin{array}{ccc}\varphi \left({\zeta }_1\mid r_{c,1},P_1\right)& \cdots & \varphi \left({\zeta }_1\mid r_{c,N_x},P_{N_x}\right)\\ ⋮& \ddots & ⋮\\ \varphi \left({\zeta }_{N_x{N}_y}\mid r_{c,1},P_1\right)& \cdots & \varphi \left({\zeta }_{N_x{N}_y}\mid r_{c,N_x},P_{N_x}\right)\end{array}\right).\hfill \end{array}$$(20)

The (N_xN_y, N_x) matrix M is called the design matrix.

This model of the spectrogram is designed to deconvolve the spectrum S₁ (λ) from the PSF. In principle, we can choose to sample it with PSF kernels separated by arbitrary distances. However, if the PSF is correctly sampled by the pixel grid, it is difficult to extract the spectrum with a spatial resolution below the typical PSF width, and therefore, below the pixel size.

At this stage, we wish to extract the spectrum from a spectrogram that is potentially contaminated by higher diffraction orders, as yielded by the pipeline steps discussed above, or that is distorted by atmospheric differential refraction. The final extraction using a full forward model that takes these physical effects entirely into account is the last step of the Spectractor pipeline. This is presented in Sect. 3.8.

3.5.2 Preparation for deconvolution

To initialise the deconvolution with parameters close to the final best model, a first PSF fit is performed transversally to the dispersion axis with 1D PSF kernels on the rotated spectrogram image,

$$I_1^{(1D)}\left(r\mid A^{(1D)},r_c,P\right)={\displaystyle \underset{i=0}{\overset{N_x}{{{\displaystyle \sum }}^{\text{}}}}}{A}_i^{(1D)}{\varphi }^{(1D)}\left(r\mid r_{c,i},P_i\right).$$(21)

This procedure is done in two steps. The first step fits the 1D parameters independently for each pixel column by applying a 5σ clipping to reject field stars and other CCD defects. In the second step, a polynomial evolution of the 1D PSF parameter vector along the dispersion axis is proposed. For instance, a polynomial evolution of the y_c,i (x_c,i) positions can model the transverse ADR, and a polynomial evolution with the width of the PSF can model a defocusing effect. The polynomial coefficients are fitted using a Gauss-Newton minimisation of a χ² (Eq. (23)) for the non-linear parameters (see Appendix B), alternating with a linear resolution for the amplitude parameters as follows.

Assuming that we gather all the spectrogram pixel values in a long vector D (as Z), we can model it as

$$D=M\left(Z\mid r_c,P\right)A+ϵ,$$(22)

with ϵ the random noise vector. The χ² function to minimise is

$${\chi }^2(A\mid P)={(D-MA)}^{T}W(D-MA),$$(23)

with W the weight matrix of dimension (N_xN_y, N_xN_y), which is the inverse of the data covariance matrix (see Sect. 3.1). In most cases, this matrix is diagonal because the pixels are all considered to be independent. The minimum of Eq. (23) is reached for the set of amplitude parameters  given by

$${\widehat{A}}^{(1D)}={\left(M^{T}WM\right)}^{-1}{M}^{T}WD.$$(24)

The covariance matrix associated with the  coefficients is C^(1D) = (M^T WM)⁻¹. At the end of this process, we obtain a first guess of the r_c and P parameters and a first guess of the amplitudes Â^(1D) with their uncertainties σ_Α^{(1 D)}, which form what we call a transverse cross-spectrum.

The result of this extraction is illustrated on simulated data in Fig. 10 left. For this simulation, we chose to ignore the second-order diffraction, and we used a wide PSF kernel with strong chromatic variations to enhance the visibility of the residuals for the 1D extraction: (γ₀,γ₁,γ₂) = (10,2,5) and (α₀,α₁,α₂) = (2, 0, 0). At first glance, the model in the top panel appears to be accurate, but the residuals show that this 1D model failed to capture the true evolution of the PSF shape with wavelength. As expected, the residuals are less dramatic with a thinner PSF, but remain measurable. This shows that a 2D PSF extraction method has to be used.

At this stage, we obtained a spectrum close to a cross-dispersion spectrum as in Horne (1986) and Robertson (1986), but informed with a fitted PSF model with a smooth polynomial wavelength evolution.

3.5.3 PSF deconvolution

Having provided ourselves with this first 1D estimate of the spectrum, we can resort to a more accurate 2D PSF modelling. When using 2D PSF kernels, the latter linear regression method enters the category of the deconvolution problems or inverse problems. Because the spectral amplitude information is mixed and diluted at a scale below the typical size of the PSF, the computation of the A vector sampled at the pixel scale inverting Eq. (19) using 2D PSF kernels may lead to results that are far from the reality, while giving an apparently good fit to data (low χ²). When this is done, a common symptom is the alternation of positive and negative values in Â, or at least with large variations, demonstrating that the problem is ill-posed. As we know by the physics that well-sampled spectra are rather continuous and differentiable functions, we enforce a regularisation method to smooth the resulting  vector.

A first fit to the data, without any prior on A, using a 1D transverse PSF model fitted independently to each column of data along the dispersion axis thus yields a vector Â^{(1 D)} that contains most of the spectral flux, especially in the smooth parts of the spectrum, but lacks precision in the rapidly varying parts. The most visible effect of the PSF is to smooth the absorption lines, and more generally, to deform the spectral information where the spectral energy density evolves rapidly (e.g. in the blue part of the spectrum), while conserving the total flux. It nonetheless provides useful information that we use as a prior A₀ on A when performing a fit using a 2D PSF kernel with a Tikhonov regularisation.

The Tikhonov regularisation method proposes to add a regularisation quantity to the χ² and minimise a new cost function,

$$\begin{array}{l}ℰ\left(A\mid r_c,P\right)={(D-MA)}^{T}W(D-MA)\hfill \\ +r{\left(A-A_0\right)}^{T}Q\left(A-A_0\right)\hfill \\ ={\chi }^2\left(A\mid r_c,P\right)+r{\chi }_{\text{pen }}^2\left(A\mid A_0\right),\hfill \end{array}$$(25) $$A_0={\widehat{A}}^{(1D)}$$(26)

where Q is the weight matrix. The last term favours  to be close to prior vector A₀, with a positive regularisation hyper-parameter r.

Minimising this ℰ(A|r_c, P) function is still a linear regression for the A parameters, whose optimal value now is

$$\widehat{A}={\left(M^{T}WM+rQ\right)}^{-1}\left(M^{T}WD+rQ{A}_0\right).$$(27)

The covariance matrix associated with  directly is

$$C={\left(M^{T}WM+rQ\right)}^{-1}.$$(28)

We tested different Q matrices on CTIO simulations, and the matrix that gives the most satisfying results when  is compared to the true amplitudes uses the Laplacian operator L,

$$L=\left(\begin{array}{ccccccc}-1& 1& 0& 0& \cdots & 0& 0\\ 1& -2& 1& 0& \cdots & 0& 0\\ 0& 1& -2& 1& \cdots & 0& 0\\ ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& \\ 0& 0& 0& 0& \cdots & -2& 1\\ 0& 0& 0& 0& \cdots & 1& -1\end{array}\right)$$(29) $$Q=L^T{U}^{T}UL,$$(30)

with U^TU the matrix proposed in Eq. (C.5),

$$U^{T}U=\left(\begin{array}{cccc}1/{\sigma }_{A_{1D}^{(1)}}^2& 0& \cdots & 0\\ 0& 1/{\sigma }_{A_{1D}^{(2)}}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1/{\sigma }_{A_{1D}^{\left(N_{{}_X}\right)}}^2\end{array}\right).$$(31)

The total variation regularisation is known to be able to retrieve very sharp features (e.g. steps or edges) when deconvolving an image. The Laplacian regularisation cannot do the same, but discontinuities are not expected in the physical spectra we observe. Moreover, the norm-2 regularisation offers analytical solution, while a norm-1 regularisation needs to conduct an iterative minimisation process (more details in Appendix C).

For the 1D transverse estimate, the deconvolution fit is performed using a Gauss-Newton minimisation of ℰ(A|P) for the non-linear P parameters (see Appendix B) and a linear resolution for the A parameters (see Sect. 3.5.2). The Gauss-Newton minimisation is repeated three times with a clipping rejection of bad pixels to remove field stars or CCD defects that can pull the final parameters in undesired directions. At the end of this process, we obtain the measured ${\widehat{r}}_c$ and $\widehat{P}$ parameters, and a measurement of the amplitudes Â(r) with their covariance matrix C(r). The vector Â(r) forms the searched spectrum, but it might be contaminated by the second diffraction order at this stage. To accelerate the computation, regions farther away than twice the PSF FWHM from the centroids are masked and set to zero with null weights (the grey regions in Fig. 10). For this process, a default value of the r regularisation parameter was chosen, without fully exploring how it could be optimised.

The PSF deconvolution problem was solved in another way in Ryan et al. (2018) for slitless spectroscopy. Instead of using a PSF model, the mathematical model was inverted using multiple exposures of the same spectrum, ideally taken at different orientations, dispersion directions, or dithered positions. In this way, the number of data leads to an invertible problem, but a damping hyper-parameter equivalent to r still needs to be introduced in Eq. (25). Ryan et al. (2018) tuned this hyper-parameter manually and determined the breaking point of a specific L-curve (see Hansen 1992). This method is mathematically similar to the one we used (described in Sect. 3.5.4) because it leads to an equilibrium between information from the prior and data.

3.5.4 Optimisation of the regularisation

The correct amount of information for which the data can be searched might be thought hard to determine. It may seem hopeless to find information at a spatial scale below the typical width of the convolution kernel with a unique spectrogram.

Therefore, we searched for the way to determine the optimal hyper-parameter r.

The regularisation parameter r is optimised via the study of the resolution operator

$$R=I-rCQ,$$(32)

with I the identity matrix. A fruitful interpretation of the R operator is given in Hansen (2010) with

$$\begin{array}{c}\mathrm{Tr}I=\mathrm{Tr}R+\mathrm{Tr}(rCQ)\iff \\ [\#\text{parameters}]=[\#\text{parameters\hspace{0.17em}resolved\hspace{0.17em}by\hspace{0.17em}data}]\\ +[\#\text{parameters\hspace{0.17em}resolved\hspace{0.17em}by\hspace{0.17em}prior}],\end{array}$$(33)

where the trace of the resolution operator gives the effective number of degrees of freedom that can be extracted from the data for a given amount of prior information. We set

$$N_{\text{dof }}=\mathrm{Tr}R.$$(34)

The optimal r parameter is found by minimising the following G(r) function, which behaves like a reduced χ²,

$$G(r)=\frac{{\chi }^2\left(\widehat{A}\mid {\widehat{r}}_c,\widehat{P}\right)}{{\left(N_x{N}_y-N_{\text{dof}}\right)}^2}.$$(35)

This method, known as general cross-validation (GCV), is extensively presented for example in Golub et al. (1979); Wahba (1990); Hansen (2010). It is demonstrated that the minimum of G(r) corresponds to the minimum of the distance |MÂ – MA_truth|², where A_truth is the true amplitude vector hidden in the data.

We tested different ways to implement the optimisation of the regularisation process in Spectractor. The most efficient result (in terms of rapidity and bias of the final result) is obtained by setting a default reasonable regularisation parameter for the fitting procedure of the amplitude  and PSF ${\widehat{r}}_c$, $\widehat{P}$ parameters, and finally to find the minimum of the G(r) function. We observed in the simulations that regardless of the r hyper-parameter that is chosen at the beginning, the process reconstructs an unbiased spectrum at the end of the ℰ(A|r_c, P) minimisation. The level of regularisation of the solution can thus be set a posteriori by finding the optimal r after minimising the G(r) function.

The result of a 2D deconvolution process is presented in the right panel of Fig. 10 for a simulation with a wide PSF kernel, but without second diffraction order to stress the benefit of the deconvolution. The residual map between the best-fitting model and the simulated spectrogram shows that the 1D transverse fit cannot extract the spectrum from the spectrogram image correctly (Fig. 10 left), while the 2D deconvolution ends with a quasi-unstructured residual map that respects the expected Gaussian distribution (Fig. 10 right). Because our model is informed by a PSF model, we are able to extract spectra with a single exposure, which involves rather small matrices compared with Ryan et al. (2018). This allows us to tune the r hyper-parameter automatically in a few seconds with a standard laptop.

If the spectrogram is not fully contained in the sensor area, the spectrum exhibits a discontinuity that causes the norm-2 regularisation to fail (the second derivative from the Laplacian operator is undefined). For a given instrumental PSF, it is therefore better to use a more dispersive grating to feed the deconvolution with more data and increase the wavelength resolution, but the spectrogram must fit inside the sensor area for regularisation techniques to be use.

3.5.5 The spectrophotometric uncertainty principle

The regularity of the deconvolved solution depends on the hyper-parameter r. The optimal r parameter is chosen as the minimum of the G(r) function represented in the top panel of Fig. 11 for a simple simulation with n_PSF = 0, γ_ο = 5, α₀ = 3 (and the same characteristics as in Table 1),

$$\varphi \left(x,y\mid r_c,P\right)=\frac{{\alpha }_0-1}{\pi {\gamma }_0^2}{\left[1+{\left(\frac{x-x_c}{{\gamma }_0}\right)}^2+{\left(\frac{y-y_c}{{\gamma }_0}\right)}^2\right]}^{-{\alpha }_0}.$$(36)

The second panel displays the ${\chi }^2\left(\widehat{A}(r)\mid {\widehat{r}}_c,\widehat{P}\right)$ function, which shows that the optimum A(r) solution does not minimise the agreement with the data (minimum χ²), but makes a compromise with a regularisation scheme (modelled by the ${\chi }_{\text{pen}}^2\left(A\mid A_0\right)$) penalty term). The lower panel illustrates that the effective number of amplitude parameters fitted by data with the optimum régularisation hyper-parameter is approximately 180. For this simulation, about 680 amplitude parameters were fitted in a spectrogram built with a constant PSF FWHM of about 5.5 pixels. Intuitively, we can conjecture that an optimum relation must exist between the typical width of the PSF kernel and the amount of information that can be searched for in data,

$$\begin{array}{c}[\text{PSF\hspace{0.17em}width}]\times [\#\text{effective\hspace{0.17em}degrees\hspace{0.17em}of\hspace{0.17em}freedom}]\\ \approx [\#\text{parameters}].\end{array}$$(37)

If too many parameters are searched for, the problem becomes ill-posed. If too few parameters are searched for, it should be compensated for with a wide PSF kernel. We thus postulate that we should have a spectrophotometric uncertainty principle of the type

$${\sigma }_{\text{PSF}}\times \frac{N_{\text{dof }}}{N_x}\gtrsim h,$$(38)

where the optimum is reached at equality, and h has the same units as σ_PSF the width of the PSF kernel. This formula gives the minimum number of the degrees of freedom required to describe data given a PSF width.

We tested the deconvolution and regularisation process on a large number of simulated spectrograms with a constant width (n_PSF = 0), but without second diffraction order. For the sake of simplicity, we tested this without any loss of generality. We tried a Gaussian PSF kernel and a Moffat kernel with two different exponents (α₀ = 2 and 3) and various γ₀ values. The results are summarised in Fig. 12. For the three models, σ_PSF was chosen to be half of the PSF FWHM. The figure shows that for any PSF kernel, the measure of the number of degrees of freedom N_dof/N_x scales as the inverse of the PSF width. They show a definite trend for the product σ_PSF × N_dof/N_x that appears to be asymptotically constant and equal to the number h close to 0.8 pixel when the PSF size is significantly greater than a few pixels. This h value varies with the S/N of the spectrogram, but for a situation, it locks the relation between σ_PSF and N_dof. This flatness of the relation shows that our procedure agrees when the extraction of information at the scale of the PSF kernel is considered optimal.

It is also noteworthy that this equation could be exploited to accelerate the computation of the PSF cubes: Instead of computing a PSF kernel for each pixel column i, we could compute it for each N_dof/N_x pixel. Because the computation times were not an issue in this paper, we left this investigation for a future project in which computation speed counts.

3.6 Wavelength calibration

Despite the astigmatism of the system, to first approximation, the slitless spectrograph obeys the usual grating formula (Eq. (1); see e.g. Murty 1962; Hall 1966; Schroeder & Inc 2000). Using the notations of Fig. 1, the grating formula can be inverted to find the relation between the u coordinate along the dispersion axis and λ.

First, we assume that the true zeroth-order position is at u₀ along the dispersion axis, but that a misfit of its centroid (see Sect. 3.2) can shift the position by a quantity $\delta u_0^{(\text{fit})}$.

The ADR also slightly spreads the zeroth-order image along the local constant azimuth line in a deterministic way depending on the airmass and the spectrum of the source. It also depends on the parallactic and camera angles, the atmosphere temperature, pressure, and humidity. It is incorporated in the wavelength calibration process as a wavelength-dependent shift δu^(ADR)(λ) of the PSF centroid position along the dispersion axis with respect to a reference wavelength λ_ref: δu^(ADR)(λ_ref) = 0.

We model this effect using the NIST metrology toolbox¹¹ recommendation of using a modified version of the Edlén equation (Edlén 1966) by Birch and Downs (Birch & Downs 1993, 1994) (see Appendix D).

The distance d(λ) between an abscissa of the spectrogram and the zeroth order then reads

$$\delta u(\lambda )=\delta u_0^{(\text{fit})}+\delta u^{(\text{ADR})}(\lambda ),$$(39) $$u_0=D_{\text{CCD}}\tan {\theta }_0\text{, }$$(40) $$d(\lambda )=u(\lambda )-u_0-\delta u(\lambda ),$$(41)

so that

$$\begin{array}{r}\hfill d\left(\lambda \mid D_{\text{CCD}},\delta u_0^{(\text{fit})}\right)=D_{\text{CCD}}[\tan \left(\arcsin \left(p{N}_{\text{eff}}\lambda +\sin {\theta }_0\right)\right)\\ \hfill -\tan {\theta }_0]-\delta u^{(\text{ADR})}(\lambda )-\delta u_0^{(\text{fit})}.\end{array}$$(42)

The bijection between the position on the CCD and the wavelength is thus parametrised with two unknown parameters, D_CCD and $\delta u_0^{(\text{fit})}$, that need to be fitted.

As a starting point, we compute a first wavelength array λ₀ from the array of distances d to the zeroth order along the dispersion axis, assuming δu₀ = 0 and given a prior value of D_CCD. To obtain a wavelength array given D_CCD and δu₀, Eq. (42) is inverted as

$$\lambda =\frac{1}{p{N}_{\text{eff}}}\left[\sin \arctan \left(\frac{d+\delta u^{(\text{ADR})}(\lambda )+\delta u_0^{(\text{fit})}}{D_{\text{CCD}}}\right)-\sin {\theta }_0\right].$$(43)

To remove the ambiguity with ADR, which also depends on wavelength, we iterate this computation five times starting from λ₀ and updating À. We verified that it is enough to converge toward a stable wavelength solution. In these steps, we associated a wavelength array λ with the amplitude array A. From this calibration, λ_ref is computed as the mean wavelength weighted by the spectrum A itself, in order to associate the maximum amplitude of the zeroth order with its brighter wavelength. With the fit of $\delta u_0^{(\text{fit})}$ and this setting of λ_ref, we ensure that we are not sensitive to the slight dispersion of the zeroth order itself.

The parameters D_CCD and $\delta u_0^{(\text{fit})}$ are fit on data using the most prominent absorption (or emission) lines on the observed stellar spectrum (typically, the hydrogen lines Hα and Ηβ, and the dioxygen lines at 762.1 nm and 686.7 nm). They are locally fitted with a polynomial background plus a Gaussian profile of unknown height, centroid, and width. A partial χ² quantity is computed for each spectroscopic line and added into a global χ².

A penalty defined as the squared distance between the fitted Gaussian centroids and the tabulated values for the detected lines, weighted by the squared S/N, is then added to the χ².

Finally, the full χ² and its penalty are normalised to the number of detected lines. This normalisation of the global χ² is necessary to avoid solutions that favour a lower number of detected lines, while the penalty gives more weights to the well-detected lines and anchors them on their tabulated values. The two parameters $\delta u_0^{(\text{fit})}$ and D_CCD are varied to minimise the penalised global χ² and find the best solution for the wavelength calibration.

The result of this process is illustrated in Fig. 13. At the top, the global χ² is represented for the wavelength calibration of the planetary nebula PNG321.0+3.9 observed at CTIO.

The sharp steps at high or low D_CCD values reflect situations in which some emission lines are detected or are not detected, which emphasizes the need to normalise the global χ² by the number of detected lines. The smoothness around the minimum is due to the penalty term on $\delta u_0^{(\text{fit})}$. In the calibrated spectrum, we can observe many emission lines that have been detected by the algorithm, and a good alignment between the tabulated values (represented by the vertical lines), the extrema of the Gaussian profiles, and those of the data curve. A summary of the detected lines with high S/N is presented in Table 2, reporting their fitted wavelength and their equivalent width (EQW).

3.7 Flux calibration

With the fitted wavelength solution, the spectrum amplitude  can be converted from ADU units into flux densities in erg s⁻¹ cm⁻² nm⁻¹, assuming that the telescope collecting area S_T, the exposure time τ, and the CCD gain G_CCD (in e⁻ ADU⁻¹) are known,

$$S_1(\lambda )=\widehat{A}\times \frac{hc{G}_{\text{CCD}}}{{𝒮}_T\tau \lambda \delta \lambda },$$(44)

where δλ is the local variation in λ within one pixel.

The end product of this pipeline is thus a background-subtracted spectrum, calibrated in wavelength and flux, that is the product of three quantities: the object SED, the instrumental transmission, and the atmospheric transmission, which might be contaminated by a second diffraction order because this latter has not yet been taken into account.

3.8 Full forward model of the spectrogram

At the end of the previous steps, we also have first estimates of the two instrumental model functions θ₁(r, λ) and Δ_p(λ), and of geometric parameters such as the zeroth-order position r₀ and the dispersion angle α. With all these ingredients, we can implement a full forward model of the data that also takes the atmospheric differential refraction (ADR) and the superposition with the second diffraction order (last stage of Fig. 3) into account.

In practice, we enrich the forward model described in the steps above with the knowledge of the ADR physics (see Sect. 3.6 and Appendix D) and with the knowledge of the second-order to first-order transmission ratio r_2/1(λ) of the spectrograph disperser. The ADR model replaces a polynomial approach to predict Δ_p(λ). In other words, it is used to predict the trace of the spectrogram on the sensor without free parameters, as long as airmass, outside pressure, outside temperature, and humidity are given. With this, two free parameters remain for the spectrogram trace on the sensor to be fully constrained: the dispersion axis angle α and δy^(fit), which compensates for a misfit of the zeroth-order centroid along the y-axis. The ratio r_2/1 (λ) can be measured on an optical test bench (see Sect. 5.1) or on on-sky data. In the full forward model, we use a new design matrix MM that is defined as

$$\tilde{M}=M\left(Z\mid r_c,P\right)+A_2{R}_{2/1}M\left(Z\mid r_c^{p=2},P^{p=2}\right),$$(45)

where A₂ is a safety normalisation parameter, R_2/1 is the transmission ratio vector computed for a given wavelength calibration, $r_c^{p=2}$ are the centroid positions of the second diffraction order PSF kernels, and P^p=2 are their shape parameters. The vector $r_c^{p=2}$ is computed using the grating formula (1) and the ADR model. P^p=2 can be fitted independently of P^p=1 but we chose to assume that the PSF shape depends more on the spectrograph defocusing towards the infrared than on the atmospheric chromatic seeing. The PSF shape parameters for the second diffraction order are thus considered the same as for the first diffraction order at the same distance of the zeroth order. We chose to set the P^p=2 vector accordingly¹². Therefore, the full forward model now includes both first- and second-order spectrogram models as

$$I_1\left(Z\mid A,r_c,P\right)+I_2\left(Z\mid A,r_c^{p=2},P^{p=2}\right)=\tilde{M}A.$$(46)

We implemented a two-step iterative method that alternates the wavelength calibration described in Sect. 3.6 and the full forward model described just above (by combining a linear fit for the A spectrum amplitudes and a Gauss-Newton descent for the non-linear parameters P) with the same r hyper-parameter that was fitted (in Sect. 3.5.4). The A parameter vector determined with PSF 2D deconvolution previously (Sect. 3.5.3) is seeded in the forward-model fit as a new prior A₀. A 20σ clipping is performed to reject the field stars and their concomitant spectrograms as well as other sensor defects. This procedure ensures that all the forward-model parameters are fitted again together on data, using the more complete model including A, P, D_CCD, $\delta u_0^{(\text{fit})}$, and A₂. Their values replaced all those that were fitted previously. The residual map obtained in this way is flatter than before, with an even smaller final χ², and consequently, with a better accuracy of all the fitted parameters.

The two main products of this step are a first-order diffraction spectrum S₁ (λ) separated from the second-order spectrum, and the second-order spectrum S₂(λ) with

$$S_2(\lambda )=r_{2/1}(\lambda )S_1(\lambda ),$$(47)

in erg s⁻¹ cm⁻² nm⁻¹ (following Eq. (44)).

At this point, we consider the second diffraction order not as a nuisance, but as a useful signal. With a strong bending due to ADR (e.g. with the dispersion axis orthogonal to the zenith direction), it can be detached from the first diffraction order on purpose to maximise the statistical power of the exposure.

In summary, a full forward model takes advantage of the higher diffraction orders as a redundant piece of data to fit all parameters, especially in the bluer part, where absorption lines are twice wider in pixels than in the first-order spectrum. This is a key advantage of the forward approach compared to the direct approach.

3.9 Validation on simulations

To test the full forward model, we simulated a spectrogram with a second diffraction order and a sharp Moffat PSF kernel to increase the spectral resolution (see the PSF parameter values in Table 3), whose shape evolves as a second-order polynomial function. Figure 14 compares the simulated data with the fitted spectrogram model and focuses on some atmospheric absorption lines: The residuals follow the expected Gaussian distribution again, even in the fast-varying regions of the spectrum.

In Fig. 15, we show that the process recovered the true spectrum injected in the simulation within the estimated uncertainties (diagonal elements from the C matrix from Eq. (28)). The agreement is excellent at all wavelengths, even around the fast-varying absorption lines. The right panel of the figure also shows the FWHM of the true PSF. The reconstructed PSF displays the same wavelength-dependent PSF FWHM as the simulated one. It is also remarkable that while the cross-spectrum issued from the transverse 1D fit described in Sect. 3.5.2 failed to recover the true spectrum and the true PSF profile because of the presence of the second diffraction order (orange curves), it still proved to be an important seed for the regularisation process because only its regularity is used because of the Laplacian operator L.

The recovered parameters are compared with the simulation values in Table 3. They are fitted together with their uncertainties in the full forward-model minimisation, which provides their full covariance matrix (Fig. 16), while the other parameters such as the star centroid are just estimated on data.

The regularisation quantities are given in Fig. 17. For a PSF FWHM between 2 to 4 pixels, we again obtain N_dof ≈ 300 out of ≈ 700 parameters. This confirms the rule of thumb given by the spectrophotometric uncertainty principle. For spectrograms such as those presented in the previous figures, the end-to-end pipeline takes 2 min on a standard laptop.

The Spectractor implementation was tested on many simulations and recovered the simulation parameters within the estimated uncertainties (68% confidence interval) for sets of parameters that were not too extreme (smooth wavelength dependence on the PSF, or the PSF kernel sampled over a few pixels). We also evaluated the extraction bias b between the true spectrum given in the simulation $S_1^{\text{truth }}(\lambda )$ and the extracted spectrum

$$b=\frac{{\displaystyle \int \text{d}}\lambda \left(S_1(\lambda )-S_1^{\text{truth }}(\lambda )\right)}{{\displaystyle \int \text{d}}\lambda S_1^{\text{truth }}(\lambda )}.$$(48)

The extraction bias was evaluated with many simulations in different cases in terms of S/N, resolution, and geometry. For the first case, the variation in S/N was simulated by multiplying the simulated spectrum by an arbitrary grey factor A₁, keeping the image background at the same level. The S/N of the simulation presented in Fig. 15 and Table 3 corresponds to A₁ = 1. We found no significant bias in the A₁ > 1 regime (Fig. 18), but a small bias appears for spectrograms with a lower S/N at the percent level. This is because a Gaussian model was used in the evaluation of the uncertainty map (see Sect. 3.1) while a Poisson distribution is more accurate at low S/N. The variation in spectrogram resolution was simulated by changing the PSF width γ₀, and a small negative bias appeared for low-resolution spectra (large γ₀). These spectra exhibit wider and shallower absorption lines than the true spectrum, which leads to these negative b values, but only at the sub-percent level. However, except for the absorption lines, the overall spectrum shape from blue to red is recovered perfectly. Finally, geometry variations were also simulated using different dispersion axis angles α, but we found no bias. In conclusion, the most important condition for extracting unbiased spectra from slitless spectrophotometry is a sufficient S/N, closely followed by a sufficiently fine spectral resolution. The adequacy of these parameters is easy to test a priori with a forward-model simulation.

For completeness, we represent the reconstructed number of degrees of freedom for all these simulations in Fig. 18. As expected, for spectrograms with a very low S/N, this number is close to zero. In that case, the extracted spectra are close to the A₀ = Â^(1D) prior spectra. When the signal increases, then N_dof strongly increases until it saturates because of pixelisation. Conversely, this number decreases when the PSF width increases because the spectrogram has a lower spectral resolution.

4 Spectrum extraction on data

The success of the spectrum extraction on data mostly depends on the model of the wavelength-dependent PSF of the telescope. If the PSF model correctly represents the reality, the residuals after the full forward-model fit of the spectrogram converges towards the expected Gaussian distribution. Otherwise, the extracted spectrum is distorted when the PSF is too far from reality.

This is illustrated in Fig. 19. In the left panels, a blazed Thorlabs grating with 300 lines mm⁻¹ was chosen to observe CALSPEC star HD111980 at CTIO on 2017 May 30. Directly placed in the filter wheel at D_CCD ≈ 55 mm from the CCD, this grating presents a rather strong defocusing that is poorly modelled by our default circular Moffat PSF, even with a fourth-order polynomial evolution with wavelength.

The building of an appropriate model for this highly defocused PSF is left for future work and is presented here as an illustration that the prior understanding of the telescope can be greatly worthwhile in a forward-fitting approach. These extractions used a sigma-clipping procedure with a 20σ threshold in order to reject only the field stars and not the poorly modelled spectrogram pixels.

However, the same PSF kernel as was used to treat the same star observed 5 min later, but with an amplitude hologram optimised to correctly focus the spectrogram at all wavelength (Moniez et al. 2021), leads to residuals between ±5σ (right figures), mostly dominated by a field star that contaminates the spectrogram at about 530 nm.

The parameters of interest for these two extractions are summarised in Table 4. We realised a posteriori that the S/N was not sufficient to fit A₂, and decided to keep it fixed at 1. The way in which the ratio r_2/1(λ) was obtained for these exposures is explained in Sect. 5.

The calibrated spectra given by Spectractor at the end of the extraction process are given in Fig. 20. Because of the strong defocusing, the spectrum from the Thorlabs grating presents broadened absorption lines, while the amplitude hologram has sharper absorption lines. The second displays a better spectral resolution that is limited by the atmospheric seeing, which argues in favour of either controlling the spectrograph PSF at the hardware level, or of being able to accurately model a defocused PSF kernel. At this point, it seems easier to adjust the spectrograph to obtain a simple PSF model than to guess the complexity of the chromatic PSF with on-sky data. Based on these spectra or spectrograms, we show in Sect. 5 how to measure the atmospheric transmission or the instrumental transmission via forward modelling.

5 A path toward measuring atmospheric transmission

One of the main objectives for building a spectrophotometric instrument and its analysis pipeline is to be able to accurately measure on-site atmospheric transmission so as to improve the photometric calibration of other telescope on the same site. For instance, the aim of the Auxiliary Telescope at Cerro Pachón is to measure the atmospheric transmission to correct the photometry of the LSST survey.

In order to discuss the capabilities of Spectractor in measuring atmospheric quantities, we first recall that its main output is a first diffraction order spectrum,

$$S_1(\lambda )=T_{\text{inst },1}(\lambda )T_{\text{atm}}\left(\lambda \mid P_a\right)S_{*}(\lambda ),$$(49)

To be able to obtain the atmospheric transmission T_atm(λ|P_a), we need to know the stellar SED and the instrumental transmission, and to have an accurate full forward model for the spectrograph.

Because the most accurate PSF model is achieved with the amplitude hologram at CTIO because of its focusing properties, we expect better results from its analysis than from the data acquired with the Thorlabs grating.

In order to inform our forward model, we need to know the disperser on-sky first-order transmission and their r_2/1(λ) ratio. While we show below how to use on-sky data to this end, we also secured the Thorlabs grating to bring it back on an optical bench at the Laboratoire de Physique Nucléaire et de Hautes Énergies (LPNHE) and measure its transmission. This was not possible for the prototype hologram used at CTIO, and we had to recover its transmission with on-sky data.

While in theory, the forward modelling of the atmospheric transmission could be based on a perfect a prior knowledge of the instrument and simply needs to fit each star exposure, in practice, all this is slightly more complex. We therefore decided to show that our pipeline can be used with intermediary steps in order to gain increasingly better understanding of the data, up to the point where the full forward modelling of the atmospheric parameters becomes possible. Because of the limited telescope time and observations, this first paper relies on a limited amount of data and aims at presenting the algorithms and procedures. Clearly for accurate results, many more data are required.

The different steps that we undertook and detail below are as follows:

• Sect. 5.1: laboratory measurement of the blazed grating transmission as a function of λ;

• Sect. 5.2.1: inference of the CTIO 0.9 m telescope transmission using data taken with the blazed grating during a stable night;

• Sect. 5.2.2: inference of the amplitude hologram transmission during the same stable night with the CTIO 0.9m telescope transmission;

• Sects. 5.2.3, 5.3: measurement of T_atm(λ|P_a) with data using the amplitude hologram with Τ_inst,1(λ), S_*(λ), and a Moffat PSF kernel.

5.1 Disperser transmission measurement

Our blazed Thorlabs 300 lines mm⁻¹ grating was studied 2 yr after the CTIO campaign on the LPNHE optical test bench. The optical bench is made of parabolic off-axis mirrors to simulate a monochromatic f / D = 18 telescope beam. A monochromator is used to select an accurately known narrow wavelength interval, and the light, after passing through the grating mounted on an xyz support, is collected on a CCD device.

The transmissions for orders 0, 1, and 2 were measured using aperture photometry at many wavelengths. Because the laboratory bench had been designed to measure filter transmissions, it could unfortunately not be adjusted to allow a S/N that was high enough in the regions below 450 nm and above 1000 nm. For these wavelength ranges, we therefore used the spreadsheet of the grating manufacturer.

To measure the r_2/1(λ) ratio, we used the optical bench measurement above 450 nm and the CTIO on-sky measurement from Moniez et al. (2021) below.

In order to measure wavelengths bluer than < 400 nm, we extrapolated the r_2/1(λ) function with an exponential model C exp [−(λ − λ₀)/τ] with three free parameters C, λ₀, and τ fitted on the laboratory data.

The first-order efficiency curve and the r_2/1(λ) curve for the blazed Thorlabs 300 lines mm⁻¹ grating are represented in Fig. 21 and were used in Spectractor when spectra taken with this grating were measured (e.g. in Fig. 19 left).

5.2 Analysis of a photometric night

To measure the CTIO 0.9 m telescope transmission, we made use of a set of spectra acquired at different airmasses during a night with stable photometric condition. The multiplicity of the airmass conditions and the hypothesis that the atmospheric transmission spectrum only varies with the quantity of atmosphere between the source and the observer allowed us to factorise the atmospheric transmission as an airmass-dependent term and the instrumental transmission term.

In order to distinguish between the average spectrum of the atmospheric transmission and the instrumental transmission spectrum, we performed the fit over the available N_s spectra by simulating the atmospheric transmission with Libradtran¹³ (Mayer & Kylling 2005; Emde et al. 2016) jointly with N_i arbitrary coefficients to sample the T_inst,1(λ) curve.

At CTIO on 2017 May 30, the night presented very stable conditions according to the in situ meteorological measurements of temperature, pressure, humidity, and also a stable seeing around 0.8″. We analysed the spectra of CALPSEC star HD111980 acquired with the Thorlabs and holographic dispersers under the hypothesis that this night was photometric. The observations cover an airmass range from 1 to 2 (see Fig. 22).

For each of the dispersers, N_s spectra were acquired and extracted. They were averaged in 3 nm bins, for which the instrumental transmission is assumed to be very smooth at this scale. The main atmospheric and hydrogen absorption lines were masked in this process.

5.2.1 Analysis of a photometric night in which the instrumental transmission was extracted

We present here how we estimated the telescope transmission from a photometric night data set and the knowledge of the laboratory measurement of the Thorlabs grating. Basically, we fitted a telescope transmission model and one atmospheric model on the collection of spectra extracted from the data collected with this disperser.

From 300 nm to 1100 nm, we simultaneously fitted the N_s = 20 good blazed grating spectra observed with the model from Eq. (9) for p = 1,

$$S_1(\lambda )=T_{\text{inst },1}(\lambda )T_{\text{atm }}\left(\lambda \mid P_a\right)S_{*}(\lambda ),$$(50)

where S_*(λ) is the binned SED of the CALSPEC star, and T_inst,1 (λ) is a vector of N_i = 250 free linear amplitude parameters.

The Libradtran atmospheric transmission simulation T_atm(λ) uses the in situ pressure, temperature, and airmass given in each exposure metadata. In addition, three common parameters P_a were fitted:

• the precipitable water vapour (PWV; in mm);

• the ozone quantity (in dobson db);

• the vertical aerosol optical depth (VAOD).

In order to account for a possible small grey variation in the atmospheric transmission, each spectrum was weighted by a grey factor $A_1^{(n)}$, with their average constrained to one,

$${\langle A_1^{(n)}\rangle }_n=1.$$(51)

The χ² we need to minimise thus reads

$${\chi }^2={\displaystyle \sum _{n=1}^{N_s}{\left[D_n-A_1^{(n)}{S}_1(\lambda )\right]}^T}{C}_n^{-1}\left[D_n-A_1^{(n)}{S}_1(\lambda )\right],$$(52)

where D_n is the data vector for spectrum number n, and C_n is its covariance matrix estimated by the Spectractor extraction pipeline. The $A_1^{(n)}$ and P_a parameters were fitted jointly via a Gauss-Newton descent and come with their covariance matrix, while the T_inst,1(λ) linear parameters were computed analytically via the usual algebra at each descent step. Because the spectrum of the star is assumed to be known, no regularisation is needed. Because the instrumental transmission is assumed to be smooth, the descent was repeated with a 5σ clipping to remove outliers.

This procedure has been tested on simulations, and we verified that it recovered the injected parameters for instrumental transmission, grey factors, and atmospheric quantities within the uncertainty ranges.

The results obtained on the CTIO data are presented in Fig. 23. Approximately 15% of the 5000 spectral data points are masked, either because they are close to a spectral line, or because they are 5σ outliers. The residuals are structured below the 2σ level in the red part of the spectra, either because of an incorrect PSF model for the redder wavelengths due to defocusing, or because of PWV variations in the atmosphere, as hinted by the spectra, vertically ordered in time.

The best T_inst,1(λ) solution that we extracted is presented in Fig. 24. The black points represent the raw fitted T_inst,1 (λ) vector, and the red curve is smoothed with a Savitzky-Golay filter of order 1 and window size 17. The error bars result from the combination of raw uncertainties from the fit plus the difference between the smoothed curve and the scattered raw black points.

This leads to larger uncertainties for the instrumental throughput where the spectra were masked, around the main absorption lines.

The transmission curve presents the expected decreases due to the loss of efficiency of the CCD. Given the measurement of the blazed Thorlabs grating at the laboratory (Fig. 21), we extracted the CTIO 0.9 m instrumental throughput from the T_inst,1 (λ) smoothed curve.

This fills the lack in a priori knowledge of the telescope throughput to inform our forward model, and it is used in the following analysis. We obtained the telescope throughput by dividing the fitted instrumental transmission by the first-order efficiency of the blazed Thorlabs grating. The atmospheric transmission results are detailed in Sect. 5.2.3.

A more accurate estimate of the instrumental transmission would need more data, both to inform a better PSF model, and to constrain the atmospheric transmission variations better.

Because the available data were limited, our goal was to illustrate that a forward-model approach can be adjusted to gain more information about the different components of the model.

We find it also noteworthy that the procedure is symmetric with respect to atmospheric transmission and telescope transmission: The need for the Libradtran model as a priori to constrain the atmospheric transmission shape could be replaced by the a priori measurement of the telescope transmission.

5.2.2 Analysis of a photometric night in which amplitude hologram transmissions were obtained

The next step needed to further inform our forward model in order to use the best-quality data to constrain the atmospheric transmission was to estimate the holographic disperser transmission. If this transmission is known, the data gathered with this disperser can be used, and we can take advantage of the fact that its PSF can be fairly well modelled with a Moffat.

In order to obtain the hologram transmission, we used the same procedure as described for the Thorlabs disperser above on data that were collected during the same photometric night, but with the holographic disperser. However, the ratio r_2/1(λ) is still a prior information needed for the full forward model. For the holographic disperser, we built r_2/1(λ) using the interpolated on-sky data presented in Moniez et al. (2021 Fig. 21) alone.

The N_s = 27 spectra are presented in Fig. 22, and the results are presented in Fig. 25. Approximately 0.5% of the 6723 data points are masked. The residuals are below 3σ. A deep absorption feature is visible around the water-absorption band at 950 nm.

As in the previous section, from the T_inst,1(λ) best fit (see Fig. 26), we deduced the transmission of the first diffraction order for the holographic disperser, using the CTIO 0.9 m telescope transmission curve determined previously. We are well aware of the systematic errors present in these results, and we stress that they are presented here to illustrate that the forward model approach we implemented can be used when a priori information about crucial components of the model is lacking.

5.2.3 Analysis of a photometric night in which atmospheric parameters were extracted

In addition to the instrumental transmissions of both dispersers, the procedures above also yield the parameters describing the mean atmospheric transmission of the night. Under the assumption that the night was photometric, these results are presented in Table 5.

The rather low value of the reduced χ² for the amplitude hologram illustrates the focusing properties of this disperser, which allow us to describe its PSF quite accurately with a simple 2D Moffat. Quantities obtained from the blazed Thorlabs grating data show lower statistical uncertainties than amplitude hologram data because their S/N is much higher (because its transmission is much higher). However, they certainly show higher unevaluated PSF systematics than the hologram measurements. The difference between the two estimates of the atmospheric transmission in Table 5 leads to variations in synthesised broad-band magnitudes for the LSST filters of about 8 mmag in the u, g, and r filters, 3 mmag in i, 1.5 mmag in z, and 4 mmag in y filter for various standard CALSPEC SEDs and supernovae at redshift 0. The millimagnitude accuracy on atmospheric transmission can thus be reached provided that the accuracy of the atmospheric parameters reaches below the difference shown in Table 5: We found that PWV must be fitted with an accuracy better than ≈ 0.05 mm to obtain a milli-magnitude accuracy in y band. For VAOD, uncertainties of about 0.001 are required for the u, g, and r bands. For ozone, a 10 db precision is enough to obtain milli-magnitude precision in r band.

Furthermore, the ozone and VAOD parameters we fitted are similar to the estimates of the global meteorological network MERRA-2¹⁴ (Gelaro et al. 2017) for the CTIO site during that night. The MERRA-2 PWV value ranges from 4 to 5 mm during the night of 2017 May 30. As MERRA-2 averages atmospheric quantities in 60 km wide cells, it can be expected that quantities with large local variations such as water vapour could differ from on-site measurements. This is even more true for CTIO, which is located at the top of a Chilean mountain. On the other hand, a high-atmosphere quantity would be expected to depart less between on-site and satellite measurements. We report the MERRA-2 values and compare them to what we extract with our forward model to illustrate that with a detailed knowledge of the telescope, the challenging problem of on-site atmospheric transmission measurement can be solved. While the quoted error bars only propagate the statistical uncertainties and are probably dominated by systematics, the tentative concordance between the parameters measured by MERRA-2 and our forward-model results supports the algorithm we developed.

For completeness, we also present the evolution of the grey parameters $A_1^{(n)}$ through the night in Fig. 27 for the holographic disperser data, together with the final correlation matrix of the fitted parameters in Fig. 28. The variation in the 27 $A_1^{(n)}$ factors is lower than 1%. This supports the first-order approximation of the night as being photometric and again shows that the procedure is able to improve our understanding of the data. It also offers venues to improve the model.

Finally, we note that the correlation matrix shows that the VAOD aerosol parameter is particularly strongly correlated to the spectrum amplitudes. This is expected because this quantity is mostly determined by the spectrum slope for λ ≈ 400 nm. Any systematic on the amplitude of the spectrum therefore directly affects the estimate of aerosols.

5.3 Atmospheric forward-model approach

After illustrating that the forward-model approach can be used to measure the telescope and disperser transmissions, which yields a set of stellar and atmospheric transmission spectra, we can proceed one step further. Assuming that our measurement of these crucial components of the forward model had been done with enough data to be accurate, we might skip the part in which the spectrum is extracted and directly fit the atmospheric parameters on the raw spectrogram.

At the cost of having access to the transmissions described above, when we model the spectrum S₁ (λ) as the product of a known instrumental transmission T_inst,1(λ), a Libradtran atmospheric model T_atm(λ|P_a), and the SED S_* (λ) of a known CALSPEC star, we can describe any observed spectrogram as

$$I\left(Z\mid A,r_c,P\right)=\tilde{M}\left(Z\mid r_c,P\right)A$$(53) $$A(\lambda )=A_1{T}_{\text{inst },1}(\lambda )T_{\text{atm }}\left(\lambda \mid P_a\right)S_{*}(\lambda ),$$(54)

with A₁ a grey factor.

As before, the parameters and their covariance matrix were estimated via a Gauss-Newton descent by minimising a χ² calculated over a single spectrogram. The fitted parameters were A₁, A₂, δy^(fit), Pa, D_CCD, α, and all the polynomial coefficients that model the wavelength dependence of the PSF kernel. Each spectrogram was fitted independently. We call this a spectrogram fit. As a comparison, and as a way to assess the quality of the stellar spectrum forward model, we also directly fit S₁(λ) and the atmospheric transmission on the stellar spectra extracted at this step. We call this a spectrum fit.

5.3.1 Qualification on simulations

The direct extraction of atmospheric parameters from a spectrogram was tested on simulations. The parameters we chose to simulate were the extracted parameters found from fitting all data spectrograms of the photometric night of 2017 May 30, involving the amplitude hologram. With these parameters, we simulated spectrograms of a CALPSEC star spanning an airmass ≈ 1 to ≈ 2 with the same seeing and atmospheric conditions as that of our data. We arbitrarily fixed the unknown P_a parameters to 300 db for ozone, 0.03 for VAOD, and 5 mm for PWV.

The result of the spectrogram fit is very similar to what was presented in Fig. 10. For comparison, we present the result of the spectrum fit of the extracted spectrum from one of the simulated spectrograms in Fig. 29. The extracted spectrum (red points), the best-fitting spectrum model (blue), and the true spectrum (green) all agree within the quoted uncertainties.

In addition, the recovered atmospheric values are compatible with the true injected values within the uncertainties, with a strong correlation between the grey parameter A₁ and the aerosols, as seen before on real data. The nightly behaviour is presented in Fig. 30. All values agree with the true values for both methods and correctly account for the variable simulated conditions.

We again note the strong correlation between the VAOD and the A₁ parameter. As mentioned before, this is an expected behaviour because aerosols specifically affect the spectrum slope in the blue and the global spectrum amplitude.

In addition to validating the spectrogram fit, theses results also show that the forward-model process and all the pipeline steps presented above do not bias the measurement of the atmospheric parameters.

5.3.2 Data analysis

The individual fits of the spectrograms and spectra extracted from CTIO data are presented in Figs. 31, 32, and 33.

The atmospheric forward modelling fits the data at the 5σ uncertainty level, but the PSF model imprints structured residuals similarly to what happens in the full forward-model case (Fig. 19). This effect is visible throughout the spectrogram and inside the dioxygen absorption line.

The spectrum presented in Fig. 32 is globally well fitted by the S₁(λ) model. The fit residuals around the main dioxygen line for data and for the simulation are compatible with the instrumental throughput uncertainties. All spectrograms and spectra of the night show the same residual patterns. The two methods yield very similar values for the atmospheric parameters. The values of the spectrogram fit are smooth in time, with a visible correlation between VAOD and A₁ parameters, while the spectrum fit values are shifted and more scattered, probably due to the higher sensitivity of this simpler procedure to outliers such as the field-star contamination of the spectra around 530 nm.

We recall that in the spectrogram fit, the raw spectrogram data are directly fitted with a model that contains the instrumental transmission for diffraction orders 1 and 2, a Libradtran atmospheric model, and models for the dispersion relation and PSF kernel.

At this point, the smoothness of the atmospheric parameter curves and the reasonable values that we obtained (low ozone and a few millimetres of precipitable water vapour) are as close to reality as can be expected with the quality and size of the data set at hand.

We again acknowledge that these atmospheric results are affected by systematic uncertainties and choices (e.g. the circularity of the PSF model, the PSF size of the second diffraction order, or the blazed grating transmission model) that affect the absolute value of the quoted parameters. Unfortunately, we do not have enough data to estimate systematics and proceed. We leave an analysis of the atmospheric transmission for a future paper, for example, based on the high-quality data set promised by AuxTel, which is dedicated to measuring atmospheric transmission at the Rubin Observatory site. With its mirror with a diameter of 1.2 m, its high-quality CCD sensor, and its fast readout electronics, the spectra of CALSPEC standards or stars with similar magnitudes can be acquired at a rate of approximately one or two per minute. Its overall observation strategy is not yet defined, but it has to infer the atmospheric transmission for each LSST pointing, presumably observing a grid of stars (CALSPEC stars or others) to infer directional atmospheric transmissions and interpolate them for the LSST pointing.

6 Summary and conclusions

Slitless spectrophotometry with forward modelling opens a path towards the acquisition of spectra with imaging telescopes that can be simply transformed into spectrographs by inserting a disperser on the light path.

We demonstrated on simulations that building a forward model of a spectrogram allows for accurate spectrophotometry, with a spectral resolution that only depends on the width of the PSF along the dispersion axis. The key of the process is a regularisation algorithm, which is fed with as much as prior information as possible (regularity of the searched spectrum, PSF parametrisation, ADR, and grating efficiency). The two key functions of the model are the dispersion function Δ_p(λ) and the PSF model ϕ(r,λ), together with the knowledge of the r_p/1(λ) ratio of diffraction-order transmissions.

We exemplified that this procedure functions on real data, with tentatively very promising results. We are aware of the limits of the data set in our possession, and we exemplified that the forward-model procedure can be used to improve our knowledge of the data, and by doing so, to inform the forward model.

We can also summarise some of the important lessons learned while implementing the Spectractor pipeline as follows:

• Forward modelling provides a modular approach in which each brick is a physical or empirical model that can be changed or improved, depending on the data particularities and S/N. The residuals indicate how the model can be improved (data rules);

• When it is implemented, forward modelling can easily simulate data sets to test new algorithms;

• The second diffraction order is not a contamination, but a signal that helps recover the blue part of the first-order spectrum. It should be taken advantage of whenever possible. This in particular means that we need to rethink the common wisdom of spectroscopy by increasing the efficiency of the grating in the second order, and use a field rotator (if available) in order to more easily separate the different diffraction orders on the sensor through ADR;

• The accurate knowledge of the PSF is thus crucial and requires dedicated data and an analysis that need to be carefully budgeted for. The PSF width sets the spectral resolution and the number of degrees of freedom that can be extracted from the data, and thus, decreasing the width is crucial.

Because the main scientific driver for the development of Spectractor is the measurement of on-site atmospheric transmission, we pushed our analysis to that point. We showed that our procedure allows us to measure the on-sky telescope transmission to the direct extraction of atmospheric parameters from spectrograms.

While the atmospheric parameters are dominated by systematic uncertainties, in particular, by our partial knowledge of the instrumental transmission and of the PSF, the comparison with satellite data shows a promising tentative agreement. Work that goes beyond this paper, which was devoted to presenting the spectrophotometry method, is a more intensive study of on-site atmosphere transmission. This requires in particular access to many more data and a specific detailed analysis to obtain accurate instrumental transmissions and PSF model.

We would finally like to acknowledge that many elements of the forward model can and will be improved as new data become available. In particular, the background might be estimated directly in the forward model, including other diffraction orders, by modelling the contamination of the field-star spectrogram and integrating the chromatic flat-fielding in the model to account for pixel efficiencies, and so. If the flat screen in the observatory can be illuminated with a system of several LEDs or lasers, a cube of wavelength-dependent relative transmission might be obtained that could be directly included in the forward model. Instead of dividing the exposures by one of them, each layer of the PSF cube ϕ_p(r, λ) could then be multiplied by the flat corresponding to the correct wavelength before the spectrogram model is integrated over λ. Using the forward model, we can thus solve the exact problem of deflation in spectrophotometry.

These ideas are worth implementing in the algorithm if required by data. On the other hand, many hardware solutions can also be implemented to increase the a priori knowledge of the instrument and to greatly improve the forward-model analysis. For instance we showed that holographic dispersers such as those presented in Moniez et al. (2021) improve the focusing of the spectrogram on the sensor on the whole visible and near-infrared range, which facilitates the PSF modelling. Its narrow width also allows a better spectral resolution. Another improvement would be using a collimated beam projector (CBP) (Coughlin et al. 2016; Souverin et al. 2022) to measure the telescope transmission at the per mill level, and monitor its evolution with time.

In conclusion, we have presented the theoretical tools, together with a detailed implementation example to add spectro-photometric ability to an imager by the insertion of a disperser on the light path. This comes at some computational cost, which is readily available today, but it also requires either a priori knowledge of the instrument or dedicated data and an analysis to bring the model to the required level of accuracy.

Acknowledgements

We are grateful to the CTIO technical staff members Hernan Tirado and Manuel Hernandez for their help during our tests with the CTIO 0.9 m telescope. We also thank Mélanie Chevance for her participation to the observations and Augustin Guyonnet for fruitful advice for the CTIO image reduction. The cost of the observations have been shared by the IJCLab (IN2P3-CNRS) and the Department of Physics and Harvard-Smithsonian Center for Astrophysics, Harvard University. F.B. is part of the FP2M federation (CNRS FR 2036) and of the project Labex MME-DII (ANR11-LBX-0023-01). This paper has undergone internal review in the LSST Dark Energy Science Collaboration. The internal reviewers were Marc Betoule, Andres Plazas-Malagon and David Rubin. J. Neveu is the primary author of the paper and of the Spectractor software, leading the analysis toward measurement of atmosphere transmission measurement. V. Brémaud implemented the atmospheric differential refraction in the forward model, and wrote the pipeline to determine the CTIO telescope throughput. F. Barret brought the mathematical frame for the regularisation procedure. S. Bongard contributed with general discussions on the spectrum extraction and atmospheric physics. Y. Copin developed the theoretical framework of slitless spectrophotometry and of the forward modelling. S. Dagoret-Campagne and M. Moniez contributed with general discussions on spectrum extraction, the data taking and analysis of CTIO data. L. Le Guillou built the specific system to measure the disperser transmission on the LPNHE optical bench and measured the dispersers. P. Antilogus, C. Juramy and E. Sepulveda built and maintained the LPNHE optical bench. The DESC acknowledges ongoing support from the Institut National de Physique Nucléaire et de Physique des Particules in France; the Science & Technology Facilities Council in the United Kingdom; and the Department of Energy, the National Science Foundation, and the LSST Corporation in the United States. DESC uses resources of the IN2P3 Computing Center (CC-IN2P3-Lyon/Villeurbanne - France) funded by the Centre National de la Recherche Scientifique; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded by UK BEIS National E-infrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration. This work was performed in part under DOE Contract DE-AC02-76SF00515.

Appendix A Astrometry

It is crucial in many regards to accurately anchor the wavelength calibration. For atmospheric transmission measurement, a small shift of about 1 nm can significantly bias the estimate of the aerosol parameters. Unfortunately, a shift like this can occur because of the poor determination of the position of the zeroth order of the spectrum, which is usually saturated with long bleeding spikes. It is difficult to localise it accurately and might not be robust enough to achieve a centroid determination precision that is better than the pixel scale for every image in every circumstances.

We thus found it useful to explore how the field stars might be used to set a precise astrometry using the astrometry.net library¹⁵.

The field star centroids were first extracted from the image using the IRAFStarFinder method from the photutils library (Bradley et al. 2020), using a 5σ clipping above the threshold.

The astrometry.net solve-field function was then called: Patches of stars were compared to known asterisms to obtain the precise location of the image on the sky as well as the transformation between image coordinates and sky coordinates in the form of a World Coordinate System (WCS) description.

This procedure may not yield subpixel precision for the centroid of the target star, in particular, if it has a high proper motion. To improve the precision, we compared the first source catalogue, whose positions were converted into sky coordinates, with the star positions from the Gaia DR2 catalogue, corrected for their proper motion.

The difference between the 50 brightest field-star coordinates in the two catalogues was subtracted from the WCS solution in order to lock it on the Gaia catalogue.

We then removed the star with the largest distance to the Gaia catalogue, reapplied astrometry.net followed by the Gaia catalogue centring and repeated this operation ten times. The astrometric solution in which the distance between the image stars and the Gaia stars was smallest was kept. This procedure minimises the effect of stars with poorly reconstructed centroids (due to saturation effects) or high proper motions. It ends with a scatter that is evaluated at ≈ 0.15″ RMS (Figure 7).

Appendix B Gauss-Newton minimisation algorithm

The gradient descent to minimise a χ² using the Gauss-Newton algorithm works as follows. We consider a data set with N data points gathered into a vector D with their uncertainties (correlated or uncorrelated), represented by a matrix W. We wish to model these data with a model m(P) depending on a set of parameters P. For a parameter vector P, the χ² is defined as

$${\chi }^2(P)={(M(P)-D)}^{T}W(M(P)-D)=R^T(P)WR(P),$$(B.1)

where M(P) is the vector of the model-predicted values for the N data points. The vector R(P) is the residual vector.

In order to find the set of parameters $\widehat{P}$ that minimises the χ² function, we searched for the zero of the χ² gradient that verifies ${\nabla }_P{\chi }^2(\widehat{P})=0$. The algorithm we used is the iterative multi-dimensional Gauss-Newton method, which we describe hereafter.

We started the minimisation with a first-guess value for the parameters P₀. A Taylor expansion at first order of the ∇_Pχ² function canbe performedaround the starting point P₀ andgives

$${\nabla }_P{\chi }^2\left(P_1\right)\stackrel{P\approx P_0}{\approx }2J_0^{T}W{R}_0+2J_0^{T}W{J}_0\delta P_1+\cdots ,$$(B.2)

with δP₁ = P₁ – P₀ and J₀ = ∇_PM(P₀) the Jacobian matrix of the model evaluated at P₀. For a linear model, that is, a model that can be written as $m(P\mid Z)={\displaystyle {\sum }_{i=1}^N{P}_{i}f(Z)}$, the ∇_Pχ²(P) is exactly equal to its first-order Taylor expansion.

The zero of the function is then approached by solving the equation ∇_Pχ²(P₁) = 0,

$${\nabla }_P{\chi }^2\left(P_1\right)=0\Rightarrow P_1=P_0-{\left(J_0^{T}W{J}_0\right)}^{-1}{J}_0^{T}W{R}_0.$$(B.3)

Because the approximation comes from the Taylor expansion and because the numerical accuracy of the Jacobian matrix computation is finite, it is unlikely that the P₁ that is found exactly cancels the χ² gradient.

We then searched for the α value that minimises the χ² function along the line parametrised by the vector α₁δP₁ where α₁ is a real number. The P₁ value solution then reads

$$P_1=P_0-{\widehat{\alpha }}_1{\left(J_0^{T}W{J}_0\right)}^{-1}{J}_0^{T}W{R}_0.$$(B.4)

The process was iterated K times,

$$P_{k+1}=P_k-{\widehat{\alpha }}_{k+1}{\left(J_k^{T}W{J}_k\right)}^{-1}{J}_k^{T}W{R}_k,$$(B.5)

until a convergence criterion is reached. For example when the value of χ²(P_k) or P_k evolution with k drops below a certain threshold.

The best-fitting model is then considered to be parametrised by the kth vector: $\widehat{P}\approx P_k$. The covariance matrix of the P_k parameters is obtained as the Hessian matrix at the minimum χ²,

$$C(\widehat{P})={\left(J_k^{T}W{J}_k\right)}^{-1}.$$(B.6)

In Spectractor, we also implemented the possibility to limit the P search within given bounds (e.g. we can impose that the amplitudes are all positive). We found that these bounds help the algorithm to converge.

Appendix C Second-order penalisation

Regularisation techniques involving the total variation are often used in image analysis for denoising or deconvolution while recovering sharp edges. One way to justify the use of a penalisation with the discretised Laplacian is to understand that it entails an automatic bound on the total variation. We show this below.

We recall some notations. The complete cost to minimise is

$$\begin{array}{l}ℰ\left(A\mid r_c,P\right)={(D-MA)}^{T}W(D-MA)\hfill \\ +r{\left(A-A_0\right)}^{T}Q\left(A-A_0\right)\hfill \\ ={\chi }^2\left(A\mid r_c,P\right)+r{\chi }_{\text{pen}}^2\left(A\mid A_0\right),\hfill \end{array}$$(C.1) $$A_0=A^{(1D)}.$$(C.2)

D is the data vector, M is the design matrix, and A is the amplitude vector that gives the spectrogram and that we wish to obtain. W is the inverse of the covariance matrix, so that when we have all the true parameters for A, r_c, and P, then χ²(A|r_c, P) is the sum of the squares of the residuals, and thus is the realisation of a random variable following a χ² law with N_xN_y degrees of freedom, hence the name of the cost, χ².

The penalisation term ${\chi }_{\text{pen}}^2\left(A\mid A_0\right)={\left(A-A_0\right)}^{T}Q\left(A-A_0\right)$ is also a quadratic term, and for Q = L^TU^TUL with the Laplacian operator L = −D^TD,

$$L=\left(\begin{array}{ccccccc}-1& 1& 0& 0& \cdots & 0& 0\\ 1& -2& 1& 0& \cdots & 0& 0\\ 0& 1& -2& 1& \cdots & 0& 0\\ ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& \\ 0& 0& 0& 0& \cdots & -2& 1\\ 0& 0& 0& 0& \cdots & 1& -1\end{array}\right)$$(C.3) $$D=\left(\begin{array}{ccccccc}1& -1& 0& 0& \cdots & 0& 0\\ 0& 1& -1& 0& \cdots & 0& 0\\ 0& 0& 1& -1& \cdots & 0& 0\\ ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& \\ 0& 0& 0& 0& \cdots & 1& -1\\ 0& 0& 0& 0& \cdots & 0& 1\end{array}\right),$$(C.4)

and U is such that

$$U=\left(\begin{array}{cccc}1/{\sigma }_{A_{1D}^{(1)}}& 0& \cdots & 0\\ 0& 1/{\sigma }_{A_{1D}^{(2)}}& \cdots & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& 0& \cdots & 1/{\sigma }_{A_{1D}^{\left(N_x\right)}}\end{array}\right)$$(C.5)

we obtain ${\chi }_{\text{pen}}^2\left(A\mid A_0\right)={\left(UL\left(A-A_0\right)\right)}^T\left(UL\left(A-A_0\right)\right)$, and it is the quadratic norm (or the squared Euclidean norm) of the vector UL(A − A₀), denoted ${\Vert UL\left(A-A_0\right)\Vert }_2^2$.

When we interpret A as the discretisation of a continuous spectrogram a by setting $A^{(i)}=a\left(\frac{i}{N_x}\right)$, and when σ is a function such that $\sigma \left(\frac{i}{N_x}\right)={\sigma }_{A_{1D}^{(i)}}$, a continuous analogue of this term would be a term of the form

$$\parallel {\left(a-a_0\right)}^{\prime \prime }{\parallel }_{2,\sigma }^2={\displaystyle {\int }_0^1{\left[{\left(a-a_0\right)}^{\prime \prime }(x)\right]}^2}\frac{\text{d}x}{{\sigma }^2(x)}$$

where a, a₀ would be functions whose A, A₀ are discretisations. As a simple consequence,

$${\displaystyle \underset{N_x\to \infty }{\lim }}{N}_x^3{\chi }_{\text{pen}}^2\left(A\mid A_0\right)=\parallel {\left(a-a_0\right)}^{\prime \prime }{\parallel }_{2,\sigma }^2$$(C.6)

because $N_{x}D\approx -\frac{d}{dx}$.

The total variation distance is defined as the (weighted) norm-1 of the gradient operator as a functional term,

$$\parallel {\left(a-a_0\right)}^{\prime }{\parallel }_{1,\sigma }={\displaystyle {\int }_0^1\left|{\left(a-a_0\right)}^{\prime }(x)\right|\frac{\text{d}x}{\sigma (x)}}.$$

We also note that

$${\displaystyle \underset{N_x\to \infty }{\lim }}{\displaystyle {\displaystyle \sum }_{i=1}^{N_x}}\left|UD\left(A-A_0\right)\right|=\parallel {\left(a-a_0\right)}^{\prime }{\parallel }_{1,\sigma }.$$(C.7)

However, by a simple argument, we can show that the 2-norm of the second derivative controls the 1-norm of the first derivative, so that minimising the former means that we also minimise the latter. In order to prove it, let f(x) = a(x) – a₀(x), we obtain

$$f^{\prime }(x)=f^{\prime }(0)+{\displaystyle {\int }_0^x{f}^{″}(s)\text{d}s=f^{\prime }(0)+}{\displaystyle {\int }_0^x{f}^{″}(s)\sigma (s)\frac{\text{d}s}{\sigma (s)}}$$(C.8) $$⩽f^{\prime }(0)+{\left({\displaystyle {\int }_0^x{\left(f^{″}(s)\right)}^2\frac{\text{d}s}{\sigma {(s)}^2}}\right)}^{1/2}{\left({\displaystyle {\int }_0^x{\sigma }^2(s)\text{d}s}\right)}^{1/2}.$$(C.9)

Thus, under the supplementary constraint that $a^{\prime }(0)=a_0^{\prime }(0)$, we have

$$\parallel {\left(a-a_0\right)}^{\prime }{\parallel }_{1,\sigma }⩽{\left({\displaystyle {\int }_0^1{\sigma }^2(s)\text{d}s}\right)}^{1/2}\parallel {\left(a-a_0\right)}^{\prime \prime }{\parallel }_{2,\sigma }.$$(C.10)

The discrete analogue asymptotically in N_x → +∞ is

$${\displaystyle {\displaystyle \sum }_{i=1}^{N_x}}{\left|UD\left(A-A_0\right)\right|}_{1,\sigma }⩽N_x\sqrt{\text{Tr}\left(U^{-2}\right){\chi }_{\text{pen}}^2\left(A\mid A_0\right)}.$$(C.11)

This shows that regularisation using the weighted quadratic second-order derivative automatically ensures an upper bound of the weighted total variation norm. This means that while it is computationally much faster, regularising by the weighted quadratic norm of the second-order derivative ensures an upper bound on the regularisation via the weighted total variation norm. Because the usual advantages of the weighted total variation norm (no assumption of a second-order derivative or research of a sparse minimiser) are not important here, the choice of the weighted quadratic norm of the second-order derivative as a loss function is completely pertinent.

Appendix D Atmospheric differential refraction

The ADR mostly depends on the pressure, the temperature, and the airmass and only loosely on the atmospheric humidity.

In our wavelength-calibration process for CALSPEC stars, the absorption line that weights most in the fit is the main dioxygen line at 762.1 nm. When the ADR is not correctly modelled and taken into account in the wavelength calibration, shifts of the absorption line minima towards the blue part of the spectrum can be observed throughout the night while the airmass of the star changes.

This is illustrated in the left panels of Figures D.1 and D.2. In the right panels, the ADR effect is included in the wavelength-calibration process through a wavelength-dependent shift of the zeroth-order centroid $\delta u_0^{(\text{ADR})}(\lambda )$. This procedure absorbs most of the line shifts when the dispersion axis is not orthogonal to the zenith direction.

For completeness, we show in Figure D.3 the angle conventions that are used in Spectractor to correctly compute the zenith direction in the image.

References

---

All Tables

All Figures