Towards reliable uncertainties in IR interferometry: The bootstrap for correlated statistical & systematic errors

TL;DR

This paper introduces a bootstrap-based method to accurately estimate uncertainties in IR interferometry by accounting for correlated errors and systematics, improving measurement reliability for stellar diameters.

Contribution

The authors develop a bootstrap approach that incorporates correlated statistical errors and systematics, enhancing the accuracy of interferometric measurements.

Findings

Correlations and systematics significantly affect diameter estimates.

Including these factors increases uncertainties by a factor of four.

The method provides more reliable error estimates for under-resolved stars.

Abstract

We propose a method to overcome the usual limitation of current data processing techniques in optical and infrared long-baseline interferometry: most reduction pipelines assume uncorrelated statistical errors and ignore systematics. We use the bootstrap method to sample the multivariate probability density function of the interferometric observables. It allows us to determine the correlations between statistical error terms and their deviation from a Gaussian distribution. In addition, we introduce systematics as an additional, highly correlated error term whose magnitude is chosen to fit the data dispersion. We have applied the method to obtain accurate measurements of stellar diameters for under-resolved stars, i.e. smaller than the angular resolution of the interferometer. We show that taking correlations and systematics has a significant impact on both the diameter estimate and its uncertainty. The robustness of our diameter determination comes at a price: we obtain 4 times larger uncertainties, of a few percent for most stars in our sample.

Tables

Table 1. Table 1: Observing log sorted by science target, calibrator, and observing night. Calibrator characteristics are H 𝐻 H magnitude, spectral type, uniform disc diameter ϑ cal subscript italic-ϑ cal \vartheta_{\mathrm{cal}} , and relative weight in the transfer function calculation, with 1.0 for calibrators very close in time and space (e.g. SCI-CAL-SCI block). This weight does not include the lesser impact of large diameter, large diameter uncertainties, and large visibility dispersion under bad conditions that arises from error propagation. Sky conditions in the direction of the target in H 𝐻 H band are the seeing and the atmospheric phase coherence time τ 0 subscript 𝜏 0 \tau_{0} . Science and calibrators observed close to, or above, the nominal limiting magnitude of the instrument are listed as ≈ \approx  lim and > >  lim. We have fitted H lim = 6.7 − 2.5 log 10 ( seeing in  V ) H_{\text{lim}}=6.7-2.5\log_{10}(\text{seeing in $V$)} to the conservative estimates from the ESO Call for Proposals, albeit it is better when τ 0 subscript 𝜏 0 \tau_{0} is large (e.g. GJ 1061).

science target		calibrator					night	conditions
name	$H$	name	$H$	sp. type	${\vartheta }_{\mathrm{cal}}$	weight	night	seeing	${\tau }_0$	ins. limit
	[mag]		[mag]		[mas]	[MJD]		[arcsec]	[ms]	sci.	cal.
GJ 1	4.73	HD 1434	3.75	K2III	$0.941\pm 0.013$	0.96	56596	$0.68\pm 0.04$	$6.8\pm 0.4$
						0.86	56849	$1.36\pm 0.17$	$3.3\pm 0.4$
						0.92	56850	—	$17.8\pm 0.0$
		HD 190	5.38	G8III	$0.431\pm 0.085$	0.64	56849	$1.13\pm 0.08$	$3.8\pm 0.3$
						0.76	56850	—	$17.8\pm 0.0$
						0.98	56889	$0.75\pm 0.13$	$22.4\pm 3.1$
		HD 214623	4.05	K2/3III	$0.776\pm 0.011$	0.22	56597	$0.92\pm 0.05$	$5.8\pm 0.4$
		HD 215709	5.61	F6V	$0.313\pm 0.062$	0.28	56889	$0.57\pm 0.05$	$28.3\pm 2.6$
		HD 216988	5.17	K0III	$0.468\pm 0.093$	0.22	56597	$0.93\pm 0.04$	$5.4\pm 0.2$
		HD 217681	5.64	K0III	$0.362\pm 0.072$	0.62	56889	$0.58\pm 0.05$	$28.7\pm 2.2$
		HD 224936	4.66	K1III	$0.635\pm 0.126$	0.53	56596	$0.68\pm 0.10$	$6.6\pm 0.9$
						0.95	56890	$0.81\pm 0.12$	$14.7\pm 1.9$
		HD 224949	4.92	K0III	$0.556\pm 0.110$	0.93	56596	$0.75\pm 0.05$	$6.0\pm 0.4$
						0.90	56597	$0.91\pm 0.05$	$5.6\pm 0.3$
						0.84	56849	$1.15\pm 0.07$	$3.7\pm 0.2$
						0.84	56850	—	$17.9\pm 0.1$
						0.97	56889	$0.85\pm 0.17$	$20.9\pm 3.5$
		HD 225101	6.23	F2V	$0.234\pm 0.047$	0.95	56597	$0.95\pm 0.05$	$5.4\pm 0.3$
		HD 32	4.80	G8III/IV	$0.538\pm 0.107$	0.91	56596	$0.80\pm 0.13$	$5.8\pm 0.9$
						0.90	56597	$0.86\pm 0.07$	$5.9\pm 0.4$
						0.89	56890	$0.95\pm 0.17$	$12.9\pm 2.0$
		HD 902	3.77	K3/4III	$0.967\pm 0.013$	0.78	56596	$0.62\pm 0.05$	$7.4\pm 0.6$
						0.93	56849	$1.51\pm 0.09$	$3.1\pm 0.1$
						0.88	56850	—	$17.9\pm 0.0$
GJ 54.1	6.75	HD 22000	7.13	G8III	$0.187\pm 0.037$	0.54	56889	$0.97\pm 0.22$	$21.3\pm 3.8$	$\approx$ lim	$>$ lim
		HD 5248	7.28	K0V	$0.163\pm 0.032$	0.92	56889	$0.84\pm 0.05$	$20.2\pm 1.3$	$\approx$ lim	$>$ lim
		HD 5911	6.51	G3V	$0.216\pm 0.043$	0.89	56597	$1.03\pm 0.17$	$5.8\pm 0.7$	$>$ lim	$\approx$ lim
		HD 5932	6.99	F0V	$0.156\pm 0.031$	0.92	56890	$0.82\pm 0.20$	$13.6\pm 2.8$	$\approx$ lim	$\approx$ lim
		HD 6022	7.15	F2V	$0.154\pm 0.031$	0.94	56889	$0.87\pm 0.03$	$20.2\pm 1.1$	$\approx$ lim	$>$ lim
		HD 6482	3.81	K0III	$0.836\pm 0.012$	0.84	56596	$0.92\pm 0.08$	$6.3\pm 1.0$	$\approx$ lim
		HD 6720	6.11	G8V	$0.288\pm 0.057$	0.84	56542	—	$5.9\pm 0.1$	$\approx$ lim
						0.96	56596	$0.98\pm 0.12$	$4.9\pm 0.6$	$\approx$ lim
		HD 7257	6.58	F3V	$0.210\pm 0.042$	0.84	56542	—	$6.3\pm 0.3$	$\approx$ lim	$\approx$ lim
						0.86	56597	$0.86\pm 0.14$	$6.1\pm 0.8$	$\approx$ lim
		HD 7495	6.26	F6V	$0.235\pm 0.047$	0.97	56596	$0.71\pm 0.06$	$7.4\pm 0.5$
						0.76	56597	$0.90\pm 0.05$	$5.4\pm 0.3$	$\approx$ lim
		HD 8406	6.50	G3V	$0.232\pm 0.046$	0.73	56596	$0.63\pm 0.04$	$6.9\pm 0.4$
						0.91	56597	$0.97\pm 0.05$	$6.6\pm 0.8$	$\approx$ lim	$\approx$ lim
GJ 86	4.25	HD 10939	5.03	A1V	$0.329\pm 0.066$	0.98	56253	$0.63\pm 0.03$	$10.6\pm 0.4$
		HD 12619	3.52	K5III	$1.139\pm 0.016$	0.85	56210	$0.68\pm 0.14$	$7.9\pm 1.4$
						0.98	56252	$0.84\pm 0.22$	$6.8\pm 1.8$
		HD 13666	3.21	K2/3III	$1.007\pm 0.014$	0.95	56252	$0.84\pm 0.23$	$6.5\pm 1.9$
		HD 15371	4.66	B7IV	$0.402\pm 0.080$	0.98	56253	$0.70\pm 0.05$	$8.9\pm 0.6$
		HD 15520	4.61	K1III	$0.624\pm 0.124$	0.85	56210	$0.74\pm 0.12$	$7.4\pm 1.3$
		HD 20640	3.09	K2III	$1.260\pm 0.016$	0.92	56252	$0.92\pm 0.10$	$5.9\pm 0.6$
		HD 21011	3.93	K0III	$0.688\pm 0.009$	0.52	56253	$0.88\pm 0.21$	$7.3\pm 1.8$
		HD 25038	4.29	K2III	$0.812\pm 0.011$	0.27	56253	$0.61\pm 0.04$	$9.3\pm 0.6$
		HD 26934	3.60	K3III	$0.984\pm 0.013$	0.66	56253	$0.74\pm 0.12$	$8.0\pm 1.2$
		HD 22663	2.14	K1III	$1.890\pm 0.022$	0.68	56252	$1.45\pm 0.44$	$3.8\pm 1.1$
		HD 28776	3.63	K0II	$0.837\pm 0.166$	0.90	56253	$0.57\pm 0.10$	$10.9\pm 1.7$

