© ESO 2020

1. Introduction

With the first detection and characterisation of an exoplanet by the Very Large Telescope Interferometer (VLTI) instrument GRAVITY (HR 8799 e, Gravity Collaboration 2019), infrared interferometry has proven to be a powerful technique for high-contrast imaging at high angular resolution. Although initially designed for observations of the galactic centre (Bartko et al. 2009), GRAVITY’s dual-feed mode combined with the recently installed integrated optics beam combiner (Perraut et al. 2018) enable spectroscopy and micro-arcsecond astrometry of exoplanets with a wide range of angular separations (Gravity Collaboration 2019).

More recently, Gravity Collaboration (2020) have used GRAVITY observations of β Pic b in order to derive reliable estimates for the mass and the C/O ratio of the young giant planet using forward modelling and free retrieval of its atmosphere. In the future, infrared interferometry will be a promising opportunity for studying giant planet formation (e.g. with Hi-5, Defrère et al. 2018) and potentially even characterising terrestrial exoplanets from space (e.g. with a formation-flying nulling interferometer, Léger et al. 1996; Mennesson & Mariotti 1997; Kammerer & Quanz 2018; Quanz et al. 2018, 2019). However, significant improvements are required on the technical side (e.g. kernel nulling, Martinache & Ireland 2018), in addition to on the data reduction side in order to achieve these ambitious goals.

Because they use the dual-feed mode of GRAVITY, the aforementioned observations are not conducted anywhere close to the diffraction limit of the interferometer, but rather the diffraction limit of a single telescope. Detecting a companion within the interferometer’s diffraction limit (a few λ/b_max), where λ is the observing wavelength and b_max is the longest baseline of the interferometer, is limited by systematic errors. While such systematic errors that are introduced by instrumental and atmospheric effects have been studied intensively (e.g. imperfect fibre coupling, Kotani et al. 2003; instrument vibrations, Le Bouquin et al. 2011; differential atmospheric piston; Colavita 1999), correlations are also introduced by the data reduction and the calibration. For instance, a systematic error might be introduced similarly to all complex visibilities measured on the science target if the instrumental transfer function obtained from the calibrator target is affected by unknown variability (Perrin 2003) and if the closure phases measured over telescope triplets of closing triangles are not mathematically independent (Monnier 2007). Nevertheless, most data reduction pipelines (e.g. the PIONIER data reduction pipeline, Le Bouquin et al. 2011; the GRAVITY data reduction pipeline, Lapeyrere et al. 2014) and model fitting routines (e.g. LITpro¹, Tallon-Bosc et al. 2008; CANDID², Gallenne et al. 2015) assume statistically independent observables. However, in order to robustly detect faint companions, or place upper limits on their brightness, a solid understanding and description of the systematic errors is inevitable.

While Lachaume et al. (2019) proposed to use the bootstrapping method (i.e. sampling with replacement, Efron & Tibshirani 1986) in order to obtain the multivariate probability density function of the squared visibility amplitudes and the closure phases, Gravity Collaboration (2020) extracted the covariances of the complex visibilities directly from the data. Although the bootstrapping method is computationally expensive, it enables estimating the systematic errors not only between the different spectral channels, baselines and triangles, but also between different observations. This enables accounting for correlations introduced by sky rotation or the calibration method, but is only applicable at a higher level when the structure of the observing sequence is known.

In this paper, we follow a similar approach to Gravity Collaboration (2020) by extracting the correlations between the squared visibility amplitudes and the closure phases directly from the data. Then, we develop an empirical model for these correlations which can be fitted to the correlations extracted from single GRAVITY pipeline products, even if only a small number of measurements is available. This enables the attainment of a systematic error estimate for every GRAVITY data set and could ultimately be included in the GRAVITY data reduction pipeline (Lapeyrere et al. 2014).

2. Methods

In Sect. 2.1 we show how we extract the correlations from individual GRAVITY pipeline products and describe their nature. In Sect. 2.2 we introduce our empirical model for these correlations and in Sect. 2.4 we present the model fitting routines with the aid of which we show the improvements that come from using our empirical correlation model. The Python code that we developed in the scope of this paper is publicly available on GitHub³.

2.1. Correlations extracted from GRAVITY data

In order to extract the correlations between the different spectral channels, baselines and triangles from GRAVITY data we use the P2VM-reduced files from the GRAVITY data reduction pipeline. These files are intermediate pipeline products which contain the individual measurements (detector read-outs) before they are averaged together. Having access to the individual measurements enables extracting the correlations from the (complex) coherent flux VIS_mbλ which is stored in the P2VM-reduced file as a data cube of shape M × B × Λ, where m = 1…M is the number of individual measurements, b = 1…B is the number of baselines and λ = 1…Λ is the number of spectral channels. From the coherent flux, we compute the squared visibility amplitudes

$$\begin{array}{c}\hfill {\mathrm{VIS}2}_{mb\lambda }=\frac{|{\mathrm{VIS}}_{mb\lambda }{|}^2}{{\mathrm{F}1\mathrm{F}2}_{mb\lambda }},\end{array}$$(1)

where F1F2_mbλ is the product of the total fluxes, and the closure phases

$$\begin{array}{c}\hfill {\mathrm{T}3}_{mt\lambda }=\mathit{K}·\angle {\mathrm{VIS}}_{mb\lambda },\end{array}$$(2)

where t = 1…T is the number of triangles, K is a stack of M matrices K which encode how the four unique triangles can be formed from the six unique baselines of the VLTI, that is

$$\begin{array}{c}\hfill \mathit{k}=\left(\begin{array}{cccccc}1& -1& 0& 1& 0& 0\\ 1& 0& -1& 0& 1& 0\\ 0& 1& -1& 0& 0& 1\\ 0& 0& 0& 1& -1& 1\end{array}\right),\end{array}$$(3)

and ∠ denotes the argument of a complex number (i.e. the phase).

Then, we compute the sample covariance of the squared visibility amplitudes and the closure phases according to

$$\begin{array}{c}\hfill {({\mathrm{\Sigma }}_X)}_{\mathit{ij}}=\frac{1}{M-1}{\displaystyle \sum _{m=1}^M(X_{\mathit{mi}}-{\overline{X}}_i)(X_{\mathit{mj}}-{\overline{X}}_j),}\end{array}$$(4)

where X is VIS2/T3 and i and j run over 1…BΛ/1…TΛ so that we obtain covariance matrices of shape (BΛ)  ×  (BΛ)/(TΛ)  ×  (TΛ) that contain the covariances between the different spectral channels and baselines/triangles. $\overline{X}$ denotes the mean of X over the individual measurements, that is

$$\begin{array}{c}\hfill {\overline{X}}_i=\frac{1}{M}{\displaystyle \sum _{m=1}^M{X}_{\mathit{mi}}.}\end{array}$$(5)

The correlations between the VIS2 and the T3 then follow by dividing the covariances by the standard deviation ${\sigma }_i=\sqrt{{\mathrm{\Sigma }}_{\mathit{ii}}}$ of the corresponding observables, that is

$$\begin{array}{c}\hfill {(C_X)}_{\mathit{ij}}=\frac{{({\mathrm{\Sigma }}_X)}_{\mathit{ij}}}{{({\sigma }_X)}_i{({\sigma }_X)}_j}.\end{array}$$(6)

The diagonal of the covariance matrix Σ equals the square of the standard deviation and that the diagonal of the correlation matrix C equals one by definition.

For developing an empirical correlation model we use data taken with GRAVITY (Gravity Collaboration 2017) at the Very Large Telescope Interferometer (VLTI) during technical time (programme 60.A-9801(U)). GRAVITY operates in the K-band (2.0–2.4 μm) and combines the light from either the four Unit Telescopes (UTs) or the four Auxiliary Telescopes (ATs) of the VLTI in order to perform interferometric imaging and astrometry by phase referencing⁴. The data used here was taken with the four UTs on the object HD 82383 (ζ Ant B) in single-field medium resolution (R = λ/Δλ ≈ 500) mode. This object is relatively bright (K-band magnitude of 6.698, cf. SIMBAD⁵), hence the short integration time of 0.85 ms for the fringe tracker and 1 s for the science camera. It has a companion at an angular separation of ∼8 arcsec (cf. WDS⁶) which is well beyond the interferometric field-of-view. By choosing a bright target with a short exposure time we make sure that there is a sufficient number of frames to compute the sample covariance (cf. Eq. (4)). The short exposure time of the fringe tracker (much less than the atmospheric coherence time t₀, which is typically ∼20 ms in the K-band, Kellerer & Tokovinin 2007) and spatially filtered nature of the GRAVITY beam combiner means that the fringe tracker data is less affected by systematic errors than many other beam combiners. We extract the correlations from a single P2VM-reduced file (GRAVI.2019-03-29T02-01-37.193_singlecalp2vmred.fits) in order to demonstrate the direct applicability of our method to the GRAVITY data reduction pipeline. Correlations extracted from other P2VM-reduced files of the same program can be found in Fig. A.1.

Figure 1 shows the correlations of the VIS2 (left panel) and the T3 (right panel) for the GRAVITY fringe tracker. There are six different baselines and four different triangles with five spectral channels each, so 30 observables for the VIS2 and 20 observables for the T3 in total. Correlations within the same baseline/triangle are highlighted with red squares and correlations between baselines having a telescope in common are highlighted with orange squares.

Fig. 1. Correlations of the VIS2 (left panel) and the T3 (right panel) for the GRAVITY fringe tracker, extracted from a single P2VM-reduced file. The axes run over the different baselines/triangles, with each individual baseline/triangle comprising five spectral channels. Correlations within the same baseline/triangle are highlighted with red squares and correlations between baselines having a telescope in common are highlighted with orange squares. We note that the correlations are computed from 46 592 individual measurements. Below each panel, the variance of the data and the names of the telescopes forming each baseline/triangle are shown.

The most dominant correlations of the VIS2 are between different spectral channels within the same baseline, with neighbouring spectral channels being affected most strongly. We suspect that these correlations are predominantly of both atmospheric or instrumental origin, since all five spectral channels follow the same optical path through the atmosphere and up to the dispersive element behind the beam combiner and before the science camera. Also, the five spectral channels do not correspond to individual pixels on the detector of the fringe tracker. In fact, the wavelengths of the five spectral channels lie somewhere between the wavelengths corresponding to the pixels on the detector of the fringe tracker, so that the values recorded by two neighbouring pixels on the detector need to be interpolated in order to find the values for the five spectral channels of the fringe tracker. This could explain the strong correlations between neighbouring spectral channels (one pixel above or below the diagonal) observed for the VIS2, but also for the T3. Furthermore, there are significant correlations between baselines having a telescope in common. Their strength is roughly half the strength of the correlations within the same baseline, which makes sense if the correlations are introduced by atmospheric or instrumental effects and affect each of the four individual beams of the interferometer separately. Also, baselines having no telescope in common are essentially uncorrelated. Hence, we conclude that most of the correlations of the VIS2 are caused by atmospheric or instrumental effects.

For the T3, we observe similar correlations between neighbouring spectral channels as for the VIS2. This makes sense since the closure phases are built from a linear combination (encoded in the matrix K) of the phase of the complex visibilities, whose absolute square are the squared visibility amplitudes. Moreover, there are significant correlations of ∼ ± 1/3 between the same spectral channels on different triangles. These are caused by the fact that each set of two different triangles has exactly one of their three baselines in common, that is each column of the matrix K has exactly two non-zero entries. If the common baseline is shared between the different triangles in parallel direction (i.e. the two entries in the corresponding column of the matrix K have the same sign), the correlation is +1/3, otherwise it is −1/3. This structure with the side-diagonals being ±1/3 can also be explained by assuming uncorrelated visibility phases (i.e. a diagonal correlation matrix

$$\begin{array}{c}\hfill {\mathit{C}}_{\angle \mathrm{VIS}}=\left(\begin{array}{cccc}1& 0& \cdots & 0\\ 0& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & 0& 1\end{array}\right)\end{array}$$(7)

of shape B × Λ for the visibility phases) and performing a basis transform

$$\begin{array}{c}\hfill {\mathit{C}}_{\mathrm{T}3}=\mathit{T}·{\mathit{C}}_{\angle \mathrm{VIS}}·{\mathit{T}}^T,\end{array}$$(8)

where T represents a matrix of shape (TΛ) × (BΛ) which maps the vector of visibility phases to the vector of closure phases and can be trivially obtained from the matrix K. Also, the observed correlations between neighbouring spectral channels on different triangles (pixels next to the side-diagonals) are naturally explained by this basis transform given the correlations of the VIS2 observed between neighbouring spectral channels on the same baseline.

Figure 2 shows the correlations of the VIS2 (left panel) and the T3 (right panel) for the GRAVITY science camera. There are six different baselines and four different triangles with 210 spectral channels each, so 1260 observables for the VIS2 and 840 observables for the T3 in total.

Fig. 2. Same as Fig. 1, but for the GRAVITY science camera. Each individual baseline/triangle comprises 210 spectral channels. We note that the correlations are computed from 32 individual measurements.

Due to the much smaller number of individual measurements if compared to the fringe tracker the correlations of the science camera are more dominated by noise. Nevertheless, we observe strong positive correlations between different spectral channels within the same baseline (i.e. inside the red squares) and significant positive correlations between baselines having a telescope in common (i.e. inside the orange squares) for the VIS2, similar to the correlations observed for the fringe tracker. Although the atmospheric turbulence and the optical elements (i.e. mirrors, delay lines, optical fibres, beam combiner) seen by the science camera are similar to those seen by the fringe tracker, the exposure time of the science camera is much longer than both the atmospheric coherence time t₀ and the fringe tracker inverse 3 dB bandwidth (Lacour et al. 2019), which means that the VIS2 correlations are expected to be decreased by a term proportional to the square of the fringe tracking error and the closure phase random errors are expected to be proportional to the cube of the fringe tracking error (Ireland 2013). Therefore, since there still are significant correlations for the science camera, they must be introduced by the (correlated) fringe tracker, forwarding the correlations shown in Fig. 1 to the science camera. For the T3, we again observe significant correlations of ∼ ± 1/3 between the same spectral channels on different triangles. On top of this, there are also weak positive correlations between different spectral channels on the same triangle (i.e. inside the red squares) and between different spectral channels on different triangles whose sign depends on whether the corresponding triangles share a baseline in parallel or anti-parallel direction. Again, these correlations are naturally explained by the basis transform T given the correlations observed for the VIS2 of the science camera.

2.2. Empirical model for the correlations

An empirical VIS2/T3 sample covariance with fewer frames than the product of the number of baselines/triangles and spectral channels is necessarily singular. It takes a number of frames much greater than this to estimate a sample covariance matrix with a condition number approaching that of the true sample covariance. For this reason, we choose to develop an analytic model for the covariance matrix Σ of the VIS2 and the T3. This model can then be fitted to the (potentially under-conditioned) covariance extracted from an arbitrary GRAVITY data set and can be used for model fitting based on log-likelihood maximisation. Most model fitting routines (e.g. LITpro, Tallon-Bosc et al. 2008; CANDID, Gallenne et al. 2015) are based on χ² minimisation, which is equivalent to log-likelihood maximisation, where

$$\begin{array}{c}\hfill {\chi }^2=R^T·{\mathbf{\Sigma }}^{-1}·R\end{array}$$(9)

and R = D − M is the residual between data and model (cf. Sect. 2.4).

Our approach is to model the correlation matrices C_VIS2 of the VIS2 and C_T3 of the T3 which have the relatively simple structure observed in Figs. 1 and 2. Moreover, the observed structure of the correlations is consistent for different data sets with different exposure times (1 s with the UTs for programme 60.A-9801(U) and 10 s with the ATs for programme 0101.C-0907(B), cf. also Figs. A.1 and A.2). Then, we compute

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}=C_{\mathit{ij}}{\sigma }_i{\sigma }_j,\end{array}$$(10)

where σ denotes the standard deviation of the data which can be obtained from the VIS2ERR and the T3PHIERR columns of the OIFITS files for example. We note that these standard deviations are used to build diagonal covariance matrices in LITpro and CANDID which assume uncorrelated data only. Of course, assuming uncorrelated data is a simplification and we discuss the problems that arise from this in Sect. 3.2.

A very important point is that Eq. (10) only holds if the errors on the VIS2 (σ_VIS2) and the T3 (σ_T3) are reliably estimated by the GRAVITY data reduction pipeline. The pipeline manual⁷ explains that the uncertainties are computed by bootstrapping over ∼10 independent samples, so that the final error on the mean measurement is estimated from the observed statistics at a slightly higher temporal frequency. There is no re-scaling or accounting for systematics in this process. In case there are less than five frames available, Monte-Carlo realisations of the theoretical photon and detector noise are added to the samples, which leads to less realistic uncertainties. However, our data sets consist of 32 frames exposures for programme 60.A-9801(U) and 20 frames exposures for programme 0101.C-0907(B), respectively. While we understand that the use of the pipeline uncertainties is a limitation and that an incorrect noise model can reduce the detection sensitivity or yield false positives (cf. e.g. Sect. 3 of Delisle et al. 2020), we also note that investigating and quantifying the credibility of these uncertainties is beyond the scope of this work.

Our models for the correlation matrices equal one on the diagonal according to the definition of a correlation matrix (cf. Eq. (6), that is every observable is 100% correlated with itself) and have one free parameter which can be determined by fitting the model to the correlations extracted from the P2VM-reduced files. For the correlation matrix of the VIS2 C_VIS2, the free parameter x represents the correlations between spectral channels within the same baseline and between baselines having a telescope in common. There are correlations of x between different spectral channels within the same baseline, correlations of x/2 between baselines having a telescope in common, and no correlations between baselines having no telescope in common (cf. left panel of Fig. 2), that is

$$\begin{array}{cc}& {\mathit{C}}_{\mathrm{VIS}2}=\left(\begin{array}{ccccccc}{\mathit{X}}_1& {\mathit{X}}_2& \cdots & \cdots & \cdots & {\mathit{X}}_2& 0\\ {\mathit{X}}_2& \ddots & \ddots & & & {\mathit{X}}_2\\ ⋮& \ddots & \ddots & {\mathit{X}}_2& 0& & ⋮\\ ⋮& & {\mathit{X}}_2& \ddots & {\mathit{X}}_2& & ⋮\\ ⋮& & 0& {\mathit{X}}_2& \ddots & \ddots & ⋮\\ {\mathit{X}}_2& & & \ddots & \ddots & {\mathit{X}}_2\\ 0& {\mathit{X}}_2& \cdots & \cdots & \cdots & {\mathit{X}}_2& {\mathit{X}}_1\end{array}\right),\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\mathit{X}}_1=\left(\begin{array}{cccc}1& x& \cdots & x\\ x& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & x\\ x& \cdots & x& 1\end{array}\right),\hfill \end{array}$$(12)

$$\begin{array}{cc}& {\mathit{X}}_2=\left(\begin{array}{ccc}x/2& \cdots & x/2\\ ⋮& \ddots & ⋮\\ x/2& \cdots & x/2\end{array}\right).\hfill \end{array}$$(13)

The correlation matrix is a block matrix consisting of B × B blocks, where each individual block is a Λ × Λ matrix. For the correlation matrix of the T3 C_T3, the free parameter y represents the correlations between spectral channels within the same triangle. Moreover, as illustrated by the basis transform T, this naturally leads to correlations of ±1/3 between the same spectral channel of different triangles and ±y/3 between different spectral channels of different triangles (cf. right panel of Fig. 2), that is

$$\begin{array}{cc}& {\mathit{C}}_{\mathrm{T}3}=\left(\begin{array}{cccc}{\mathit{Y}}_1& {\mathit{Y}}_2& \cdots & {\mathit{Y}}_2\\ {\mathit{Y}}_2& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & {\mathit{Y}}_2\\ {\mathit{Y}}_2& \cdots & {\mathit{Y}}_2& {\mathit{Y}}_1\end{array}\right),\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\mathit{Y}}_1=\left(\begin{array}{cccc}1& y& \cdots & y\\ y& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & y\\ y& \cdots & y& 1\end{array}\right),\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\mathit{Y}}_2=\left(\begin{array}{cccc}\pm 1/3& \pm y/3& \cdots & \pm y/3\\ \pm y/3& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \pm y/3\\ \pm y/3& \cdots & \pm y/3& \pm 1/3\end{array}\right).\hfill \end{array}$$(16)

The correlation matrix is a block matrix consisting of T × T blocks, where each individual block is a Λ × Λ matrix. The sign is positive if the two triangles share a baseline in parallel direction and negative if they share a baseline in anti-parallel direction.

We fit the previously described model to the correlations of the VIS2 and the T3 which we extracted from the single P2VM-redcued file of GRAVITY introduced in Sect. 2.1. Figure 3 shows the extracted and the model correlations (top panels) and the extracted and the model covariances (bottom panels) for the VIS2. The free parameter x takes a value of ∼3.2e−1. Figure 4 shows the same for the T3 and the free parameter y takes a value of ∼7.4e−2.

Fig. 3. Correlations of the VIS2 for the GRAVITY science camera, extracted from a single P2VM-reduced file (upper left panel) and our one-parameter model fitted to them (upper right panel). The bottom panels show the corresponding covariances obtained by multiplying the correlation Cij with the product of the standard deviations σiσj. Correlations/covariances within the same baseline are highlighted with red squares and correlations/covariances between baselines having a telescope in common are highlighted with orange squares.

Fig. 4. Same as Fig. 3, but showing the correlations/covariances of the T3 and our one-parameter model fitted to them for the GRAVITY science camera.

2.3. Simulated and real data

In order to demonstrate the improvement that comes from taking into account the correlations between the data we perform companion injection and recovery tests with simulated and real data. Therefore, we use GRAVITY data of ζ Ant B from the technical time programme 60.A-9801(U) and of HIP 78183 from the normal programme 0101.C-0907(B), PI M. J. Ireland, listed in Table 1. Both objects were observed in single-field medium resolution mode, but the former one with the four UTs and the latter one with the four ATs (medium configuration D0-G2-J3-K0).

From the file GRAVI.2019-03-29T02-01-37.193_singlecalp2vmred.fits we have already extracted the covariances and correlations and fitted our empirical models to them (cf. Sect. 2.2). For the companion injection and recovery tests with simulated data, we simply use these models and the uv-tracks u and v of the files belonging to programme 60.A-9801(U) listed in Table 1 in order to obtain a realistic uv-coverage over ∼20 min (cf. Fig. 5). We simulate the complex visibility of a uniform disc with an unresolved companion according to

$$\begin{array}{cc}& {\mathrm{VIS}}_{\mathrm{bin}}=\frac{{\mathrm{VIS}}_{\mathrm{ud}}+fexp(-2\pi i(\frac{{\mathrm{\Delta }}_{\mathrm{RA}}u}{\lambda }+\frac{{\mathrm{\Delta }}_{\mathrm{Dec}}v}{\lambda }))}{1+f},\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\mathrm{VIS}}_{\mathrm{ud}}=\frac{2J_1(\pi \theta b)}{\pi \theta b},\hfill \end{array}$$(18)

Fig. 5. Fourier u- and v-baselines of our simulated data extracted from three on-sky observations of GRAVITY using the four UTs over ∼20 min.

where 0 ≤ f ≤ 1 is the relative flux of the companion, Δ_RA and Δ_Dec are the on-sky separation in the direction of the celestial north and the celestial east between the companion and its host star, λ is the observing wavelength, J₁ is the first order Bessel function of first kind, θ is the angular diameter of the uniform disc and $b=\sqrt{u^2+v^2}$ is the distance between the two telescopes observing the object. The squared visibility amplitudes and the closure phases follow according to

$$\begin{array}{cc}& {\mathrm{VIS}2}_{\mathrm{bin}}={|{\mathrm{VIS}}_{\mathrm{bin}}|}^2,\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\mathrm{T}3}_{\mathrm{bin}}=\mathit{k}·\angle {\mathrm{VIS}}_{\mathrm{bin}}.\hfill \end{array}$$(20)

Then, we add correlated noise to the VIS2_bin and the T3_bin by drawing from a multivariate normal distribution with covariance Σ_VIS2, fit and Σ_T3, fit, which we obtain from our correlation model C_VIS2, fit and C_T3, fit (cf. Sect. 2.2) and assuming a standard deviation of σ_VIS2 = 0.01 and σ_T3 = 1 deg.

For the companion injection and recovery tests with real data, we extract the correlations of the visibility amplitudes VISAMP (instead of the squared visibility amplitudes VIS2) and the closure phases T3 from the P2VM-reduced files belonging to programme 0101.C-0907(B) listed in Table 1, fit our empirical models to them and compute the covariances using Eq. (10) and the errors from the corresponding final GRAVITY pipeline products (the “singlesciviscalibrated” files). Using the VISAMP instead of the VIS2 can yield better results in some cases where the normalisation of the VIS2 is not done properly by the GRAVITY data reduction pipeline. From the final GRAVITY pipeline products, we also extract the VISAMP and the T3 and inject an unresolved companion according to

$$\begin{array}{cc}& {\mathrm{VISAMP}}_{\mathrm{inj}}=\mathrm{VISAMP}·|{\mathrm{VIS}}_{\mathrm{bin}}|,\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\mathrm{T}3}_{\mathrm{inj}}=\mathrm{T}3+\mathit{k}·\angle {\mathrm{VIS}}_{\mathrm{bin}},\hfill \end{array}$$(22)

where we set VIS_ud to one. The VISAMP are simply the square root of the VIS2, so that our correlation model and fitting routines can be equivalently applied in the high signal-to-noise regime.

2.4. Model fitting

We search for faint companions in the data by fitting the model for a uniform disc with an unresolved companion (cf. Eq. (17)) to it. We maximise the log-likelihood of the model by minimising its χ² in order to find the best fit parameters p_fit of the model, that is

$$\begin{array}{c}\hfill p_{\mathrm{fit}}={\mathrm{argmin}}_p({\chi }^2)={\mathrm{argmin}}_p(R^T·{\mathbf{\Sigma }}^{-1}·R),\end{array}$$(23)

where R = D − M is the residual between data and model and p = (f, Δ_RA, Δ_Dec, θ) is the four-dimensional parameter vector of the model.

In order to find the global minimum of the χ² within a given range of companion separations, we first find a prior for the uniform disc diameter θ₀ by fitting the corresponding model (cf. Eq. (18)) to the data. Then, we perform a set of minimisations with priors on a Δ_RA − Δ_Dec grid, the uniform disc diameter θ₀, and a small relative flux f₀ = 1e−3. This is necessary since the χ² hyper-surface is bumpy (i.e. has many local extrema) if projected onto the Δ_RA − Δ_Dec surface and the BFGS algorithm which is used to minimise the χ² converges on local minima. The bumpiness is a result of the sparse uv-coverage of a long-baseline optical interferometer which causes the sensitivity to vary substantially over the FOV.

The above method relies on the covariance matrix Σ being invertible. This is not the case for a sample covariance that is estimated from a small number of frames, which is usually singular, and is the reason why we develop an empirical covariance model. However, our empirical model for the covariances of the closure phases is also singular, since the fourth triangle can be written as a linear combination of the other three. There are multiple solutions to this problem, and for simplicity we decide to completely ignore the data recorded on the fourth triangle since it theoretically is redundant anyway⁸. There are more sophisticated methods to keep the data recorded on the fourth triangle, such as the “jackknife” method (i.e. averaging over four model fits using data recorded on different sets of three triangles), projection into a sub-space that preserves the information in the covariance matrix (Blackburn et al. 2020), and the approach from Kulkarni (1989) which is adding a small numerical value ϵ ≪ 1 to the diagonal of the covariances of the closure phases, that is

$$\begin{array}{c}\hfill {\mathbf{\Sigma }}_{\mathrm{T}3,\mathrm{invertible}}={\mathbf{\Sigma }}_{\mathrm{T}3,\mathrm{fit}}+ϵ·\mathrm{id},\end{array}$$(24)

where id is the identity matrix, so that the covariance matrix becomes numerically invertible.

Finally, in order to determine the statistical significance of a detected companion, we compute the probability that the binary model is preferred over the uniform disc model according to

$$\begin{array}{c}\hfill P=1-{\mathrm{CDF}}_{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}}}\left(\frac{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}}{\chi }_{\mathrm{red},\mathrm{ud}}^2}{{\chi }_{\mathrm{red},\mathrm{bin}}^2}\right),\end{array}$$(25)

where CDF_Nd.o.f is the χ² cumulative distribution function with N_d.o.f degrees of freedom, ${\chi }_{\text{red},\text{ud}}^2$ is the reduced χ² of the best fit uniform disc model and ${\chi }_{\text{red},\text{bin}}^2$ is the reduced χ² of the best fit binary model (cf. Gallenne et al. 2015).

If the host star is essentially unresolved (i.e. θbλ ≪ 1) and the companion is at high contrast (i.e. f ≪ 1) one can linearise the VIS2_bin and the T3_bin as a function of the relative flux of the companion f according to

$$\begin{array}{cc}& {\mathrm{VIS}2}_{\mathrm{bin}}\propto 1+f,\hfill \end{array}$$(26)

$$\begin{array}{cc}& {\mathrm{T}3}_{\mathrm{bin}}\propto f.\hfill \end{array}$$(27)

A more detailed derivation of this relationship can be found in Appendix B. Let D be the data, M_ref a reference binary model which is normalised to the relative flux of the companion f_ref, and Σ the covariances between the data, that is

$$\begin{array}{cc}& D=\left(\begin{array}{c}\mathrm{VIS}2-1\\ \mathrm{T}3\end{array}\right),\hfill \end{array}$$(28)

$$\begin{array}{cc}& M_{\mathrm{ref}}=\left(\begin{array}{c}({\mathrm{VIS}2}_{\mathrm{bin},\mathrm{ref}}-1)/f_{\mathrm{ref}}\\ {\mathrm{T}3}_{\mathrm{bin},\mathrm{ref}}/f_{\mathrm{ref}}\end{array}\right),\hfill \end{array}$$(29)

$$\begin{array}{cc}& \mathbf{\Sigma }=\left(\begin{array}{cc}{\mathbf{\Sigma }}_{\mathrm{VIS}2}& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Sigma }}_{\mathrm{T}3}\end{array}\right),\hfill \end{array}$$(30)

where VIS2_bin, ref and T3_bin, ref are the binary model VIS2 and T3 evaluated at a reference relative flux f_ref = 1e−3. Then, the best fit relative flux f_fit and its uncertainty σ_ffit follow according to

$$\begin{array}{cc}& f_{\mathrm{fit}}=\frac{M_{\mathrm{ref}}^T·{\mathbf{\Sigma }}^{-1}·D}{M_{\mathrm{ref}}^T·{\mathbf{\Sigma }}^{-1}·M_{\mathrm{ref}}},\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\sigma }_{f_{\mathrm{fit}}}=\frac{1}{\sqrt{M_{\mathrm{ref}}^T·{\mathbf{\Sigma }}^{-1}·M_{\mathrm{ref}}}},\hfill \end{array}$$(32)

(cf. Le Bouquin & Absil 2012; Kammerer et al. 2019). Equation (31) can be computed on a Δ_RA − Δ_Dec grid, and the best fit parameters p_fit follow from the grid position which minimises

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red}}^2=\frac{{\chi }^2}{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}}}=\frac{R^T·{\mathbf{\Sigma }}^{-1}·R}{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}}},\end{array}$$(33)

where N_d.o.f is the number of the degrees of freedom. This grid search technique is commonly used in order to find the global minimum of the (f, Δ_RA, Δ_Dec, θ) parameter space and its corresponding ${\chi }_{\text{red}}^2$ (e.g. Absil et al. 2011; Gallenne et al. 2015). However, the statistical structure of the grid is complex due to redundancy and periodicity in sensitivity originating from the very limited uv-coverage of a sparse interferometer such as the VLTI (cf. Fig. 5). Therefore, the detection significance is derived using Eq. (25) which yields the probability that the binary model is preferred over the uniform disc model (without any companion).

3. Results

We evaluate the impact of our full covariance model by performing model fitting and companion injection and recovery tests with both simulated and real GRAVITY data. In Sect. 3.1, we simulate data without any astronomical object, that is correlated noise only, and use model fitting to determine the fundamental detection limits when assuming uncorrelated data (i.e. a diagonal covariance) and correlated data (i.e. our full covariance model). In Sects. 3.2 and 3.3 we inject companions with different relative fluxes and separations into simulated and real GRAVITY data and try to recover them, again assuming both uncorrelated and correlated data.

3.1. Model fitting to correlated noise

In order to compare the fundamental detection limits when assuming uncorrelated and correlated data, we simulate 100 GRAVITY data sets of an unresolved host star without any companion (i.e. θ = f = 0 in Eq. (17)) affected by realistic correlated noise (Σ_VIS2, fit and Σ_T3, fit, cf. Sect. 2.3). Then, we use Eq. (31) in order to compute the best fit relative flux f_fit on a Δ_RA − Δ_Dec grid for each of the 100 simulated data sets, first assuming uncorrelated data (i.e. a diagonal covariance diag(Σ_VIS2, fit) and diag(Σ_T3, fit)) and then assuming correlated data (i.e. our full covariance model) in Eq. (31). Since no companion was injected into the data, the best fit flux ratios of these grids represent the fundamental contrast floor. Any companion with a higher contrast (i.e. smaller flux) would not be distinguishable from the noise. By computing an azimuthal average of these grids, we obtain a fundamental 1 − σ contrast curve (i.e. best fit relative flux vs. angular separation curve).

The mean of the 100 azimuthal averages obtained for each of these two scenarios (model fitting assuming uncorrelated data in blue and correlated data in orange) is shown in the left panel of Fig. 6. The contrast floor remains roughly constant at a contrast of ∼6e−4 outward of an angular separation of ∼5 mas for the scenario assuming uncorrelated data. This is because at a contrast of ∼6e−4 one is dominated by the systematic (i.e. the correlated) errors. However, for the scenario assuming correlated data, the fundamental 1 − σ detection limit continues to decrease with increasing angular separation. At an angular separation of ∼10 mas it is already a factor of four better than the limit assuming uncorrelated data (cf. dashed black curve). We note that such a behaviour has already been observed by Ireland (2013) for orthogonal kernel phases and statistically independent kernel phases, which are obtained by projecting the orthogonal kernel phases into an eigenspace with zero covariances. Its reason is that at small angular separations, the detection limits rely on the average of the VIS2 and the T3 over the spectral channels, while with increasing angular separations the VIS2 and T3 vary within the spectral bands and the impact of the correlations is growing. A flat uncorrelated contrast curve (as a function of angular separation) is further consistent with previous works on interferometric observables assuming uncorrelated data (e.g. Absil et al. 2011; Gallenne et al. 2015).

Fig. 6. Left panel: contrast curve (i.e. azimuthal average of the best fit relative flux) for simulated data of an unresolved host star without any companion affected by correlated errors, computed with model fitting assuming uncorrelated data (blue curve) and correlated data (orange curve). Both curves show the mean contrast curve over 100 simulated data sets and the shaded region highlights its standard deviation. The dashed black line shows the ratio of the blue and the orange curve, representing the improvement (i.e. the factor by which the detection limits improve) when using our correlated error model instead of the classical uncorrelated one. Right panel: same, but for the real GRAVITY data introduced in Sect. 2.3.

Furthermore, when using our full covariance model, the contrast floor is also more stable for different representations of the noise (highlighted by the shaded regions in the left panel of Fig. 6 which show the standard deviation of the contrast curves over the 100 simulated data sets) meaning that the derived detection limits can be regarded more robust (i.e. independent of the exact representation of the noise which is a random component). Hence, if an observer is only working with a small number of data sets, they will still be able to derive universally valid detection limits.

The right panel of Fig. 6 shows the same plot, but for model fitting to the real GRAVITY data consisting of the three files belonging to programme 0101.C-0907(B) listed in Table 1. Since there is only one real GRAVITY data set we cannot compute or show any standard deviation. The plot looks similar, except for the fundamental detection limits being about two orders of magnitude worse and the ratio between the two scenarios being a lot smaller due to much weaker correlations being present in the real data if compared to the simulated data (cf. Fig. A.1).

3.2. Injection and recovery tests (simulated data)

As a next step, we perform companion injection and recovery tests with simulated data, in order to compare the empirical detection limits when assuming uncorrelated and correlated data. Therefore, we simulate GRAVITY data sets (affected by correlated noise) of a 1 mas uniform disc (the host star) and inject companions with a range of relative fluxes and at different positions around the host star, that is

$$\begin{array}{cc}& f\in [{10}^{-4},{10}^{-3.75},{10}^{-3.5},\dots ,{10}^{-1.5}],\hfill \end{array}$$(34)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{RA}}\in [-30,-25,-20,\dots ,30]\mathrm{mas},\hfill \end{array}$$(35)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{Dec}}\in [-30,-25,-20,\dots ,30]\mathrm{mas},\hfill \end{array}$$(36)

using Eq. (17). Then, we perform model fitting with priors on a Δ_RA − Δ_Dec grid in order to find the global minimum of the χ² (cf. Eq. (23)). We note that this method is similar to how CANDID searches for companions for example. Similar to before (cf. Sect. 3.1), we perform the model fitting once assuming uncorrelated data (i.e. a diagonal covariance diag(Σ_VIS2, fit) and diag(Σ_T3, fit)) and once assuming correlated data (i.e. our full covariance model). We classify an injected companion as recovered if the best fit relative flux f_fit differs by no more than 10% from the injected one f_inj and the best fit position (Δ_RA, fit, Δ_Dec, fit) differs by no more than one resolution element of the interferometer from the injected one (Δ_RA, inj, Δ_Dec, inj), that is

$$\begin{array}{cc}& |f_{\mathrm{fit}}-f_{\mathrm{inj}}|/f_{\mathrm{inj}}<0.1,\hfill \end{array}$$(37)

$$\begin{array}{cc}& \sqrt{{({\mathrm{\Delta }}_{\mathrm{RA},\mathrm{fit}}-{\mathrm{\Delta }}_{\mathrm{RA},\mathrm{inj}})}^2+{({\mathrm{\Delta }}_{\mathrm{Dec},\mathrm{fit}}-{\mathrm{\Delta }}_{\mathrm{Dec},\mathrm{inj}})}^2}<\frac{{\lambda }_{\mathrm{mean}}}{2b_{\max }},\hfill \end{array}$$(38)

where λ_mean is the mean of the observed wavelength range (∼2.2 μm for GRAVITY) and b_max is the longest baseline of the interferometer (∼130 m for observations with the VLTI UTs).

The left panel of Fig. 7 shows the fraction of recovered companions as a function of the relative flux of the companion for model fitting assuming uncorrelated data (blue points) and correlated data (orange points). These values are summed over all positions around the host star with 5 mas ≤ ρ ≤ 45 mas, where $\rho =\sqrt{{\mathrm{\Delta }}_{\mathrm{RA},\mathrm{inj}}^2+{\mathrm{\Delta }}_{\mathrm{Dec},\mathrm{inj}}^2}$ is the angular separation, so to avoid any significant influence from the 1 mas uniform disc (the host star). Although one of the findings in this paper is that the contrast curve is not flat outward a few λ/b_max when accounting for the data correlations (cf. Fig. 6), it is still a reasonable simplification to sum over positions with different angular separations.

Fig. 7. Left panel: fraction of correctly recovered companions as a function of the relative flux of the injected companion from injection and recovery tests with simulated GRAVITY data (cf. Sect. 3.2), assuming uncorrelated data (blue points) and correlated data (orange points) for the model fitting. The blue and the orange curves are logistic growth functions fitted to the data points. Right panel: same, but from injection and recovery tests with real GRAVITY data (cf. Sect. 3.3).

The overplotted blue and orange curves are logistic growth functions

$$\begin{array}{c}\hfill g(x)=\frac{L}{1+e^{-k(x-x_0)}},\end{array}$$(39)

fitted to the data points (in log-space), where L = 1 is the upper growth barrier, k is the logistic growth rate and x₀ is the midpoint. The fact that the orange curve is shifted towards the left compared to the blue curve means that the detection limits are fainter when assuming correlated data instead of uncorrelated data. This could be expected since the data is affected by correlated noise and correctly accounting for these correlations in the model fitting should lead to fainter detection limits. For any given fraction of detections (i.e. any point on the y-axis), the ratio of the contrasts (i.e. the x-values) of the blue and the orange curve gives the improvement that comes from our full covariance model over the conventional diagonal covariance model. This ratio varies between ∼3 and ∼4, depending on the fraction of detections, with an average value of ∼3.5. This means that the detection limits improve by a factor of ∼3.5 when assuming correlated data instead of uncorrelated data.

In order to check the validity of our detection criterion (cf. Eq. (25)) we count the number of companions in different categories, which we represent in a confusion matrix

$$\begin{array}{c}\hfill \zeta =\left(\begin{array}{cc}\#\mathrm{of}\mathrm{true}\mathrm{positives}& \#\mathrm{of}\mathrm{false}\mathrm{negatives}\\ \#\mathrm{of}\mathrm{false}\mathrm{positives}& \#\mathrm{of}\mathrm{true}\mathrm{negatives}\end{array}\right),\end{array}$$(40)

where positive/negative refers to a detection being classified as significant/insignificant according to our detection criterion. Hence, in an ideal world, there would only be true positives or true negatives and ζ would be a diagonal matrix. Obviously, this is not the case in the real world where the data is affected by noise. If we choose an optimistic detection criterion (e.g. the significance needs to be above 1 − σ) there will be many false positives (i.e. many detections that are classified as significant, but which are no true companions) and if we choose a pessimistic detection criterion (e.g. the significance needs to be above 5 − σ) there will be many false negatives (i.e. many detections that are classified as insignificant, but which are true companions).

The confusion matrices from our companion injection and recovery tests for a 3 − σ detection criterion are

$$\begin{array}{cc}& {\zeta }_{\mathrm{diag}}=\left(\begin{array}{cc}674& 3\\ 517& 654\end{array}\right),\hfill \end{array}$$(41)

$$\begin{array}{cc}& {\zeta }_{\mathrm{full}}=\left(\begin{array}{cc}826& 213\\ 12& 797\end{array}\right),\hfill \end{array}$$(42)

for assuming uncorrelated data (ζ_diag) and correlated data (ζ_full). In the former case, there is a large fraction of false detections being classified as significant (517/1171 ≈ 44%), whereas in the latter case this fraction (12/809 ≈ 1%) is roughly consistent with a 3 − σ result. The number of true detections is higher when assuming correlated data (1039 = 826 + 213) than when assuming uncorrelated data (677 = 674 + 3) because the detection limits are fainter, and although the number of false negatives (i.e. true companions being classified as insignificant) is a lot higher, the number of true positives is still higher when assuming correlated data. In summary, when using our full covariance model there are less detections above 3 − σ significance than when using the conventional diagonal covariance model (838 = 826 + 12 vs. 1191 = 674 + 517), but the number of true positives (i.e. true companions being classified as significant) is still higher and significant detections are much more reliable since there are almost no false positives. Accounting for the correlations is therefore clearly preferred over ignoring them.

Before proceeding to the injection and recovery tests with real data we also assess the robustness of our correlation model with respect to errors in the model parameters x and y. Therefore, we repeat the injection and recovery tests with wrong correlation and covariance matrices where x and y are only 50% and 25% of their true values respectively. We find that the number of false positives or false negatives increases slightly, but not significantly. This was expected since the detection of asymmetric structure (such as a companion) is governed by the T3 whose correlations are dominated by the correlations of ±1/3 originating from shared baselines among different triangles. For scenarios where the VIS2 have a larger impact on the model (e.g. measuring stellar diameters) we expect that errors in the model parameters, especially x, have a more significant impact.

3.3. Injection and recovery tests (real data)

In the previous section it is obvious that our full covariance model would outperform the conventional diagonal covariance model, since we simulated data affected by correlated noise. Therefore, the crucial next step is to validate our methods with real GRAVITY data sets. For this purpose, we extract the correlations of the VISAMP and the T3 from the files belonging to programme 0101.C-0907(B) listed in Table 1, fit our empirical models to them and use the VISAMP and the T3 data from the corresponding final GRAVITY pipeline products as noise model.

Since the real data is affected by bright speckles arising from an imperfect calibration, for which our correlation model does not account, we subtract the theoretical VISAMP and T3 of the best fit companion from the data before performing the injection and recovery tests. This also helps us to enter the medium-contrast regime (f_fit ⪅ 10%) where the linearisation of the binary model (see Appendix B) holds. The parameters of the subtracted best fit companion are

$$\begin{array}{c}\hfill p_{\mathrm{sub}}=(f,{\mathrm{\Delta }}_{\mathrm{RA}},{\mathrm{\Delta }}_{\mathrm{Dec}},\theta )=(0.0383,0.20,-5.46,0),\end{array}$$(43)

and the corresponding detection map is shown in Fig. C.1. The parameters were obtained assuming correlated errors. An extension of our correlation model to inter-observation correlations, for instance arising from the calibration process, is left for future work.

Then, we compute the covariances from the correlations, the VISAMPERR, and the T3ERR from the final GRAVITY pipeline products using Eq. (10) and inject companions with

$$\begin{array}{cc}& f\in [{10}^{-3},{10}^{-2.75},{10}^{-2.5},\dots ,{10}^{-0.5}],\hfill \end{array}$$(44)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{RA}}\in [-30,-25,-20,\dots ,30]\mathrm{mas},\hfill \end{array}$$(45)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{Dec}}\in [-30,-25,-20,\dots ,30]\mathrm{mas},\hfill \end{array}$$(46)

using Eqs. (21) and (22). In order to obtain empirical detection limits when assuming uncorrelated and correlated data, we then repeat the model fitting described in the previous Section.

The fraction of correctly recovered companions as a function of the relative flux of the companion for both scenarios (uncorrelated noise: blue points and correlated noise: orange points) is shown in the right panel of Fig. 7, again overplotted with logistic growth functions fitted to the data points (cf. Sect. 3.2). The plot looks similar to the one from the injection and recovery tests with simulated data and confirms the applicability of our full covariance model to real GRAVITY data. Of course, the empirical detection limits are about one to two orders of magnitude worse and the improvement that comes from our full covariance model (i.e. the lateral shift of the orange curve with respect to the blue curve) is only a factor of ∼2 (consistent with the right panel of Fig. 6 which also shows an improvement by a factor of ∼2) due to weaker correlations being present in the data used for the injection and recovery tests with real data. In summary, our full covariance model still brings a singificant improvement over the convential diagonal covariance model.

4. Conclusions

Correlated noise is placing fundamental detection limits on interferometric data. From on-sky VLTI/GRAVITY data, we extract and illustrate the correlations present in the data and develop an empirical model in order to describe them. This empirical model is sufficiently simple for it to be fitted to the correlations extracted from a single GRAVITY data product and could therefore be directly integrated into the GRAVITY data reduction pipeline and made available to the community as part of the OIFITS 2 file (which has a well-defined standard for providing covariance matrices Duvert et al. 2017).

Then, we evaluate the impact of our full covariance model by performing model fitting and companion injection and recovery tests with both simulated and real GRAVITY data. Our methods are based on χ² = R^T ⋅ Σ⁻¹ ⋅ R minimisation, where we compare the scenarios assuming uncorrelated data (i.e. a diagonal covariance matrix diag(Σ)) and correlated data (i.e. a full covariance matrix Σ following from our empirical correlation model). We show that accounting for the correlations that we find to be present in GRAVITY data could yield to an improvement in the detection limits by a factor of up to ∼3.5 over ignoring them. Moreover, the obtained detection limits (and therefore also potential detections) can be regarded more robust in the former case. We also highlight the problems which arise from ignoring the correlations, as it is done in model fitting pipelines such as LITpro (Tallon-Bosc et al. 2008) and CANDID (Gallenne et al. 2015) so far, and discuss that conventional detection criteria based on χ² statistics are strongly biased towards false positives (i.e. detections that are no true companions).

The empirical correlation model presented in this paper is a simple one-parameter model derived from GRAVITY data, but is arguably also applicable (with small modifications) to other instruments such as VLTI/PIONIER for example. It only treats the correlations between the different observables, but not yet between different frames or targets (such as the science and the calibrator target). We choose this approach in order to enable a simple implementation into existing data reduction and model fitting pipelines. Especially with the increasing availability of computing power, the use of full covariance matrices for describing the correlated noise in interferometric data should become a standard. Collaborations around future instruments should provide estimated data covariances as part of the official data reduction pipelines.

In the future, we aim to compare our empirical correlation model with the data covariances derived from bootstrapping (e.g. Lachaume et al. 2019) and extend our model in order to account for correlations between different frames and targets. Finally, we will re-analyse several marginal detections of companions around Cepheid stars from Gallenne et al. (2013, 2014, 2015) by properly accounting for the correlated noise.

Acknowledgments

The authors would like to thank Jean-Baptiste Le Bouquin for helpful feedback on the GRAVITY data reduction pipeline. M. J. I. acknowledges support from the ESO visitor programme. This work has been partly supported by the Australian Research Council’s Discovery Projects (DP190101477). The manuscript was also substantially improved following helpful comments from an anonymous referee.

References

Appendix A: Collection of correlations

Fig. A.1. Same as Fig. 2, but extracted from three other P2VM-reduced files from programme 60.A-9801(U). From top to bottom: GRAVI.2019-03-29T01-46-28.155_singlecalp2vmred.fits, GRAVI.2019-03-29T01-57-13.182_singlecalp2vmred.fits, GRAVI.2019-03-29T01-59-31.188_singlecalp2vmred.fits.

Fig. A.2. Same as Fig. 2, but showing the correlations of the VISAMP instead of the VIS2, extracted from the three P2VM-reduced files used for the injection and recovery tests with real data in Sect. 3.3 (programme 0101.C-0907(B)). From top to bottom: GRAVI.2018-04-18T08-08-19.739_singlescip2vmred.fits, GRAVI.2018-04-18T08-12-10.749_singlescip2vmred.fits, GRAVI.2018-04-18T08-20-04.769_singlescip2vmred.fits.

Appendix B: Linearised model

If the host star is essentially unresolved (i.e. θbλ ≪ 1) and the companion is at high contrast (i.e. f ≪ 1) one can linearise the VIS2_bin and the T3_bin as a function of the relative flux of the companion.

Consider the complex visibility of the binary model VIS_bin in the case where θ → 0 ⇔ VIS_ud → 1 and f ≪ 1, then the VIS2 of this model is

$$\begin{array}{cc}& |{\mathrm{VIS}}_{\mathrm{bin}}{|}^2\approx \frac{1}{{(1+f)}^2}[{(1+fcos(x))}^2+f^2{sin}^2(x)]\hfill \end{array}$$(B.1)

$$\begin{array}{cc}& =\frac{1}{{(1+f)}^2}[1+2fcos(x)+f^2]\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& =(1-2f+\mathcal{O}(f^2))[1+2fcos(x)+f^2]\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& =1+2fcos(x)-2f+f^2-4f^{2}cos(x)+\mathcal{O}(f^3)\hfill \end{array}$$(B.4)

$$\begin{array}{cc}& =1+f(2cos(x)-2)+\mathcal{O}(f^2)\hfill \end{array}$$(B.5)

and the phase (or argument) of this model (and therefore any linear combination of phases such as the T3) is

$$\begin{array}{cc}& \angle {\mathrm{VIS}}_{\mathrm{bin}}=arctan\frac{\mathrm{Im}{\mathrm{VIS}}_{\mathrm{bin}}}{\mathrm{Re}{\mathrm{VIS}}_{\mathrm{bin}}}\hfill \end{array}$$(B.6)

$$\begin{array}{cc}& \approx \frac{\mathrm{Im}{\mathrm{VIS}}_{\mathrm{bin}}}{\mathrm{Re}{\mathrm{VIS}}_{\mathrm{bin}}}\hfill \end{array}$$(B.7)

$$\begin{array}{cc}& =\frac{-fsin(x)}{1+fcos(x)}\hfill \end{array}$$(B.8)

$$\begin{array}{cc}& \approx f\frac{-sin(x)}{1}\hfill \end{array}$$(B.9)

for x = −2πi(Δ_RAu/λ + Δ_Decv/λ). Hence, in the high-contrast regime, one has VIS2_bin ∝ 1 + f and T3_bin ∝ f.

Appendix C: Detection map for real data

Fig. C.1. Detection map for the GRAVITY data used for the injection and recovery tests in Sect. 3.3. The host star is located in the centre of the map and the cyan circle highlights the position of the best fit companion. Its parameters and reduced chi-squared are shown at the top and the bottom of the map. North is up and east is left.

All Tables

All Figures

Fig. 1. Correlations of the VIS2 (left panel) and the T3 (right panel) for the GRAVITY fringe tracker, extracted from a single P2VM-reduced file. The axes run over the different baselines/triangles, with each individual baseline/triangle comprising five spectral channels. Correlations within the same baseline/triangle are highlighted with red squares and correlations between baselines having a telescope in common are highlighted with orange squares. We note that the correlations are computed from 46 592 individual measurements. Below each panel, the variance of the data and the names of the telescopes forming each baseline/triangle are shown.

In the text

	Fig. 2. Same as Fig. 1, but for the GRAVITY science camera. Each individual baseline/triangle comprises 210 spectral channels. We note that the correlations are computed from 32 individual measurements.
In the text

Fig. 3. Correlations of the VIS2 for the GRAVITY science camera, extracted from a single P2VM-reduced file (upper left panel) and our one-parameter model fitted to them (upper right panel). The bottom panels show the corresponding covariances obtained by multiplying the correlation Cij with the product of the standard deviations σiσj. Correlations/covariances within the same baseline are highlighted with red squares and correlations/covariances between baselines having a telescope in common are highlighted with orange squares.

In the text

	Fig. 4. Same as Fig. 3, but showing the correlations/covariances of the T3 and our one-parameter model fitted to them for the GRAVITY science camera.
In the text

	Fig. 5. Fourier u- and v-baselines of our simulated data extracted from three on-sky observations of GRAVITY using the four UTs over ∼20 min.
In the text

Fig. 6. Left panel: contrast curve (i.e. azimuthal average of the best fit relative flux) for simulated data of an unresolved host star without any companion affected by correlated errors, computed with model fitting assuming uncorrelated data (blue curve) and correlated data (orange curve). Both curves show the mean contrast curve over 100 simulated data sets and the shaded region highlights its standard deviation. The dashed black line shows the ratio of the blue and the orange curve, representing the improvement (i.e. the factor by which the detection limits improve) when using our correlated error model instead of the classical uncorrelated one. Right panel: same, but for the real GRAVITY data introduced in Sect. 2.3.

In the text

Fig. 7. Left panel: fraction of correctly recovered companions as a function of the relative flux of the injected companion from injection and recovery tests with simulated GRAVITY data (cf. Sect. 3.2), assuming uncorrelated data (blue points) and correlated data (orange points) for the model fitting. The blue and the orange curves are logistic growth functions fitted to the data points. Right panel: same, but from injection and recovery tests with real GRAVITY data (cf. Sect. 3.3).

In the text

	Fig. A.1. Same as Fig. 2, but extracted from three other P2VM-reduced files from programme 60.A-9801(U). From top to bottom: GRAVI.2019-03-29T01-46-28.155_singlecalp2vmred.fits, GRAVI.2019-03-29T01-57-13.182_singlecalp2vmred.fits, GRAVI.2019-03-29T01-59-31.188_singlecalp2vmred.fits.
In the text

Fig. A.2. Same as Fig. 2, but showing the correlations of the VISAMP instead of the VIS2, extracted from the three P2VM-reduced files used for the injection and recovery tests with real data in Sect. 3.3 (programme 0101.C-0907(B)). From top to bottom: GRAVI.2018-04-18T08-08-19.739_singlescip2vmred.fits, GRAVI.2018-04-18T08-12-10.749_singlescip2vmred.fits, GRAVI.2018-04-18T08-20-04.769_singlescip2vmred.fits.

In the text

	Fig. C.1. Detection map for the GRAVITY data used for the injection and recovery tests in Sect. 3.3. The host star is located in the centre of the map and the cyan circle highlights the position of the best fit companion. Its parameters and reduced chi-squared are shown at the top and the bottom of the map. North is up and east is left.
In the text