5.3.6 External calibration

Author(s): Paolo Montegriffo

After the internal calibration has put all the observations onto a common photometric system and the reference fluxes have been computed for all sources, the external calibration effectively provides a means to relate this to other photometric systems and allows for a physical interpretation of the data. This is done providing a response function and an absolute zero point for each pass band. The strategy adopted for the pass band modelling is based on a forward model approach: a response model whose shape depends on a small number of adjustable parameters is fine tuned till the predicted integrated magnitudes fit with the internally calibrated ones for a set of flux calibrators (the SPSS). The model has been designed to reproduce the nominal response curves when all parameters are set to zero. Two different models have been developed: a basic one where the nominal curve is multiplied by a low order polynomial to model possible large scale differences between actual and nominal response (hereafter the Analytic model), and a second one where the nominal curve is summed up with a linear combination of special basis functions derived with Principal Component Analysis (hereafter PCA) to model residuals with the actual filter curve (hereafter the Hybrid model).

Nominal Passband Model

The nominal photonic pass band is modelled as the product of the following quantities:

RNom(λ)={T0(λ)ρatt(λ)Q(λ)forGT0(λ)ρatt(λ)Q(λ)Tp(λ)forXP (5.8)

where

  1. 1.

    T0(λ) is the telescope (mirrors) reflectivity;

  2. 2.

    ρatt(λ) is the attenuation due to rugosity (small-scale variations in smoothness of the surface) and molecular contamination of the mirrors;

  3. 3.

    Q(λ) is the CCD QE;

  4. 4.

    Tp(λ) is the prism (fused silica) transmittance curve (including filter coating on their surface) for the BP/RP case.

All these quantities were initially measured by Airbus DS during on-ground laboratory test campaigns. However, the filter wavelength cut-ons and cut-offs in the prism transmittance curves change over the focal plane due to uneven thickness of the prism coating producing bandwidth non uniformity: the expected variation is of about 3 nm in the position of the BP cut-off and of about 2 nm in the position of the RP cut-on (hereafter for convenience only the generic term cut-off will be used for both instruments). Moreover, combining measurements taken all over the focal plane should result in a mean instrument passband with smeared cut-offs. To properly take into account these effects the prism transmittance curve is modelled as follows:

Tp(λ)=Tp(λ)×Rcut(λ;μ,σ) (5.9)

where Tp(λ) is the tabulated transmittance curve provided by Airbus DS truncated right before and after the cut-off wavelength for the BP and RP case respectively:

Tp(λ)BP={Tp(λ)forλ<647nmTp(647nm)forλ647nm (5.10)
Tp(λ)RP={Tp(λ)forλ>650nmTp(650nm)forλ650nm (5.11)

and Rcut(λ;μ,σ) is a parametric curve based on the Gauss error function and its complementary function:

Rcut(λ;μ,σ)={12erfc((λ-λco-μ)σco+|σ|)forBP12(1+erf((λ-λco-μ)σco+|σ|))forRP (5.12)

where λco represents the nominal cut-off wavelength (659.6 nm for BP, 640.7 nm for RP), and σco represents the nominal cut-off slope (4.71 for BP, 3.26 for RP). When the two parameters μ and σ are set to zero the nominal transmissivity curve is reproduced to a relative accuracy better than 0.5%.

Response Analytic Model

The Response Analytic Model represents the mean system photon-counting response function and is the product of the Nominal Passband model with a polynomial function of low degree:

S(λ)=RNom(λ;μ,σ)×i=0degpiλi (5.13)

This model has been designed to correct for large scale discrepancies between the actual and the nominal response curves, such as those produced by the telescope contamination effect.

Response Hybrid Model

The Response Hybrid Model is an alternative model for the mean system photon-counting response function, and consists of the sum of the Nominal Passband model and a linear combination of special basis functions computed with PCA on a large set of theoretical response curves (hereafter the training set):

S(λ)=RNom(λ;μ,σ)+i=0nRriRi(λ) (5.14)

The training set was built from the extensive set of measurements made by Airbus DS of all the response curve components seen in Eq. (5.8) and available in the DPAC Calibration Data Base (CDB). A MonteCarlo-like procedure was run to apply the relative errors to all measured curves and generate a set of 10,000 different response curves for each pass band. The nominal passband model is then subtracted from these curves, and residuals are arranged in a matrix which is decomposed with Singular Value Decomposition (SVD) into the product of three arrays B×D×BT, where D is a diagonal matrix whose elements (singular values) are in non-increasing order, and B is an orthonormal matrix whose columns form the orthonormal basis set Ri(λ) used in Eq. (5.14). The ordering of singular values ensures that the best approximation of the residuals is achieved for any given number of parameters nR.

The integrated photometry Instrument Model

Given a source with an energy flux distribution per wavelength unit fλ(λ) in units of W nm-1 m-2, the corresponding predicted mean flux can be modelled as:

=PA(fλ(λ)λ109c)S(λ)dλ (5.15)

where PA is the telescope pupil area in unit of m2, is the Planck constant, c is the velocity of light in vacuum, and S(λ) is the mean system photon-counting response function model; in practice, the quantity between brackets represents the source photon flux distribution in units of photons s-1 nm-1 m-2 where λ is expressed in nm.

Given a set of N SPSS, the optimal set of parameters of S(λ) model to represent the actual photometric system pass band is found by minimising the cost function defined as:

cost=i=1NLδ(ai) (5.16)

where Lδ(ai) is the Huber loss function defined as:

Lδ(ai)={12ai2for|ai|δδ(|ai|-12δ)otherwise (5.17)

and the variable ai refers to the weighted residuals for the ith SPSS, i.e., the difference between observed and predicted mean flux normalised by the corresponding standard deviation:

ai=Ii¯-iσi (5.18)

When the instrument model implies the fit of the response cut-off parameters, as in the BP and RP case, since the model is not linear with respect to the fitting parameters, the minimisation of the cost function is achieved with the Differential Evolution algorithm described by Storn and Price (1997); in all other cases this is obtained by standard linear least squares regression. The usage of the Huber loss function should give intrinsically more robustness against outliers with respect to a traditional χ-squared loss based algorithm; however, all solvers implement a sigma rejection scheme to under-weight possible outliers in the data.

Solving for the model parameters in the BP/RP case

It can be demonstrated that the integrated magnitudes of the SPSS do not carry enough information to determine the filter cut-off shape with sufficient accuracy, so prior knowledge was needed to constrain the possible variation of the cut-off parameters: the Gaia mean BP and RP spectra could in principle provide such knowledge, however at the time of Gaia DR2 data analysis this kind of information was not yet available, therefore the fitting algorithm was set to allow only for a maximum cut-off shift of 3 nm around the nominal position, according to the Airbus DS findings about the bandwidth variation, as already mentioned in Section 5.3.6. Moreover, the difference between actual and nominal cut-off shapes was expected to be a 2nd order correction to the nominal passbands and, in order to avoid as much as possible local minima in the complex cost function space, an iterative fitting scheme was implemented, where:

  • the cut-off parameters are kept fixed while the overall response shape parameters (pi coefficients for Eq. (5.13) or ri coefficients for Eq. (5.14)) are fitted with a linear solver;

  • the response shape parameters are kept fixed while the cut-off parameters (μ and σ of Eq. (5.8)) are evaluated with the Differential Evolution solver;

The two steps were repeated till convergence was reached.

The G, GBP and GRP magnitude scales in the VEGAMAG system

The weighted mean flux I¯ provided by the internal calibration can be used to define an instrumental magnitude:

Ginstr=-2.5logI¯ (5.19)

The Gaia magnitude G scale is defined by adding a zero point, G0, to the instrumental magnitude, as:

GGinstr+G0 (5.20)

Traditionally, the value of G0 can be determined by the comparison between the synthetic photometry computed on a set of calibrators with known SEDs and their corresponding instrumental magnitudes: however, since the current passband model has been optimized to reproduce the measured SPSS mean fluxes, the zeropoint computation is simpler. The specific algorithm to compute synthetic photometry depends on the adopted magnitude system: in the Gaia DR2 case, in accordance with Gaia DR1, the photometry is computed in the VEGAMAG system, this being the de-facto standard used internally by the Gaia data processing during the years of preparation of the mission.

The VEGAMAG system is defined so that Vega’s colours are all zero: this is equivalent to normalizing all the observed fluxes to the flux of Vega. The reference spectrum chosen for the Gaia photometric system is a Kurucz/ATLAS9 model (Buser and Kurucz 1992) for the Vega spectrum (CDROM 19) with Teff = 9550 K, log g = 3.95 dex and [M/H]=-0.5 dex, and vmicro = 2 km s-1. Its visual magnitude is assumed to be V0 = 0.0 mag.

Given a source with an energy flux distribution fλ(λ), the integrated energy flux per unit wavelength interval fλ, measured by a photon-counting detector with response S(λ), can be expressed as:

fλ=fλ(λ)S(λ)λ𝑑λ (5.21)

The synthetic magnitude in the VEGAMAG system is given by:

G = -2.5logfλ(λ)S(λ)λ𝑑λfλVega(λ)S(λ)λ𝑑λ+GVega (5.22)

Since all Vega colours by definition are zero, GVega=0.0 mag. Isolating the zero point G0 from Eq. (5.20) and substituting G with the expression in Eq. (5.22), one can set:

G0 = -2.5log(fλ(λ)S(λ)λ𝑑λI¯) (5.23)
+2.5log(fλVega(λ)S(λ)λ𝑑λ)

Using Eq. (5.15) one obtains:

fλ(λ)S(λ)λ𝑑λI¯I¯×109cPA (5.24)

Given the passband model optimization described in Section 5.3.6, the expectation value for I¯ is 1 and Eq. (5.23) becomes:

G0=+2.5log(PA(fλVega(λ)λ109c)S(λ)𝑑λ) (5.25)

The error on magnitudes associated to G0 can be computed by propagating the errors on S(λ) parameters into the above equation. The final uncertainty in the externally calibrated magnitude (G) in Eq. (5.20)) is

σG=(1.0857σI¯I¯)2+(σG0)2 (5.26)

where σI¯ is given by the internal calibration.

Zero point computation in the AB scale

The AB system (Oke and Gunn 1983) is defined in such a way that an object with constant flux per unit frequency interval has zero colour:

ABν=-2.5logfν-48.60 (5.27)

where fν is the flux in erg cm-2 s-1 Hz-1. The constant term is defined to set AB = 0 mag for a source s0 with fν,s0=3.63110-20 erg s-1 Hz-1 cm-2.

Gaia fluxes are expressed in SI units (W m-2 Hz-1) and hence the value of the constant becomes -56.10 mag. According to Bessell and Murphy (2012), Eq. 5.27 can be generalised to be used with broad photometric bands:

AB=-2.5logfν-56.10 (5.28)

with

fν=10-9fλ(λ)S(λ)λ𝑑λS(λ)c𝑑λ/λ (5.29)

Using Eq. 5.28 to describe Gaia external magnitudes in the AB system, GAB, and isolating the zero point from Eq. 5.20 one gets:

G0,AB=GAB-Ginstr=-2.5log(fνIs¯)-56.10 (5.30)

The quantity Q=fνI¯ can be estimated by taking the weighted mean of all the values measured of the SPSS:

Qs=1NSPSSwsqss=1NSPSSws (5.31)

where qs=fν,sIs¯ and ws=1σqs2
being

σqs=|qs|(σfλ,sfλ,s)2+(σIs¯Is¯)2 (5.32)

where σfλ,s is computed by propagating the error of the SPSS SEDs ground-based measurements with Eq. 5.21, while σIs¯ is given by the internal calibration.

The standard uncertainty on the weighted mean Q is computed by using an approximation formula given by Cochran (1977):

σQ2 = NSPSS(NSPSS-1)(s=1NSPSSws)2[s=1NSPSS(wsqs-w¯Q)2 (5.33)
-2Qs=1NSPSS(ws-w¯)(wsqs-w¯Q)+Q2s=1NSPSS(ws-w¯)2]

Summarizing the AB scenario, the externally calibrated magnitude (GAB) and its uncertainty (σGAB) in the AB scale for a source with an internally calibrated mean flux Is¯ is given by:

GAB = -2.5logIs¯-2.5log(Q)-56.10 (5.34)
σGAB = (1.0857σIs¯Is¯)2+(1.0857σQQ)2 (5.35)

Results

The data exported from internal calibration which matched the SPSS IDs consisted of 93 records out of the 94 sources present in the V1 SPSS list (one star was missed by the cross-match identification process because of its high proper motion): a set of 64 SPSS were extracted from the gold sources (i.e. sources with colour in the GBP-GRP =0.0-4.0 range, where the time link-calibration could be applied, see Section 5.4.3) set while the remaining 29 were recovered from the silver set sources (i.e. sources with colour outside the GBP-GRP =0.0-4.0 range, where an ad hoc calibration had to be applied, see Section 5.4.3). The distribution of collected data as a function of the nominal G magnitude and GBP-GRP colour is shown in Figure 5.7.

Figure 5.7: Left: Synthetic G magnitude distribution of the SPSS exported by the internal calibration pipeline: yellow and grey colours represent the distribution of gold and silver sources respectively, while blue indicates the missing SPSS (not identified by the cross-match, see text). Right: same data as in left panel, plotted against the synthetic GBP-GRP colour.

To provide a visual comparison between the Gaia DR2 photometric system properties and the pre-launch nominal one, we show in Figure 5.8 the difference between synthetic and instrumental magnitudes as a function of the GBP-GRP colour for the 93 SPSS: while the G magnitudes are comparable with internally calibrated photometry, the GBP data exhibit a complex behaviour and clearly show a rather strong colour term for sources bluer than GBP-GRP -0.2 mag. This colour term is presumably due to the combination of internal calibrators colour selection with time link calibrations to correct for telescopes transmissivity variations due to contamination. The GRP case shows a clear colour term with redder sources being fainter than expected.

Figure 5.8: Comparison between synthetic and instrumental magnitudes as a function of source GBP-GRP colour for G (top), GBP (centre) and GRP (bottom) obtained by adopting the nominal passbands. The point colours encode the G magnitude of the sources (red dots showing fainter sources); triangles represent the silver sources. The black horizontal line represents the mean value of the plotted values.

Figure 5.9 shows the calibrated Gaia pass bands for the G, GBP and GRP photometric systems; each curve is compared against the nominal counterpart and the ±1σ level is represented as a shaded region. The G and GBP curves have been modelled by the Response Hybrid Model using respectively 4 and 3 basis functions. The GBP curve exhibits a pronounced bump in the blue side with respect to the nominal curve, that is needed to correct for the colour term left into the silver sources data by the Time Link Calibration. The GRP case was more problematic to optimize, because the output curve resulting from the Hybrid Response Model had an implausible shape independently of the number of basis function used: this is an evidence that the actual pass band shape is significantly different from the nominal one and the PCA based representation fails to give an acceptable result. For this reason a 1st degree polynomial function was used for the Analytic Response Model of the GRP.

Figure 5.9: Calibrated Gaia pass bands for G (top), GBP (centre) and GRP (bottom) Gaia DR2 photometric systems. Shaded areas represent the ±1σ level. Black lines represent nominal response curves while the coloured ones are the calibrated curves.(Evans et al. 2018)

Figure 5.10: Left: Comparison between synthetic and instrumental magnitudes as a function of source GBP-GRP colour for G (top), GBP (centre) and GRP (bottom). The point colours encode the G magnitude of the sources (red dots showing fainter sources); triangles represent the silver sources. Right: Residuals with respect to the fitted zeropoints as a function of source magnitude; the point colours encode the GBP-GRP colour of the sources (red dots showing reddest sources).(Evans et al. 2018)

Once the calibrated pass bands are used to recompute synthetic photometry, the comparison with SPSS mean fluxes shows that the colour terms are effectively removed (Figure 5.10 top panels). The solid black lines represent the derived external calibration zero points which are given in Table 5.2.

Table 5.2: Photometric zeropoints in the VEGAMAG system.
Band ZP error
G 25.6884 0.0018
GBP 25.3514 0.0014
GRP 24.7619 0.0019

These passbands were released for internal DPAC use in July 2017 and adopted in the photometric processing to calibrate the magnitudes to be released in Gaia DR2. These passbands are available to be downloaded at the following url: https://www.cosmos.esa.int/documents/29201/1645651/GaiaDR2_Passbands_ZeroPoints.zip.

A following analysis showed that the BP/RP spectra were better described by taking into account an offset in the RP cutoff wavelength position. A new set of revised passbands was then derived in October 2017 to introduce this and a few other minor changes. These updated passbands can be downloaded at the following url: https://www.cosmos.esa.int/documents/29201/1645651/GaiaDR2_Revised_Passbands_ZeroPoints.zip.

The left panels of Figure 5.10 show residuals with respect to the adopted zeropoint as function of G, GBP and GRP magnitudes respectively to investigate for possible systematics left by the internal calibration. The only appreciable effect is visible in the G case, where sources with G 12 mag seem to have slightly higher residuals with respect to fainter sources.

Finally, Table 5.3 gives the zeropoints for the three passbands in the AB system: these values can be used to convert Gaia photometry to the AB system.

Table 5.3: Photometric zeropoints in the AB system.
Band ZP error
G 25.7934 0.0018
GBP 25.3806 0.0014
GRP 25.1161 0.0019