SPectro-Interferometric Data Analysis Software Tool

# Spidast: a new modular software to process spectro-interferometric measurements††thanks: Applied on observations made with ESO telescopes at the Paranal Observatory under Belgian VISA Guaranteed Time programme ID 083.D-029(A/B), 084.D-0131(A/B), 086.D-0067(A/B/C)

P. Cruzalèbes, Y. Rabbia, A. Jorissen, A. Spang, S. Sacuto, E. Pasquato, A. Chiavassa, O. Chesneau and P. Fréville
Laboratoire Lagrange, UMR 7293, Université de Nice-Sophia Antipolis, CNRS, Observatoire de la Côte d’Azur, Bd de l’Observatoire, B.P. 4229,
F-06304 Nice cedex 4, France
Institut d’Astronomie et d’Astrophysique, Université Libre de Bruxelles, Campus Plaine C.P. 226, Bd du Triomphe, B-1050 Bruxelles, Belgium
Department of Astronomy and Space Physics, Uppsala Astronomical Observatory, Box 515, S-751 20 Uppsala, Sweden
Institute of Astronomy, University of Vienna, Türkenschanzstraße 17, A-1180 Vienna, Austria
Club d’Astronomie Copernic, Maison des Associations, 642 rue des Batteries, F-83600 Fréjus, France
E-mail: pierre.cruzalebes@oca.eu
Version 2013 April 04
###### Abstract

Extracting stellar fundamental parameters from SPectro-Interferometric (SPI) data requires reliable estimates of observables and with robust uncertainties (visibility, triple product, phase closure). A number of fine calibration procedures is necessary throughout the reduction process. Testing departures from centro-symmetry of brightness distributions is a useful complement. Developing a set of automatic routines, called SPIDAST (made available to the community) to reduce, calibrate and interpret raw data sets of instantaneous spectro-interferograms at the spectral channel level, we complement (and in some respects improve) the ones contained in the amdlib Data Reduction Software. Our new software SPIDAST is designed to work in an automatic mode, free from subjective choices, while being versatile enough to suit various processing strategies. SPIDAST performs the following automated operations: weighting of non-aberrant SPI data (visibility, triple product), fine spectral calibration (sub-pixel level), accurate and robust determinations of stellar diameters for calibrator sources (and their uncertainties as well), correction for the degradations of the interferometer response in visibility and triple product, calculation of the Centro-Symmetry Parameter (CSP) from the calibrated triple product, fit of parametric chromatic models on SPI observables, to extract model parameters. SPIDAST is currently applied to the scientific study of 18 cool giant and supergiant stars, observed with the VLTI/AMBER facility at medium resolution in the K band. Because part of their calibrators have no diameter in the current catalogs, SPIDAST provides new determinations of the angular diameters of all calibrators. Comparison of SPIDAST final calibrated observables with amdlib determinations shows good agreement, under good and poor seeing conditions.

###### keywords:
methods: data analysis – techniques: interferometric – stars: late-type
pagerange: SPIDAST: a new modular software to process spectro-interferometric measurementsthanks: Applied on observations made with ESO telescopes at the Paranal Observatory under Belgian VISA Guaranteed Time programme ID 083.D-029(A/B), 084.D-0131(A/B), 086.D-0067(A/B/C)Referencespubyear: tbd

## 1 Introduction

The power of optical-infrared interferometry to obtain information about the astronomical source morphology (including the angular size) is now well-established. To properly determine the source properties, the quality of the measurements is an issue, which is still the subject of active research. As with any other measuring apparatus, the absolute calibration of the instrument (including atmosphere) requires careful attention. Derived from the measurement of the mutual degree of coherence of the incident radiation field, on spatial frequencies sampled by the aperture-array configuration, the final interferometric observables are non-linear mixes of noisy quantities, and of parameter estimates with their own uncertainties.

In this paper, we propose to revisit and extend the existing data processing and calibration methods, in the aim to obtain reliable estimates and robust uncertainties for calibrated measurements of visibility and complex triple product. The careful reduction process, described in the present paper, has been elaborated for the scientific study of a sample of 18 bright cool giant and supergiant stars (see Table 1). Measurements were obtained with the VLTI/AMBER facility at medium spectral resolution (=1500) in the band, using triplets of 1.80-m auxiliary telescopes. Observations have been conducted during 15 observing nights between May 2009 and December 2010 (under Belgian VISA Guaranteed Time). The aim of the present work is to test the data reduction and calibration software, on data with various qualities, that we started to develop in 2006.

## 2 Defining the basic interferometric observables

The coherent flux is provided by the measurement of the instantaneous observables delivered by the AMBER Data Reduction Software . This quantity traces the sine-like molulated component superimposed on the continuum component, in the observed intensity distribution. It is computed using a linear fit of the individual interferograms in the detector plane, for each spectral channel (see Tatulli et al. 2007, for details).

At this point, we define the other interferometric observables:

1. the squared flux, which is the product of the photometric fluxes and , associated with the baseline

 f2ij=pipj; (1)
2. the squared visibility, which is the ratio of the squared modulus of the coherent flux to the squared flux

 v2ij=∣∣cij∣∣2f2ij; (2)
3. the complex bispectrum (Weigelt 1977), which is the product of the 3 coherent fluxes contributing to the baseline triplet

 b123=c12c23c31; (3)
4. the triple product, which is the ratio of the bispectrum to the product of the 3 fluxes contributing to the baseline triplet

 t123=b123f12f23f31. (4)
5. the closure phase, which is the argument of the complex bispectrum (or triple product)

 ψ123=argb123=argt123, (5)

## 3 Reducing the raw data

To process the temporal series of frames (composing the “exposures”) produced by the AMBER instrument, we use the data reduction software amdlib, which consists of a core library of C functions plus a high-level interface in the form of a Yorick plugin (a high-level language environment, similar to IDL or Matlab, which is in the public domain).

The standard procedure to extract the basic interferometric observables from the AMBER raw data includes a frame-selection step, which provides a set of “good” frames according to a given criterion, or a sequence of several selection criteria (Malbet et al. 2011). The criterion used the most is based on the Signal to Noise Ratio (SNR) of the fringe-contrast in each frame. Two selection methods are provided in standard by amdlib:

