Issue 
A&A
Volume 684, April 2024



Article Number  A21  
Number of page(s)  29  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202347422  
Published online  01 April 2024 
Slitless spectrophotometry with forward modelling: Principles and application to measuring atmospheric transmission^{★}
^{1}
Sorbonne Université, CNRS, Université de Paris, LPNHE,
75252
Paris Cedex 05, France
^{2}
Université ParisSaclay, CNRS, IJCLab,
91405
Orsay, France
email: jeremy.neveu@universiteparissaclay.fr
^{3}
MODAL’X, UPL, Univ. Paris Nanterre, CNRS,
92000
Nanterre, France
^{4}
Univ. Lyon, Univ. Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822,
69622
Villeurbanne, France
Received:
10
July
2023
Accepted:
20
November
2023
Context. In the next decade, many optical surveys will aim to answer the question of the nature of dark energy by measuring its equationofstate parameter at the per mill level. This requires trusting the photometric calibration of the survey with a precision never reached so far on many sources of systematic uncertainties. The measurement of the onsite atmospheric transmission for each exposure, or for each season or for the full survey on average, can help reach the per mill precision for the magnitudes.
Aims. This work aims at proving the ability to use slitless spectroscopy for standardstar spectrophotometry and its use to monitor onsite atmospheric transmission as needed, for example, by the Vera C. Rubin Observatory Legacy Survey of Space and Time supernova cosmology program. We fully deal with the case of a disperser in the filter wheel, which is the configuration chosen in the Rubin Auxiliary Telescope.
Methods. The theoretical basis of slitless spectrophotometry is at the heart of our forwardmodel approach to extract spectroscopic information from slitless data. We developed a publicly available software called Spectractor, which implements each ingredient of the model and finally performs a fit of a spectrogram model directly on image data to obtain the spectrum.
Results. We show through simulations that our model allows us to understand the structure of spectrophotometric exposures. We also demonstrate its use on real data by solving specific issues and illustrating that our procedure allows the improvement of the model describing the data. Finally, we discuss how this approach can be used to directly extract atmospheric transmission parameters from the data and thus provide the base for onsite atmosphere monitoring. We show the efficiency of the procedure in simulations and test it on the limited available data set.
Key words: atmospheric effects / instrumentation: spectrographs / techniques: spectroscopic / cosmology: observations
The software Spectractor is available at https://github.com/LSSTDESC/Spectractor
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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, highredshift supernovae were fainter in red bands than what can be inferred from lowredshift 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 forwardmodelling 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 datataking 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 fluxweighted 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 nearinfrared (nearIR) 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 perthousand 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 (AuxTel) (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 InterAmerican 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 onsky 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 forwardmodelling analysis. For example, the field of view can be small and may contain only one star, or it can be crowded, with many sourcedispersed images (socalled spectrograms hereafter) superimposed 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 chargecoupled device (CCD) sensor. In the following, the positions in the sensor plane are parametrised with the coordinates r = (x, y), and the zaxis 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_{0} 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 xaxis. 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}, 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 pointlike 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),
The angles are those of the projection in the plane perpendicular to the grating lines (see Fig. 1); θ_{0} 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^{1} 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 zerothorder image. In the following, we assume that the focus has been made on the zeroth order unless otherwise specified.
Fig. 1 Geometry of a simple slitless spectrograph in the plane orthogonal to the CCD (top) and parallel to the CCD (bottom). A convergent beam with incident angle θ_{0} is focused on a CCD at position r_{0} = (x_{0},y_{0}), but also passed through a grating with N_{eff} grooves per millimetre at a distance D_{CCD} above the CCD. The beam is deflected at an angle θ(λ) along the mean dispersion axis u, which forms an angle α with respect to the CCD x axis, but is focused somewhere above the sensor. Atmospheric refraction adds a supplementary dispersion along and transversally to the mean dispersion axis. 
2.2 Theoretical model of a spectrogram
A theoretically perfect image of a light source can be modelled by its spatiospectral flux density C(r, λ). For a point source, we can separate the spectral and spatial distributions as
with S_{*}.(λ) spectral energy density (SED) of the astrophysical object, and δ the Dirac distribution. The observed spectral and spatial distribution is then
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_{0}(r, λ) = Α δ(λ – λ_{0}) × δ (r – r_{0}), with A the received source flux at λ_{0}. 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 ϕ_{0}(r, λ) whose spatial integral is normalised to one. Therefore, by definition of the PSF, the image of a monochromatic point source centred on r_{0} on the CCD can be described as a convolution product,
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 pointsource 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 dispersedimaging 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_{0},y_{c,p}(λ) – y_{0}} on the CCD with respect to the zerothorder position (x_{0}, y_{0}), for a diffraction order p and a wavelength λ. This quantity can be computed by applying the classical grating formula (1). With a pointlike monochromatic source, the image recorded on the CCD is modelled as
with Δ_{0}(λ_{0}) = r_{0} and Α_{p} the flux density at wavelength λ_{0} 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_{*}(λ_{0}) 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
with
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 unstack 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.
Fig. 2 General illustration of the spectrograph model and of the purpose of the spectrum extraction, (a) Illustration of the image formation for a monochromatic star observed through a slitless spectrograph, (b) Same for a polychromatic star, (c) Real acquired image: The detector generally is a blackandwhite sensor, with noise, a structured and chromatic background, and field star zeroth and firstorder spectrograms. The goal of Spectractor is to extract the colour information from this image to obtain the stellar spectrum. 
Fig. 3 Overview of the Spectractor pipeline. The green ellipses represent external inputs needed for the spectrum extraction, and blue ellipses stand for the Spectractor products. The bottom stage represents the full forwardmodelling method for extracting the spectrum from the raw spectrogram. 
2.3 Spectractor
We call Spectractor^{2},^{3} the computer suite we wrote to analyse the future AuxTel images and the images obtained on the Cerro Tololo InterAmerican 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.
Zerothorder centroid. The main inputs of Spectractor are a preprocessed 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 derotated 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^{4} 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^{−1} cm^{−2} nm^{−1} 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 preprocessed^{5} but not rotated exposure using a 2D PSF model, and a model for the atmospheric differential refraction (ADR)^{6}. 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 forwardmodel 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 MayJune 2017 at the Cerro Tololo InterAmerican Observatory (CTIO) 0.9 m Cassegrain telescope (f/13.7, scale at focal plane 60 μm″^{−1}). This telescope is equipped with a cooled Tek2K CCD device of 2048 × 2046 pixels, read by four amplifiers^{7}. Two filter wheels are installed. The first wheel in the light path was used to insert broadband 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^{−1}) ref. GT5003, and an amplitude holographic optical element (around 350 lines mm^{−1}) 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 150200 (Moniez et al. 2021).
Figure 4 shows an example of the data we obtained: The dispersion axes are nearly horizontal along the xaxis 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 pixeltopixel 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 onsite 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^{8} (Bohlin et al. 2014, 2020) star HD111980 with an amplitude hologram. The main characteristics of these data are summarised in Table 1.
Fig. 4 Processed exposure of CALSPEC star HD111980 observed at CTIO with an amplitudehologram grating of 350 lines mm^{−1} on 2017 May 30. The greyscale gives the exposure intensity in ADU s^{−1}. The yellow circle indicates the zerothorder position of HD111980. 
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 realdata image presented in Table 1, but with variations in the unknown parameters such as the PSF model, the amount of secondorder 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,
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 xaxis. 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
Parametrisation with Legendre polynomials has the advantage to give an equal weight to all polynomial coefficients during χ^{2} 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 Forwardmodel extraction of a spectrum
In this section, we follow and describe in detail the steps of the Spectractor pipeline in order to obtain the firstorder spectrum S_{1} (λ) 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.
Fig. 5 Zerothorder images from different celestial objects, observed in different seeing conditions with several dispersers. The red cross shows the fitted centroid using the Spectractor method from Sect. 3.2. All images have saturated pixels. 
3.1 Uncertainty evaluation
When it is not given, we must start to build the uncertainty map of the exposure when the exposure is preprocessed. The uncertainties on the pixel values are estimated using the CCD gain G_{CCD}(x, y) (in electrons ADU^{−1}) and its readout noise σ_{ro} (in electrons). Other sources of noise, such as those coming from the flatfielding 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
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 signaltonoise 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 Zerothorder 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 zerothorder centroid (x_{0}, y_{0}) 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 astrometric 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 zerothorder 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 zerothorder image, as close as needed so that the targeted star is the brightest object. Then the cropped image is projected onto the x and (yaxis: The maximum of the two projections sets a new approximation of the zerothorder position. From there, saturated pixels are detected, and a null weight is associated with them. A 2D secondorder 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 nonsaturated 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^{9} (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 onehalf 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 zerothorder shift δu_{0} used during the wavelengthcalibration process to account for a mistake in the evaluation of the zerothorder centroid.
Fig. 6 Illustration of the fitting process to find the zerothorder centroid for the planetary nebula PNG321.0+3.9. Left: zerothorder data image of the target. Middle: bestfitting 2D circular Moffat model. Right: residuals. 
Fig. 7 Difference between the fitted centroids of the target stars and the centroid recovered using the WCS estimate locked on the Gaia catalogue (blue and orange lines) on the x and yaxis at CTIO during the night of 2017 May 30, as a function of the exposure date. The histogram of the differences is presented on the right. 
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 xaxis, 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 forwardmodel 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 xaxis 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
The two eigenvalues of the Hessian matrix Η are calculated as
with . 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 xaxis can be computed analytically. For instance, we have for λ_
with the trigonometric formula tan 2α = 2 tan α/ (l – tan^{2} α). 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 xaxis. 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.
Fig. 8 Disperser rotation angle estimation. Left: Angle α on the spectrogram image for the 5% pixels with the highest λ_ values found in the Hessian matrix, and the dispersion axis as the median of the angles (black line, shifted upward for clarity). The masked pixel with low λ_ values are indicated (light grey). Right: Histogram of the selected angles from the left panel and its median (vertical black line). 
Fig. 9 Background estimation on a CTIO image (in arbitrary units). Top: laterals bands to the spectrogram (masked behind the rectangular grey region) where the background is estimated. The grey patches indicate masked sources. Middle: fitted background using the SExtractor method with the evaluation boxes in red (here final size is 20 × 20 pixels). Bottom: residuals normalised to their uncertainties. Right: distribution of the normalised residuals in units of σ. 
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 nonflat, 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 (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 . 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^{10} algorithm for background extraction, which is roughly based on the bilinear interpolation of the median value inside the boxes, after a sigmaclipping 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 backgroundsubtracted 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 firstorder spectrum S_{1}(λ) and estimating the wavelengthdependent PSF. The extraction of the spectral information S_{1}(λ) 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 crossdispersion direction, possibly with weights, to form what we call a crossdispersion spectrum. Horne (1986) and Robertson (1986) presented an optimal method to achieve this. However, this method leads to distorted spectra in case of a wavelengthdependent 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 Firstorder 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 xaxis,
with r = (x, y) the vector of the pixel coordinates, A the amplitude parameter vector along the dispersion axis of the spectrogram, and φ(rr_{c,i}, P_{i}) the 2D PSF kernel whose integral is normalised to one. This kernel depends nonlinearly 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 Δ_{1}(λ), and the vector A is the flux density S_{1}(λ) integrated within the pixels. When we index all the N_{x}N_{y} spectrogram pixels as a long vector Z = (ζ_{1},…, ζN_{x}N_{y}), the Eq. (18) takes a matricial form,
with
The (N_{x}N_{y}, N_{x}) matrix M is called the design matrix.
This model of the spectrogram is designed to deconvolve the spectrum S_{1} (λ) 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,
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 GaussNewton minimisation of a χ^{2} (Eq. (23)) for the nonlinear 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
with ϵ the random noise vector. The χ^{2} function to minimise is
with W the weight matrix of dimension (N_{x}N_{y}, N_{x}N_{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
The covariance matrix associated with the Â coefficients is C^{(1D)} = (M^{T} WM)^{−1}. 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 crossspectrum.
The result of this extraction is illustrated on simulated data in Fig. 10 left. For this simulation, we chose to ignore the secondorder diffraction, and we used a wide PSF kernel with strong chromatic variations to enhance the visibility of the residuals for the 1D extraction: (γ_{0},γ_{1},γ_{2}) = (10,2,5) and (α_{0},α_{1},α_{2}) = (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 crossdispersion spectrum as in Horne (1986) and Robertson (1986), but informed with a fitted PSF model with a smooth polynomial wavelength evolution.
Fig. 10 Results from the deconvolution of a simulated spectrogram of the CALSPEC star HD111980 with a Moffat PSF kernel, without secondorder diffraction contamination, (γ_{0},γ_{1},γ_{2}) = (10,2,5) and (α_{0}, α_{1}, α_{2}) = (2, 0, 0). Left: simulated spectrogram (top), best,fitting spectrogram model (middle) and residuals in units of σ (bottom) after the rotation process and the 1D transverse fit. Right: same simulated spectrogram with its original rotation (top), bestfitting spectrogram model (middle), and residuals in units of σ (bottom) after the deconvolution process. All colour maps are normalised by the maximum of the simulated spectrogram. The grey areas designate masked pixels. 
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 χ^{2}). 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 illposed. As we know by the physics that wellsampled 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_{0} 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 χ^{2} and minimise a new cost function,
where Q is the weight matrix. The last term favours Â to be close to prior vector A_{0}, with a positive regularisation hyperparameter r.
Minimising this ℰ(Ar_{c}, P) function is still a linear regression for the A parameters, whose optimal value now is
The covariance matrix associated with Â directly is
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,
with U^{T}U the matrix proposed in Eq. (C.5),
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 norm2 regularisation offers analytical solution, while a norm1 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 GaussNewton minimisation of ℰ(AP) for the nonlinear P parameters (see Appendix B) and a linear resolution for the A parameters (see Sect. 3.5.2). The GaussNewton 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 and 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 hyperparameter equivalent to r still needs to be introduced in Eq. (25). Ryan et al. (2018) tuned this hyperparameter manually and determined the breaking point of a specific Lcurve (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 hyperparameter r.
The regularisation parameter r is optimised via the study of the resolution operator
with I the identity matrix. A fruitful interpretation of the R operator is given in Hansen (2010) with
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
The optimal r parameter is found by minimising the following G(r) function, which behaves like a reduced χ^{2},
This method, known as general crossvalidation (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}^{2}, 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 , parameters, and finally to find the minimum of the G(r) function. We observed in the simulations that regardless of the r hyperparameter that is chosen at the beginning, the process reconstructs an unbiased spectrum at the end of the ℰ(Ar_{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 bestfitting 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 quasiunstructured 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 hyperparameter 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 norm2 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.
Fig. 11 Optimisation of the r hyperparameter for a simple simulation with n_{PSF} = 0, γ_{0} = 5, α_{0} = 3 (and the same characteristics as in Table 1). Top: G(r) function (blue) and the minimum position r = 0.0467 (vertical black line). Middle: χ^{2} evolution with r. Bottom: trace of the resolution matrix Tr R. The intersection with the black line gives the effective number of parameters fitted by data, here, 180. 
3.5.5 The spectrophotometric uncertainty principle
The regularity of the deconvolved solution depends on the hyperparameter 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, α_{0} = 3 (and the same characteristics as in Table 1),
The second panel displays the function, which shows that the optimum A(r) solution does not minimise the agreement with the data (minimum χ^{2}), but makes a compromise with a regularisation scheme (modelled by the ) penalty term). The lower panel illustrates that the effective number of amplitude parameters fitted by data with the optimum régularisation hyperparameter 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,
If too many parameters are searched for, the problem becomes illposed. 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
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 (α_{0} = 2 and 3) and various γ_{0} 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.
Fig. 12 Product of the PSF FWHM with N_{dof}/N_{x} as a function of the PSF FWHM for three different PSF models: a Gaussian kernel (orange stars), a Moffat kernel with α_{0} = 2 (blue points), and a Moffat kernel with α_{0} = 3 (green points). 
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 zerothorder position is at u_{0} along the dispersion axis, but that a misfit of its centroid (see Sect. 3.2) can shift the position by a quantity .
The ADR also slightly spreads the zerothorder 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 wavelengthdependent 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^{11} 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
so that
The bijection between the position on the CCD and the wavelength is thus parametrised with two unknown parameters, D_{CCD} and , that need to be fitted.
As a starting point, we compute a first wavelength array λ_{0} from the array of distances d to the zeroth order along the dispersion axis, assuming δu_{0} = 0 and given a prior value of D_{CCD}. To obtain a wavelength array given D_{CCD} and δu_{0}, Eq. (42) is inverted as
To remove the ambiguity with ADR, which also depends on wavelength, we iterate this computation five times starting from λ_{0} 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 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 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 χ^{2} quantity is computed for each spectroscopic line and added into a global χ^{2}.
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 χ^{2}.
Finally, the full χ^{2} and its penalty are normalised to the number of detected lines. This normalisation of the global χ^{2} is necessary to avoid solutions that favour a lower number of detected lines, while the penalty gives more weights to the welldetected lines and anchors them on their tabulated values. The two parameters and D_{CCD} are varied to minimise the penalised global χ^{2} and find the best solution for the wavelength calibration.
The result of this process is illustrated in Fig. 13. At the top, the global χ^{2} 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 χ^{2} by the number of detected lines. The smoothness around the minimum is due to the penalty term on . 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).
Fig. 13 Wavelengthcalibration process on a planetary nebula spectrum. Top: global function and its minimum (black circle) for the planetary nebula PNG321.0+3.9 observed at CTIO. The sharp steps at the top of the plot arise when certain lines are detected. Bottom: calibrated spectrum of PNG321.0+3.9. The vertical lines indicate emission or absorption lines that are detected, positioned at their tabulated values. The dashed blue lines show the fitted background (whose degree depends on its length), and the plain blue lines show the Gaussian profiles fitted on data. 
Emission lines detected in the spectrum of planetary nebula PNG321.0+3.9 with a S/N above 10.
3.7 Flux calibration
With the fitted wavelength solution, the spectrum amplitude Â can be converted from ADU units into flux densities in erg s^{−1} cm^{−2} nm^{−1}, assuming that the telescope collecting area S_{T}, the exposure time τ, and the CCD gain G_{CCD} (in e^{−} ADU^{−1}) are known,
where δλ is the local variation in λ within one pixel.
The end product of this pipeline is thus a backgroundsubtracted 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 θ_{1}(r, λ) and Δ_{p}(λ), and of geometric parameters such as the zerothorder position r_{0} 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 secondorder to firstorder 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 zerothorder centroid along the yaxis. The ratio r_{2/1} (λ) can be measured on an optical test bench (see Sect. 5.1) or on onsky data. In the full forward model, we use a new design matrix MM that is defined as
where A_{2} is a safety normalisation parameter, R_{2/1} is the transmission ratio vector computed for a given wavelength calibration, are the centroid positions of the second diffraction order PSF kernels, and P^{p=2} are their shape parameters. The vector 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^{12}. Therefore, the full forward model now includes both first and secondorder spectrogram models as
We implemented a twostep 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 GaussNewton descent for the nonlinear parameters P) with the same r hyperparameter 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 forwardmodel fit as a new prior A_{0}. 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 forwardmodel parameters are fitted again together on data, using the more complete model including A, P, D_{CCD}, , and A_{2}. 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 χ^{2}, and consequently, with a better accuracy of all the fitted parameters.
The two main products of this step are a firstorder diffraction spectrum S_{1} (λ) separated from the secondorder spectrum, and the secondorder spectrum S_{2}(λ) with
in erg s^{−1} cm^{−2} nm^{−1} (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 firstorder 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 secondorder 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 fastvarying 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 fastvarying absorption lines. The right panel of the figure also shows the FWHM of the true PSF. The reconstructed PSF displays the same wavelengthdependent PSF FWHM as the simulated one. It is also remarkable that while the crossspectrum 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 forwardmodel 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 endtoend 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 and the extracted spectrum
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_{1}, 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} = 1. We found no significant bias in the A_{1} > 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 γ_{0}, and a small negative bias appeared for lowresolution spectra (large γ_{0}). These spectra exhibit wider and shallower absorption lines than the true spectrum, which leads to these negative b values, but only at the subpercent 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 forwardmodel 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_{0} = Â^{(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.
Fig. 14 Same as Fig. 10, but with the addition of a secondorder diffraction spectrogram and after the full forwardmodel fitting procedure. 
Parameters of a spectrogram simulation and Spectractor estimations, recovered with a full forwardmodel approach.
Fig. 15 Results from the full forwardmodel fit of a simulated spectrogram of the CALSPEC star HD111980 with a Moffat PSF kernel. A second diffraction order was simulated. Top left: spectrum output from the full forwardmodel process (orange) compared with the true spectrum injected in the simulation (blue) and the intermediate 1D transverse fit used as the prior vector A_{0} (green). Bottom left: residuals between the true spectrum and the 1D and 2D fits normalised by their estimated uncertainties. Top right: FWHM of the true PSF (blue, right below the green curve) compared with the fitted PSF during the full forwardmodel procedure (green) or the 1D transverse fit (orange). Bottom right: difference between the true PSF FWHM injected in the spectrogram simulation and the 1D and 2D fits. All blue curves are covered by the orange curves. 
Fig. 16 Correlation matrix of the full forwardmodel fitting of a spectrogram simulation at the end of the GaussNewton descent. 
Fig. 18 Extraction bias b (blue points) as a function of the S/N materialised by the amplitude A_{1} factor (top), of the PSF width γ_{0} (middle), and of the rotation angle α (bottom). For completeness, the effective number of degrees of freedom N_{dof} is represented with red diamonds for each simulation. 
4 Spectrum extraction on data
The success of the spectrum extraction on data mostly depends on the model of the wavelengthdependent PSF of the telescope. If the PSF model correctly represents the reality, the residuals after the full forwardmodel 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^{−1} 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 fourthorder 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 forwardfitting approach. These extractions used a sigmaclipping 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_{2}, 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 onsky 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.
Fig. 19 Same as in Fig. 14, but for CTIO data: CALSPEC star HD111980 observed on 2017 May 30 with a blazed Thorlabs grating 300 lines mm^{−1} (left), for which the PSF is out of focus and deviates from a Moffat model, and an amplitude hologram of 350 lines mm^{−1} (right), better focused on the CCD. 
Fig. 20 Calibrated spectra of CALSPEC star HD111980 observed at CTIO on 2017 May 30 with a blazed Thorlabs grating 300 lines mm^{−1} (left) and an amplitude hologram of 350 lines mn^{−1} (right). The CALSPEC SEDs are given for comparison (scaled for convenience). The two dispersers do not have the same transmission curves, which explains the different shapes of the spectra. The vertical lines indicate emission or absorption lines that are detected, positioned at their tabulated values. Locally, the dashed blue lines show the fitted continuum, and the plain blue lines are the Gaussian profiles fitted on absorption lines. 
Parameters of the CTIO exposures and Spectractor estimates.
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 onsite 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,
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 onsky firstorder transmission and their r_{2/1}(λ) ratio. While we show below how to use onsky 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 onsky 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^{−1} grating was studied 2 yr after the CTIO campaign on the LPNHE optical test bench. The optical bench is made of parabolic offaxis 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 onsky 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 [−(λ − λ_{0})/τ] with three free parameters C, λ_{0}, and τ fitted on the laboratory data.
The firstorder efficiency curve and the r_{2/1}(λ) curve for the blazed Thorlabs 300 lines mm^{−1} 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 airmassdependent 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^{13} (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.
Fig. 21 First diffractionorder transmission of the blazed Thorlabs grating with 300 lines mm^{−1} (blue and black) and the ratio of transmissions of the order 2 over the order 1 (red, blue, and orange). 
Fig. 22 Spectra from the CALSPEC star HD111980 acquired during the night of 2017 May 30 at CTIO, which was assumed to be photometric, with the blazed Thorlabs grating with 300 lines mm^{−1} (top) and an amplitude hologram with 350 lines mm^{−1} (bottom). The curves are coloured according to the acquisition airmass z. In the amplitude hologram, a field star creates a spike around 530 nm. 
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,
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 , with their average constrained to one,
The χ^{2} we need to minimise thus reads
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 and P_{a} parameters were fitted jointly via a GaussNewton 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 SavitzkyGolay 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 firstorder 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 forwardmodel 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.
Fig. 23 Multispectra fit of a CTIO photometric night using the blazed grating with 300 lines mm^{−1}. Top: D_{n} data spectra binned in 3 nm intervals, indexed vertically by their index n and with a coloured amplitude. Masked regions are shown in grey. Middle: bestfitting spectrum models indexed vertically by their index n. Bottom: residual map. 
Fig. 24 Measured T_{inst,1}(λ) curve for the CTIO 0.9 m telescope equipped with a Thorlabs 300 lines mm^{−1} blazed grating (black points) from a photometric night. The red curve is a smoothing using a Savitzky–Golay filter of order 1 and a window size 17. 
Fig. 25 Multispectra fit of a CTIO photometric night using the amplitude hologram with 350 lines mm^{−1}. Top: D_{n} data spectra binned in 3 nm intervals, indexed vertically by their index n and with a coloured amplitude. Masked regions are shown in grey. Middle: bestfitting spectrum models indexed vertically by their index n. Bottom: residual map. 
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 bestquality 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 onsky 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 waterabsorption 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.
Fig. 26 Measured T_{inst,1}(λ) curve for the CTIO 0.9 m telescope equipped with an amplitude holographic grating of about 350 lines mm^{−1} (black points) from a photometric night. The red curve shows a smoothing using a Savitzky–Golay filter of order 1 and a window size 11. 
Atmospheric parameters P_{a} fitted for the same photometric night of 2017 May 30, observing CALSPEC star HD111980 with two different dispersers.
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 χ^{2} 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 broadband 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 millimagnitude 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 millimagnitude precision in r band.
Furthermore, the ozone and VAOD parameters we fitted are similar to the estimates of the global meteorological network MERRA2^{14} (Gelaro et al. 2017) for the CTIO site during that night. The MERRA2 PWV value ranges from 4 to 5 mm during the night of 2017 May 30. As MERRA2 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 onsite measurements. This is even more true for CTIO, which is located at the top of a Chilean mountain. On the other hand, a highatmosphere quantity would be expected to depart less between onsite and satellite measurements. We report the MERRA2 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 onsite 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 MERRA2 and our forwardmodel results supports the algorithm we developed.
For completeness, we also present the evolution of the grey parameters 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 factors is lower than 1%. This supports the firstorder 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.
Fig. 27 Measured grey absorption factors for the CTIO night of 2017 May 30, in which CALSPEC star HD111980 was observed, for each spectrum (ordered in time). 
5.3 Atmospheric forwardmodel approach
After illustrating that the forwardmodel 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_{1} (λ) 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
with A_{1} a grey factor.
As before, the parameters and their covariance matrix were estimated via a GaussNewton descent by minimising a χ^{2} calculated over a single spectrogram. The fitted parameters were A_{1}, A_{2}, δ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_{1}(λ) and the atmospheric transmission on the stellar spectra extracted at this step. We call this a spectrum fit.
Fig. 28 Correlation matrix for the multispectrum parameters fitted for the CTIO night of 2017 May 30, in which CALSPEC star HD111980 was observed. 
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 bestfitting 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_{1} 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_{1} 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 forwardmodel 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 forwardmodel 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_{1}(λ) 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_{1} 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 fieldstar 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 highquality 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 highquality 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.
Fig. 29 Results from the S_{1}(λ) model fit on the spectrum of CALSPEC star HD111980, extracted from a simulated spectrogram. Left: spectrum data (red) compared with the bestfitting model (blue) and the model uncertainties (light blue band) due to the CTIO telescope transmission uncertainties, the true injected spectrum in the simulation (green), and the residuals (bottom). Middle: zoom around the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. Right: correlation matrix of the fitted parameters. 
Fig. 30 Evolution of the atmospheric parameters in the simulation of the night of 2017 May 30, in which CALSPEC star HD111980 was observed with an amplitude hologram. The blue points show the spectrum fits, while the orange points show the spectrogram fits. The green line gives the true values injected in the simulations. 
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 diffractionorder 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 forwardmodel 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 firstorder 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 onsite atmospheric transmission, we pushed our analysis to that point. We showed that our procedure allows us to measure the onsky 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 onsite 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 fieldstar spectrogram and integrating the chromatic flatfielding 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 wavelengthdependent 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 forwardmodel 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 nearinfrared 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 spectrophotometric 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.
Fig. 31 Results from the atmospheric forward model of CALSPEC star HD111980 with a Moffat PSF kernel (the shape parameters evolve as a fourthorder polynomial function), observed with an amplitude hologram with 350 lines mm^{−1}. Left: spectrogram data (top), bestfitting spectrogram model (middle), and residuals in units of σ (bottom) Middle: zoom around the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. All colour maps are normalised by the maximum of the simulated spectrogram. Right: correlation matrix of the fitted parameters. 
Fig. 32 Results from the S_{1}(λ) model fit on the spectrum of CALSPEC star HD111980. Left: spectrum data (red) compared with the bestfitting model (blue, mostly behind the green curve) and the model uncertainties (light blue band) due to the CTIO telescope transmission uncertainties, and the residuals (bottom). Middle: zoom on the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. Right: correlation matrix of the fitted parameters. 
Fig. 33 Evolution of the atmospheric parameters as a function of the exposure time during the night of 2017 May 30, in which CALSPEC star HD111980 was observed with an amplitude hologram. The blue points show the spectrum fits, and the orange points show the spectrogram fits. 
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 (IN2P3CNRS) and the Department of Physics and HarvardSmithsonian Center for Astrophysics, Harvard University. F.B. is part of the FP2M federation (CNRS FR 2036) and of the project Labex MMEDII (ANR11LBX002301). This paper has undergone internal review in the LSST Dark Energy Science Collaboration. The internal reviewers were Marc Betoule, Andres PlazasMalagon 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. DagoretCampagne 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 (CCIN2P3Lyon/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. DEAC0205CH11231; STFC DiRAC HPC Facilities, funded by UK BEIS National Einfrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration. This work was performed in part under DOE Contract DEAC0276SF00515.
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^{15}.
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 solvefield 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 fieldstar 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 GaussNewton minimisation algorithm
The gradient descent to minimise a χ^{2} using the GaussNewton 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 χ^{2} is defined as
where M(P) is the vector of the modelpredicted values for the N data points. The vector R(P) is the residual vector.
In order to find the set of parameters that minimises the χ^{2} function, we searched for the zero of the χ^{2} gradient that verifies . The algorithm we used is the iterative multidimensional GaussNewton method, which we describe hereafter.
We started the minimisation with a firstguess value for the parameters P_{0}. A Taylor expansion at first order of the ∇_{P}χ^{2} function canbe performedaround the starting point P_{0} andgives
with δP_{1} = P_{1} – P_{0} and J_{0} = ∇_{P}M(P_{0}) the Jacobian matrix of the model evaluated at P_{0}. For a linear model, that is, a model that can be written as , the ∇_{P}χ^{2}(P) is exactly equal to its firstorder Taylor expansion.
The zero of the function is then approached by solving the equation ∇_{P}χ^{2}(P_{1}) = 0,
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_{1} that is found exactly cancels the χ^{2} gradient.
We then searched for the α value that minimises the χ^{2} function along the line parametrised by the vector α_{1}δP_{1} where α_{1} is a real number. The P_{1} value solution then reads
The process was iterated K times,
until a convergence criterion is reached. For example when the value of χ^{2}(P_{k}) or P_{k} evolution with k drops below a certain threshold.
The bestfitting model is then considered to be parametrised by the kth vector: . The covariance matrix of the P_{k} parameters is obtained as the Hessian matrix at the minimum χ^{2},
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 Secondorder 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
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 χ^{2}(Ar_{c}, P) is the sum of the squares of the residuals, and thus is the realisation of a random variable following a χ^{2} law with N_{x}N_{y} degrees of freedom, hence the name of the cost, χ^{2}.
The penalisation term is also a quadratic term, and for Q = L^{T}U^{T}UL with the Laplacian operator L = −D^{T}D,
and U is such that
we obtain , and it is the quadratic norm (or the squared Euclidean norm) of the vector UL(A − A_{0}), denoted .
When we interpret A as the discretisation of a continuous spectrogram a by setting , and when σ is a function such that , a continuous analogue of this term would be a term of the form
where a, a_{0} would be functions whose A, A_{0} are discretisations. As a simple consequence,
because .
The total variation distance is defined as the (weighted) norm1 of the gradient operator as a functional term,
We also note that
However, by a simple argument, we can show that the 2norm of the second derivative controls the 1norm 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_{0}(x), we obtain
Thus, under the supplementary constraint that , we have
The discrete analogue asymptotically in N_{x} → +∞ is
This shows that regularisation using the weighted quadratic secondorder 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 secondorder 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 secondorder derivative or research of a sparse minimiser) are not important here, the choice of the weighted quadratic norm of the secondorder 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 wavelengthcalibration 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 wavelengthcalibration process through a wavelengthdependent shift of the zerothorder centroid . This procedure absorbs most of the line shifts when the dispersion axis is not orthogonal to the zenith direction.
Fig. D.1 Absorption line Hβ at 486.3 nm of CALSPEC star HD111980 observed during the night of 2017 May 30 while it goes from an airmass of 1 to an airmass of approximately 2. Left: Without modelling the ADR effect in the wavelengthcalibration process. Right: With the ADR effect in the wavelengthcalibration process. The profiles are better aligned. 
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.
Fig. D.2 Same as Figure D.1, but for the Fe absorption line at 430.8 nm. 
Fig. D.3 Angle conventions for the prediction of the ADR effect in the spectrogram. Top: Celestial coordinates in right ascension α_{d} (RA) and declination δ (DEC) system, with q the parallactic angle that sets the direction of the local zenith positively eastward. Bottom: Angle conventions in the Spectractor image, with ϕ_{c} the northwest axis angle with respect to the horizontal axis of the camera. 
References
 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Marriner, J., Regnault, N., et al. 2013, A&A, 552, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Antier, S., Bertin, E., et al. 2023, A&A, 670, A119 [Google Scholar]
 Birch, K. P., & Downs, M. J. 1993, Metrologia, 30, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Birch, K. P., & Downs, M. J. 1994, Metrologia, 31, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Bohlin, R. C., Gordon, K. D., & Tremblay, P.E. 2014, PASP, 126, 711 [Google Scholar]
 Bohlin, R. C., Hubeny, I., & Rauch, T. 2020, AJ, 160, 21 [Google Scholar]
 Bolton, A. S., & Schlegel, D. J. 2010, PASP, 122, 248 [NASA ADS] [Google Scholar]
 Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, https://doi.org/18.5281/zenodo.4844744 [Google Scholar]
 Burke, D. L., Axelrod, T., Blondin, S., et al. 2010, ApJ, 720, 811 [Google Scholar]
 Burke, D. L., Rykoff, E. S., Allam, S., et al. 2017, AJ, 155, 41 [Google Scholar]
 Coughlin, M., Abbott, T. M. C., Brannon, K., et al. 2016, in Observatory Operations: Strategies, Processes, and Systems VI, eds. A. B. Peck, R. L. Seaman, & C. R. Benn (SPIE) [Google Scholar]
 Edlén, B. 1966, Metrologia, 2, 71 [CrossRef] [Google Scholar]
 Emde, C., BurasSchnell, R., Kylling, A., et al. 2016, Geosci. Model Dev., 9, 1647 [Google Scholar]
 Gelaro, R., McCarty, W., Suárez, M. J., et al. 2017, J. Climate, 30, 5419 [Google Scholar]
 Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215 [Google Scholar]
 Hall, J. T. 1966, Appl. Opt., 5, 1051 [Google Scholar]
 Hansen, P. 1992, SIAM Rev., 34, 561 [Google Scholar]
 Hansen, P. 2010, Discrete Inverse Problems: Insight and Algorithms, Fundamentals of Algorithms (Philadelphia: Society for Industrial and Applied Mathematics) [Google Scholar]
 Hazenberg, F. 2019, Theses, Sorbonne Université, France [Google Scholar]
 Horne, K. 1986, PASP, 98, 609 [Google Scholar]
 Houston, J. M., & Rice, J. P. 2006, Metrologia, 43, S31 [Google Scholar]
 Ingraham, P., Clements, A. W., Ribeiro, T., et al. 2020, SPIE Conf. Ser., 11452, 114520U [Google Scholar]
 Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782 [Google Scholar]
 Li, M., Li, G., Lv, K., et al. 2019, MNRAS, 484, 2403 [Google Scholar]
 Mayer, B., & Kylling, A. 2005, Atmos. Chem. Phys., 5, 1855 [Google Scholar]
 Moniez, M., Neveu, J., DagoretCampagne, S., et al. 2021, MNRAS, 506, 5589 [Google Scholar]
 Murty, M. V. R. K. 1962, J. Opt. Soc. Am., 52, 768 [Google Scholar]
 Neveu, J., Brémaud, V., DagoretCampagne, S., & FisherLevine Merlin. 2021, Astrophysics Source Code Librairy [record ascl:2104.004] [Google Scholar]
 Outini, M., & Copin, Y. 2020, A&A, 633, A43 [Google Scholar]
 Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 [Google Scholar]
 Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [Google Scholar]
 Robertson, J. G. 1986, PASP, 98, 1220 [Google Scholar]
 Rubin, D., Aldering, G., Antilogus, P., et al. 2022, Uniform Recalibration of Common Spectrophotometry Standard Stars onto the CALSPEC System using the SuperNova Integral Field Spectrograph, (UC Berkeley) [Google Scholar]
 Ryan, R. E., Casertano, S., & Pirzkal, N. 2018, PASP, 130, 034501 [Google Scholar]
 Schroeder, D., & Inc, E. I. 2000, Astronomical Optics, Electronics & Electrical (Elsevier Science) [Google Scholar]
 Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Souverin, T., Neveu, J., Betoule, M., et al. 2022, Measurement of telescope transmission using a Collimated Beam Projector, 56th Rencontres de Moriond on Cosmology, Jan 2022, La Thuile, Italy. [Google Scholar]
 The Nearby Supernova Factory (Buton, C., et al.) 2013, A&A, 549, A8 [Google Scholar]
 Wahba, G. 1990, Spline Models for Observational Data, CBMSNSF Regional Conference Series in Applied Mathematics (Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM)) [Google Scholar]
Another choice could have been to assume that the spectrograph does not suffer from defocus, and thus arguing that the PSF shape parameters for the second diffraction order are the same than that of the first diffraction order at the same wavelength whatever the distance to the zeroth order. For CTIO images, our first choice leads to better fits to data.
All Tables
Emission lines detected in the spectrum of planetary nebula PNG321.0+3.9 with a S/N above 10.
Parameters of a spectrogram simulation and Spectractor estimations, recovered with a full forwardmodel approach.
Atmospheric parameters P_{a} fitted for the same photometric night of 2017 May 30, observing CALSPEC star HD111980 with two different dispersers.
All Figures
Fig. 1 Geometry of a simple slitless spectrograph in the plane orthogonal to the CCD (top) and parallel to the CCD (bottom). A convergent beam with incident angle θ_{0} is focused on a CCD at position r_{0} = (x_{0},y_{0}), but also passed through a grating with N_{eff} grooves per millimetre at a distance D_{CCD} above the CCD. The beam is deflected at an angle θ(λ) along the mean dispersion axis u, which forms an angle α with respect to the CCD x axis, but is focused somewhere above the sensor. Atmospheric refraction adds a supplementary dispersion along and transversally to the mean dispersion axis. 

In the text 
Fig. 2 General illustration of the spectrograph model and of the purpose of the spectrum extraction, (a) Illustration of the image formation for a monochromatic star observed through a slitless spectrograph, (b) Same for a polychromatic star, (c) Real acquired image: The detector generally is a blackandwhite sensor, with noise, a structured and chromatic background, and field star zeroth and firstorder spectrograms. The goal of Spectractor is to extract the colour information from this image to obtain the stellar spectrum. 

In the text 
Fig. 3 Overview of the Spectractor pipeline. The green ellipses represent external inputs needed for the spectrum extraction, and blue ellipses stand for the Spectractor products. The bottom stage represents the full forwardmodelling method for extracting the spectrum from the raw spectrogram. 

In the text 
Fig. 4 Processed exposure of CALSPEC star HD111980 observed at CTIO with an amplitudehologram grating of 350 lines mm^{−1} on 2017 May 30. The greyscale gives the exposure intensity in ADU s^{−1}. The yellow circle indicates the zerothorder position of HD111980. 

In the text 
Fig. 5 Zerothorder images from different celestial objects, observed in different seeing conditions with several dispersers. The red cross shows the fitted centroid using the Spectractor method from Sect. 3.2. All images have saturated pixels. 

In the text 
Fig. 6 Illustration of the fitting process to find the zerothorder centroid for the planetary nebula PNG321.0+3.9. Left: zerothorder data image of the target. Middle: bestfitting 2D circular Moffat model. Right: residuals. 

In the text 
Fig. 7 Difference between the fitted centroids of the target stars and the centroid recovered using the WCS estimate locked on the Gaia catalogue (blue and orange lines) on the x and yaxis at CTIO during the night of 2017 May 30, as a function of the exposure date. The histogram of the differences is presented on the right. 

In the text 
Fig. 8 Disperser rotation angle estimation. Left: Angle α on the spectrogram image for the 5% pixels with the highest λ_ values found in the Hessian matrix, and the dispersion axis as the median of the angles (black line, shifted upward for clarity). The masked pixel with low λ_ values are indicated (light grey). Right: Histogram of the selected angles from the left panel and its median (vertical black line). 

In the text 
Fig. 9 Background estimation on a CTIO image (in arbitrary units). Top: laterals bands to the spectrogram (masked behind the rectangular grey region) where the background is estimated. The grey patches indicate masked sources. Middle: fitted background using the SExtractor method with the evaluation boxes in red (here final size is 20 × 20 pixels). Bottom: residuals normalised to their uncertainties. Right: distribution of the normalised residuals in units of σ. 

In the text 
Fig. 10 Results from the deconvolution of a simulated spectrogram of the CALSPEC star HD111980 with a Moffat PSF kernel, without secondorder diffraction contamination, (γ_{0},γ_{1},γ_{2}) = (10,2,5) and (α_{0}, α_{1}, α_{2}) = (2, 0, 0). Left: simulated spectrogram (top), best,fitting spectrogram model (middle) and residuals in units of σ (bottom) after the rotation process and the 1D transverse fit. Right: same simulated spectrogram with its original rotation (top), bestfitting spectrogram model (middle), and residuals in units of σ (bottom) after the deconvolution process. All colour maps are normalised by the maximum of the simulated spectrogram. The grey areas designate masked pixels. 

In the text 
Fig. 11 Optimisation of the r hyperparameter for a simple simulation with n_{PSF} = 0, γ_{0} = 5, α_{0} = 3 (and the same characteristics as in Table 1). Top: G(r) function (blue) and the minimum position r = 0.0467 (vertical black line). Middle: χ^{2} evolution with r. Bottom: trace of the resolution matrix Tr R. The intersection with the black line gives the effective number of parameters fitted by data, here, 180. 

In the text 
Fig. 12 Product of the PSF FWHM with N_{dof}/N_{x} as a function of the PSF FWHM for three different PSF models: a Gaussian kernel (orange stars), a Moffat kernel with α_{0} = 2 (blue points), and a Moffat kernel with α_{0} = 3 (green points). 

In the text 
Fig. 13 Wavelengthcalibration process on a planetary nebula spectrum. Top: global function and its minimum (black circle) for the planetary nebula PNG321.0+3.9 observed at CTIO. The sharp steps at the top of the plot arise when certain lines are detected. Bottom: calibrated spectrum of PNG321.0+3.9. The vertical lines indicate emission or absorption lines that are detected, positioned at their tabulated values. The dashed blue lines show the fitted background (whose degree depends on its length), and the plain blue lines show the Gaussian profiles fitted on data. 

In the text 
Fig. 14 Same as Fig. 10, but with the addition of a secondorder diffraction spectrogram and after the full forwardmodel fitting procedure. 

In the text 
Fig. 15 Results from the full forwardmodel fit of a simulated spectrogram of the CALSPEC star HD111980 with a Moffat PSF kernel. A second diffraction order was simulated. Top left: spectrum output from the full forwardmodel process (orange) compared with the true spectrum injected in the simulation (blue) and the intermediate 1D transverse fit used as the prior vector A_{0} (green). Bottom left: residuals between the true spectrum and the 1D and 2D fits normalised by their estimated uncertainties. Top right: FWHM of the true PSF (blue, right below the green curve) compared with the fitted PSF during the full forwardmodel procedure (green) or the 1D transverse fit (orange). Bottom right: difference between the true PSF FWHM injected in the spectrogram simulation and the 1D and 2D fits. All blue curves are covered by the orange curves. 

In the text 
Fig. 16 Correlation matrix of the full forwardmodel fitting of a spectrogram simulation at the end of the GaussNewton descent. 

In the text 
Fig. 17 Same as Fig. 11, but for the simulation illustrated in Fig. 15. 

In the text 
Fig. 18 Extraction bias b (blue points) as a function of the S/N materialised by the amplitude A_{1} factor (top), of the PSF width γ_{0} (middle), and of the rotation angle α (bottom). For completeness, the effective number of degrees of freedom N_{dof} is represented with red diamonds for each simulation. 

In the text 
Fig. 19 Same as in Fig. 14, but for CTIO data: CALSPEC star HD111980 observed on 2017 May 30 with a blazed Thorlabs grating 300 lines mm^{−1} (left), for which the PSF is out of focus and deviates from a Moffat model, and an amplitude hologram of 350 lines mm^{−1} (right), better focused on the CCD. 

In the text 
Fig. 20 Calibrated spectra of CALSPEC star HD111980 observed at CTIO on 2017 May 30 with a blazed Thorlabs grating 300 lines mm^{−1} (left) and an amplitude hologram of 350 lines mn^{−1} (right). The CALSPEC SEDs are given for comparison (scaled for convenience). The two dispersers do not have the same transmission curves, which explains the different shapes of the spectra. The vertical lines indicate emission or absorption lines that are detected, positioned at their tabulated values. Locally, the dashed blue lines show the fitted continuum, and the plain blue lines are the Gaussian profiles fitted on absorption lines. 

In the text 
Fig. 21 First diffractionorder transmission of the blazed Thorlabs grating with 300 lines mm^{−1} (blue and black) and the ratio of transmissions of the order 2 over the order 1 (red, blue, and orange). 

In the text 
Fig. 22 Spectra from the CALSPEC star HD111980 acquired during the night of 2017 May 30 at CTIO, which was assumed to be photometric, with the blazed Thorlabs grating with 300 lines mm^{−1} (top) and an amplitude hologram with 350 lines mm^{−1} (bottom). The curves are coloured according to the acquisition airmass z. In the amplitude hologram, a field star creates a spike around 530 nm. 

In the text 
Fig. 23 Multispectra fit of a CTIO photometric night using the blazed grating with 300 lines mm^{−1}. Top: D_{n} data spectra binned in 3 nm intervals, indexed vertically by their index n and with a coloured amplitude. Masked regions are shown in grey. Middle: bestfitting spectrum models indexed vertically by their index n. Bottom: residual map. 

In the text 
Fig. 24 Measured T_{inst,1}(λ) curve for the CTIO 0.9 m telescope equipped with a Thorlabs 300 lines mm^{−1} blazed grating (black points) from a photometric night. The red curve is a smoothing using a Savitzky–Golay filter of order 1 and a window size 17. 

In the text 
Fig. 25 Multispectra fit of a CTIO photometric night using the amplitude hologram with 350 lines mm^{−1}. Top: D_{n} data spectra binned in 3 nm intervals, indexed vertically by their index n and with a coloured amplitude. Masked regions are shown in grey. Middle: bestfitting spectrum models indexed vertically by their index n. Bottom: residual map. 

In the text 
Fig. 26 Measured T_{inst,1}(λ) curve for the CTIO 0.9 m telescope equipped with an amplitude holographic grating of about 350 lines mm^{−1} (black points) from a photometric night. The red curve shows a smoothing using a Savitzky–Golay filter of order 1 and a window size 11. 

In the text 
Fig. 27 Measured grey absorption factors for the CTIO night of 2017 May 30, in which CALSPEC star HD111980 was observed, for each spectrum (ordered in time). 

In the text 
Fig. 28 Correlation matrix for the multispectrum parameters fitted for the CTIO night of 2017 May 30, in which CALSPEC star HD111980 was observed. 

In the text 
Fig. 29 Results from the S_{1}(λ) model fit on the spectrum of CALSPEC star HD111980, extracted from a simulated spectrogram. Left: spectrum data (red) compared with the bestfitting model (blue) and the model uncertainties (light blue band) due to the CTIO telescope transmission uncertainties, the true injected spectrum in the simulation (green), and the residuals (bottom). Middle: zoom around the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. Right: correlation matrix of the fitted parameters. 

In the text 
Fig. 30 Evolution of the atmospheric parameters in the simulation of the night of 2017 May 30, in which CALSPEC star HD111980 was observed with an amplitude hologram. The blue points show the spectrum fits, while the orange points show the spectrogram fits. The green line gives the true values injected in the simulations. 

In the text 
Fig. 31 Results from the atmospheric forward model of CALSPEC star HD111980 with a Moffat PSF kernel (the shape parameters evolve as a fourthorder polynomial function), observed with an amplitude hologram with 350 lines mm^{−1}. Left: spectrogram data (top), bestfitting spectrogram model (middle), and residuals in units of σ (bottom) Middle: zoom around the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. All colour maps are normalised by the maximum of the simulated spectrogram. Right: correlation matrix of the fitted parameters. 

In the text 
Fig. 32 Results from the S_{1}(λ) model fit on the spectrum of CALSPEC star HD111980. Left: spectrum data (red) compared with the bestfitting model (blue, mostly behind the green curve) and the model uncertainties (light blue band) due to the CTIO telescope transmission uncertainties, and the residuals (bottom). Middle: zoom on the dioxygen line at 762 nm and zoom on the H_{2}O absorption band around 950 nm. Right: correlation matrix of the fitted parameters. 

In the text 
Fig. 33 Evolution of the atmospheric parameters as a function of the exposure time during the night of 2017 May 30, in which CALSPEC star HD111980 was observed with an amplitude hologram. The blue points show the spectrum fits, and the orange points show the spectrogram fits. 

In the text 
Fig. D.1 Absorption line Hβ at 486.3 nm of CALSPEC star HD111980 observed during the night of 2017 May 30 while it goes from an airmass of 1 to an airmass of approximately 2. Left: Without modelling the ADR effect in the wavelengthcalibration process. Right: With the ADR effect in the wavelengthcalibration process. The profiles are better aligned. 

In the text 
Fig. D.2 Same as Figure D.1, but for the Fe absorption line at 430.8 nm. 

In the text 
Fig. D.3 Angle conventions for the prediction of the ADR effect in the spectrogram. Top: Celestial coordinates in right ascension α_{d} (RA) and declination δ (DEC) system, with q the parallactic angle that sets the direction of the local zenith positively eastward. Bottom: Angle conventions in the Spectractor image, with ϕ_{c} the northwest axis angle with respect to the horizontal axis of the camera. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.