© The Authors 2025

1 Introduction

A major challenge in astronomical observations is to accurately transform the instrumental photometry into a calibrated flux in a well-defined passband. Such calibration is typically required for either measuring the flux of the source (in physical units), or to make a comparison between the flux measured by different telescopes. Ground-based surveys regularly provide an absolute calibration of their photometry to the 1-2% level (e.g., Padmanabhan et al. 2008, Ofek et al. 2012a,b, Masci et al. 2019, Magnier et al. 2020). However, breaking this ∼1% barrier is difficult. Furthermore, the photometry based on standard methods is not provided in a well defined band. One limiting factor is the intrinsic instrumental difference between astronomical facilities. Even when the same filter is used by different telescopes, they do not necessarily share the same optics, coatings, and detectors, which may modify slightly the effective transmission. Furthermore, when comparing the flux taken with the same filter using two telescopes, their effective transmission may be different due to e.g., local atmosphere. Furthermore, the variable atmosphere constantly changes the effective transmission, resulting in flux variations even if the source is constant. All these complications impose the historical limits on the possibility of achieving absolute calibrations with an accuracy of less than a few percent (see Stubbs & Tonry 2006 for a review).

The classical photometric calibration method was to observe a set of standard stars at different air mass, and to use this to measure the system photometric zero point as a function of air mass (e.g., using Landolt stars; Landolt 1992). The main limitation of this approach is the simplistic assumption that the zero point depends only on air mass (i.e., that atmospheric transmission is otherwise uniform). A better approach would be to have a set of standard stars with a density high enough such that some will be present in any random observation (e.g., Skrutskie et al. 2006; Ofek 2008; Ofek et al. 2012b; Magnier et al. 2020).

A powerful approach for photometric calibration adopted by several sky surveys is the übercalibration method (Padmanabhan et al. 2008). In this approach the photometry is first calibrated internally by using stars that appear in two, somewhat overlapping, images defining a self-consistent system for instrumental magnitudes. After that, an external calibration is performed using a set of photometric standard stars that serve as calibrators with accurate measurements of their spectral energy distribution (SED). A recent application of this calibration approach is the one used for the Gaia mission (Gaia Collaboration 2016), for which the photometric G band has been defined (Carrasco et al. 2016). Indeed in the case of Gaia , its photometry, as well as spectrophotometric measurements are consistently calibrated to the level of a few mmag (Carrasco et al. 2016; De Angeli et al. 2023). However, for ground-based observations, the temporal variations in the effective passbands are limiting the accuracy of this method.

A conceptual solution to this problem is to somehow measure the atmospheric and system transmission for each image obtained by the telescope. Previous works have discussed the importance of having a precise description of the time-dependent optical transmission of the atmosphere by e.g., realtime measures on site (e.g., Burke et al. 2010, Burke et al. 2018, also Stubbs et al. 2007 for a comprehensive overview). Burke et al. (2010) shows how auxiliary spectroscopic measurements of a set of calibrator stars during the night allow to characterize the atmospheric absorption of different atmospheric constituents and their variation. However, auxiliary spectroscopic measurements on site require dedicated instruments to monitor the fields observed by the main survey (e.g., the Auxiliary calibration telescope of the Vera C. Rubin Observatory, LSST Science Collaboration 2009).

The method presented in this work consists of fitting the system transmission for each image by combining parametric models for the instrumental transmission and for the atmospheric transmission, fitted using a large number of stellar calibrators across the field of view. The method aims to obtain high accuracy and stability of the photometry over time, and to break the 1% barrier of the current state-of-the-art calibrations for groundbased facilities. Critical to this work is the existence of accurately calibrated low-resolution spectra of a large number of stars over the entire sky. Therefore, this work became feasible only after the appearance of the Gaia catalog of ∼220 million low-resolution spectra (De Angeli et al. 2023), that provide a sufficient number of diverse calibrators in wide-field images. A similar approach was suggested by Montegriffo et al. 2023, where the instrument model for the Gaia BP and RP spectrophotometers was derived by fitting correction terms to the nominal, pre-launch passband of the instrument. In our approach, we also make use of synthetic photometry to match the predicted and observed counts in our detector for each source, following the technique outlined in Bessell & Murphy 2012. Following the description of the methodology, we explore the effects of atmospheric transmission variations on the photometry. We demonstrate the suggested method using data obtained by the Large Array Survey Telescope (LAST; Ofek et al. 2023a, Ben-Ami et al. 2023) and show that it can provide stable and consistent absolute photometric calibration to a level of 3-5 mmag. Throughout this work, the term absolute calibration refers to the process of determining a solution for the system’s effective transmission. Since the Gaia spectra used here are calibrated against the CALSPEC scale (Altavilla et al. 2021; Montegriffo et al. 2023), the resulting calibration is inherently tied to the CALSPEC standards (Bohlin et al. 2014).

In Sect. 2, we present our proposed solution to this problem. In Sect. 3, we dissect the transmission model into all of its components. In Sect. 4, we describe the approach to the model optimization based on synthetic photometry, and in Sect. 5 we evaluate the sensitivity of our model to variations of the atmospheric parameters. In Sect. 6 we apply the method on real data from LAST, while in Sect. 7 we discuss the results. In Sect. 8 we describe the code that we used to implement the method and in Sect. 9 we summarize the work and discuss future efforts. In Appendix E we suggest a standard for FITS headers to store the photometric transmission parameters.

2 Description of the methodology

The problem of deriving the throughput as a function of the wavelength λ, given a set of calibrators with a known spectrum, can be generally formulated as follows: $$F_{inst,j}=\sum _i^N{\alpha }_i{T}_i(\lambda )S_j(\lambda ).$$(1)

Here F_{inst, j} is the observed instrumental flux of the j-th calibrator, S_j(λ) is its known spectral flux at the wavelength λ, T_i(λ) is a set of N positive-definite function-basis that span the spectral space. Finally, αi are the free parameters we would like to find, where the system transmission is given by $$T(\lambda )=\sum _i^N{\alpha }_i{T}_i(\lambda ).$$(2)

The flux and spectra in this problem can be generally written in physical (calibrated) units. In this work, we convert them to instrumental photon counts only when required (e.g., Sect. 4). A necessary condition for solving Equation (1) is that the spectra of the calibrators S_j(λ) will be diverse enough such that finding α_i will be possible.

A straightforward approach to solve for α_i is to assume the T_i(λ) are top-hat functions, and in this case, the problem can be solved using linear least square or positive linear least square methods. To validate the feasibility of this approach we have performed the following simulations. We selected 1000 random Gaia spectra, generated a transmission curve, and calculated the synthetic flux of the Gaia spectra using this transmission. Next, we added noise to the synthetic fluxes (to mimic F_{inst, j}) and attempted to solve Equation (1). Assuming T_i is a set of tophat functions we were able to recover α_i and the input T(λ) only when the photometric relative errors were unrealistically small (≲10⁻³). This indicates that the problem requires some regularization. A reasonable regularization is to assume that the transmission function we are trying to recover is smooth, or composed of physical components. For example, when we redid the simulations using T_i which are polynomials or other functions (see below) that problem became solvable even in the presence of realistic relative noise (e.g., ∼0.1).