1. retaining a given percentage of frames sorted according to their SNR values (called the “percentage” method),

2. retaining the frames with SNR values higher than a given threshold (called the “threshold” method).

With both methods, the choice of the optimal value (assumed not biasing the final results) of the selection criterion (in percentage of threshold) is done a posteriori, from a sequence of processings using different values of the selection criterion (see Millour et al. 2007, for details). Two drawbacks of this are: the need for a several-steps procedure (to be performed manually), and the risk for biasing the final results if the selection criterion is ill-determined.

In our method, that we have implemented in a specific Yorick script added to amdlib, we apply the automatic procedure:

1. remove only the aberrant measurements (in squared visibility and triple product), for each spectral channel independently. Using this specific script results in amounts of useful data larger than the amount obtained with the standard procedure, without increasing biasing effect (only the lower and upper tails of the input data histograms are rejected);

2. assign a weight to each individual measurement. After many preliminary tests, we decided to choose weights based on the SNR of the coherent flux, and on the piston deviation;

3. compute the temporal weighted average of the squared visibility and triple product for each spectral channel, with uncertainties derived from the bootstrap technique invented by Efron (1979).

These automatic operations mentioned above, needing neither several-step selection nor manual choice of the optimal value of the selection criterion, are described in more details in the following subsections.

### 3.1 Rejecting the aberrant data

To remove the aberrant measurements, we use a combination of physical and statistical basic criteria:

1. strictly positive sum and product of the fluxes in the photometric channels, since negative flux measurements result from shortcomings in the determination of the continuum in each spectrum;

2. coherent-flux SNR higher than unity, since inclusion of data with low SNR values reduces the reliability of the estimators of the interferometric observables (Millour et al. 2008);

3. squared visibility greater than and lower than (“box-plot filtering”), where denotes the first quartile, the third quartile, and the interquartile range (Hoaglin et al. 1983). With triple product data, this third criterion applies to the squared visibility for each baseline, and box-plot filtering on is added to the list of criteria.

Based on the removal of aberrant measurements at the spectral channel level, rather than on the selection of the “best” frames, our approach provides an amount of useful data larger than the amount obtained with the standard frame-selection procedure. For data shown in Fig. 1, obtained for a “good” seeing on the calibrator target  Lib, the number of removed input data varies only slightly from one spectral channel to the other (typically less than 3%). Besides, various amount of data are rejected, according to the seeing conditions.

Using the seeing parameters (seeing angle) and (coherence time), the percentages of rejected data, respective to three baselines and two seeing conditions are given in the following lines (where “good seeing” refers to =0.5 , =13 , and “poor seeing” refers to =1.7 , =1.3 ):

• shortest baseline A0–D0 (32 ): 2% (good seeing), to 9% (poor seeing);

• intermediate baseline D0–H0 (64 ): 2% to 16%;

• longest baseline A0–H0 (96 ): 2% to 25%;

• in addition, for the baseline triplet A0–D0–H0, the rejection of data amounts 7% and 50%.

Figure 2 shows the variation of the final squared visibility (at 2.2 μ) produced with amdlib, w.r.t. the percentage of selected frames, based on the SNR and pertaining to the calibrator target  Lib, observed with the shortest baseline A0–D0 (32 ) in good and poor seeing conditions (top panel and bottom panel respectively). The squared visibility produced with our method is shown for comparison.

In this example, the position of the plateau of each curve indicates a value of the selection threshold of 50%, giving the a posteriori determination with amdlib. There we find that our automatic method gives a squared-visibility value close to the 50% value, without a need for the several-steps procedure, in good and poor seeing conditions. Besides, Fig. 3 shows that the standard amdlib frame-selection method keeps the extremal V values in poor seeing condition, which may bias the result of the calculation of the average over the frames contained in each exposure. On the contrary, the box-plot filtering used with our method rejects only the lower and upper tails of the V histogram, thus leading to estimates of the temporal average, more reliable than with the frame-selection method.

Next subsection describes the weighting scheme used to average the non-aberrant data over each exposure. Note that selecting the frames with the standard amdlib procedure is equivalent to assigning unity weights to the selected frames, and null weights to the others.

### 3.2 Computing the raw observables

To reduce the instrumental effects on the data, at sthe pectral channel level, we compute in each exposure (300 frames in our example) the weighted average of the squared visibility, and of the complex triple product. For the visibility, as weight associated with each frame we use the ratio of the SNR to the relative excursion of the piston, where the relative excursion of the piston itself is the ratio of the piston excursion (absolute value) for a given frame, to the average over the whole exposure. For the triple product, each individual weight is given by the geometrical mean of the three weights associated with the three baselines. The final uncertainties are derived from the variance of the distribution of the weighted means, obtained by random sampling with replacement of the original series of data (direct-bootstrap method) (Efron & Tibshirani 1993). Note that our specific script computes the real and imaginary parts of the triple product, from which we extract the calibrated closure phase (see Sect. 4.2).

Figure 4 shows an example of squared visibility and triple product produced by the reduction script (before calibration) from one single exposure of 300 input frames, obtained with the calibrator target  Lib, for the three baselines A0–D0 (32 ), D0–H0 (64 ), and A0–H0 (96 ).

In complement to these interferometric observables, we compute two other quantities used further in the calibration procedure (see Sect. 4.2.3):

1. the weight of each exposure, defined as the product of two ratios: one is the average of the SNR to its variance; the other is the average of the inverse of the piston excursion (absolute value) to its variance;

2. the final number of non-aberrant data, for each spectral channel.

## 4 Calibrating the data

The calibration process is the key point to obtain accurate estimates of the “true” (intrinsic) observables of the science targets. To correct the measurements for environmental and instrumental instabilities, we observe reference stars (the calibrators), with known angular diameters and independently well-defined brightness distributions.

To calibrate the wavelength-dependent measurements obtained with VLTI/AMBER, we have developed a library of IDL functions, included in the SPIDAST (SPectro-Interferometric Data Analysis Software Tool) modular software suite (Cruzalèbes et al. 2008, 2010), which allows one to:

