Test of the temporal stability of the differential photometry in MUSE

Edmund Christian Herenz (ESO Fellow),
Fuyan Bian (ESO)

Rationale

The MUSE integral field spectrograph records spectra over its field of view on 24 unit spectrographs. Spectrophotometric calibration of the science data uses a single response function computed from standard star observations. Differences in the response between different spectrographs are accounted for by spectrally slowly varying multiplicative corrections that are calculated from twilight flat field observations and internal illumination lamp exposures. (For more information see the documentation for the MUSE data reduction system and Weilbacher et al. 2020.)

Scientific applications that aim at monitoring spectral variability with MUSE require differential spectrophotometry to be temporally stable. This means that for a time series of observations the ratios of measured flux-densities across the field of view stay constant between sources without variability. Examples for such applications for MUSE are tests for stellar-variability (Giesers et al. 2018) or quasar reverberation mapping (as suggested in Selman et al. 2019). So far quantitative tests on the stability of the differential spectrophotometry of MUSE have been performed on multi-epoch observations of globular clusters (Giesers et al. 2018; see also Appendix A.1 in Giesers 2020). According to these analyses the average spectrophotometric variations across the field of view are ≲ 0.05 mag for high signal-to-noise ratio (SNR ≳ 20) spectra. Moreover, Weilbacher et al. (2015, Sect. 3.3) state that the relative flux calibration is accurate within 5 %.

We here present an independent test of the differential spectrophotometric accuracy for MUSE wide-field mode observations. For this purpose we use regular standard star observations that are taken as part of the MUSE calibration plan. For our experiment we exploit that some standard star fields contain multiple stars. In the present analysis we focus on stars in the field around the spectrophotometric standard LTT 3218.

Our experiment is thought as a first step to quantify and monitor the temporal stability of the differential spectrophotometry. In order to facilitate further analyses, we provide the Python 3 code used in the present analyses. We also provide an a link to the files that store all the extracted spectra in the MUSE field of view around LTT 3218.

Method

Data retrival

We first need to retrieve the observations of the standard star LTT 3218 from the ESO Science Archive Facility (SAF). For this purpose we query the SAF with the the python3 library astroquery (Ginsburg et al. 2019), that in the used version (0.4.2.dev6684) provides a front-end to the SAF's MUSE Raw Data Query Form.

In order to properly reduce these standard star observations, we require the associated calibration data being downloaded as well. For normal science observations, the SAF automatically provides these associations with the CalSelector tool. However, no associated calibrations are offered when the SAF is queried to download a standard star exposure. As a workaround, we queried the SAF for science observations taken in the same night and in the same instrument mode as the standard, since this delivers the required calibration data for reduction of the standard itself. This means that some standard star observations could not be used in the present analysis, as in some nights no science observations were taken in an instrument mode for which the standard was present. A table in the appendix provides an overview of the available data for LTT 3218 in 2018.

For each of the useful observations we queried the SAF to download first OBJECT exposure taken in this mode, and in the archive request the CalSelector was set to provide us with the associated raw calibrations. This operation was carried out with the script query_standard_and_calib.py.

Data Reduction

We automated the reduction of the downloaded data sets by making use of the association tree that is delivered as an XML file by the SAF for each OBJECT exposure. In particular the sub-tree association category="STD" in those XML files contains all the needed associations for reducing a standard star exposure to a datacube. To parse the XML tree programmatically we use Python 3's Element Tree XML API. Furthermore, the calls to the MUSE pipeline (here version 2.8.3 was used) recipes are made via Python-CPL Library (Streicher and Weilbacher 2012, Streicher 2016). The script automating the reduction of the raw standard star data into datacubes ready for analysis is called reduce_automat_v2.py.

For the reduction all adjustable parameters of the pipeline are left at default values. In the final muse_scipost recipe of the reduction chain a response curve needs to be provided. As we are interested in testing the stability of relative photometry, we here did not create a response curve from the reduced standard-star observation, but we used the response curve from the static pipeline calibrations. We also resampled all cubes to the same world coordinate system 3D grid, which was created from the first reduction.

Extraction of stellar spectra

Before extracting stellar spectra we need to know the positions of the stars in the field of view. To build this catalogue we used a median image from three 50 Å band images around the central wavelength plane of the datacubes as detection image. This approach is also used to detect stars in the astrometric calibration routine of the MUSE pipeline (Sect. 3.12 in Weilbacher et al. 2020). We then apply the DAOFIND algorithm (Stetson 1982) to the detection image (using the implementation from the Python3 photutils package; Bradley et al. 2021). We opted for a threshold of 6σ and a Gaussian convolution kernel FWHM of 5px. Moreover, we masked regions around the standard star itself, as here often multiple false detections happened due to diffraction spikes.