There are many options for choosing a base of T_i, which are likely to work. Here, we decided to use a set of functions which are physically motivated, but to add freedom via low-order perturbing functions. In addition, atmospheric absorption can be attributed to several known mechanisms with some free parameters. In the next section we describe the basis T_i used in this work. The proposed basis contains a large number of free parameters. However, by investigating the effect of each parameter we select a subset of parameters that have the largest effect on the results and describe the transmission to acceptable accuracy.

3 The transmission model basis

The model for the overall transmission includes a series of different components such as the instrumental transmittance of the telescope optical elements (T_tel), the quantum efficiency (QE) of the camera (T_det), the atmospheric transmission (T_atm), and a component that introduces additional perturbations in the model (T_pert). The model can be described as follows: $$T(\lambda )=T_{tel}(\lambda )T_{det}(\lambda )T_{atm}(\lambda )T_{pert}(\lambda ).$$(3)

Each of these components, is composed of several independent sub-components that we present in details below. Some of these components may be known and other may have some free parameters.

3.1 Telescope and detector

The first two terms of Equation (3), T_tel(λ) and T_det(λ), describe the known transmission of the different instrumental components of our system. The first one describes the components of the telescope itself, and can be in turn described as the factorization of the transmissions of the different components (TC_i) within the optical path in the telescope: $$T_{tel}(\lambda )=\prod _i^N{T}_{C_i}.$$(4)

As for T_det (λ) is usually the product of the quantum efficiency (T_QE (λ)) of the camera installed on the telescope, as typically obtained by laboratory measurements and the transmission of the filter T_F (λ) used for observations $$T_{det}(\lambda )=T_{QE}(\lambda )T_F(\lambda ).$$(5)

3.2 Perturbations

In order to introduce some freedom in the shape of the instrumental transmissions measured in the laboratory, one can consider adding perturbations to the model with a number of free parameters. This is very important since even small changes in the transmission can lead to measurable differences in the flux estimation of sources (see Appendix A). One possible way to add such perturbations, is to use a base of low-order (N) polynomial functions described by the general form $$T_{pert}(\lambda )=\sum _{n=0}^N{c}_n{P}_n(\lambda ),$$(6)

where c_n are the fitted coefficients and P_n is the n-th order polynomial. A similar perturbation approach is applied in Montegriffo et al. 2023, where a basis of Legendre polynomials is used and the perturbation is the exponential function of the linear combination of the basis, to guarantee the non-negativity of the perturbations.

Another approach, which we use in this work, is to replace the measurement of T_QE(λ) of the camera with an analytic function with a similar shape, but with a number of free parameters that allow adjusting its shape. For example, in this work we replace T_det with a skewed Gaussian function multiplied by a polynomial basis $$T_{det}=f(\lambda ;A,\mu ,\sigma ,\gamma )exp(\sum _{n=0}^N{c}_n{P}_n(\lambda )),$$(7)

where the polynomial basis is treated as argument of the exponential function in order to have non-negative values for the transmission model. The skewed Gaussian function f (λ; A, μ, σ, γ) is expressed as $$f(\lambda ;A,\mu ,\sigma ,\gamma )=\frac{A}{\sigma \sqrt{2\pi }}{e}^{[-(\lambda -\mu {)}^2/2{\sigma }^2]}\{1+erf[\frac{\gamma (\lambda -\mu )}{\sigma \sqrt{2}}]\},$$(8)

where A is the amplitude, μ is the center, σ is the equivalent width and γ is the skeweness. Such modeling allows to introduce modifications in the shape of the total effective transmission introduced by the optics and other external factors. This method is useful, for example, when the actual QE curve of the detector is similar to a skewed Gaussian function. In Appendix B, we explore the effects of this choice on the results.

3.3 Atmospheric transmission

The contribution of the atmosphere to the overall transmission is, in most cases, the major factor that limits the accuracy of flux measurements to the 1-2% level. The atmospheric effects are of two different kinds: scattering and absorption. The former includes Rayleigh scattering of light by molecules of gases in the atmosphere, and Mie scattering by aerosol particles. All remaining processes are due to absorption of light by atmospheric constituents (gases of different species). In this work, we use analytical and empirical parameterizations for the atmospheric components and the relative air mass of the constituents from the SMARTS v2.9.8 model (see Gueymard 2019 and references therein).

The total transmission due to light traveling throughout the atmosphere can be computed as the product of the transmissions from each of the N (=5) processes considered: $$T_{atm}(\lambda )=\prod _{i=1}^{i=N}{T}_{atm,i}(\lambda ).$$(9)

Each T_atm,i is usually described by the analytical law $$T_{atm,i}(\lambda )=e^{-m_i{\tau }_i(\lambda ;\hat{\gamma })},$$(10)

where m_i is the constituent’s air mass over the column at a fixed zenith angle θ_z given by $$m_X=\frac{1}{\cos ({\theta }_z)+{\eta }_0\ast {\theta }_z^{{\eta }_1}({\eta }_2-{\theta }_z{)}^{{\eta }_3}}.$$(11)

Here X defines the component, and η₁, η₂ and η₃ are parameter-izations in Gueymard 2019 to estimate the effective airmass of each atmospheric component. We list in Appendix D the parameters from SMART v2.9.8. The term τi(λ; γ̂) is the constituent optical depth at a specific wavelength λ per unit air mass, and γ̂ is a set of additional parameters. In our transmission model, we consider a total of five atmospheric components: Rayleigh scattering, aerosol scattering, ozone absorption, water vapour absorption, and absorption from a uniform mixture of gases. The respective contribution of these components is discussed below and it is shown in Figure 1.

Rayleigh scattering: describes the scattering of light by gas molecules. The optical depth of this process can be described by a fourth degree polynomial of the wavelength λ and is parameterized by average empirical results in the literature (SMART v2.9.8 uses results from Bodhaine et al. 1999, Bucholtz 1995 and Tomasi et al. 2005) $${\tau }_{R,\lambda }=\frac{p}{p_0}(117.3405{\lambda }^4-1.5107{\lambda }^2+0.017535-0.00087743{\lambda }^{-2}{)}^{-1},$$(12)

where p₀ = 1013.25 mbar is the standard pressure, p the groundlevel pressure (in mbar), and λ is the wavelength measured in μm. The effective airmass of the component determining the Rayleigh scattering is calculated from Equation (11) using the values of the η parameters in Table D.1. The only free parameter for this component is the surface pressure p. The top-left panel of Figure 1 shows the transmission profile from Rayleigh scattering at different air mass, using p = 965 mbar¹.

Aerosol scattering: describes the interaction with aerosol mixtures in the atmosphere. The analytical formulation for the optical depth of this component (also known as Ångstrom’s law) is $${\tau }_{a,\lambda }={\tau }_{a500}(2\lambda {)}^{-\alpha },$$(13)