1. link each spectral channel to a wavelength value. We use a specific method based on the correlation of calibrator measured spectra with synthetic templates given by the MARCS model (Gustafsson et al. 2008);

2. measure and correct for the degradation of the spatial coherence. As calibrators, we use stars with brightness distribution of Limb-Darkened (LD) discs, for which angular diameters are given by fits of synthetic spectra on the wide-band spectrophotometric measurements in the infrared (Cruzalèbes et al. 2010).

### 4.1 Computing the spectral shifts

To reach a high precision in the angular diameter estimation using model fitting, spectral calibration is a critical point (see e.g., Wittkowski et al. 2008; Domiciano de Souza et al. 2008; Stefl et al. 2011). Given the lack of any internal instrumental module for wavelength calibration in the optical setup of AMBER (Robbe-Dubois et al. 2007; Petrov et al. 2007), amdlib provides calibrated wavelength tables, computed from a theoretical polynomial dispersion law (Mérand et al. 2010), but with coefficients still badly known, leading to wavelength shifts up to 10  in the detector plane, at medium spectral resolution (Malbet et al. 2011). To correct for this drawback, amdlib shifts the wavelength table using the correlation of the measured spectrum with a template table containing the telluric lines.

To improve the precision of the wavelength calibration, observers usually compute the coefficients of the dispersion law by identifying visually some prominent telluric lines in the measured spectrum (e.g., Kraus et al. 2009; Weigelt et al. 2011; Ohnaka et al. 2011). We propose an alternative automatic approach, based on the cross-correlation of the calibrator measured spectra with synthetic templates produced by the MARCS + Turbospectrum codes (Alvarez & Plez 1998; Plez 2012), including the telluric lines as well. This method was previously suggested by Plez (2003) for Gaia, and by Decin et al. (2004) or Decin & Eriksson (2007) to calibrate the SPITZER spectrograms. It is particularly suitable for stellar spectra showing easily identifiable spectral features, as it is the case with our sample of cool giant calibrators, showing strong CO bands in the 2.126–2.474 μ spectral range (Martí-Vidal et al. 2011). Figure 5 shows the synthetic spectrum produced by the MARCS + Turbospectrum code, for the calibrator  Lib. Some reference spectral lines are shown for identification of the corresponding stellar atmospheric elements.

For each exposure, obtained with each calibrator, and with the working instrumental setting, our automatic spectral calibration process contains the following steps:

1. apply the heliocentric and systemic radial velocity corrections, and multiply the synthetic spectrum by the atmospheric and instrumental transmittance profiles. Atmospheric transmission data for the southern sites (CTIO, Chile) are given by the USAF atmospheric code PLEXUS (Cohen, Wheaton & Megeath 2003);

2. remove the continuum parts of the raw and synthetic spectra, and normalize the resulting spectra. As a first approximation, we estimate the continuum by decreasing the spectral resolution to =40. The main drawback of this method is to produce apparent “pseudo-continua” that are lower than the real continuum levels, even in almost line-free regions (Rix et al. 2004). Since our goal is only to perform spectral calibration by correlation, this drawback has no effect on the final result;

3. divide the spectrum in contiguous sub-windows of the same size, and compute the wavelength shift by phase-correlation (Vera & Torres 2008) in each sub-window. If the spectrum is divided in sub-windows, we apply a polynomial law of degree to correct for the wavelength shifts. Sub-pixel precision is obtained by embedding the cross-power spectrum in the middle of a two-times larger array of zeros, before computing the phase-correlation function, by inverse Fourier Transform of the final cross-power spectrum (Guizar-Sicairos, Thurman & Fienup 2008);

4. finds the position of the barycenter of the correlation peak for each sub-window, computes the coefficients of the interpolation law (assumed to be polynomial), and corrects the global wavelength table;

5. calculates the residuals of the wavelength shifts, using once again the phase-correlation method on the wavelength-corrected spectrum.

For each observing night, the final wavelength table, associated with the working instrumental setting, is obtained by ensemble average of the corrected tables (of all exposures with all calibrators measured with the same instrumental setting), rejecting the tables showing wavelength-shift residuals higher than 8  (5 times the nominal spectral resolution). Figure 6 shows the final continuum-corrected spectrum, given by the spectral-calibration procedure, with the calibrator  Lib. The model spectrum drawn with the dashed line is produced by the MARCS + Turbospectrum code.

To show the difference between the wavelength-calibration provided by SPIDAST and by amdlib, we plot in Fig. 7 the corrections from the initial polynomial dispersion law of AMBER, given by the two softwares. While amdlib computes only a unique correction value applied over the whole spectrum, our method computes a polynomial correction (here parabolic from three sub-windows), with sub-pixel deviations with respect to the polynomial law, caused by the final ensemble averaging of the corrected wavelength tables. The choice of the degree of the polynomial correction provided by SPIDAST allows the accuracy of the wavelength correction to be adapted to the scientific programme, reaching a sub-pixel precision if necessary.

### 4.2 Calibrating spectro-interferometric data

In order to assess and correct for the measurement defects on the science targets, proper interferometric calibration (not yet fully supported with amdlib) is needed. The Instrument Response Function (IRF), also called “transfer function” (Perrin 2003; van Belle & van Belle 2005; Boden 2007; Cruzalèbes et al. 2010), is derived from observations of calibrators as concomitant as possible to the ones of the science targets.

The standard calibration method is based on the assumption that a given “raw” (measured) quantity of any observable , is equal to the product of the “true” (intrinsic) quantity , multiplied by a global degradation factor (in other words the IRF). Such a factor cannot practically be modellised with sufficient reliability (numerous parameters and some of them remaining unknown).

To calibrate the quantity related to the science target, a second assumption applies, as soon as appropriate conditions (described later on) are satisfied and can be then determined by using the relation

 qsciQsci=qcalQcal=RQ, (6)