Above procedure, implemented in the Python3 script star_cats_for_cubes.py, was applied to all reduced datacubes. This script also creates a plot that shows the positions of the detected stars overlaid onto the detection image. Using those images we selected visually the highest quality exposure with a stellar catalogue that did not contain any erroneous detections. This exposure - catalogue combination (2019-03-10 / WFM-AO-N) defines the reference catalogue and the plot for this exposure is shown below.

2019-03-10 WFM-AO-N detection image with DAOFIND detected sources indicated. This catalogue defines our reference frame and the IDs of the stars are used to identify the extracted spectra in all exposures.

Sifting through all detection images we sometimes noted significant (∼ 1'') offsets between the positions of all stars relative to the fixed pixel coordinates. This is because the absolute astrometry of a MUSE observation is determined relative to the telescope guide star, which is chosen at the liking of the telescope operator when acquiring a target, and the astrometric accuracy varies between different guide stars. This effect is demonstrated in the following plot.

Offset between coordinates of detected stars in the reference catalogue (blue points) and a catalogue from a different observation (2019-05-08 / WFM-AO-N; orange points). Note that the latter catalogue also contains a few erroneous detections. The green crosses show the transformation of the reference catalogue coordinates, obtained with the method described below, that account for the offsets relative reference catalogue. (This figure is created by the script test_coordtrafo.py.)

To account for these offsets between reference catalogue coordinates (x_r, y_r) and coordinates in the other catalogues (x_t, y_t) before extracting the spectra in each catalogue we compute a linear transformation of the coordinates (i.e., x_t = a × x_r + b × y_r + c and y_t = d × x_r + e × y_r + d) for each catalogue relative to the reference catalogue. This is achieved with the catalogue matching algorithm quad-kd of Heyl (2013). (The transformations neccesary here are typically characterised by a ≈ e ≈ 1 and b ≈ d ≈ 0, i.e. pure translations.)

Equipped with the coordinates of the stars for each observation we then extracted the spectra for each star in 5 different apertures of diameter (0.6, 1.0, 1.4, 1.8, 2.2)'' or, equivalently, (3, 5, 7, 9, 11) pixels. For each spectral extraction the variances are propagated accordingly from the STAT extension of the datacube. To illustrate the quality of the extracted spectra we show below a plot of four extracted spectra.

Extracted spectra from the reference stars 25, 6, 22, and 15 (from top left to bottm right). The x-axes are the indexes of the spectral layers and the y-axis is flux in arbitrary units. Note that spectral energy distribution of the extracted is not real, as an arbitrary response curve was used for extracting the spectra. (Plots of the individual spectra from the spectral extractions are done with the script plot_refspec.py.)

The routine performing the coordinate transformation with respect to the reference catalogue and the spectral extractions for the individual exposures is given in the Python3 script extract_spectra.py.

(For the extraction of the reference spectra the coordinate transformation is not required. For those spectra extraction is handled by the script extract_reference_spectra.py.)

Results

In the following we present a two simple tests that characterise the temporal stability of the differential photometry. We focus here on results in the I_Bessel band (taken from the ADPS; Moro and Munari 200). Our results are based on the spectra extracted in the 1.8'' diameter aperture, as this aperture provided the optimal SNR for the extracted spectra for the majority of observations. (We provide a table with the SNR measurements in I_Bessel in the appendix.)

First, we show the timeline of the stability for the first four pairs of the highest SN stars (6/4, 4/22, 22/2, and 2/15; we do not use the star 25 as it is close to the edge off the field of view, and therefore subject to different flux loss due to differences in atmospheric seeing).

Timeline visualising the deviations of the flux density ratio (F1/F2) between two stars in the IBessel band relative to the mean <F1/F2> over all displayed measurements from 2019: (F1/F2)/<F1/F_2> - 1. We show here the four pairs with the highest SNR. (The date range 2019-06-01 -- 2019-11-01 has been omitted from the plot since in this range LTT 3218 was not observed. The individual panels of the figure are created with the script spec_spat_diff_iband_test.py.)

For the four pairs with the highest SNRs standard deviations from the mean are always less then 10 % (standard deviations are 0.05, 0.04, 0.04, and 0.03 for the pair 6/4, 4/22, 22/2, and 2/15, respectively). Moreover, throughout the period of 2019 no strong systematic trends are apparent. However, on a shorter timescale (from 2019-04-01 to 2019-05-30) a gentle systematic trends could be envisioned for the 6/4 pair. Interestingly, while those stars are very close to each other, but their fluxes are recorded by different unit spectrographs (channels): star 6 falls on channel 5, while star 4 is recorded on channel 4. Nevertheless, over the whole period this trend is statistically not significant (Spearman ρ=0.125, p=0.32).

As a second test, we investigated how the standard deviation over the whole time-series depends on the SNR of the stars used. This is shown in the following figure. Since our analysis is based on pairs of stars, we plot on the x-axis the geometric mean of the two stars involved. (As above, spectra from star 25 are excluded from the analysis.)

Standard deviation of the relative flux-density variations over all 2019 observations as a function of SNR. (This figure is created with the script spec_spat_diff_sn_test.py.)

Concluding Remarks

Above results present a first test of the temporal stability of the differential photometry using standard star exposures. The main difference of the analysis here to the analysis in Giesers (2020) is, that the here used exposures do not involve any dither pattern or rotation. Thus, the observational dataset analysed here is more susceptible to differences between the response in individual channels.

We here provide all code that was used to perform the present analyses. This not only ensures reproducibility of the results presented so far, but it also enables future work to build upon our work here.

References

Bradley, L. et al. 2021, astropy/photutils: 1.1.0, Zenodo.
Giesers, B. 2020, PhD Thesis, Georg-August-Universität Göttingen
Giesers, B. et al. 2018, MNRAS 475, L15
Ginsburg A. et al. 2019, AJ 157, 98
Heyl, J.S. 2013, MNRAS 433, 935 / http://ascl.net/2107.022
MUSE-DRP: MUSE Data Reduction Pipeline, ascl:1610.004
Moro, D. and Munari, U. 2000, A&AS 147, 361 Selman, F.J. 2019, AN 341, 26
Stetson, P.B. 1982, PASP 99, 191
Streicher O. 2016, ascl:1612.001
Streicher, O. and Weilbacher P. 2012, ASPC 461, 853
Weilbacher et al. 2020, A&A 641, A28
Weilbacher, P. et al. 2015, A&A 582, A114

Appendix

Table of available data for LTT 3218 in 2018

The table below (standard_overview_table_2018.ascii in this repository) is created with the script query_standard_overview_table.py; adjust the variable eso_username if you re-use the scripts in this repository. The first column of table provides the dates for which the standard LTT 3218 was observed in at least one instrument mode and the following columns indicate whether a science target (SCI; observation for which calibration data can be retrieved with our workaround) or not (NOSCI; observation not usable with our workaround) is available.

DATE	WFM_NOAO_N	WFM_NOAO_E	WFM_AO_N	WFM_AO_E
2018-01-03	SCI
2018-01-04	SCI
2018-01-22	SCI	NOSCI
2018-01-27	NOSCI
2018-01-29			NOSCI
2018-02-12	SCI	NOSCI	SCI	NOSCI
2018-02-20	SCI	NOSCI
2018-03-08	SCI
2018-03-09	SCI	SCI	NOSCI	SCI
2018-03-10				SCI
2018-03-11	SCI	NOSCI	NOSCI	SCI
2018-03-12	NOSCI	NOSCI	SCI	SCI
2018-03-19	SCI	NOSCI	SCI	SCI
2018-03-20	SCI		SCI
2018-03-24	SCI
2018-03-31			NOSCI
2018-04-03			SCI
2018-04-07	SCI	SCI	SCI	NOSCI
2018-04-08	SCI	NOSCI	SCI	NOSCI
2018-04-09	SCI	SCI	SCI	NOSCI
2018-04-11	SCI	NOSCI	SCI	SCI
2018-04-12	SCI	NOSCI	SCI	NOSCI
2018-04-13	SCI	NOSCI	SCI	NOSCI
2018-04-15	SCI	SCI	SCI	NOSCI
2018-04-16	SCI	NOSCI	SCI	NOSCI
2018-04-17			SCI
2018-05-06	SCI	SCI	SCI	NOSCI
2018-05-10	NOSCI	SCI	SCI	NOSCI
2018-05-11	SCI	NOSCI	SCI	SCI
2018-05-12	SCI	NOSCI	SCI	SCI
2018-05-13	SCI	NOSCI	NOSCI	NOSCI
2018-05-14	SCI	SCI	SCI	SCI
2018-10-16	SCI
2018-11-28	SCI
2018-11-29	SCI		SCI
2018-12-02	SCI	SCI	SCI	SCI
2018-12-03	SCI		SCI
2018-12-04	SCI	SCI	SCI	NOSCI
2018-12-05	SCI	SCI
2018-12-06	NOSCI	NOSCI	SCI	NOSCI
2018-12-07			SCI	NOSCI
2018-12-08	SCI		SCI
2018-12-09			SCI
2018-12-10	SCI		SCI
2018-12-11	NOSCI		SCI
2018-12-12	SCI	NOSCI	SCI	NOSCI

SNR of the used stars in I_Bessel band

The following table lists the obtained I~Bessel SNR for the stars by their IDs. Spectra have been extracted in apertures of 1.8'' diameter.

id	SNR
1	369
2	1037
3	134
4	4576
5	490
6	5961
7	166
8	476
9	265
10	203
11	597
12	161
13	555
14	714
15	1031
16	197
17	539
18	483
19	372
20	238
21	116
22	1353
23	140
24	314
25	1092

The SNR values are computed with the script sn_refspecs.py. A machine readable table of the SNR values for all extraction apertures is provided in sn_refspecs_table.ascii.