Table 2. Table 4: Tuning of the parameters in the interpolation of the transfer function for all setups (Sect.  3.1.2 , Eqs.  3a ,  3b ). MJD: Modified Julian Day; #: Number of spectral setups; DIT: total integration time during a scan; N fowler subscript 𝑁 fowler N_{\text{fowler}} : number of Fowler (i.e. non-destructive) reads of the infrared detector per scan position; N opd subscript 𝑁 opd N_{\text{opd}} : number of scan positions; τ t subscript 𝜏 𝑡 \tau_{t} : the time-scale to determine the weight of calibrators taken close in time; τ α subscript 𝜏 𝛼 \tau_{\alpha} : the alt+az angular distance scale to determine the weight of calibrators close by in the sky; ε 𝜀 \varepsilon : the minimum relative error in the calibrator visibility considered when determining the weight of a calibrator. The first setup is a good example of an unstable night when only calibrators very close in time to the science observation ( ≲ 20 less-than-or-similar-to absent 20 \lesssim 20  min and 15 deg) have a significant weight. Many nights are stable enough and free of alt+az polarisation effects (e.g. MJD 56384–56850).

MJD	#	DIT	$N_{\text{fowler}}$$\times$$N_{\text{opd}}$	${\tau }_t$	${\tau }_{\alpha }$	$\epsilon$
$56000+$		[s]		[h]	[deg]	[%]
210	1	0.39	4$\times$512	0.3	15	2
	1	0.60	1$\times$1024	0.8	20	1
252	4	all	all	0.8	$+\infty$	1
253	1	0.23	1$\times$512	0.8	20	1
325–326	5	all	all	0.8	$+\infty$	1
384–850	21	all	all	0.8	20	1
889	3	all	all	0.8	$+\infty$	1
890	3	all	all	0.8	20	1

Table 3. Table 5: The different uncertainty models (Sect.  4 ). The first four models used uncorrelated statistical errors, the three following ones correlated statistical errors, and the last ones include correlated systematic errors. Data uncertainties can be obtained from pndrs error propagation, the bootstrap variances, the bootstrap variances plus a noisy correlation matrix representing uncorrelated data (“corr. noise”), or the bootstrap covariances. Baseline-dependent statistical errors of uncertain origin may be “ignored” or “fit” on each baseline so that baseline χ r 2 = 1 subscript superscript 𝜒 2 r 1 \chi^{2}_{\mathrm{r}}=1 . Systematic errors can either be modelled with “correlation” matrix or one can just perform a “rescaling” of all statistical errors to a global χ 2 = 1 superscript 𝜒 2 1 \chi^{2}=1 . For the sake of comparison we also include the incorrect correlation matrix (“PPP corr.”) that may lead to Peelle’s Pertinent Puzzle. Fitted visibilities can either be the mean of the bootstraps (“mean”) or all bootstraps (“bootstraps”). In the bootstrap case, the diameter value and uncertainty are obtained from the median and 1- σ 𝜎 \sigma confidence interval of the obtained distribution. Mathematical expression for uncertainties : either an error vector or a covariance matrix.

Model	Data uncertainties	Baseline-dependent	Systematic	Fitted	Mathematical expression
		statistical errors	errors	visibilities	used for the uncertainties
VAR/PROP	propagation	ignored	rescaling	mean	Eq. (4b)	$\Delta V_d^2$
VAR	variances	ignored	rescaling	mean	Eq. (4d)	$\Delta V_d^2{}_{}{}^{\prime }$
VAR/BS	variances	ignored	rescaling	bootstraps	Eq. (4d)	$\Delta V_d^2{}_{}{}^{\prime }$
VAR/NOISE	variances	ignored	rescaling	mean	Eq. (9b)	${\tilde{\Sigma }}_{d{d}^{\prime }}$
	+ corr. noise
COV	covariances	ignored	rescaling	mean	Eq. (4c)	${\Sigma }_{d{d}^{\prime }}$
COV/BL	covariances	fit	rescaling	mean	Eq. (4c)	${\hat{\Sigma }}_{d{d}^{\prime }}$ with ${\sigma }^{\text{sys}}={\sigma }^{\text{wl}}=0$
COV/BL+BS	covariances	fit	rescaling	bootstraps	Eq. (4c)	${\hat{\Sigma }}_{d{d}^{\prime }}$ with ${\sigma }^{\text{sys}}={\sigma }^{\text{wl}}=0$
SYS/PPP	covariances	fit	PPP corr.	mean	Eq. (4c)	${\hat{\Sigma }}_{d{d}^{\prime }}$ with ${\overline{V}}_d^2=V_d^2$
SYS	covariances	fit	correlation	mean	Eq. (4c)	${\hat{\Sigma }}_{d{d}^{\prime }}$
SYS/BS	covariances	fit	correlation	bootstraps	Eq. (4c)	${\hat{\Sigma }}_{d{d}^{\prime }}$

Table 4. Table 6: Uniform disc diameters of the sample using different model fitting routines of Sect.  4 : ϑ prop subscript italic-ϑ prop \vartheta_{\text{prop}} (uncorrelated least squares using the error propagation of the reduction pipeline), ϑ var subscript italic-ϑ var \vartheta_{\text{var}} (uncorrelated least squares using observed variances), ϑ cov subscript italic-ϑ cov \vartheta_{\text{cov}} (least squares using observed covariances), ϑ bl subscript italic-ϑ bl \vartheta_{\text{{bl}}} (least squares using observed covariance and fitting baseline-dependent statistical errors), and ϑ sys subscript italic-ϑ sys \vartheta_{\text{sys}} (correlated least squares with baseline systematics). For each model, the following relative errors on the visibilities are given: σ st superscript 𝜎 st \sigma^{\text{st}} is obtained from the data, σ b bl superscript subscript 𝜎 𝑏 bl \smash{\sigma_{b}^{\text{bl}}} is derived by modelling an additional baseline-dependent statistical (i.e. uncorrelated) error term, σ sys superscript 𝜎 sys \sigma^{\text{sys}} is obtained by fitting the strength of the systematics (fully correlated on the same baseline). The reduced χ r 2 subscript superscript 𝜒 2 r \chi^{2}_{\mathrm{r}} of the fits is also given. For GJ 581, the reduced chi squared of the COV/BL model is below one, so systematics adjust to zero and σ sys = 0 superscript 𝜎 sys 0 \sigma^{\text{sys}}=0 for SYS.

Star	VAR/PROP			VAR			COV			COV/BL				SYS
	$\vartheta$	${\chi }_r^2$	${\sigma }^{\text{st}}$	$\vartheta$	${\chi }_r^2$	${\sigma }^{\text{st}}$	$\vartheta$	${\chi }_r^2$	${\sigma }^{\text{st}}$	$\vartheta$	${\chi }_r^2$	${\sigma }^{\text{st}}$	${\sigma }_b^{\text{bl}}$	$\vartheta$	${\chi }_r^2$	${\sigma }^{\text{st}}$	${\sigma }_b^{\text{bl}}$	${\sigma }^{\text{sys}}$
	[mas]		[$\%V^2$]	[mas]		[$\%V^2$]	mas		[$\%V^2$]	[mas]		[$\%V^2$]		[mas]		[$\%V^2$]