with self explanatory notations, stating that the degradation factor is identical for both the calibration source and the science target. This assertion is all the more valid that the required conditions are fullfilled. The relation here above yields an empirical determination of , via , relying on commonplace statistical tools. The main condition to satisfy is that the IRF is stable enough, and the standard approach to avoid too large variations of the IRF is to perform sequences, where observation of calibrator and science targets are interlaced and repeated several times, rapidly sampling the IRF. So, the actual stability requirement is that the IRF might be slowly variable over a given “calibrator-science-calibrator” observation sequence. Besides, the calibrator and the science target should be close enough regarding time of observation, angular separation, and brightness (in the spectral domain of work), so that adjustments of the whole interferometer would not be significantly modified from one source to the other. In these conditions, only the atmospheric turbulence might cause substantial degradations of the measurements. For this latter case, servo-loop systems are to be used to reduce turbulence effects.

Usually, calibrators are found in dedicated catalogs and the corresponding sources should meet the desired requirements. However, the number of calibrators is necessarily limited and somewhat depends on the type of science targets, so that the selected calibrators might stand out from the ideal conditions: point-like source, small angular separation, brightness matching.

In practice, the limited choice of calibrators makes it necessary to accept “partially” resolved targets (van Belle & van Belle 2005), angular separations counted in degrees, and brightness mismatch amounting some units of magnitudes. For example, a bright calibrator is rarely found in a close angular vicinity of a given bright science target. Such a situation nevertheless remains convenient for the calibration procedure, since switching between science targets and calibrators is fast enough and does not affect the configuration of the interferometer, and its response as well.

This procedure starts with a first selection of calibrators via the SearchCal tool, created by the JMMC working group “Catalogue of calibration sources”, providing rough estimates for angular diameters, with comparatively large uncertainties. Following this first selection, the sources are searched in the calibrator catalogues of Bordé et al. (2002) and Mérand et al. (2005), in order to find angular diameters of better quality. Since some calibrators have not a diameter in the published catalogues, it is necessary to make our own determinations of those diameters (see Subsect. 4.2.1). Moreover, in order to control these determinations and to build a homogeneous set of diameters, we apply our own method to recalculate all the diameters for the selected calibrators. The good agreement found between our determinations and the ones of the Bordé’s Catalogue of Calibrator Stars for Long-Baseline Stellar Interferometry (CCSL) (as shown in Table 2, and in Fig. 8), attests to the reliability of our specific determinations of angular diameters.

#### 4.2.1 Determining calibrator angular diameters

To derive the IRF for each spectro-interferometric (SPI) observable, we firstly need to get a reliable estimate of the angular diameter of each calibrator, and its associated error. Fitting the stellar models of Kurucz (1979) on the absolutely calibrated spectro-photometric templates of Cohen et al. (1999), the CCSL contains 374 calibrators, with limb-darkened angular diameters ranging from 1 to 3 . Since our selected calibrators might have no angular-diameter estimate (one third of the calibrators of our dataset have no value found in the catalogues), or poorly estimated, we have added to SPIDAST various routines for determination of the angular diameter from indirect methods, presented in detail in Cruzalèbes et al. (2010). These routines compute the angular diameter:

1. by combining the linear diameter, deduced from absolute luminosity and effective temperature from the Morgan-Keenan (MK) spectral type, with the parallax;

2. by using experimental laws based on the interstellar-corrected colour index (surface brightness method);

3. by scaling synthetic spectra on broadband photometric measurements (infrared flux method);

4. by fitting synthetic spectra on infrared spectro-photometric measurements (see Sect. 6).

Table 2 gives the angular diameter values of our calibrators found in the CCSL, as well as our final estimates, derived from the fit of MARCS + Turbospectrum synthetic spectra on photometric or spectro-photometric measurements. In the last column, we also give the estimates found in the JMMC Stellar Diameter Catalogue (JSDC) of Lafrasse et al. (2010), with accuracies between 7 and 10% (Delfosse 2004; Bonneau et al. 2006).

Figure 8 shows the comparison of our estimates with the CCSL values. The mean difference between our estimates and the CCSL is 5%. The mean difference for the three targets with the available calibrated templates of Cohen (i.e.  Ori,  Scl, and  Cet) is lower: 3.5%. This shows the satisfactory agreement between our results and the CCSL, which confirms the reliability of our approach for determination of the calibrator angular diameters, implemented in SPIDAST.

#### 4.2.2 Modelling the calibrator visibility and triple product

Once we determined the angular diameter of each calibrator, and its associated uncertainty, we derive the IRF for each SPI observable, computing the ratio of the observable measured on the calibrator, to the synthetic observable, assumed to represent the “true” calibrator observable. Each synthetic SPI observable is derived from the synthetic coherent flux

1. the synthetic flux (Eq. 1) is given by the synthetic coherent flux computed for the null baseline (visibility unity);

2. the squared synthetic visibility (Eq. 2) by the ratio of the squared synthetic coherent flux to the synthetic squared flux;

3. the synthetic bispectrum (Eq. 3) by the product of the synthetic coherent flux for the three baselines composing each triplet;

4. and the synthetic triple product (Eq. 4) by the ratio of the synthetic bispectrum to the cubed synthetic flux.

Applying the van-Cittert-Zernike theorem, the model coherent flux of each calibrator, assumed to emit a luminous intensity with circular symmetry, is derived from the Hankel transform of the source intrinsic spectral radiance (in )

 Cmodelλ=πϕ24∫∞0Lλ(r)J0(πrϕBijλ)rdr, (7)

where is the stellar angular diameter, the impact parameter, i.e. the distance from the centre of the stellar disc ( at disc center), the length of the projected baseline of the (i,j)-pair of apertures, and J the Bessel function of order zero.

Since compact photospheres are known to deviate from simple uniform discs (see Hajian et al. 1998, and references therein), radially decreasing from their photometric centre, we use the MARCS + Turbospectrum codes to produce models of Centre-to-Limb Variation (CLV) profiles, with input parameters derived from the MK spectral type. Figure 9 shows the map of typical CLV profiles given by the model, where the reduced radial coordinate is given by = (Young 2003).