where λ is measured in μm, τ_a500 is the aerosol optical depth (AOD) at 500 nm and α is the Ångstrom exponent that depends on the average grain size of the aerosol mixture. Therefore, this component has two free parameters. In the typical conditions of rural aerosol model τ_a500 ≈ 0.084, while in desert-dust aerosol conditions it gets up to τa₅₀₀ ≈ 0.54. In desert-like conditions, α is constant through the whole wavelength range of interest in our work, and has typical value ≈0.6 in moderate presence of dust (Gueymard 2019). The upper-central panel of Figure 1 shows the transmission profile from aerosol scattering at different air mass, using τ_a500 = 0.084 and α = 0.6. The effective airmass of the component determining the aerosol scattering is calculated from Equation (11) using the values of the η parameters in Table D.1.

Ozone absorption: is the absorption due to the total O₃ column in the atmosphere. Absorption from O₃ suppresses severely the transmission in the near-UV band below about 360 nm, and its optical depth is described by the Beer-Lambert-Bouguer law $${\tau }_{o,\lambda }=u_o{A}_{o,\lambda },$$(14)

where u₀ is the length of the total O₃ column expressed in Dobson Units (1 DU = 10 μm-atm) and A_oλ is the known spectral absorption of O₃. Tabulated values for A_oλ are taken from SMARTS v2.9.8 (Gueymard 2019). The parameter u₀ is the only parameter that can be assumed free in this component. The upper-right panel of Figure 1 shows the transmission profile from O₃ absorption at different air mass, assuming u₀ = 300 DU. Other than the severe absorption below ∼360nm mentioned above, the transmission profile shows another absorption feature in the Chappuis band (500 nm < λ < 700 nm) that accounts for up to ∼8% attenuation at 600 nm for moderate air mass. The effective airmass of the O₃ column is calculated from Equation (11) using the values of the η parameters in Table D.1.

Water vapor absorption: is the absorption due to the water vapor in the atmosphere. Unlike the previous components, its transmission has a different functional dependence with the air mass, and it can be described as (Gueymard 2001): $$T_{v,\lambda }=e^{-[(m_w{p}_w{)}^{1.05}{f}_v^{n(\lambda )}{B}_v(p)A_{v,\lambda }{]}^{c(\lambda )}}.$$(15)

Here m_w is the air mass for the water vapor (calculated with Equation (11), using the values of η for H₂O in Table D.1), and p_w the total amount of precipitable water (in cm). The scaling factor f^n(λ) compensates pressure inhomogeneities along the water vapor column following the Curtis-Gordon approximation (Pierluissi et al. 1989) and B_v(p) is a correction term dependent on the atmospheric pressure, that accounts for the different absorption intensity with respect to the distance from the band center. Lastly, c and n are wavelength-dependent exponents, which are tabulated in Gueymard 2019. Apart from using the parameterizations given in Gueymard (2001) and the effective water vapor air mass of the observation, the only free parameters in this component are the amount of precipitable water p_w and the atmospheric pressure p. The value of p_w depends strongly on the site location and conditions, with typical values between 1-2 cm.

Uniform gas mixture absorption: is the absorption due to a uniform mixture of gases following the Standard US atmosphere model (National Geophysical Data Center 1992), to include remaining contributions not considered for the other components. Each gaseous component contributes to the total transmittance with the functional form given by Gueymard (2001): $$T_{g,\lambda }=e^{-(m_g{u}_g{A}_{g,\lambda }{)}^a}.$$(16)

Here m_g is the gas mixture air mass, calculated with Equation (11) and the values of the parameters η from Table D.1 of all the gaseous components from O₂ to NH₃ weighted by their respective abundance in the Standard US atmosphere model. The term u_g is the altitude-dependent gaseous pathlength and A_g,λ is the spectral absorption coefficient. The tabuled values of u_g and A_g,λ, as well as the values of the exponent a for different wavelength ranges, are the same as implemented in the SMARTS v2.9.8 model (Gueymard 2019). The free parameters of this component are the surface pressure and temperature.

The transmittance of each atmospheric component included in our model is shown in the subpanels of Figure 1. Each transmittance is presented at zenith angles (θ_z) of 0°, 15°, 30°, 45° and 60°. The bottom-right panel of Figure 1 shows the resulting transmittance after multiplying all the single components.

Fig. 1 Transmittance of each atmospheric component at different zenith angles. The panels in the top row show from left to right the transmittance for Rayleigh scattering, aerosol scattering and ozone absorption. In the bottom panels, left to right panels are the transmittance of the absorption by water vapor, absorption by the uniform mixture of gases, and the product of all atmospheric components.

4 Fitting procedure

We outline here the method for the optimization of the overall transmission model described in the previous sections. This is obtained by comparing the synthetic photometry with the modeled transmission of sources with calibrated spectra with the instrumental photometry (photon counts) measured from the images. An overview of synthetic photometry can be found in Bessell & Murphy 2012, while here we describe the elements we adopt in our procedure.

Fluxes measured in instrumental photometry are usually reported in total number of photons within a fixed photometric aperture or point-spread function (PSF) fitting. We define here as N_e the rate of measured electrons (in e⁻/s) obtained by dividing the integral photon flux by the image exposure time. For a source with spectral flux F(λ) and a system with transmission T(λ), the electron rate is given by $$N_{e,syn}=\frac{A_g}{hc}{\int }_{\lambda }F(\lambda )\lambda T(\lambda )d\lambda ,$$(17)

where A_g is the geometrical collecting area of the telescope, h is the Planck constant and c the speed of light. In the integral term of Equation (17), the only term containing free parameters is T(λ), defined in the most general form in Equation (3). The optimization to find the set of best-fit parameters ($\hat{\varrho }$) for the transmission T($\lambda ;\hat{\varrho }$) consists then in minimizing the residuals between the measured instrumental electron rate and the predicted electron rate from synthetic photometry: $$R=N_{e,obs}-N_{e,syn}.$$(18)

So far we assumed that the transmission is constant across the field of view. This is generally not accurate, as differences can be caused by factors like non-uniform illumination or response of the sensor, or residuals from dark and flat field corrections (e.g., Padmanabhan et al. 2008; Ofek et al. 2012a). Another kind of non-uniformity across the field of view can be introduced by the presence of cloud structures. This is known as gray extinction (as it is wavelength-independent) and can introduce significant spatial structures in the zero-point distribution across the field of view (e.g., Burke et al. 2014). We therefore add a polynomial term to Equation (18) to correct for these non-uniformities $$P(x,y)=\sum _{i=1}^{n_x}\sum _{j=1}^{n_y}{k}_{i,j}{P}_i(x)P_j(y),$$(19)

where P_i(x) and P_j(y) are a basis of low-order (n) polynomials, x and y the pixel coordinates of the sensor, and k_i,j the fitted coefficients. It is convenient to minimize these residuals in logarithmic space, therefore the field corrections P(x,y) can be expressed in magnitudes and the formulation of the residuals to be minimized becomes $$\overline{R}=2.5lo{g}_{10}(\frac{A_{g}C}{hc{N}_{e,obs}}{\int }_{\lambda }F(\lambda )\lambda T(\lambda ;\hat{\varrho })d\lambda )+P(x,y;\hat{k}),$$(20)