GL1	0.780(01)	2.52	3.1	0.782(01)	3.43	2.4	0.797(03)	3.44	2.4	0.796(02)	1.08	2.4	2.0	0.794(05)	1.00	2.4	2.0	1.0
GL54.1	0.400(07)	3.33	6.6	0.407(06)	3.31	6.7	0.408(10)	2.95	6.7	0.402(09)	1.63	6.7	4.3	0.397(35)	1.00	6.7	4.3	4.3
GL86	0.662(02)	2.12	4.4	0.648(02)	3.10	4.1	0.683(05)	3.09	4.1	0.678(04)	1.27	4.1	3.8	0.683(15)	1.00	4.1	3.8	2.6
GL229	0.865(03)	5.96	4.9	0.860(03)	6.95	3.5	0.888(05)	3.26	3.5	0.870(04)	1.27	3.5	3.3	0.840(12)	1.00	3.5	3.3	1.6
GL273	0.743(02)	1.92	4.9	0.744(02)	2.87	4.3	0.760(05)	2.80	4.3	0.752(04)	1.14	4.3	2.2	0.763(10)	1.00	4.3	2.2	1.6
GL370	0.522(02)	2.04	2.2	0.518(03)	2.73	2.0	0.536(05)	2.77	2.0	0.532(04)	1.15	2.0	1.4	0.530(12)	1.00	2.0	1.4	1.6
GL406	0.521(04)	2.02	6.5	0.522(04)	2.03	6.5	0.539(08)	2.18	6.5	0.540(07)	1.17	6.5	2.8	0.562(20)	1.00	6.5	2.8	2.7
GJ433	0.390(10)	2.95	4.5	0.363(11)	3.93	3.7	0.426(15)	2.91	3.7	0.422(13)	1.28	3.7	2.5	0.457(30)	1.00	3.7	2.5	2.6
GL447	0.513(07)	1.63	8.1	0.519(08)	2.73	4.8	0.539(16)	2.97	4.8	0.533(12)	1.29	4.8	6.1	0.524(29)	1.00	4.8	6.1	2.3
GL551	1.059(03)	4.04	3.2	1.051(03)	4.62	3.1	1.038(05)	3.29	3.1	1.068(04)	1.17	3.1	3.0	1.066(07)	1.00	3.1	3.0	1.1
GJ581	0.467(04)	1.43	3.7	0.471(05)	2.61	2.7	0.460(08)	2.15	2.7	0.464(07)	0.97	2.7	2.3	0.464(07)	0.97	2.7	2.3	0.0
GJ628	0.634(02)	1.96	1.8	0.636(02)	1.85	1.8	0.645(07)	2.23	1.8	0.646(05)	1.13	1.8	0.9	0.644(14)	1.00	1.8	0.9	1.6
GJ667C	0.390(08)	1.90	4.8	0.371(08)	3.01	3.6	0.401(15)	2.86	3.6	0.411(11)	1.05	3.6	3.1	0.406(14)	1.00	3.6	3.1	0.7
GJ674	0.705(06)	1.10	3.7	0.704(06)	1.21	3.7	0.695(18)	1.70	3.7	0.700(15)	1.07	3.7	1.4	0.720(37)	1.00	3.7	1.4	2.6
GJ729	0.601(03)	1.89	3.6	0.602(03)	1.81	3.7	0.620(08)	2.42	3.7	0.615(06)	1.27	3.7	0.9	0.625(20)	1.00	3.7	0.9	2.4
GL785	0.681(09)	6.64	2.1	0.630(10)	21.5	1.7	0.619(11)	7.04	1.7	0.720(08)	1.82	1.7	2.4	0.677(21)	1.00	1.7	2.4	2.0
GJ832	0.807(04)	2.36	2.5	0.801(04)	3.12	2.2	0.803(05)	2.21	2.2	0.798(05)	1.29	2.2	1.2	0.794(10)	1.00	2.2	1.2	1.0
GJ876	0.680(02)	2.14	4.5	0.683(02)	2.49	4.2	0.692(04)	2.69	4.2	0.691(03)	1.12	4.2	3.2	0.686(09)	1.00	4.2	3.2	1.7
GL887	1.298(01)	4.01	2.9	1.297(01)	5.08	2.5	1.302(02)	4.27	2.5	1.300(02)	1.16	2.5	3.7	1.297(04)	1.00	2.5	3.7	1.4
GL1061	0.292(14)	1.93	3.2	0.313(12)	2.27	2.9	0.332(18)	1.92	2.9	0.332(16)	1.16	2.9	1.9	0.344(34)	1.00	2.9	1.9	1.8

Equations

$$V_{db}^{2\text{raw}}$$

$$\Delta V_d^{2\text{raw}}$$

$$V_{cb}^{2\text{ud}}$$

$$\Delta V_c^{2\text{ud}}$$

$$T_{cb}$$

$$\Delta T_c$$

$$T_{db}$$

$$w_k$$

$$V_{db}^2$$

$$\Delta V_d^2$$

$${\Sigma }_{d{d}^{\prime }}$$

$${\Delta V_d^2}^{\prime }$$

$$\Delta V^2\approx 2(1-V^2){\sigma }^{\text{wl}}.$$

$${\hat{V}}_d^2=E_d^{\text{sys}}(V_d^2+E_d^{\text{bl}})+E_d^{\text{wl}}.$$

$${\mathcal{V}}^{\mathrm{sys}}=\mathrm{Var}(E_d^{\text{bl}})$$

$${\mathcal{V}}_b^{\mathrm{bl}}=\mathrm{Var}(E_d^{\text{sys}}{V}_d^2)$$

$${\mathcal{V}}^{\mathrm{wl}}=\mathrm{Var}(E_d^{\text{wl}})$$

$$(\Delta {\hat{V}}_d^2{)}^2$$

$${\hat{\Sigma }}_{d{d}^{\prime }}$$

$${\varrho }_{d{d}^{\prime }}^{(b)}$$

$${\tilde{\Sigma }}_{d{d}^{\prime }}^{(b)}$$

Full text

We have applied the method to obtain accurate measurements of stellar diameters for under-resolved stars, i.e. smaller than the angular resolution of the interferometer. We show that taking correlations and systematics has a significant impact on both the diameter estimate and its uncertainty. The robustness of our diameter determination comes at a price: we obtain 4 times larger uncertainties, of a few percent for most stars in our sample.

1 Introduction

Long-baseline interferometry consists in recombining the light from several telescopes to measure interference fringes on an astronomical object (Lawson, 2000). The key observable is the visibility, a complex number containing the contrast and phase of fringes obtained on a telescope pair. It was first introduced in radioastronomy (Bracewell, 1958) to generalise Michelson’s term (a synonym for contrast in Michelson & Pease, 1921). In the ideal case, the visibility is the Fourier transform of the object’s image taken at a spatial frequency related to the telescope separation (Bracewell, 1958; Honsberger, 1975; Labeyrie, 1975). In the infrared and optical, though, the atmospheric turbulence shifts the fringes in milliseconds, so that only the square visibility amplitude and partial phase information can be retrieved (Roddier & Lena, 1984), for instance via the closure phase, the sum of phases over a telescope triplet (Roddier, 1986; Cornwell, 1987). A robust determination of the uncertainties on these observables is paramount to ensure confidence on the physical parameters derived from model fitting.

“Statistical” errors—deviations of expected mean zero—on interferometric observables are produced by fast-varying intrinsic, atmospheric, and instrumental effects. In addition to the detector and photon noises, a few sources of errors that have been observed at the Very Large Telescope Interferometer (VLTI) are the differential atmospheric piston (Colavita, 1999; Esposito et al., 2000), imperfect fibre injection (Kotani et al., 2003), mechanical vibrations (Le Bouquin et al., 2011, Sect. 5.6), background fluctuations (Absil et al., 2004), detector efficiency variations due to cooling cycle (Absil et al., 2004), and 50 Hz electronic noise from the power grid (Absil et al., 2004). The relatively large number of sources of uncertainties made it difficult to obtain a reliable assessment of the precision of interferometric measurements. In particular, the theoretical estimate using photon and detector noises is unrealistically low, so most processing software tools need heuristics to provide a better one, for instance using the dispersion of a given data set. One of the common pitfalls is the assumption that the statistical errors on visibility measurements are uncorrelated (Meimon, 2005) and follow a Gaussian distribution. In particular, most public data processing software tools do not determine these correlations, in particular those for the VLTI instruments MIDI111MID-infrared Interferometric instrument (Hummel & Percheron, 2006), AMBER222Astronomical Multi-BEam combineR (Millour et al., 2008), PIONIER333Precision Integrated-Optics Near-infrared Imaging ExpeRiment (Le Bouquin et al., 2011), and GRAVITY (ESO GRAVITY pipeline team, 2018). The same happens with popular model-fitting tools (e.g. Litpro, see Tallon-Bosc et al., 2008) or image reconstruction programmes (e.g. MIRA, see Thiébaut, 2008). Also, the Optical Interferometric FITS (OIFITS v. 1, see Pauls et al., 2005) format did not provide a codified way to document correlations in its first and most used version. Only very recently has an update to the standard given specifications for a covariance matrix (OIFITS v. 2, see Duvert et al., 2017).