With “partially” resolved calibrators, the deviation from the UD model increases with the baseline length. Figure 10 shows the squared visibility and the real part of the triple product produced by the MARCS and the UD models, for the calibrator  Ceti (=3.320.02 ), deduced from measurements obtained with the D0–I1–G1 baseline triplet (G1–I1=46 , D0–G1=69 , and D0–I1=79 ). In this example, the differences in visibility and triple product between the UD model and the LD-MARCS model are larger than the uncertainty of the synthetic observables (derived by propagating the uncertainty on the angular diameter), hence the justification to use the LD model instead of the UD model.

#### 4.2.3 Correcting for the degradations of the IRF

The main source of non-stationary instabilities which affects the fringe formation is the atmospheric-phase turbulence (Roddier 1981). This issue is partially solved thanks to the use of fringe-tracking servo-loops scanning the fringes with cycle periods smaller than the seeing coherence time, defined as the time over which the phase fluctuations remain coherent (Davis & Tango 1996; Kellerer & Tokovinin 2007). Thus, the residual instabilities are expected to cause only slow drifts of the IRF between observations of the science targets and their associated calibrators. In order to estimate the response at the time when each science target is measured, linear interpolation between the successive measurements on its calibrator(s) is legitimate (Perrin 2003), provided that the “calibrator-science-calibrator” observation sequence, repeated several times, rapidly samples the evaluation of the IRF (Berger, Dumas & Kaüfer 2011).

To perform the calibration of the squared visibility and triple product data, SPIDAST applies the following procedure, working at spectral channel level:

1. Compute the measured IRF for each exposure, given by the ratio of the calibrator measured observable to the calibrator synthetic observable (at the same spectral resolution, for the same baseline or baseline triplet). Values and uncertainties are derived from the 4th-order approximation formulae of Winzer (2000). For triple product data, the approximation formulae apply separately to the real and imaginary parts of the complex ratio. Since the calibrator angular-diameter uncertainties are smaller than 10% (see Table 2), the uncertainty of a model observable is simply obtained from its first partial derivative w.r.t. the angular diameter, multiplied by the angular diameter uncertainty.

2. Average the IRF measurements over each OB. As weight associated with each exposure, we use the ratio of the associated weight (defined in Sect. 3.2) to the variance of the calibrator “raw” observable.

3. Determine the IRF at the mean time of each exposure obtained with the science target, from linear interpolation, or polynomial fit, on the averaged IRF measurements. Figure 11 shows the temporal variation of the measured and interpolated response in squared visibility and triple product, obtained during one observing night, with three successive calibrators.

4. Divide the observable measured on the science target for each exposure of a given OB, by each interpolated value of the IRF obtained for the same OB (using Winzer’s formulae). Applying a method similar to that used with the getCal tool of the Cal. Inst. Tech. (2008), a weight is associated with each calibrated ratio, which combines information on the observations of the calibrator and of the science target (science–calibrator angular separation and observation-time delay, coherence time, number of non-aberrant data).

5. Compute the weighted average over each OB of the ratios used in the calibration, which gives the final calibrated SPI observable (see Subsect. 4.2.4). The argument of the final calibrated triple product gives the final calibrated closure phase.

Although the calibration procedure used by SPIDAST is based on the standard calibration method, great care has been taken to provide reliable uncertainty estimates, thanks to the use of weights tracing the data quality at different levels of the processing. Thus, SPIDAST produces large uncertainties on final calibrated observables obtained with input data of poor quality. Estimating reliable error bars is of crucial importance to get calibrated data, usable for the fitting process of parametric models, even on data of uneven quality (see Sect. 6).

#### 4.2.4 Comparing the final calibrated observables with amdlib

Although the calibration routine provided with the amdlib package is not yet officially validated (Malbet et al. 2011), we compute the differences in squared visibility and closure phase between SPIDAST and amdlib, and compare them with the associated uncertainties.

Figure 12 shows the results obtained on the science target  Oph, with the baseline triplet A0–D0–H0, under good and poor seeing conditions. The histograms of the uncertainties show that the uncertainties in squared visibility estimated with amdlib are smaller than with SPIDAST (whatever the seeing), while it is the contrary for the uncertainties in closure phase.

To explain the disagreement between the uncertainties in squared visibility, we remind that SPIDAST takes into account instrumental and environmental conditions in the calculation of the uncertainties, while amdlib uses only unweighted variance. Thus, we consider as underestimated the uncertainties in squared visibility produced by amdlib.

Regarding the closure phase, it is more difficult to explain why the uncertainties from amdlib are larger than the ones from SPIDAST. Indeed, no details are found in Malbet et al. (2011), expected to describe the amdlib calculation. SPIDAST computes the uncertainties on the real and imaginary parts of the calibrated triple product, and derives the uncertainties on the final closure phase from them. A deep analysis based on the source code of the calibration routine provided by amdlib reveals that amdlib computes the uncertainties in the final closure phase, adding the uncertainties in closure phase of the science target and of the estimated transfer function. Since the closure phase is less stable than the triple product, the final uncertainties in closure phase handled with amdlib are larger than with SPIDAST.

To evaluate the agreement between the two procedures, we compare the differences between the final calibrated SPI observables they produce, to the associated uncertainties. Since the median differences are smaller than the median uncertainties (Fig. 12), we conclude that the agreement between the two softwares is satisfactory, SPIDAST giving more reliable uncertainties than amdlib.

## 5 Measuring the deviation from circular symmetry

The closure phase couples the phases of the Fourier transform of the source brightness distribution at the three spatial frequencies , , and , probed by the three baselines in the following way: . In the case of a centro-symmetrical brightness distribution on-axis, the phase of the Fourier transform is uniformly zero and is naturally equal to zero. If this source is off-axis, the phase is linear w.r.t. the spatial frequency, so that , and, subsequently, also here is equal to zero. Thus, a non-zero phase closure is a clear indication, at least qualitatively, for a deviation from circular symmetry.

In Sect. 5.2, we introduce a new parameter (called the centro-symmetry parameter), based on the triple product, more sensitive to deviation from centro-symmetry than the closure phase.

### 5.1 Global closure phase

The integration of the real and imaginary parts of the triple product , over the observation spectral band , leads to the global closure phase, defined by Ragland et al. (2006) and Tatebe et al. (2006) as

 ψband=atan ∫λmaxλminITtrue(λ)dλ∫λmaxλminRTtrue(λ)dλ. (8)