where $\hat{\varrho }$ is the set of free parameters in the transmission model, k̂ is a set of free parameters in the spatial polynomial corrections (Equation (19)) and C is a normalization factor that corrects for the effective collecting area of the system. R is now in magnitude space.

5 Sensitivity to atmospheric parameters

The atmospheric component of our transmission model has a total of six free parameters (p,τ_a500, α, u₀, p_w and T) that were described in Sect. 3.3. Since some of these parameters are degenerate, and most of them are not expected to vary significantly outside of a physically-motivated range, it is useful to investigate the sensitivity of the photometry to the variation of the free parameters. Hence, we assume a mock broadband transmission model and we calculate the variation of a putative zero-point of the system (hence the integrated transmission) for variations of these parameters within typical ranges. In the panels of Figure 2 we show the variation (in magnitudes) of the zero point for selected bands by varying each time a single atmospheric parameter while keeping the rest of the parameters at their nominal values. We make this exercise on the LAST passband (introduced in Sect. 6) and on the u,g,r,i,z SDSS filters (Doi et al. 2010). For most of the parameters, there is a linear dependency to the induced variation in the zero point. Two exceptions are the parameters α and p_w , which introduce non-linear variations on large ranges.

For the LAST band, the parameter α from the aerosol component introduces variations <1% even at the edges of the range considered, and p_w is supposed to vary between 1-2 cm on an average day, which will still introduce variations <1% in the most extreme cases. Moreover, the behavior of both parameters can be considered fairly linear across small intervals within the fiducial range. The plots show that differences between the real and assumed/measured value for the other atmospheric parameters will introduce an error on the zero-point of <1% in the LAST band, with the exception of τ_a500 that can lead to errors of the order of ∼10%. From this exercise, we can identify two parameters with contributions to the photometric error in the LAST band that could exceed the 5 mmag level which are the aerosol τ_a500 and the amount of precipitable water p_w. Therefore, in principle one can consider keeping these parameters free in the model when applied to LAST data.