The assumption of uncorrelated Gaussian measurement errors could not be further from the truth. For instance, closure phases are not independent (Monnier, 2007). Also, several random effects impact the different spectral channels of a same observation in the same way, such as the blurring of fringes due to the turbulent atmosphere (Lawson, 2000, Sect. 7.5 “Atmospheric Biases”). Finally, all observations of an observing sequence are impacted in the same way by the errors on the calibrators, virtually leading to correlations between all data points (Perrin, 2003), even collected in different runs at different facilities. The assumption that the fringe contrasts or phases follow a Gaussian distribution is not confirmed by experience either (Schutz et al., 2014, in the case of AMBER). Also, when deriving the instrumental transfer function, a weighted average of square visibilities and closure/differential phases is obtained, using a few calibrators observed close to the science targets. In most cases this average does not follow Gaussian distribution, as Perrin (2003) note.

For these reasons, Perrin (2003) proposed an analytic formalism to propagate the non-Gaussian correlated uncertainties of the square visibility amplitudes in an approximate, yet relatively accurate way. Several authors have applied these results to the FLUOR444Fiber Linked Unit for Optical Recombination instrument (at IOTA555Infrared and Optical Telescope Array, then CHARA666Center for High Angular Resolution Array, see for instance Perrin et al., 2004; Absil et al., 2006; Berger et al., 2006) to our knowledge the only one for which correlations have been regularly determined.

In addition to the statistical errors that one can infer from the data and/or noise modelling, there are “systematic” errors that typically plague interferometric data and impact all the data of a given set of observations in a similar and poorly understood way. Part of the biases are removed either theoretically (e.g. group-delay dispersion, see Zyvagin et al., 2003) or calibrated out (e.g. polarisation, Haguenauer et al., 2000) by the measurement of the “instrumental visibility” on stars of known geometry, ideally unresolved ones, with the underlying assumption that it varies slowly enough to be interpolated to science observations with sufficient precision (Hanbury Brown et al., 1974; Perrin, 2003).

However, complete removal does not happen. Colavita et al. (2003) measured $\approx 5\%$ systematic errors on the calibrated squared visibility amplitudes at the Keck Interferometer by observing binaries of known orbital parameters. More recently, high-precision diameter measurements using sufficiently well resolved stars with CHARA (White et al., 2018; Karovicova et al., 2018) were shown to significantly differ (several $\sigma$ and up to 15%) from values previously obtained from under-resolved interferometric observations, leading the authors to conclude that they were plagued with undiagnosed systematics. At VLTI, Le Bouquin et al. (2009) and Kervella et al. (2004a) also tried to identify the origins of the large systematics, discarding the uncertainty on central wavelength (calibration errors are reported to be 0.35 to 0.50% at PIONIER at VLTI and PAVO777Precision Astronomical Visible Observations at CHARA, respectively by Huber et al., 2012; Gallenne et al., 2018), atmospheric jitter and injection efficiency due to seeing, and instrumental variations during the observation ($\kappa$ matrix). Possible sources are a differential polarisation effect (several percents, Le Bouquin et al., 2012) and the way the bias is removed in Fourier space during the PIONIER data processing (see Le Bouquin et al., 2009, in the case of FINITO888Fringe-tracking Instrument of NIce and TOrino). We also stress that calibrators may be an additional source of systematics: for instance, an unsuspected binary with a flux ratio of 1:100 would likely go undetected in closure phase with PIONIER ($\lesssim 2.3$ deg), yet could account for a bias in the squared visibility amplitudes of up to 2%.

While some authors deal with the uncertainty on the calibrators’ diameters as a systematic error (e.g. Creevey et al., 2015; Perraut et al., 2016, but we include them in the correlated “statistical” errors for practical reasons), few take into account other systematics, despite all the evidence. Huber et al. (2012) introduce systematic errors of the order of one percent in order to account for errors in the spectral calibration, but they are still much smaller than the observed and unexplained errors.

In this paper, we present a model that make use of a technique that is relatively easy to implement into existing pipelines, is more accurate than the classical error propagation, and requires no additional analytic developments. Known as the bootstrap method (Efron & Tibshirani, 1993), it consists in randomly selecting interferograms to feed the data reduction software. Repeating it enough times, we generate a sampling of the multivariate probability density function of the square visibility amplitudes and closure phases. Bootstrapping was originally introduced in interferometry by Kervella et al. (2004b) in order to determine the statistical errors for VINCI999VLT Interferometer commissioning instrument at the VLTI. In addition, we treat the systematic errors as additional correlated term whose magnitude is left as a free parameter.

In Sect. 2, we present the data set, acquired for the companion paper by Rabus et al. (2019), that has led us to undertake this work. Sect. 3 details our modelling of the correlated statistical errors and systematic errors. We then show (Sect. 4) how the estimate uniform disc diameters with PIONIER at the VLTI is significantly impacted by the level of detail in the error modelling. We summarise and conclude in Sect. 5.

2 Data set

In a companion paper by Rabus et al. (2019), we have needed to obtain accurate stellar diameters of marginally resolved M dwarfs in order to calibrate the mass-radius relation down to the fully convective regime. While the present paper focuses on the data processing method that leads to reliable uncertainties, we find useful to show our main findings on actual data. Our sample consists of 20 under-resolved late-type stars of the solar neighbourhood, 13 M dwarfs of the original programme by Rabus et al. plus 7 backup targets acquired during the observing campaign.

Table 1 summarises the science and calibrator observations that we have carried out. The observation strategy was to observe each science target with different calibrators and on different nights. Each science observation was bracketed by calibrator observations close in time ($\approx 10$ min) and altitude + azimuth (alt + az) positions (a few degrees when possible). Calibrators were chosen so that as to minimise uncertainties on the calibration of the transfer function, i.e. either unresolved or resolved with a small diameter uncertainty. We performed the selection with the searchCal tool (Chelli et al., 2016) provided by the Jean-Marie Mariotti Center (JMMC). All calibrators have indirect (i.e. non-interferometric) diameter determinations, so they are immune to the biases investigated in this paper. We used the Calibrator stars for 200 m baseline interferometry by Mérand et al. (2005, hereafter MER05, using the method: absolute spectro-photometric calibration) for 15 large ($\sim 1$ mas) K0III-K5III calibrators with a typical precision of the order of 1–2% on diameter and $\lesssim 1$% in visibility calibration, the Catalogue of Calibrator Stars for LBSI by Bordé, P. et al. (2002, using the same method) for one K1III star with a similar level of precision, and the JMMC Stellar Diameters Catalogue (JSDC, Bourges et al., 2017, using the method: photometric calibration) for 65 smaller ($\lesssim 0.5$ mas) calibrators with a typical precision of 10–20% in diameter and $\lesssim 2$% in visibility. Swihart et al. (2017) have showed that the indirect spectro-photometric calibration method for large calibrators is not biased, as their catalogue is consistent both with interferometric measurements and catalogue by MER05. For most of our science stars (15 of 20) we observed small calibrators from JSDC or ones smaller than the target from MER05. For 4 of our science targets (GJ 1, GJ 54.1, GJ 86, GJ 370) we used one or more calibrators from MER05 that are more resolved than the target, but we also included several smaller ones from JSDC to mitigate the possible impact of a large, unexpected error in a calibrator’s diameter: GJ 1 and GJ 54.1, main targets of Rabus et al. (2019), have been observed together with 8 and 9 different JSDC calibrators, respectively, in addition to the 3 and 1 from MER05. While the bracketing calibrators have the largest impact on the calibration of a given science observation, other calibrators taken for other targets on the same night and with the same instrumental setup contribute to some extent to the transfer function, typically if they were taken within an hour of the science target and relatively close in the sky. For any given target and observing night, Table 1 lists all relevant calibrators with their relative weight (see Eq. 3b in Sect. 3.1.2 for the weighing as a function of distance in time and position) as well as sky conditions.

Table 4 gives an overview of the instrumental setups used and the fine-tuning of calibration parameters.

3 Data processing

We call “data set”, with index $d$ ($d$ in $1\ldots N_{\text{data}}$), a set of ${N_{\text{int}}}\sim 10^{2}$ interferograms $\mathcal{I}_{di{}}$, indexed by $i$ ($1\ldots{N_{\text{int}}}$), taken in quick succession with a single telescope pair in a single spectral channel. An interferogram is the temporal scan as function of optical path difference (OPD). Each data set of a science target will result in exactly one calibrated squared visibility amplitude $V^{2}_{d{}}$, which we will refer to as “visibility”. The interferograms of a data set are assumed to share the projected baseline length $u_{d}=\boldsymbol{B}_{d}\cdot{\boldsymbol{\hat{r}}_{d}}/\lambda_{d}$, where $\lambda_{d}$ is the effective wavelength of the spectral channel, $\boldsymbol{B}_{d}$ is the mean separation between the telescopes, and ${\boldsymbol{\hat{r}}_{d}}$ the radial unit vector representing the target’s location in the horizontal coordinate system (elevation and azimuth).

An “observation” consists of several data sets taken at the same time, with different values of $u_{d}$, for several baselines are used at the same time, sometimes also different wavelengths. In the present case, we have used PIONIER with six telescope pairs and one to three spectral channels, so each observations consist of 6 or 18 data sets. In a single telescope pointing, five observations are usually performed in a row (for a total of $5$–10 min). Over the observing runs, we have acquired data for a few dozens of pointing positions per scientific target, collecting of the order of one thousand data sets per star.

A “setup” is a unique combination of configuration of the telescope array (stations used), instrumental setup (spectral dispersion, readout mode, scanning speed) and observing night. Because of the span in stellar brightness among the sources and the varying observing conditions, a few setups are used each night. In a setup, several calibrator stars and one or more scientific targets are observed. All data taken with the same setup can show some level of correlated systematics due to the wavelength calibration error.

We shall call “baseline” a unique combination of a telescope pair and a setup. In particular, we will consider that data sets taken with the same stations with different spectral configurations or on different nights originate from different baselines as they are likely to have different systematic error terms.

The PIONIER data reduction software (hereafter pndrs, see Sect. 5 of Le Bouquin et al., 2011), that we have modified, determines uncalibrated visibilities, computes the instrumental transfer function for calibrators of known diameters, interpolates it for scientific targets taken with the same setup, and derives the calibrated visibilities for these sources.

Our treatment of uncertainties proceeds in four steps: we determine the “statistical” errors that can be inferred from the noise in the data and the uncertainty on the diameters of the calibrators (Sect. 3.1), we model an additional error term to account for the dispersion of the reduced visibilities at each baseline (Sect. 3.2), introduce a highly correlated, systematic error term to account for the discrepancy between the reduced visibilities of different baselines (Sect. 3.3), and model the wavelength calibration error as an additional systematic error term (Sect. 3.4). Sect. 3.5 gives the resulting covariance matrix.

3.1 “Statistical” errors

3.1.1 Raw visibilities

The “raw” visibility (uncalibrated visibility) is determined by correcting the fringe contrast of the interferograms from different atmospheric and instrumental effects such as the finite bandwidth and the flux imbalance between the beams. For the sake of clarity, we will assume that it is obtained separately on each interferogram and then averaged. However, the details of how pndrs computes visibility amplitudes may vary according to setup and processing mode (see Sect. 5 of Le Bouquin et al., 2011, for further details of the data processing). The bootstrap method will work as long as the resulting visibility is a function of a significant number of scans or frames. For instance, it fails for the science detector of GRAVITY with long exposures (ESO GRAVITY pipeline team, 2018), but it should work for all other VLTI instruments.

Uncertainties are determined by the bootstrap method. The bootstraps $\mathcal{B}_{di{b}}$ ($b$ in $1\ldots{N_{\text{boot}}}$) are ${N_{\text{boot}}}$ sets of ${N_{\text{int}}}$ interferograms picked at random with repeats from the original set $\mathcal{I}_{di{}}$. With the exception that the first bootstrap is the original set, i.e. $\mathcal{B}_{di{1}}=\mathcal{I}_{di{}}$, we pick ${N_{\text{int}}}({N_{\text{boot}}}-1)$ independent random numbers $r_{ib}$ uniformly in $1\ldots{N_{\text{int}}}$, so that $\mathcal{B}_{di{b}}=\mathcal{I}_{dr_{ib}}$ for $i\geq 2$. The random numbers $r_{ib}$ are the same for all data sets of the same observation, so that cross-channel and cross-baseline correlations are correctly measured.

We then obtain ${N_{\text{boot}}}$ raw visibilities $V^{2\,\text{raw}}_{d{b}}$ ($1\leq b\leq{N_{\text{boot}}}$) by averaging the visibility $\mathcal{V}^{2}(\mathcal{B}_{di{b}})$ obtained for each interferogram. Together with its uncertainty, it is given by

[TABLE]

In this work, we picked ${N_{\text{boot}}}=5\times 10^{3}$ so that we can derive a covariance matrix for a few thousands of data sets and use the (very slightly biased) sample variance to avoid numerical issues. The original version of pndrs published by Le Bouquin et al. (2011) also bootstraps the interferograms (with ${N_{\text{boot}}}\sim 10^{2}$) to determine the uncertainty, but discards them afterwards. It keeps the value $V^{2\,\text{raw}}_{d{1}}\pm\Delta V^{2\,\text{raw}}_{d{}}{}$ of Eq. (1) and propagates the errors assuming an uncorrelated multivariate Gaussian distribution. We have modified the software to keep all the bootstraps down to the final product and get an empirical sampling of the calibrated visibility distribution.

3.1.2 Instrumental transfer function

There are still instrumental effects, difficult to compute, in the raw visibilities, so that the fringe contrast is lower than expected from a theoretical point of view. To remove them, the transfer function (also known as instrumental visibility) is calculated on unresolved sources or targets of precisely known geometry, observed in the vicinity of the scientific targets. Then it is interpolated for science targets. Some care has to be taken to include the uncertainties on the calibrators’ diameters and variations of the transfer function due to changing atmospheric conditions.

If the star is a calibrator with known uniform disc diameter $\vartheta\pm\Delta\vartheta$, then ${N_{\text{boot}}}$ diameters $\vartheta_{b}$ are picked at random assuming a Gaussian distribution $\mathcal{N}(\vartheta,\Delta\vartheta^{2})$, with the exception that $\vartheta_{1}=\vartheta$. This is done once per calibrator star for all the observing runs, so that correlations from calibrator errors are correctly propagated to the final data.

For data set $c$ corresponding to a calibrator and bootstrap number $b$, the ratio of the raw visibility $V^{2\,\text{raw}}_{c{b}}$ to the theoretical uniform disc visibility $V^{2\,\text{ud}}_{c{b}}$ yields the transfer function $T_{c{b}}$:

[TABLE]

where the theoretical visibility amplitude of a uniform disc is given by the function $\mathcal{V}^{2}_{\text{disc}}(x)=|2J_{1}(x)/x|^{2}$ with $J_{1}$ the Bessel function of the first kind.

If the star is a scientific target, the reduction software interpolates the transfer function using calibrator observations $c_{1},\ldots,c_{n}$ obtained on the baseline (same telescope pair and setup). Calibrator observations close in time and/or position in the sky and with smaller error bars are given more weight. If, for bootstrap number $b$, the transfer functions $T_{c_{1}b},\ldots,T_{c_{n}b}$ are determined at time $t_{c_{1}},\ldots,t_{c_{n}}$ and horizontal coordinates $\boldsymbol{\hat{r}}_{c_{1}},\ldots,\boldsymbol{\hat{r}}_{c_{n}}$, then the estimated transfer function for a science observation $d$ at time $t_{d}$ and position ${\boldsymbol{\hat{r}}_{d}}$ is given by

[TABLE]

where $\alpha(\boldsymbol{\hat{r}}_{1},\boldsymbol{\hat{r}}_{2})$ is the alt + az difference between telescope pointing positions $\boldsymbol{\hat{r}}_{1}$ and $\boldsymbol{\hat{r}}_{2}$. $\tau_{t}$, $\tau_{\alpha}$, and $\varepsilon$ are constants within a given setup. $\tau_{t}$ is the time-scale of variations of the transfer function, typically of the order of 1 hour in the original pndrs software, but we have shortened it for a few very “agitated” nights (see Table 4). $\tau_{\alpha}$ is the angular size on the sky over which the transfer function varies. In the original software, $\tau_{\alpha}=+\infty$ so that calibrators have a weight independent of their position relative to the science target. However, in some nights, polarisation effects led us to use a finite value (see Table 4). $\varepsilon=0.01$ (also used in the original pndrs) is the minimum relative uncertainty we consider in the determination of the weight of calibrator observations, in order to avoid that a calibrator with an unexpectedly low uncertainty biases the transfer function.

The calculation of $\Delta T_{d{}}{}$ includes three terms:

the error propagation using Eq.  (3a), which includes the uncertainty on the diameter of the calibrators; 2. 2.

an interpolation uncertainty taking into account the varying atmospheric conditions between target and calibrator, which we measure by comparing the interpolated transfer function for calibrator observations with the measured one; 3. 3.

an extrapolation uncertainty for data not bracketed by calibrators (which we have avoided).

An example for the first two uncertainty terms of the transfer function is given in Sect. 3 of Nuñez et al. (2017) in the case of weights linear in time. The full expression for $\Delta T_{d{}}{}$ is not given here for it follows exactly the same steps as in the original software.

3.1.3 Calibrated visibilities

The calibrated visibilities $V^{2}_{d{b}}$, their uncertainties $\Delta V^{2}_{d{}}{}$, and their covariances $\Sigma_{d{d^{\prime}}}$ are given by:

[TABLE]

differ to some extent from the ones derived by error propagation ($\Delta V^{2}_{d{}}{}$ in Eq. 4b) in the original software. The reasons are correlations, departure from the Gaussian distribution, and non-linearity in the calculations in the reduction software (in particular the divisions).

3.2 Baseline-dependent “statistical” errors

We have noted that within a given baseline, it frequently happens that the dispersion of the data points is higher than what our carefully deduced uncertainties suggest. Because most of the stars in our programme are bona fide centro-symmetric targets, under-resolved in the $H$ band at VLTI, the uniform disc model (or any symmetric model, see Lachaume, 2003) should fit the data correctly. Thus, we expect the reduced chi squared of a least squares fit to the data of a single baseline to be close to unity. While it is the case on some baselines, it might be significantly higher ($\approx 2$) on others.

It is quite possible that the determination of the instrumental transfer function in quickly changing conditions is not perfect. In particular, the interpolation of the transfer function (Eq. 3a) would fail if conditions changed abruptly. In general, our knowledge on the variation of the transfer function is very limited so we rely on an imperfect, generic smoothing law (Eq. 3b).

We also considered that time-correlations may impact the bootstrap method. In Sect. 3.1.1, the determination assumes that the $\mathcal{V}^{2}(\mathcal{I}_{di{}})$, for $i$ in $1\cdots{N_{\text{int}}}$, are not time-correlated. Unfortunately, the data are too noisy to infer a meaningful auto-correlation function and model its impact on the result. However, we know that the observing cadence ($\sim 1$ s) is significantly slower than the typical interferometric coherence time ($\sim 100$ ms in the IR, see Perrin, 1997; di Folco et al., 2003; Glindemann, 2011) and atmospheric turbulence time $\tau_{0}$ ($\sim 10$ ms in the IR). So, we expect little correlation from the differential atmospheric piston and fibre injection variations. In the case unsuspected sources of temporal correlation did show up, we would expect to underestimate uncertainties. However, we have checked that, even with highly correlated visibilities, the impact is very small: in the case of a correlation coefficient of 0.9 between consecutive interferograms, we simulated batches of 100 measurements, which yielded a bootstrap estimate of $\Delta V^{2\,\text{raw}}_{d{}}{}$ within 10% of the correct value.

In the absence of a clear understanding of these errors, we decide to model them as an uncorrelated additive term $E^{\text{bl}}_{d}$. It has expectancy 0, standard deviation $\smash{\sigma_{b}^{\text{bl}}}V^{2}_{d{}}$, and correlations 0. For each baseline where the reduced chi squared $\chi^{2}_{\mathrm{r}}$ differs significantly from 1101010$\chi^{2}_{\mathrm{r}}$ is expected to be $1.00\pm 0.22$ with a 3-$\sigma$ confidence interval for an observation with two telescope pointings and three spectral channels, but we have found values of 2 to 4 on some baselines., the value of $\smash{\sigma_{b}^{\text{bl}}}$ is adjusted so that a least-squares fit to the data of this baseline (and only these data) yields a reduced $\chi^{2}_{\mathrm{r}}=1$.

3.3 “Systematic” errors

We have tried to avoid the commonest systematic error sources with PIONIER. Any calibrator showing hints of binarity, thus being able to skew the transfer function for baselines parallel to the binary separation, was excluded. Calibrators were chosen close (a few degrees when possible) to the science targets so that the differential polarisation effect is minimised. On nights were this effect was impacting the data processing with the default pndrs parameters, the instrumental transfer function was interpolated in a way that calibrators distant in alt+az position are filtered out (see Sect. 3.1.2 and Table 4).

In spite of these efforts, there is strong evidence of systematics in our data. Visually (see Fig. 6, in particular GJ 86, GJ 1061, GJ 1), the data of some baselines do not align with those at other baselines. This first look is confirmed by a reduced $\chi^{2}>1$ on most of our fits in spite of a detailed modelling of statistical errors (Sects. 3.1 & 3.2).

We model these systematic errors by a multiplicative term $E^{\text{sys}}_{d}$ with high correlations along the same baselines. It has expectancy 1, standard deviation $\sigma^{\text{sys}}$, and correlation coefficients $\varrho\approx 1$ for all observations $d$ of the same baseline, and 0 for observations on different baselines. The value of $\sigma^{\text{sys}}$ is adjusted so that a least squares fit to the full data set yields a reduced $\chi^{2}=1$.

3.4 Wavelength calibration errors

Gallenne et al. (2018) showed that the PIONIER wavelength calibration has a relative uncertainty of $\approx 0.35$% by alternating observations of a well know binary with PIONIER and wavelength-stabilised second-generation VLTI instrument GRAVITY (Eisenhauer et al., 2011). The uncertainty on the central wavelength of spectral channels of a given spectral configuration, $\sigma^{\text{wl}}=\Delta\lambda/\lambda$ produces an uncertainty on the projected baseline $u_{d}$ of the same amount (0.35%). To model it exactly, we should introduce correlated uncertainties along the $x$-axis ($u_{d}\pm\sigma^{\text{wl}}u_{d}$) in addition to those on the $y$-axis ($V^{2}\pm\Delta V^{2}$). Since the model is a continuously differentiable function of baseline, we decide instead to translate this small error in baseline into a small error in visibility. For under-resolved objects, $1-V^{2}_{d{}}\propto u_{d}^{2}$ so that we obtain:

[TABLE]

We model this systematic error source as an additive term $E^{\text{wl}}_{d}$ of zero expectancy, standard deviation $2\sigma^{\text{wl}}(1-V^{2}_{d{}})$ and correlation $\approx 1$ for all observations taken with the same spectral setup on the same night. Because Gallenne et al. (2018, see their Fig. 1) show that the wavelength error varies from one night to the other, we have assumed that correlations are zero for data taken on different nights and/or different spectral configurations.

For observations taken on a single night in a single spectral setup, we can expect a diameter uncertainty of $0.35\%$. However, according to our observation strategy, most of the stars in our sample were observed on different nights, so we expect these systematics to have a lower impact on the final uncertainty determination (0.15 to 0.25%).

3.5 Final uncertainty determination

The final quantity we measure and fit is therefore

[TABLE]

The determination of the uncertainties and the covariance matrix must be done with some care. It is well known that a naïve error propagation in presence of high correlations can lead, and indeed leads as we discovered in our data, to a least-squares fit that falls far off the data and yields a parameter estimation inconsistent with a more careful analysis. The paradox, known as Peelle’s pertinent puzzle (Peelle, 1987), has an easy remedy in the case of a fit by a constant value (Neudecker et al., 2012) or a set of constant values (Neudecker et al., 2014). In their analytic derivation of the covariances, the weighted average of the measurements replaces the individual measurements. We need to generalise their results to suit our needs, since the different visibilities of the same baseline do not have the same expected value ($\mathcal{V}^{2}_{\text{disc}}$ is not constant). As Zhao & Perey (1992) already noted, there is no longer an obvious (analytic) candidate for the weighted average in the case of a non-linear model. We decided to use the most “natural” estimate: we replace the visibilities $V^{2}_{d{}}$ by their estimator $\bar{V}^{2}_{d{}}$ obtained from a least-squares model fit to the data of the given baseline and setup. In the case of a constant model, it would give the same result as Neudecker et al. (2014).

We introduce the following variances related to the baseline-dependent statistical errors $E^{\text{bl}}_{d}$, the systematic errors $E^{\text{sys}}_{d}$, and the wavelength calibration errors $E^{\text{wl}}_{d}$:

[TABLE]

The final error and correlation matrix are given by

[TABLE]

where $\mathcal{B\!L}_{d}$ and $\mathcal{S\!S}_{d}$ are the baseline and spectral setup of data set $d$. In order to prevent a (numerically) singular matrix, we use $\varrho=0.95$ instead of 1. Since errors are relatively small (a few percents) the second-order terms ($<0.1$%) arising when propagating the errors from Eq. 6 have been ignored in Eqs. 8a & 8b.

4 Impact of the uncertainty model

We assess the respective importance of the correlation and the level of detail used in the determination of uncertainties by comparing the results obtained with different error models. These error models, explained below, are briefly summarised in Table 5. The values of the uniform disc diameters obtained from these models are given in Table 6 and Fig. 1. For the most discrepant error models, the deviation between estimates is shown in Fig. 4.

In all models, a least-squares fit of the value of the diameter is performed on the calibrated visibilities $V^{2}_{d{}}$, sometimes for each bootstraps $b$ using $V^{2}_{d{b}}$. When the fits are performed on each bootstrap, the diameter estimate is the median of obtained values and the uncertainty is obtained with the 1-$\sigma$ confidence interval.

We have gathered the different models into three groups. In the first one uncorrelated errors (uncertainties $\Delta V^{2}_{d{}}{}$ or ${\Delta V^{2}_{d{}}{}}^{\prime}$) are assumed and no systematics are taken into account. They are called VAR (“variance”). The second group, uses correlated errors (covariance matrix $\Sigma_{d{d^{\prime}}}$ or $\widehat{\Sigma}_{d{d^{\prime}}}$) but no systematics, they are called COV (“covariance”). The last group of models use both correlated statistical errors and systematics, they are called SYS (“systematics”). Within these groups small differences in the error handling have been considered, they are specified after a slash.

4.1 Error propagation

The two first error models we introduce and compare take neither correlations nor systematics into account. They differ in the way the error propagation is performed from the raw to calibrated visibilities. In the first model, VAR/PROP, the standard quadratic addition of error terms along the way is performed (fit to $V^{2}_{d{}}\pm\Delta V^{2}_{d{}}{}$). The second one, VAR, determines the error from the bootstraps (fit to $V^{2}_{d{}}\pm{\Delta V^{2}_{d{}}{}}^{\prime}$). VAR/PROP and VAR are displayed as the first two points in the left panels of Fig. 1.

If errors are not strictly Gaussian, as expected from the quotient in Eq. (4a), the bootstrap method samples the real distribution of the calibrated visibilities and estimates correctly its average and confidence intervals, while the standard propagation may not.

In most cases (see Fig. 4, top-left panel), the estimates for the diameters using VAR/PROP and VAR are consistent with each other and differ by at most 2-$\sigma$. For two stars (10% of the sample), GJ 785, and GJ 86, the results are at least 3-$\sigma$ apart. We conclude that departure from a Gaussian distribution has not a significant impact in most cases.

4.2 Correlation of statistical errors

The VAR error model group assumes that errors on the visibilities are not correlated ($V^{2}_{d{}}\pm\Delta V^{2}_{d{}}{}$ is fit), while the COV model group does ($V^{2}_{d{}}$ is fit with correlation matrix $\Sigma_{d{d^{\prime}}}$). We expect that positive correlations increase the errors bars on the diameter, as the different data points become partially redundant. They can also shift its estimated value, as groups of correlated data lose their weight relative to uncorrelated data.

This is indeed what we observe for most of our stars (in the left panel of Fig. 1, 2nd and 5th points from the left labelled VAR and COV). The difference is quite significant as 40% of the stars show a discrepancy of at least 2-$\sigma$ (see Fig. 4, middle-left panel).

4.3 Baseline-dependent statistical errors

As explained in Sect. 3.2, the data at some baselines show more dispersion than the statistical errors determined during the data reduction. An additional baseline-dependent error term has been introduced to account for this discrepancy in model COV/BL.

As it can be seen in Fig. 4 (for star GJ 785), including these errors can significantly modify the diameter estimate. When a few baselines are impacted by very high noise of unknown origin, most of the time bad and fast varying atmospheric conditions, the inclusion of an additional error term decreases the weight of these baselines in the fit (for GJ 785, see two baselines around 50 and 70–80 M$\lambda$ in Fig. 6, right column, 3rd panel from the bottom, which show higher dispersion than error bars predict as well as some bias towards a less resolved star). GJ 785 is a K star with no hint of specially high activity, so it is unlikely that this $V^{2}$ dispersion is produced by a short-term variability. A relatively fast atmospheric turbulence ($\tau_{0}=5$–7 ms in $H$ with a decent seeing of 0.7–0.9) may explain part of it. Two other stars of the sample show a high, unexplained dispersion of measurements at several baselines which could be imputed to an intrinsic short-term variability: for M-type flare stars GJ 54.1 and GJ 447, $\smash{\sigma_{b}^{\text{bl}}}=4.3$% and 6.1%, respectively but, in their case, the dispersion only impacts the uncertainty, not the diameter estimate. For the remaining 85% of the sample, the unexplained dispersion is moderate ($\smash{\sigma_{b}^{\text{bl}}}=2.3\pm 0.9$%) and the COV & COV/BL diameter estimates are consistent with each other.

4.4 Systematic errors

Usually, error models assume that the uncertainties on data points are underestimated and rescale them so that a reduced chi squared of 1 is obtained. This is indeed what we do in our models that exclude systematic errors (VAR and COV model groups).

In the SYS model group, however, no rescaling of error bars occurs. Instead a relative error term $\sigma^{\text{sys}}$ is introduced, as explained in Sect. 3.3 and the fit is performed with a correlated covariance matrix ($\widehat{\Sigma}_{d{d^{\prime}}}$ in Eq. 8b). The value of $\sigma^{\text{sys}}$ is fit to a reduced chi squared of one.

The difference in estimated diameter is often significant as 25% of the stars show a discrepancy of at least 1-$\sigma$ (see e.g. Fig. 4, bottom-left panel) between the COV and SYS models.

In Fig. 1, one can also see (left panel, 6th and 9th points from the left, labelled COV/BL and SYS) that systematics tend to increase the diameter uncertainties (as positive correlation usually do). A comparison of uncertainties in the VAR and SYS models is given in Fig. 5 as a function of the resolution factor, which we define as the ratio of the UD diameter of the star to the nominal resolution of the interferometer $\lambda/B_{\mathrm{max}}$. The typical error in SYS is 4.8% for a resolution factor of 0.2 ($\approx 0.45$ mas diameter at 140 m baselines in the $H$ band) with respect to 1.2% in the standard VAR model. Since the VAR model is widely used in the literature, we conclude that diameter uncertainties may be significantly underestimated. The relative uncertainty scales approximately as the inverse of the 2.5th power of the resolution, making determinations for under-resolved objects ($\lesssim 0.4$ mas at VLTI in $H$) difficult. Since the fringe contrast loss $1-V^{2}$ scales as the square of the diameter, we expected the diameter uncertainty to scale as the inverse of the second power of the resolution factor. The additional drop in precision may be attributed to a loss of precision as the stars get fainter. In Fig. 1, we can see that 3 of the 4 most under-resolved targets have been observed very close to the sensitivity limit of PIONIER.

For reference, Fig. 1 shows another error model, SYS/PPP, that uses the naïve but erroneous correlation matrix leading to Peelle’s pertinent puzzle. Indeed, the diameter estimates are significantly off for about half of the sample.

4.5 Noise in the covariance matrix

With model VAR/NOISE we aim to assess the influence of the noise in the correlation to disentangle noise systematics from actual correlation influences. We measure the noise level in the covariances $\Sigma_{d{d^{\prime}}}$, generate noisy matrices for uncorrelated data, and perform a data fit using these matrices.

Most data are only lightly correlated, for different nights and instrument configurations are loosely correlated by the small uncertainty on the calibrators diameter. So, the histogram of the correlation matrix values features a central component, a Gaussian-like distribution around a small positive constant, corresponding to the mostly uncorrelated data. In addition it has a bump with large positive correlations corresponding to the small fraction of highly correlated data (spectral channels of the same observation, for instance). We fit the width of the central component to get the noise level.

Then, we generate ${N_{\text{boot}}}$ covariance matrices for noisy uncorrelated data

[TABLE]

where $M_{kd}^{(b)}$ are picked using independent normal Gaussian distributions and $\varepsilon$ is adjusted to the noise level. The largest dimension along $k$ that allows $\varrho_{d{d^{\prime}}}^{(b)}$ to be positive definite is used.

For each of these noisy covariance matrices, we perform a least-squares fit to all calibrated visibilities $V^{2}_{d{}}{}$ taken at baselines $u_{d}$. Let $\vartheta_{\text{corr noise}}^{(b)}$ ($b$ in $1\cdots{N_{\text{boot}}}$) be the values of the UD diameters. The UD diameter estimate $\vartheta_{\text{corr noise}}$ is obtained from their average. Its uncertainty $\Delta\vartheta_{\text{corr noise}}$ is the quadratic sum of the model uncertainty $\Delta\vartheta_{\text{var}}$ and the dispersion of $\vartheta_{\text{corr noise}}^{(b)}$.

As we can see from the data in Fig. 1 (left panels, second and fourth point), there is no significant difference between the VAR and VAR/NOISE models for any of our stars. We conclude that the impact of correlations is a real effect, not a bias introduced by noisy data.

4.6 Non-Gaussian uncertainties

In order to assess the influence of non Gaussian uncertainties, we performed least squares fit to each bootstrap and look at the distribution of the estimates for the uniform disc diameter. We expect that a distribution of visibilities with significant skewness, kurtosis, or long tails would be reflected in the distribution of diameter estimates, and, in turn, in its mean value and/or uncertainty.

For model VAR/BS, $V^{2}_{d{b}}\pm\Delta V^{2}_{d{}}{}$ ($1\leq b\leq{N_{\text{boot}}}$) is fitted, yielding a set of uniform disc diameter estimates $\vartheta_{\text{var bs}}$. The median value and 1-$\sigma$ confidence interval of $\vartheta_{\text{var bs}}$ yields the VAR/BS diameter estimate $\vartheta_{\text{var}}\pm\Delta\vartheta_{\text{var}}$. The same is performed for error models with correlated statistical errors (COV/BL+BS) and correlated statistical and systematic errors (SYS/BS).

Figure 1 (left panel) clearly shows that there is no significant difference in the diameter estimate for VAR and VAR/BS, COV/BL and COV/BL+BS, and SYS and SYS/BS.

4.7 “Bad baseline” scenario

We additionally check for the scenario that a “bad baseline”, e.g. tainted with a strong systematic, may significantly alter the value of the uniform disc estimate. Given that baseline, we fit the data that are taken at any other baseline and compare the diameter estimates with the one obtained with all baselines. The process is repeated for each baseline. The diameter estimates with one baseline removed are reported on the right side of Fig. 1.

Two examples o. “bad baselines” can be seen in Fig. 6. For GJ 876, a baseline around 75 to 80 megacycles displays very high dispersion ($V^{2}$ between 0.5 and 1.2) despite small error bars ($\sim 0.05$). This issue does not originate from a bad calibrator, because this kind of dispersion would only be seen with a well resolved, low contrast binary. That would clearly been seen (1) in closure phases which was not the case and (2) as a strong bias towards a more resolved object with most points below the fit. The inclement weather during this particular observation is a likely explanation with a seeing below average (1.05 arcsec in $H$) and a short coherence time (6 ms in $H$). Interestingly, this “bad baseline” has very little impact on the diameter estimates as GJ 876 has all its estimates (VAR, COV, COV/BL, SYS models) within less than 1% dispersion, probably because it shows little bias (the data approximately average the $V^{2}\approx 0.85$ estimate). Another example of “bad baseline” has already been mentioned for GJ 785 (less resolved around 50 and 70-80 M$\lambda$ in Fig. 6, left column, 3rd panel from the bottom), this time with a clear bias. Its origin is unknown, but it is unlikely to be a calibrator issue either: if there was a large error on the diameter of the unique, large calibrator HD 196387, there would be an unseen error on GJ 785’s diameter but very little additional dispersion between baselines. The calibrator has no measured closure phase ($\lesssim 0.02$ rad) and no detectable near-infrared excess (no evidence from JHK modelling in McDonald et al., 2012), so binarity or circumstellar material are very unlikely to account for a 5% error in the visibility calibration. The bias induced by the “bad baseline” reflects in Fig. 1, right panel: with this baseline removed (3rd point of the 6), the estimates of the VAR, COV/BL, and SYS models are consistent, but with it (left panel) they differ by several $\sigma$.

These two “bad baselines” are extreme cases that could be solved easily by removing the data from the fit. However, there are targets (GJ 86, GJ 229, GJ 447, GJ 551) where the visibility plot (Fig. 6) does not show anything obviously wrong, but the analysis of removing one baseline in the fits produces a significant difference in the VAR and/or COV/BL diameter estimates (right panel of Fig. 1).

When systematic errors are included, the impact of “bad baselines” disappears completely (right panel of Fig. 1): all diameter fits with one baseline removed are consistent with each other. The price to pay for this stability is larger error bars as we have shown in Sect. 4.4 and Fig. 5.

It appears, therefore, that observations with a few configurations over a few nights are enough to prevent a single baseline to significantly bias the diameter estimate as long as systematic errors are modelled.

4.8 Bad calibrator scenario

For several of our targets (GJ 1, GJ 54.1, GJ 86, GJ 370) one or more calibrators from MER05 more resolved than the target were used together to smaller, unresolved ones. In the bad luck scenario that a resolved calibrator’s diameter is off by several standard deviations, it would bias the transfer function at the two or three large baselines it is observed with. Fortunately, for these four stars, there is no evidence that a few baselines skew the diameter estimate. The difference in the COV estimates when one baseline removed from the fit is $<1$% (right panel of Fig. 1) and the COV and SYS estimate are within 1% of each other. We are confident that the diameter estimates of these stars are not significantly impacted by a bad calibrator.

5 Conclusion

Many astronomical objects (young stellar objects, late-type dwarfs) are too faint in the visible to be observed by optical interferometers and require kilometric baselines to be fully resolved in the infrared. For this reason, we are bound to infer the object’s properties from under-resolved observations. In the case of stellar diameters, we are fortunate enough that there are no degeneracies in the (unique) parameter estimation, but the under-resolved character has a strong impact in terms of precision. While the geometric size of a fully resolved object can be estimated within a fraction of a percent, uncertainties and systematics of 5–10% on a single observation are common when the target is under-resolved.

Our study had the main objective to partially overcome these limitations by a careful observation layout and to better quantify the remaining uncertainties. To that end, our SYS model fully takes into account correlated statistical uncertainties and systematics.

The main results are

• •

The error model has a significant impact on the estimate of the uniform disc diameter and its uncertainty. Correlations between visibilities have a strong impact (Sects. 4.2) as well as baseline systematics (Sect. 4.4). In both cases, the estimates may differ by more than three standard deviations. The additional errors atop the statistical errors determined by the data processing software are $3.3\pm 1.2$ % on the square visibilities (mean and dispersion in our 20 surveyed stars), which is in line with the generally accepted value of 5% (Colavita et al., 2003). However, the purely systematic term (highly correlated errors) is only $1.8\pm 0.9$ %.

• •

A few observations with different configurations and/or nights are usually enough to avoid a significant bias by a single baseline and instrumental configuration (Sect. 4.7) provided that systematic errors are taken into account. It confirms the usual observation strategy by, for instance, Boyajian et al. (2012), Gallenne et al. (2012), and Rabus et al. (2019) in the case of stellar diameters of under-resolved stars.

• •

Departure from a Gaussian distribution has no significant impact (Sect. 4.6) except indirectly in the error propagation (Sect. 4.1). It needs not to be modelled provided that the uncertainties on the calibrated visibilities are obtained from the bootstraps (this work) or from an analytic determination of the probability density function (Perrin, 2003).

• •

The uncertainty on the diameter is four times larger than the ones modelled with a standard least-squares fit to uncorrelated data. For a diameter of 0.45 mas in $H$ band at 140 m baselines (typical of VLTI/PIONIER), our typical uncertainty is 4.8% (SYS model) with respect to 1.2% with the usual determination (VAR/PROP and VAR models).

We have offered here a relatively easy way, albeit numerically intensive, to obtain correlations between observables, by means of bootstrapping. Even if most reduction pipelines do not propagate correlations, it is possible to run them a large number of times on randomised data (interferograms and calibrator diameters are picked at random) to obtain the multivariate probability density function of the interferometric observables. Systematics errors such as a bad calibrator not bad enough to be detected, rapidly varying atmospheric conditions between science target and calibrators, or instrumental systematics (e.g. differential polarisation) have typically cast a doubt on the robustness of parameter estimation. We have provided a method to deal with these systematics in a relatively inexpensive way: the covariance matrix of the least-squares fit is modified to include a relative systematic error term. The price to pay for a more robust estimation (accurate, i.e. non-biased) is a significantly larger uncertainty.

Therefore, we strongly recommend that future interferometric studies take into account correlated errors and take time to model systematics.

Acknowledgements

Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO IDs 090.D-0917, 091.D-0584, 092.D-0647, and 093.D-0471. R.L., M.R., A.J., and R.B. acknowledge support from CONICYT project Basal AFB-170002. A.J. acknowledges support from FONDECYT project 1171208, BASAL CATA PFB-06, and project IC120009 “Millennium Institute of Astrophysics (MAS)” of the Millenium Science Initiative, Chilean Ministry of Economy. R.B. acknowledges additional support from project IC120009 “Millenium Institute of Astrophysics (MAS)” of the Millennium Science Initiative, Chilean Ministry of Economy. This work made use of the Smithsonian/NASA Astrophysics Data System (ADS) and of the Centre de Données astronomiques de Strasbourg (CDS). This research made use of Astropy, a community-developed core Python package for Astronomy (Price-Whelan et al. 2018). We would like to thank the referee for carefully reading our manuscript and for giving constructive comments which substantially helped improving the paper.