A value of the global closure phase close to zero (less than 1 , for our sample) suggests a high degree of centro-symmetry of the brightness distribution, related to the observation spectral band. We compute the uncertainty on the global closure phase thanks to the direct-bootstrap method, using random sampling with replacement of the spectral dataset of the calibrated triple product. The values of the global closure phases in , and their associated uncertainties, are given in the third column of Table 3, for each science target.

### 5.2 Centro-Symmetry Parameter (CSP)

When the positive and negative values of along the spectral domain almost mutually compensate (as seen on Fig. 13 for TX Psc and W Ori), the spectral integration produces a nearly null global closure phase. As mentioned above, this can be taken as a hint for centro-symmetry. However, is clearly non-null in some parts of the spectrum, and this is a hint for a deviation from centro-symmetry. To rule out the contradiction, and so as to consider this latter possibility, we introduce a new estimator, that we call the Centro-Symmetry Parameter (CSP), similar to the global closure phase, but using instead the absolute value of in the numerator, and the modulus of in the denominator,

 CSP=asin ∫λmaxλmin|ITtrue(λ)|dλ∫λmaxλmin|Ttrue(λ)|dλ. (9)

As for the global closure phase, a small CSP value (less than 2 , for our sample) suggests a high degree of centro-symmetry. Significantly high CSP values, as judged from their uncertainties (bootstrap method again), require to use asymmetric brightness distribution models, for the fitting on the SPI data. The values of the global CSP in the -band, and their uncertainties, are given in the second column of Table 3.

To illustrate the difference between the global closure phase and the CSP, we plot in Fig. 13 three quantities, for 3 scientific targets showing a global closure phase close to zero,

1. in the panels on the left is shown the imaginary part of the triple product , divided by the spectral mean444defined as the integral w.r.t. the wavelength, divided by the spectral bandwidth of the real part . Note that the spectral mean of this ratio is equal to the tangent of the global closure phase, as defined in Eq. (8);

2. in the panels forming the central column is shown the absolute value of the imaginary part of the triple product , divided by the spectral mean of the modulus . Note that the spectral mean of this ratio is equal to the sine of the global CSP, as defined in Eq. (9);

3. in the panels on the right is shown the absolute value of the real part of the triple product , divided by the spectral mean of the modulus .

The three panels on the left show three typical behaviours for  :

1. uniformly close to zero ( Oph);

2. decreasing symmetrically around zero (TX Psc);

3. increasing symmetrically around zero (W Ori).

For each of these situations, the integration over the spectrum produces a nearly null result (hint for centro-symmetry). However the CSP allows to detect deviations from centro-symmetry (see Fig. 14).

If the real and imaginary parts of the triple product are wavelength-independent, for the observation spectral range, one can show that . If not, no analytical relation linking the two quantities can be derived from Eqs. (8) and (9), because of the integrals in the numerators and the denominators. Figure 14 shows the CSP in displayed versus the global closure phase in the same spectral band, for each science target, choosing the OB which gives the smallest uncertainty on the CSP. Almost all the stars fall along the diagonals of the (, CSP) diagram, which trace the relation CSP = . One (W Ori), however, is flagged as asymmetrical with the CSP indicator, but not with . This is precisely why we favour the CSP over , as discussed above. T Cet has the largest values for CSP and as well, indicating large deviation from centro-symmetry (may be due to strong asymmetries at the wavelengths of the CO bands or even suggestive of the presence of a binary companion). The second largest CSP value pertains to R Scl, recently revealed as a wide binary by Maercker et al. (2012), using the ALMA array at millimetric wavelengths. In addition, the last column of Table 3 shows that, except for T Cet, the SNR for the CSP is comparable or larger than the SNR for the global closure phase.

## 6 Fitting parametric chromatic models

In order to interpret the final calibrated SPI observables, SPIDAST provides a fitting routine, based on the modified gradient-expansion algorithm, very similar to the algorithm invented by Levenberg (1944), and improved by Marquardt (1963). The fit applies on any calibrated SPI measurements (to be chosen between visibility, flux, coherent flux, bispectrum, triple product, and closure phase), using a library of single-component or composite parametric chromatic models, characterised by the Fourier spectrum of their intensity distribution, and the associated first-order partial derivatives w.r.t. the model parameters. Figure 15 shows examples of visibility-fitting results obtained with SPIDAST, using two models (UD and LD-MARCS). This fitting routine, applied on spectro-photometric data, is used to determine the calibrator angular diameters from the MARCS + Turbospectrum synthetic spectra (see Subsect. 4.2.1).

Note that the LITpro software, developed by the JMMC working group “Model fitting”, uses a set of elementary geometrical and center-to-limb darkening functions, as well (Tallon-Bosc et al. 2008). However, it does not offer, contrary to SPIDAST, the possibility to fit stellar-disk models with synthetic tabulated radiance data (or exitance, for fits on spectro-photometric measurements), such those produced by the MARCS + Turbospectrum code.

Since the uncertainties of the final calibrated data are not normally distributed, the covariance matrix, that comes out of the chi-squared fit, cannot be used to infer the parameter uncertainties (Press et al. 2007). To compute the uncertainties, SPIDAST uses the residual-bootstrap method, described in detail in Cruzalèbes et al. (2010):

1. “synthetic” datasets are produced from random sampling with replacement of the Pearson residuals (difference between calibrated and fitted values divided by the uncertainty of the observed value), added to the initial fitted values;

2. fitting the model on these new datasets produces a set of chi-squared minima, following a probability distribution, from which we extract the boundaries of the confidence interval, with a given confidence level;

3. the parameter values associated with these chi-squared boundaries give the upper and lower limits of the parameter estimates, leading to asymmetric uncertainties.

The whole procedure allows the computing of reliable estimates of the parameter uncertainties.

## 7 Conclusion

In the present paper, we introduce the new SPIDAST software, developed since 2006 with the aim to reduce, calibrate, and interpret the visibility and triple product measurements obtained with the VLTI/AMBER facility. SPIDAST contains a whole set of modules, which can be launched separately or in an automatic batch file, and summarised herafter:

1. The raw data reduction used by SPIDAST computes the weighted average of non-aberrant data, at spectral-channel level, providing estimates of the SPI observables using an automated procedure, while the method presently used with amdlib selects the “best” frames according to a quality threshold determined a posteriori, after several trials.

2. The wavelength calibration procedure performed by SPIDAST provides spectral shifts following a polynomial law, tracing them at the channel level. This is done by computing the cross-correlation of the measured spectra of the calibrators with their synthetic spectra produced by the MARCS model, while amdlib only provides a constant spectral shift over the spectrum, from the correlation with the telluric lines.

3. For selected calibrators not included in the calibrator catalogues, SPIDAST provides several routines for estimation of the angular diameter with indirect methods, the most accurate being the fit of stellar-atmosphere model spectra given by MARCS on spectro-photometric data. The calibrator synthetic observables are derived using CLV functions produced by MARCS.

4. To obtain an accurate interferometric calibration (via an automated procedure), SPIDAST: (1) divides the calibrator raw data with the associated synthetic observables in each spectral channel, which gives the instrumental response function in squared visibility and triple product; (2) interpolates or fits the response at the time of each exposure on the science target; (3) divides the science raw data with the interpolated/fitted response, which gives the science calibrated observables for each exposure; (4) computes their weighted average over each OB. At each processing step, the uncertainties are computed thanks to the bootstrap method applied on the weighted means.

5. Using the real and imaginary parts of the calibrated triple product, SPIDAST measures the deviation from centro-symmetry of the brightness distribution of each science target in the observation spectral band, thanks to a new parameter, more sensitive to asymmetries than the global closure phase.

6. Finally, SPIDAST proposes a complete fitting tool, using a set of parametric and chromatic models, and accepting input tables of flux/intensity synthetic data.

Such a careful calibration process of SPI data is a crucial step for their trustworthy astrophysical exploitation, which is the topic of two associated papers (Cruzalèbes et al. 2013b, c). Parameter extraction using non-linear fits of source model, as well as aperture-synthesis image reconstruction, need reliable estimates of the calibrated observables, with robust uncertainties. Our reduction, calibration, and fitting routines also apply to any other spectral datasets, including spectroscopic data. We made the SPIDAST software public: the source code of any program of our software suite can be obtained by sending an e-mail to the first author of the present paper.

## Acknowledgments

The authors thank the ESO-Paranal VLTI team for supporting their AMBER observations, especially the night astronomers A. Mérand, G. Montagnier, F. Patru, J.-B. Le Bouquin, S. Rengaswamy, and W.J. de Wit, the VLTI group coordinator S. Brillant, and the telescope and instrument operators A. and J. Cortes, A. Pino, C. Herrera, D. Castex, S. Cerda, and C. Cid. They also thank the Programme National de Physique Stellaire (PNPS) for supporting part of this collaborative research. A.J. is grateful to B. Plez, K. Eriksson, and T. Masseron for their ongoing support on the use of the MARCS code. S.S. was (partly) supported by the Austrian Science Fund through FWF project P19503-N16; A.C. is post-doctoral fellow from F.R.S.-FNRS (Belgium; grant 2.4513.11); E.P. is supported by PRODEX. The data were partially reduced using the publicly available data reduction software package amdlib, kindly provided by the Jean-Marie Mariotti Center. This study used the SIMBAD and VIZIER databases at the CDS, Strasbourg (France), and NASAs ADS bibliographic services.

## References