It is important to remark that the previous considerations are for the specific case of the LAST broad passband. The plots in Figure 2 show how the atmospheric components can affect very differently each single narrower bands (see, e.g., Figure 1. The contribution of each parameter to the photometric error depends on the specific band considered, e.g., the u band is more sensitive to Rayleigh scattering and the aerosol α parameter, while the precipitable amount of water p_w dominates the contribution to the z band. Therefore we can conclude that, in our approach, the choice of the subset of free atmospheric parameters will depend on the specific band to be calibrated.

We repeat the same exercise for variations in air mass, in order to quantify the effect of uncertainties in the estimation of the airmass from the geometrical and empirical approximations used in the model. We select the same passbands used in the previous section and we calculate the variation of the zero point as a function of the air mass (Figure 3). The zero point has a linear dependency on the air mass, and for the case of the broad LAST band it varies by ∼0.3 mag per unit air mass. Therefore, even assuming an uncertainty on the total effective air mass of the order of 1%, the error in the photometric zero point is <1% in the case of LAST. Among the selected bands, the SDSS u filter shows the strongest dependency to air mass variations, with a slope of ∼0.7 mag per unit air mass.

In addition to the uncertainties on the air mass, there is also the difference in the effective air mass in the different parts of the image (compared to the air mass calculated at the image center). This effect increases with the size of the image, and it is more prominent for observations made at large zenith angles. We find that the polynomial corrections in Equation (19) are effective in correcting also this non-uniformity. Alternatively, one can estimate the position dependent air mass for each star.

Fig. 2 Change in zero point magnitude in selected bands as a function of each atmospheric parameter. The data points are obtained by freeing the atmospheric parameter of interest in the transmission model and calculating the variation of the zero point as a function of the parameter. The reference zero point for each curve is calculated for the atmospheric values in Table 2 and T = 15.8 °C.

Fig. 3 Change in zero point magnitude in selected bands as a function of the air mass. The data points are obtained by freeing the air mass while keeping each transmission fixed.

6 Tests on real data from LAST

Here we apply the method to data obtained with the Large Array Survey Telescope (LAST). We first give a brief description of the LAST system and the photometry pipeline. We then present the transmission model for the LAST system and the results of the method on a variety of images.

6.1 The LAST data and instrumental photometry

LAST is an optical telescope system located at the Weizmann Observatory in the Israeli Negev desert. It is an array of 28 cm f/2.2 Rowe-Ackermann-Schmidt telescopes designed to study the variable and transient sky (Ofek & Ben-Ami 2020,Ofek et al. 2023a, Ben-Ami et al. 2023). The full LAST array will consist of 72 telescopes (32 already operative), with each telescope providing a field of view of 3.3 deg × 2.2 deg with a plate scale of 1.25 arcsec/pixel on camera. The limiting magnitude (AB) of the system is 19.6 (21.0) in 20 s (20 × 20 s) exposures.

The application of the absolute photometric calibration method described in this section is based on the products of the LAST data pipeline described in Ofek et al. 2023b. The pipeline processes and calibrates 24 sub-divisions of each single image, and applies two different methods to extract the photometry for each source. The first method is aperture photometry with 3 different apertures (2, 4 and 6 pixel radius) and it is calculated summing all the pixels within each radius, after FFT-shifting the sources to the pixel origin. The second method is PSF-fit photometry, where the PSF is estimated for each sub-image and fit to all sources. As expected, the performance of the aperture photometry is better for bright sources (G < 16 mag) and PSF photometry has better performance at the faint end. Since our procedure uses bright calibrators in the image, our default is to calibrate the aperture photometry (with 6 pixel radius).

6.2 The instrumental transmission of LAST

After introducing a general formulation for the instrumental transmission model in Sect. 3.1, we now define each term specifically for the application to LAST data. In the case of LAST, the optical elements in each telescope are the primary mirror, the corrector plate and the field flattener lenses. Therefore, we rewrite Equation (4) as $$T_{tel}=T_{XLT}{T}_{lenses},$$(21)

where T_XLT is the mirror transmission with the Celestron Star-bright XLT coating and T_lenses the transmission of the corrector plate and field flattener. We consider the reflectivity of the primary mirror (T_XLT) and the transmission of the corrector plate (T_corr) from the StarBright XLT coating technical datasheet². The available measurements cover only a wavelength range that is smaller than the one where the LAST camera is sensitive. Therefore, these transmission curves are modelled with second-order polynomial functions and extrapolated to the desired range. Although measurements for the field flattener lens system are not available, we assume that its transmission has a similar wave dependence to that of the corrector.

The other instrumental term for the transmission is T_det from Equation (5). In the case of LAST, no filter is used and therefore $$T_{det}=T_{QE},$$(22)

where TQE is the quantum efficiency measured for the QHY600-PH CMOS camera mounted on each LAST telescope (Ofek et al. 2023a). Since the QE is dependent on the incidence angle, we expect the nominal QE to vary across the field of view. In order to allow some optimization of the effective camera QE, we fit the laboratory measurements from Ofek et al. (2023a) with the skewed Gaussian function multiplied by a perturbation term (Equations (7) and (8)), with an 8th order Legendre polynomial basis. This procedure is done once and the polynomials parameters listed in Table 1 are held fix during the calibration. We also fix the parameters A, γ and σ, but we leave μ as a free parameter. In Appendix C we show the effects of this method across the field of view when this freedom is not included in the model, and color terms are observed in the final model residuals after the calibration.

To summarize, Figure 4 shows each instrumental component of the transmission and the resulting shape T_OTA (solid red line) given by $$T_{OTA}=T_{tel}{T}_{det}.$$(23)

6.3 Atmospheric parameters for LAST

The only atmospheric parameter currently measured on each telescope mount of the LAST array simultaneously to observations and registered in the image header is the temperature T.

The reported value is used for each image calibration and kept fixed during the model optimization. The simultaneous measurement of the surface pressure p will be included in the future in the data products, meanwhile, in this application we fix it at the average value of 965 mbar (hPa) at the site location. For the Ozone column, measurements³ from stations nearby the LAST site show that u₀ varies slowly with time, with an average value of ∼300 DU and σ ∼ 30 DU. Deviations in this range, correspond to a variation in the photometric zero point <1 mmag. Therefore we choose to fix the value of u₀ to 300 DU in the model. Lastly, given the fact that the aerosol parameter α is highly correlated with the aerosol optical depth τ_a500 and its small contributions to the zero-point variations (see Figure 2), we fix it to the typical value of 0.6 indicated by Gueymard (2019) for desert-like/rural conditions. As shown in Sect. 5, among all six atmospheric parameters, only the aerosol optical depth (τ_a500) and the amount of precipitable water (p_w) can significantly contribute to the photometric error assuming a broadband transmission model. This is also the case for LAST data, which are acquired in clear filter mode. Therefore these two parameters are kept free in our model, and are fitted simultaneously. A summary of the values for the atmospheric parameters used in the calibration of LAST data, their correspondent atmospheric component and their value or status in the fit is presented in Table 2.

Fig. 4 Optical tube assembly instrumental transmission components. The black dots are the laboratory measurements from Ofek et al. 2023a of the QHY600PH QE. The dashed lines show the models of the single instrumental components, while the solid red line shows the total transmission given by the product of the single components.

6.4 Sources with Gaia XP spectra in the LAST fields

An important element required for this approach is to have, in each image, a sufficient number of stellar objects with calibrated spectra. This condition has historically been very difficult to meet, since the coverage of spectroscopy campaigns is usually limited, and their photometric calibration is often poor. For our goal of calibrating each observation independently, it is important to have a large number of calibrators spread uniformly across the sky and across a wide range of magnitudes and colors.

In the Gaia Data Release 3 (Gaia -DR3), a sample of 220 million low-resolution spectra for sources with G < 17.65 mag are available in both the BP and RP bands. These spectra cover the wavelength range from 330-1050 nm (De Angeli et al. 2023; Montegriffo et al. 2023), they are widely distributed across the sky and their coverage of colors and magnitudes make this an excellent sample to be used in our method.

For the application on data from LAST, we generally use the sample of ∼35 million externally-calibrated BP and RP (named as XP in Gaia products) spectra with G ≤ 15 and pre-sampled in the range 336-1020 nm. In the selection process of sources within each sub-image, we select only those with a probability⁴ of being stellar objects >90%. Figure 5 shows the average number of available calibrators per sub-image for each field observed in the all-sky survey footprint of LAST. The median average number of calibrators in each sub-image is ∼130, and 4% of the fields have, on average, <30 calibrators. When the number of calibrators is less than 30, we extend the magnitude threshold for calibrators to G <16, and we retrieve raw spectra for these additional calibrators. When selecting from this larger sample, we externally calibrate the raw spectra with the GaiaXPy⁵ software. We then resample them on an evenly-spaced wavelength array with Δλ = 0.2 nm, to match the sampling of the pre-calibrated spectra. The spectra of fainter stars are typically of lower quality.

As shown in Figure 4, the instrumental transmission for LAST has a percent-level sensitivity in the ranges 300-336 nm and 1020-1100 nm. Therefore, we extrapolate the spectra outside the nominal wavelength range, assigning a constant flux value to F(λ) equal to the first (last) calibrated value in the nominal wavelength range. An example of this extrapolation is shown in Figure 6, where the calibrated spectrum is shown as a continuous line and the extrapolated one as a dashed line.

Fig. 5 Average number of calibrators with Gaia XP spectrum and G < 15 in area of ≈0.3 sq. deg. The vertical dashed orange line marks the median of the distribution (∼130).

Fig. 6 Example of externally calibrated Gaia spectrum (solid blue line) and its extrapolation (dashed orange line) in the 300-1100 nm range.

6.5 Application

For the absolute calibration of LAST data, the model has a total of fourteen free parameters. These are divided as follows: one overall normalization (C, Equation (20)), a set of three parameters in the transmission model ($\hat{\varrho }$, Equation (20)), and a set of ten parameters in the polynomial field correction term (k̂, Equation (20)). $\hat{\varrho }$ contains the center of the QE model (μ) and two atmospheric parameters (p_w and τ_a500). The last set of parameters k̂ are the coefficient of the polynomial field correction terms (Equation (19)) for which we use a basis of Chebyshev polynomials up to the 4th with n_x = 0,1 for n_y ≤ 4 and n_y = 0,1 for n_x ≤ 4. The instrumental coordinates x, y ranging from 0 to 1726 pix in each sub-image are transformed into the range of -1 and 1.

For the model optimization, a particular attention goes in the estimation of the errors to weight the residuals in Equation (20). Since the uncertainty on the Gaia spectra (typically on the order of 1-10% per wavelength bin) dominates the error budget on the photometry of a single source, we estimate the photometric error derived by error propagation of the error in Gaia spectra for each bin, multiplied with a nominal model of the LAST transmission, to obtain an error on the counts rate.

The model is fitted independently on each of the 24 subdivisions of each image (i.e., in a 0.55 deg × 0.55 deg field of view). The fit is performed in two main steps: firstly, the normalization C, the center μ of the QE model and the ten k̂ coefficients of the polynomial field correction are fitted simultaneously. In the second step, these parameters are fixed to the best-fit values and the two atmospheric parameters are freed and fit simultaneously. For each step, 3 iterations are performed and a 3σ clipping is applied in order to remove stellar calibrators that have large residuals. Given the loose criteria in the selection of the stellar calibrators, this procedure helps with removing variable sources and problematic calibrators. Depending on the field, typically a fraction between 10-50% of calibrators are clipped during the iterations, keeping a minimum value of 30 calibrators. The model residuals are minimized using the Levenberg-Marquardt algorithm for least-square fitting (Levenberg 1944; Newville et al. 2015).

Once obtained the best-fit model for the system transmission, the flux measurements of the sources in the subimages are converted into magnitudes in the AB system (Oke 1974; Oke & Gunn 1983) calculating the zero-point using a flat-spectrum source with F(v) = 3631 Jy and following the prescriptions from Bessell & Murphy 2012.

7 Results based on LAST data

We apply the calibration procedure described in the previous section to a sample of LAST data. We compare the results of both single 20 s exposures and 20 × 20 s coadded images, to compare the calibration with increased overall signal-to-noise ratio in the image. The behavior of the calibration residuals throughout the image is summarized in Figure 7, were we show diagnostic plots for a LAST image (24 sub-images fitted independently) for a single 20 s exposure (top block) and for an image resulting from the coaddition of 20×20 s exposures. For the single 20 s exposure, residual distributions (top left) are very similar across all the sub-images, with an overall Gaussian-like distribution centered around 0 with σ= 0.5% in the example shown. The standard deviation of the residual distribution varies slightly depending on the field and observing conditions, but we find it to be generally <1%. The top right panel of Figure 7 shows the fitted transmissions of each sub-image. It can be noticed that the best-fit solutions are very similar to each other, with small position-dependent differences across the image. The mid panels of Figure 7 show the model residual scatter as a function of the instrumental coordinates of the sensors on the X (green dots, mid-left panel) and Y (red dots, mid-right panel) axes for all sub-images. No correlated trends are observed, proving the effectiveness of the polynomial field-correction term in the model. The bottom-left panel shows the residual scatter as a function of the Gaia BP-RP color of the calibrators, and no correlated trends are shown (see Appendix B for an example of color-dependent trends). The bottom-right panel shows the residual scatter as a function of the mean Gaia G magnitude. This shows a slow average increase towards higher magnitudes, as expected from the degrading of the quality of the spectra and photometry.

The bottom block of Figure 7 shows the same residuals distributions this time for a 20 × 20 s coadded image. As expected, the increase in the overall signal-to-noise ratio leads to an improvement of the calibration accuracy.

In Appendix C, we discuss the effect of the polynomial corrections introduced in Sect. 4 on the residuals across the sub-images and on the resulting spatial distribution of the photometric zero point.

7.1 Zero-point error from non-parametric bootstrap

Another way to quantify the accuracy of the model is to estimate the standard error on the zero-point via a non-parametric bootstrap process. To do so, we follow the procedure in Efron & Tibshirani (1994), and realize N calibrator sets for each sub-image by drawing randomly with repetitions from the original set of calibrators. Each random set of calibrators has the same size as the original one and it is used to derive a new set of best-fit parameters for the transmission model. After computing the coordinate-dependent transmission for all N bootstrap representations, the bootstrap standard error on the zero point can be calculated for each point in the image (in our case, we create a grid with 30 pixel spacing). We show in Figure 8 the distributions of the standard error evaluated for all sub-images in a single 20 s exposure (blue) and for all sub-images in a20×20 s coadded image (orange). The distribution for the single 20 s exposure has a median value of ∼5 mmag, with values ranging from ∼2 mmag and a longer tail ranging up to ∼12 mmag. This tail is likely due to the degrading at the very edges of each sub-image. The median value of ∼5 mmag is a good representation of the zeropoint error we get when we apply this test to many observations taken with different LAST units. The orange distribution shows the error on the zero point for each grid point in a 20 × 20 s coadded image. In this case, the median of the distribution is lower at ∼3 mmag, with values ranging from ∼1 mmag to ∼7 mmag.

Fig. 7 Summary of residual distributions for a single 20 s exposure (top) and for 20×20 s coadded exposures (bottom). The top-left panel shows the overall distribution of the residuals, with the best-fit standard deviation assuming a Gaussian shape. The top right panel shows the 24 transmissions of each sub-image of a single image, fitted independently. The middle panels show the residuals (for all sub-images) as a function of the instrumental coordinates of the camera sensor. The bottom panels show the residuals as a function of the Gaia color (left) and the Gaia G magnitude (right), respectively.

7.2 Stability over time

In order to quantify the stability over time of the photometric calibration, we consider a series of 20×20 s observations per night of the same field, taken for a total of 8 nights between 3-12 March 2024. Each batch of 20 images is observed continuously, and in LAST defines a standard visit for which coadded images are created. After calibrating each sub-image individually, we consider all sources within a single sub-image. We show in Figure 9 the robust standard deviation (σ′)⁶ as a function of the average AB magnitude for all sources (blue points) within the sub-image with at least 80 detections in the whole sample of observations. The points are calculated for both aperture and PSF photometry for each source, then the one with the smaller σ′ is chosen. Typically, sources with M_AB ≤ 16 have better σ′ with aperture photometry. The scatter reaches a 5-7 mmag level towards the bright end of the distribution. At MAB < 12 the sources are saturated in the image. The orange points show the σ′ versus M_AB for the same sources in coadded images (20×20 s) for the same observations. Compared to the blue points, the distribution spreads on a wider range of σ′ because it is calculated on only 8 images resulting from the coaddition of the 160 original images into groups of 20. However, the improvement in stability from these images is clear, as the σ′ reaches values as low as a few mmag on the bright end, and stay below 1% up to M_AB ∼ 16.

It is worth noting that secondary color terms might be present in the calibration, mostly coming from the reddening of sources at high air masses. In the application to LAST data we find these secondary color terms are similar between different telescopes and stable in observations form multiple nights. The color term is linear with typical values of 10-50 mmag/airmass, depending on the color of the source. This correction can be easily calculated after the main photometric calibration, by having two or more observations at different airmass.

Fig. 8 Distributions of the bootstrap standard error on the AB zero point for all LAST sub-images of a single exposure (blue) and for 20 × 20 s coadded images (orange). The dashed lines mark the medians of each distribution.

7.3 Nightly evolution of the transmission

It is interesting to have a glimpse at how the transmission varies during the night and between different nights. To do so, in the upper panel of Figure 10 we show the evolution of the overall normalization parameter C from Eq. (20), which dominates the amplitude of the transmission as function of time. We show data from the same nights of the temporal stability study in Sect. 7.2, but this time we use observations from several fields observed throughout each night during the sky survey operations. We see, as expected in the period considered (March 2024), a worsening of the conditions of the site during the night, mostly due to the increased humidity and formation of thin layers of clouds. During the fifth night we see more dramatic variations of conditions with gradual drops of C from ∼0.55 to ≤0.3 within 1-2 hours and then subsequent gradual improvements. In the lower panel we show the corresponding airmass of each observations, to remark that these trends are not correlated with airmass changes, but can be interpreted as a broad metric of the observing conditions. However, we stress on the fact that the parameter C is an overall normalization of the model, therefore its long-term variation can be due also to instrumental effects. For example, we notice average increases of this parameter up to 10-20% after regular cleaning operations of the telescope corrector plates.

Fig. 9 Robust standard deviation (σ′) as a function of the mean AB magnitude for all sources in a LAST sub-image using data from 20 single images per night (blue points) and one 20×20 s coadded image per night (orange points) for a total of 8 nights. The points are calculated for both aperture and psf photometry for each source, then the one with smaller σ′ is chosen.

Fig. 10 Upper panel: overall normalization parameter for several skysurvey observations obtained during the same nights of the temporal stability study in Sect. 7.2. Each blue point represents the parameter C (from Eq. (20)) for a coadded image. Lower panel: corresponding airmass of each observation (orange points).

8 Software

The implementation of this method is publicly available on Github² in its python implementation. A MATLAB version is under development and it will be included in the AstroPack/MAATv2⁸ Astronomy & Astrophysics Toolbox (Ofek 2014). For the atmospheric components, our code is based on the parameterizations and implementations from the SMARTS v2.9.8 model (Gueymard 2019) and its python implementation from Ruiz-Arias 2022. Gaia spectra are retrieved either through Astroquery (Ginsburg et al. 2019) or a fast local catalog matching through catsHTM (Soumagnac & Ofek 2018). Lastly, the model optimization is performed using the lmfit package (Newville et al. 2015).

At the date of this preprint, the software contains mainly libraries to apply the model on LAST data. Libraries to handle data products from other instruments will be added in the futures, as well as more generic base models to make the code more versatile to any filter shape. In Appendix E we suggest formats to report the calibration parameter in FITS headers.

9 Summary

We presented a novel approach to the absolute photometric calibration of ground-based telescopes. This method consists of deriving the transmission of the system and aims to break the 1% barrier of accuracy that is typical of the current state-of-the-art absolute calibrations. Furthermore, it is designed to over come the problem of comparing photometric data obtained using different telescopes, or the same telescope at different atmospheric conditions. In this case, a rough spectral model of the source can be used to convert observations taken with different instruments to the same passband. A key point of this method is that the transmission model can be fit independently for each image, given a large enough number of calibrators in the field of view.

From an application to LAST data we showed that, by fitting an appropriate subset of the instrumental and atmospheric parameters for each image, we achieve a median accuracy of the photometric zero point that ranges between 3-5 mmag, depending on the total exposure of the image. Moreover, stability of the zero-point accuracy is observed over long time scales. The improvement observed in 20×20 s coadded images is of a factor ∼1.5-2. We notice that this is lower than a factor of $\sqrt{n}\sim$4.5, likely because we are reaching the systematic error level in Gaia photometry of about 3 mmag. Lastly, the spatial photometric corrections that are simultaneously fitted with the transmission provide a uniform accuracy across the whole image and do not show dependencies on the spectral color of the calibrators.

Thanks to its high versatility, this approach can be applied to any ground-based facility that has a large enough field of view, with the appropriate changes on the model assumptions for the instrumental response. The image-by-image calibration approach also allows for an efficient calibration of existing and future data from other telescopes. This approach is especially useful for wide field surveys like LAST, for which the high multiplicity of the instrument array requires an accurate, image-wise calibration approach. Since one of the basic requirements for this approach is to have a fair number of Gaia spectra, this approach is more efficient for wide-field instruments. However, if the current sample of ∼220 million spectra measured by Gaia will increase in the future, it will increase the efficiency of this method also for instruments with a smaller field of view.

The method presented in this work will be adopted as the official absolute calibration for the LAST telescopes and it will be expanded to be applied also to monitor the evolution of the system throughput in space telescopes (e.g., ULTRASAT; Shvartzvald et al. 2024). With space telescopes, dedicated campaigns observing a sample of calibrators at regular time intervals will be sufficient to fit the transmission from multiple images.

Acknowledgements

We would like to thank the anonymous referee for useful comments on the manuscript. We would also like to thank G. Altavilla for useful comments on the preprint. S.G. is grateful for the support of the Koshland Family Foundation. E.O.O. is grateful for the support of grants from the Will-ner Family Leadership Institute, André Deloro Institute, Paul and Tina Gardner, The Norman E Alexander Family M Foundation ULTRASAT Data Center Fund, Israel Science Foundation, Israeli Ministry of Science, Minerva, BSF, BSF-transformative, NSF-BSF, Israel Council for Higher Education (VATAT), Sagol Weizmann-MIT, Yeda-Sela, and the Rosa and Emilio Segre Research Award. This research is supported by the Israeli Council for Higher Education (CHE) via the Weizmann Data Science Research Center, and by a research grant from the Estate of Harry Schutzman.

Appendix A Sensitivity to variations of the transmission average wavelength

We provide here an example to compare the sensitivity of small variations in the average wavelength of a passband. We consider a nominal transmission for the LAST system, described in Sect. 6.2 and we vary its average wavelength (λ₀, Koornneef et al. 1986) defined as $${\lambda }_0=\frac{\int T(\lambda )\lambda d\lambda }{\int T(\lambda )d\lambda }$$(A.1)

between 547.0 nm and 551.0 nm. We select two Gaia spectra from a blue and a red star, respectively, and we calculate the magnitude of the source as observed in each passband assuming a zero point of 25 mag. Figure A.1 shows the magnitude difference in each transmission calculated starting from the nominal passband with λ₀ = 549.05 nm. We see how these small variations in the transmission have a small effect (< 10⁻³) for the blue star, while they reflect in a difference of a few % for the red star. This example shows how even small changes in the system transmission can affect dramatically the photometric stability.

Fig. A.1 Comparison of the magnitude variation as a function of the effective wavelength of the transmission between a blue and a red star.

Appendix B The importance of fitting the transmission center μ

Here we explore the effect of fitting the shape of the analytical model used for the detector quantum efficiency. We run the calibration on the single 20 s exposure calibrated in Sect. 7 by keeping fixed the parameter μ of the skewed Gaussian model at μ = 390.1 (the best-fit value on the data from Ofek et al. 2023a). Figure B.1 shows the calibration residuals for all sources with a Gaia counterpart as function of their color. Differently to the behavior in Figure 7, in this case a correlation with the Gaia BP-RP color can be seen in the distribution. For this particular case, the correlation coefficient is ∼ −0.3, while stronger correlations can be generally found in individual subframes throughout the image.

Appendix C The importance of the zero-point field corrections and the effect of gray structures

We repeat the same exercise but this time we remove the polynomial field correction terms from the model. We show in Figure C.1 an example of calibration residuals in a LAST subimage when the field correction terms from Equation (20) are removed. A strong gradient of the residuals from the bottom left to the top right corner of the frame can be distinguished in the colormap.

Fig. B.1 Residuals vs Gaia BP-RP color for a fixed-shape QE analytical model in the calibration procedure.

Fig. C.1 Residuals as function of the detector coordinates in a LAST sub-image. In this case, the field correction terms are removed from the model. The residuals are calculated for each calibrator star (black dots) and linearly interpolated in a grid.

When present, these gradients are observed with different patterns that only depend on the specific image. Therefore these might be related to residual non-uniformity after the flat-field correction or some temporary instrumental gradient on the camera sensor. We observe non-uniformity in the residuals up to a ∼10% level in some of the most extreme cases. As already mentioned in Sect. 4, the polynomial field corrections (together with the overall model normalization), are important to correct for wavelength-independent (gray) extinction caused by clouds. However, it is hard in the current formulation of the method to disentangle the contribution of gray extinction and instrumental non-uniformity. Moreover, other effects might contribute significantly (e.g., PSF variation across the field of view). Nevertheless, from the spatial distribution of the zero-point, it is possible to distinguish cases in which instrumental non-uniformities are dominant, to cases where the field is dominated by gray structures. In the left panel of Figure C.2 we show an image in which the zero point spatial distribution is dominated by the typical vignetting pattern (likely residuals from flat-field corrections). In the right panel, we show an image in which the gray structures are likely dominant, and affect significantly the spatial distribution of the zero-point values. It is interesting to notice how the structures are consistent across multiple sub-images, which are calibrated independently. This supports the robustness and consistency of our approach. This method, utilizing a large number of LAST telescopes observing a continuous patch of the sky, may be used to measure the power spectrum of the gray extinction. This may be investigated in a future work.

Appendix D Airmass parameterizations for atmospheric components

Here we report in Table D.1 the parameterizations from the SMARTS v2.9.8 model (Gueymard 2019) used to calculate the effective airmass of each atmospheric component. The details on which component is included for each atmospheric transmission are described in Sect. 3.3.

Appendix E Suggested standard for FITS header keywords

Here we suggest a common scheme for FITS⁹ header keywords to report the values of the model parameters in each image. Assuming the more general form of the transmission model, about 30 parameters are needed (e.g., Tables 1,2).

In our scheme, each header keyword will appear in the format PT_〈PARNUM〉_〈TYPE〉, where PT defines parameters for the photometric transmission, 〈PARNUM〉 is a 3-digits running number to distinguish between different parameters, and 〈TYPE〉 can be:

• N: for keywords storing the extended unique name for the parameter (e.g., ATM_aerosol_alpha);

• V: for keywords storing the parameter value;

• E: optional, for keywords storing the parameter error;

• F: for keywords storing a boolean flag associated with the parameter that indicates if the parameter was free in the fit (1) or fixed (0).

Each parameter keyword will then appear up to four times with a different 〈TYPE〉. Note that we do not force a unique mapping between the name and 〈PARNUM〉. In addition, we suggest special keywords for the whole photometric calibration that are listed in Table E.1.

The user has the freedom to use custom names for parameters, but we want to stress on the importance to adopt the same convention in the name structure of the header keywords. We list examples for the header keywords used in this work in Table E.2.

We suggest also a keyword scheme to store the transmission after sampling it in a number of wavelength bins N_λ and to report the value of the transmission in each wavelength bin. This will also require to report four additional special header keywords (Table E.1) for the minimum wavelength (PT_LMIN), maximum wavelength (PT_LMAX) sampling width (PT_LDEL) and wavelength units (PT_LUN). Alternatively, a sampled version of the QE or the full transmission can be stored in an external file which path is stored in the special keywords PT_QE and PT_TRSM. These special keywords are listed in Table E.1. The transmission values are stored for each bin in three keywords:

• PT_〈BINNUM〉_L which contains the central wavelength of the bin;

• PT_〈BINNUM〉_T which contains the transmission value in the bin;

• PT_〈BINNUM〉_E optional, which contains the transmission error in the bin.

Also here, 〈BINNUM〉 is a running 3-digits number. These keywords will therefore have the form, e.g., for bin 1: PT_001_L, PT_001_T and PT_001_E.

Independently from the convention adopted for the header keywords, each source in the catalog results will have an additional column that is the AB zero-point calculated at the physical location of the source on the camera sensor.

Fig. C.2 Zero-point 2-D distribution for an image with low presence of cloud structures (left) and one with significant presence of cloud structures (right). The dashed lines mark the sub-images, which are calibrated independently.

References

All Tables

All Figures

Fig. 1 Transmittance of each atmospheric component at different zenith angles. The panels in the top row show from left to right the transmittance for Rayleigh scattering, aerosol scattering and ozone absorption. In the bottom panels, left to right panels are the transmittance of the absorption by water vapor, absorption by the uniform mixture of gases, and the product of all atmospheric components.

In the text

	Fig. 2 Change in zero point magnitude in selected bands as a function of each atmospheric parameter. The data points are obtained by freeing the atmospheric parameter of interest in the transmission model and calculating the variation of the zero point as a function of the parameter. The reference zero point for each curve is calculated for the atmospheric values in Table 2 and T = 15.8 °C.
In the text

	Fig. 3 Change in zero point magnitude in selected bands as a function of the air mass. The data points are obtained by freeing the air mass while keeping each transmission fixed.
In the text

	Fig. 4 Optical tube assembly instrumental transmission components. The black dots are the laboratory measurements from Ofek et al. 2023a of the QHY600PH QE. The dashed lines show the models of the single instrumental components, while the solid red line shows the total transmission given by the product of the single components.
In the text

	Fig. 5 Average number of calibrators with Gaia XP spectrum and G < 15 in area of ≈0.3 sq. deg. The vertical dashed orange line marks the median of the distribution (∼130).
In the text

	Fig. 6 Example of externally calibrated Gaia spectrum (solid blue line) and its extrapolation (dashed orange line) in the 300-1100 nm range.
In the text

Fig. 7 Summary of residual distributions for a single 20 s exposure (top) and for 20×20 s coadded exposures (bottom). The top-left panel shows the overall distribution of the residuals, with the best-fit standard deviation assuming a Gaussian shape. The top right panel shows the 24 transmissions of each sub-image of a single image, fitted independently. The middle panels show the residuals (for all sub-images) as a function of the instrumental coordinates of the camera sensor. The bottom panels show the residuals as a function of the Gaia color (left) and the Gaia G magnitude (right), respectively.

In the text

	Fig. 8 Distributions of the bootstrap standard error on the AB zero point for all LAST sub-images of a single exposure (blue) and for 20 × 20 s coadded images (orange). The dashed lines mark the medians of each distribution.
In the text

	Fig. 9 Robust standard deviation (σ′) as a function of the mean AB magnitude for all sources in a LAST sub-image using data from 20 single images per night (blue points) and one 20×20 s coadded image per night (orange points) for a total of 8 nights. The points are calculated for both aperture and psf photometry for each source, then the one with smaller σ′ is chosen.
In the text

	Fig. 10 Upper panel: overall normalization parameter for several skysurvey observations obtained during the same nights of the temporal stability study in Sect. 7.2. Each blue point represents the parameter C (from Eq. (20)) for a coadded image. Lower panel: corresponding airmass of each observation (orange points).
In the text

	Fig. A.1 Comparison of the magnitude variation as a function of the effective wavelength of the transmission between a blue and a red star.
In the text

	Fig. B.1 Residuals vs Gaia BP-RP color for a fixed-shape QE analytical model in the calibration procedure.
In the text

	Fig. C.1 Residuals as function of the detector coordinates in a LAST sub-image. In this case, the field correction terms are removed from the model. The residuals are calculated for each calibrator star (black dots) and linearly interpolated in a grid.
In the text

	Fig. C.2 Zero-point 2-D distribution for an image with low presence of cloud structures (left) and one with significant presence of cloud structures (right). The dashed lines mark the sub-images, which are calibrated independently.
In the text