• Alvarez & Plez (1998) Alvarez R., Plez B., 1998, A&A, 330, 1109
• Berger et al. (2011) Berger J.-P., Dumas C., Kaüfer A., 2011, VLT Paranal Science Operations AMBER User Manual. ESO, Garching b. München, Germany
• Boden (2007) Boden A. F., 2007, NewA, 51, 617
• Bonneau et al. (2006) Bonneau D., Clausse J.-M., Delfosse X., et al., 2006, A&A, 456, 789
• Bordé et al. (2002) Bordé P., Coudé du Foresto V., Chagnon G., Perrin G., 2002, A&A, 393, 183
• Cal. Inst. Tech. (2008) Cal. Inst. Tech. 2008, The getCal Interferometric Observation Planning Tool Suite. NASA Exoplanet Science Inst., 2.10.7 edn
• Cohen et al. (1999) Cohen M., Walker R. G., Carter B., et al., 1999, AJ, 117, 1864
• Cohen et al. (2003) Cohen M., Wheaton W. A., Megeath S. T., 2003, AJ, 126, 1090
• Cruzalèbes et al. (2013b) Cruzalèbes P., Jorissen A., Rabbia Y., et al., 2013b, in preparation
• Cruzalèbes et al. (2013c) Cruzalèbes P., Jorissen A., Rabbia Y., et al., 2013c, in preparation
• Cruzalèbes et al. (2010) Cruzalèbes P., Jorissen A., Sacuto S., Bonneau D., 2010, A&A, 515, A6
• Cruzalèbes et al. (2008) Cruzalèbes P., Spang A., Sacuto S., 2008, in Kaufer A., Kerber F., eds, The 2007 ESO Instrument Calibration Workshop Calibration of AMBER Visibilities at Low Spectral Resolution. Springer, pp 479–482
• Davis & Tango (1996) Davis J., Tango W., 1996, PASP, 108, 456
• Decin & Eriksson (2007) Decin L., Eriksson K., 2007, A&A, 472, 1041
• Decin et al. (2004) Decin L., Morris P. W., Appleton P. N., et al., 2004, ApJSS, 154, 408
• Delfosse (2004) Delfosse X., 2004, Technical report, Stellar diameter estimation from photospheric indices. JMMC
• Domiciano de Souza et al. (2008) Domiciano de Souza A., Bendjoya P., Vakili F., Millour F., Petrov R. G., 2008, A&A, 489, L5
• Efron (1979) Efron B., 1979, Ann. Statist., 7, 1
• Efron & Tibshirani (1993) Efron B., Tibshirani R., 1993, An Introduction to the Bootstrap, CRC Press edn. Vol. 57 of Monographs on Stat. and Appl. Prob., Chapman and Hall
• Guizar-Sicairos et al. (2008) Guizar-Sicairos M., Thurman S. T., Fienup J. R., 2008, Opt. Lett., 33, 156
• Gustafsson et al. (2008) Gustafsson B., Edvardsson B., Eriksson K., et al., 2008, A&A, 486, 951
• Hajian et al. (1998) Hajian A. R., Armstrong J. T., Hummel C. A., et al., 1998, ApJL, 496, 484
• Hoaglin et al. (1983) Hoaglin D. C., Mosteller F., Tukey J. W., 1983, Understanding Robust and Exploratory Data Analysis
• Kellerer & Tokovinin (2007) Kellerer A., Tokovinin A., 2007, A&A, 461, 775
• Kraus et al. (2009) Kraus S., Weigelt G., Balega Y. Y., et al., 2009, A&A, 497, 195
• Kurucz (1979) Kurucz R. L., 1979, ApJSS, 40, 1
• Lafrasse et al. (2010) Lafrasse S., Mella G., Bonneau D., et al., 2010, in W. C. Danchi et al. ed., Optical and infrared interferometry II Vol. 7734 of SPIE Conf. Ser., Building the ’JMMC Stellar Diameters Catalog’ using SearchCal. San Diego, CA
• Levenberg (1944) Levenberg K., 1944, Q. Appl. Math., 2, 164
• Maercker et al. (2012) Maercker M., Mohamed S., Vlemmings W. H. T., et al., 2012, Nature, 490, 232
• Malbet et al. (2011) Malbet F., Millour F., Absil O., Duvert G., 2011, User Manual - Issue 3.0.3 JMMC-MAN-2720-0001, AMBER Data Reduction Software, www.jmmc.fr/doc/approved/JMMC-MAN-2720-0001.pdf. JMMC
• Marquardt (1963) Marquardt D. W., 1963, Indian J. Indust. Appl. Math., 11, 431
• Martí-Vidal et al. (2011) Martí-Vidal I., Marcaide J. M., Quirrenbach A., et al., 2011, A&A, 529, A115
• Mérand et al. (2005) Mérand A., Bordé P., Coudé Du Foresto V., 2005, A&A, 433, 1155
• Mérand et al. (2010) Mérand A., Stefl S., Bourget P., et al., 2010, in W. C. Danchi et al. ed., Optical and infrared interferometry II Vol. 7734 of SPIE Conf. Ser., Perspectives for the AMBER Beam Combiner. San Diego, CA
• Millour et al. (2008) Millour F., Petrov R., Malbet F., et al., 2008, in Kaufer A., Kerber F., eds, The 2007 ESO Instrument Calibration Workshop AMBER on the VLTI: data processing and calibration issues. Springer, pp 461–470
• Millour et al. (2007) Millour F., Petrov R. G., Chesneau O., et al., 2007, A&A, 464, 107
• Ohnaka et al. (2011) Ohnaka K., Weigelt G., Millour F., et al., 2011, A&A, 529, A163
• Perrin (2003) Perrin G., 2003, A&A, 400, 1173
• Petrov et al. (2007) Petrov R. G., Malbet F., Weigelt G., et al., 2007, A&A, 464, 1
• Plez (2003) Plez B., 2003, in U. Munari ed., GAIA Spectroscopy: Science and Technology Vol. 298 of ASP Conf. Ser., Cool star atmospheres and spectra for GAIA: MARCS models. pp 189–198
• Plez (2012) Plez B., 2012, Astrophys. Source Code Libr., p. 1205.004
• Press et al. (2007) Press W. H., Teukolsky S. A., Vetterling W. T., et al., 2007, Numerical Recipes 3rd Edition - The Art of Scientific Computing. Cambridge Univ. Press
• Ragland et al. (2006) Ragland S., Traub W. A., Berger J., et al., 2006, ApJL, 652, 650
• Rix et al. (2004) Rix S. A., Pettini M., Leitherer C., et al., 2004, ApJL, 615, 98
• Robbe-Dubois et al. (2007) Robbe-Dubois S., Lagarde S., Petrov R. G., et al., 2007, A&A, 464, 13
• Roddier (1981) Roddier F., 1981, Prog. Opt., 19, 281
• Stefl et al. (2011) Stefl S., Mérand A., Le Bouquin J.-B., et al., 2011, in The Origin and Fate of the Sun: Evolution of Solar-mass Stars Observed with High Angular Resolution Problems of the AMBER wavelength calibration. to be published, Garching b. München, Germany
• Tallon-Bosc et al. (2008) Tallon-Bosc I., Tallon M., Thiébaut E., et al., 2008, in Optical and Infrared Interferometry Vol. 7013 of SPIE Conf. Ser., LITpro: a model fitting software for optical interferometry. Marseille, France, pp 70131J-1–70131J-9
• Tatebe et al. (2006) Tatebe K., Chandler A. A., Hale D. D. S., Townes C. H., 2006, ApJL, 652, 666
• Tatulli et al. (2007) Tatulli E., Millour F., Chelli A., et al., 2007, A&A, 464, 29
• van Belle & van Belle (2005) van Belle G. T., van Belle G., 2005, PASP, 117, 1263
• Vera & Torres (2008) Vera E., Torres S., 2008, in EURASIP ed., 16th ESP Conf. (EUSIPCO 2008) Subpixel accuracy analysis of phase correlation registration methods using aliased imagery. Lausanne, Switzerland
• Weigelt et al. (2011) Weigelt G., Grinin V. P., Groh J. H., et al., 2011, A&A, 527, A103
• Weigelt (1977) Weigelt G. P., 1977, Opt. Commun., 21, 55
• Winzer (2000) Winzer P., 2000, Rev. Sci. Instrum., 71, 1447
• Wittkowski et al. (2008) Wittkowski M., Boboltz D. A., Driebe T., et al., 2008, A&A, 479, L21
• Young (2003) Young J., 2003, in Perrin G., Malbet F., eds, Observing with the VLTI Vol. 6 of EAS Publ. Ser., From Visibilities to Science with Simple Models. pp 181–190
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters