© ESO, 2012

1. Introduction

Astronomical interferometers combine the light from multiple telescopes to form interference fringes. The amplitude and phase of the fringes corresponding to interference between a given pair of telescopes can be combined into a single complex number, the fringe complex amplitude, which in an ideal interferometer is proportional to a single Fourier component of the object brightness distribution.

The complex amplitude forms the fundamental observable in phase-stable interferometric experiments, for example at radio wavelengths, but at optical wavelengths the phase of the fringes is unstable on timescales of order milliseconds, because of optical pathlength perturbations introduced by the Earth’s atmosphere. As a result, observers need to combine the information from multiple short-exposure interferograms to recover the desired science data. Since the phase of the fringes is randomly fluctuating, simply averaging the complex amplitude over multiple interferograms would lead to a vanishingly small result. Instead two observables, the power spectrum and bispectrum, are calculated from each interferogram and these can be averaged over frames as both are immune (to first order) to the atmospheric and instrumental phase perturbations. The power spectrum serves to capture information about the modulus of the complex amplitude, whilst the phase of the bispectrum (also known as the “closure phase”) captures information about the phase of the complex amplitude on closed triplets of baselines.

At the low light levels typically encountered in astronomical interferometry, the effects of the noise arising in the detection process become important. This noise can usually be modelled as a mixture of Poisson noise, for example photon noise, and signal-independent Gaussian noise, for example readout noise in the detector electronics. We use the term “detection noise” to refer to the combination of all the noise processes in the detection process including the photon noise.

Detection noise produces random measurement errors which can be reduced by averaging over many interferograms, but in the case of the power spectrum and bispectrum it also produces systematic biases in the averaged estimates. There has been considerable study in the literature on these biases (see for example Goodman & Belsher 1976a,b; Sibille et al. 1979; Wirnitzer 1985; Beletic 1989; Perrin 2003). Some authors have derived formulae for the power spectrum bias which allow for both Poisson and Gaussian noise sources, but existing bispectrum bias estimates have been derived under the assumption of Poisson noise only. As a result, it is difficult to derive an unbiased estimate of the closure phase from observations where both Poisson and Gaussian noise are significant. Such conditions are increasingly present in interferometric observations made at near-infrared wavelengths so the absence of such an estimator is becoming a critical limitation to the acquisition of accurate interferometric data as we will demonstrate in Sect. 8.

Much of the existing power spectrum and bispectrum analysis relies on the direct Fourier transform of the interferogram to recover the fringe complex amplitudes. Even under noise-free conditions the use of a Fourier-based estimator will return inaccurate estimates of the fringe complex amplitude if the interferogram is, for example, non-uniformly sampled. To account for this, Tatulli & Le Bouquin (2006) introduced a pixel-to-visibility transformation matrix (P2VM) for deriving fringe amplitude estimates in circumstances where the Fourier-based estimators are inaccurate, which includes the Fourier estimators as a special case.

There has been considerable study of the variance of the power spectrum and bispectrum (Goodman & Belsher 1976a,b; Dainty & Greenaway 1979; Sibille et al. 1979; Wirnitzer 1985; Roddier 1987; Ayers et al. 1988), but only the work of Tatulli & Le Bouquin (2006) applies the more general P2VM approach, and in this case the variance is calculated under assumptions valid only for purely Gaussian noise processes. Knowing the variance, and hence the signal-to-noise ratio is important for predicting the integration times and detection limits prior to an observation run.

This paper builds on the introduction by Gordon et al. (2010) and previous work by Thorsteinsson & Buscher (2004) to present a general formalism for computing biases and variances in the presence of detection noise. Our formalism has two advantages over current treatments: (a) it can be applied when the fringe complex amplitudes are derived using any linear estimator, including the P2VM and Fourier transform as special cases and (b) it yields accurate answers for any additive mixture of Poisson noise and Gaussian noise. We derive formulae for the biases affecting both the power spectrum and bispectrum under these conditions and show how our results allow the construction of bias-free estimators. We present simulations to show that using our new bispectrum estimator can significantly enhance the accuracy of measuring the closure phase under realistic observation conditions. We also use our formalism to derive expressions for the variance of the power spectrum, and explain the method by which the reader can calculate the bispectrum variance if required.

2. The interferometric equations

2.1. The interferogram

We will start from a specific idealised example and use this to introduce our more general model which will be applicable to a wide range of realistic scenarios. In our example, the light from N_tel telescopes is combined to produce a 1-dimensional interference pattern, which consists of the superposition of a set of sinusoidal fringes at different fringe frequencies f^ij corresponding to the interference between telescopes i and j. The interference pattern is sampled to produce a set of N_pix intensity measurements  {i_p} . We shall call this discrete set of intensity measurements the “interferogram”. In our example the interferogram is given by: ${\mathit{i}}_{\mathit{p}}\mathrm{=}\frac{{\mathit{C}}^{\mathrm{00}}}{{\mathit{N}}_{\mathrm{pix}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{pix}}}\sum _{\mathit{i}\mathit{<}\mathit{j}}^{{\mathit{N}}_{\mathrm{tel}}}\mathrm{Re}\left[{\mathit{C}}^{\mathit{ij}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}\right]\mathrm{+}\mathrm{noise}\mathit{,}$(1)where C^ij denotes the complex amplitude of the fringes corresponding to interference between light from telescopes i and j, C⁰⁰ denotes an offset due to the light from all the telescopes (we shall call this the “zero-spacing flux”), and α_p is the coordinate of pixel p.

Recognising that Eq. (1) represents a sum over a set of Fourier components, we can derive an estimate c^ij of the complex amplitudes from the sampled intensity data using a discrete Fourier transform (DFT): ${\mathit{c}}^{\mathit{ij}}\mathrm{=}\sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{pix}}}{\mathit{i}}_{\mathit{p}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}\mathit{,}$(2)where we note that this equation can also be used to estimate the zero-spacing flux by setting f⁰⁰ = 0. In noise-free conditions, the estimator c^ij will be equal to the true amplitudes C_ij provided that the samples are evenly spaced in the coordinate α and the frequencies f^ij are integer multiples of the fundamental frequency given by f_fundamental = 1/ [N_pix(α₂ − α₁)] . We will call this set of conditions the “DFT conditions”.

We can rewrite the above equations in the form: ${i}\mathrm{=}{W}{C}\mathrm{+}\mathrm{noise}\mathit{,}$(3)and: ${c}\mathrm{=}{H}{i}\mathit{,}$(4)respectively where i =  [i₁,i₂,...]  is a vector containing the pixel intensities, C =  [C⁰⁰, Re {C¹²}, Im {C¹²}, Re {C¹³}, Im {C¹³}, ...] is a vector describing the true complex amplitudes, c =  [c⁰⁰,Re {c¹²}, Im {c¹²}, Re {c¹³}, Im {c¹³}, ...] is a vector describing the estimated complex amplitudes, W is a matrix known as the “visibility-to-pixel matrix” (V2PM) and H is the pixel-to-visibility matrix (P2VM) mentioned before. When the DFT conditions hold, the matrices W and H will consist of sine and cosine terms.

If we now assume that the DFT conditions do not hold, for example if the samples are not evenly spaced or there are “dead” pixels, then it is often still possible to set up a linear set of equations in the form of Eq. (3) by modifying the matrix W. For the case of any “working” beam combiner, we can always find the corresponding matrix H such that: ${H}{W}\mathrm{=}{I}\mathit{,}$(5)where I is the identity matrix, i.e. Eq. (5) is the condition that H is a pseudo-inverse of W. It is then clear that Eq. (4) can be used to derive a value for c such that c = C under noise-free conditions. The pseudo-inverse generalises the Fourier inverse to the cases in which the DFT conditions may not hold; in general the pseudo-inverse matrix will not consist of simple Fourier terms.

Equation (4) serves to define the basis of our more general model for an interferometric measurement. We assume only a discrete set of intensity measurements  {i_p} related to a set of complex fringe amplitudes  {C^ij} . We do not require a simple Fourier relationship between the amplitudes and the pixel intensities; we only require that there exists a set of estimators for  {C_ij} which are linear in  {i_p} . In other words we require only that there exists at least one set of complex weights $\mathrm{\left\{}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathrm{\right\}}$ such that the complex amplitude estimators: ${\mathit{c}}^{\mathit{ij}}\mathrm{=}\sum _{\mathit{p}}{\mathit{i}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathit{,}$(6)have the property that c^ij = C^ij in the limit where detection noise can be neglected. This generalised model can be applied to a wide variety of beam combination scenarios: the interferograms can be sampled temporally or spatially or a combination of the two, and multiple telescopes can be combined “all-in-one” into a single multiplexed fringe pattern or pairwise onto separate detectors.

Note that pseudo-inverse techniques are only presented as one possible route to deriving $\mathrm{\left\{}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathrm{\right\}}$; under circumstances where the forward matrix W is not known or not fixed (for example if the fringe pattern has not been spatially filtered) alternative techniques to derive $\mathrm{\left\{}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathrm{\right\}}$, perhaps based on knowing the statistical properties of W, may be needed. Nevertheless, the results presented here will be valid independent of the technique used to derive values for ${\mathit{H}}_{\mathit{p}}^{\mathit{ij}}$.

2.2. Interferometric observables

Under phase-unstable conditions, it is common to choose the power spectrum, S^ij and bispectrum, B^ijk as the interferometric observables: ${\mathit{S}}^{\mathit{ij}}\mathrm{=}{\left|{\mathit{C}}^{\mathit{ij}}\right|}^{\mathrm{2}}\mathit{,}$(7)and: ${\mathit{B}}^{\mathit{ijk}}\mathrm{=}{\mathit{C}}^{\mathit{ij}}{\mathit{C}}^{\mathit{jk}}{\mathit{C}}^{\mathit{ki}}\mathit{.}$(8)Typically the values of C^ij will change randomly from interferogram to interferogram due to atmospheric disturbances, but the values of S^ij and B^ijk will be more stable. However, unless otherwise specified, we do not assume any particular distribution for the values of C_ij; for example the results will apply equally well under incoherent conditions where the average value of C_ij is zero or conditions where a fringe tracker is being used and the value of C_ij is stable from frame to frame. Our aim is to determine suitable estimators, ${\mathit{S}}_{\mathrm{0}}^{\mathit{ij}}$, ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$, for the power spectrum and bispectrum respectively that do not suffer from bias in the presence of noise. For this analysis these observables are calculated on an interferogram-by-interferogram basis and then averaged. Derived quantities, such as the closure phase, are calculated using the averaged observables, rather than on an interferogram-by-interferogram basis themselves.

3. Noise model

The signal recorded by a real detector system will always have some noise component, which becomes especially important at low light levels. We model two distinct sources of noise, photon shot noise and detector read noise.

3.1. Photon noise

Photon noise is always present in a real interferogram as a result of the quantum nature of the detection process which converts the incoming photons to an electrical signal. For a given pixel p, the probability of N_p photoevents occurring is given by the Poisson distribution: ${\mathit{P}}_{\mathrm{Poisson}}\mathrm{\left(}{\mathit{N}}_{\mathit{p}}\mathrm{\right|}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{)}\mathrm{=}\sum _{\mathit{n}}^{\mathrm{\infty }}\mathit{\delta }\mathrm{(}{\mathit{N}}_{\mathit{p}}\mathrm{-}{\mathit{n}}^{\mathrm{)}}\frac{{\mathrm{\Lambda }}_{\mathit{p}}^{{\mathit{N}}_{\mathit{p}}}}{{\mathit{N}}_{\mathit{p}}\mathrm{!}}\exp \mathrm{[}\mathrm{-}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{]}\mathit{,}$(9)where Λ_p is the ideal noise-free intensity, also called the “classical intensity” or “high-light-level intensity”. Note that this equation implicitly defines the normalisation of Λ_p in units of photoevents, and we require that i_p = Λ_p under noise-free conditions, thereby defining the normalisation of i_p.

3.2. Read noise

Read noise is present in all optical and IR detectors and is often significant at low light levels. Assuming the read noise is independent of the number of photoevents, the probability of a given measurement error ϵ_p = i_p − N_p is given by the Gaussian distribution: ${\mathit{P}}_{\mathrm{Gaussian}}\mathrm{\left(}{\mathit{ϵ}}_{\mathit{p}}\mathrm{\right|}{\mathit{\sigma }}_{\mathit{p}}\mathit{,}{\mathit{N}}_{\mathit{p}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}}\exp \left[\mathrm{-}\frac{{\mathrm{\left(}{\mathit{ϵ}}_{\mathit{p}}\mathrm{\right)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}\right]$(10)where ${\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}$ is the variance of the noise on pixel p; note that the variance is not assumed to be the same on all pixels.

3.3. Other sources of noise

There are additional noise sources that follow the same statistics as photon or read noise. Detector dark current, D_p, for example is also described by the Poisson distribution and could be included in our formalism by replacing Λ_p with Λ_p + D_p in Eq. (9). For the sake of readability dark current will not be explicitly included in this paper.

3.4. Combined noise model

If the photon noise and read noise are statistically independent then the probability of obtaining the noisy data, i_p, from the underlying ideal data, Λ_p, is given by: $\mathit{P}\mathrm{\left(}{\mathit{i}}_{\mathit{p}}\mathrm{\right|}{\mathrm{\Lambda }}_{\mathit{p}}\mathit{,}{\mathit{\sigma }}_{\mathit{p}}\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{i}}_{\mathit{p}}}{\mathit{P}}_{\mathrm{Poisson}}\mathrm{\left(}{\mathit{N}}_{\mathit{p}}\mathrm{\right|}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{)}{\mathit{P}}_{\mathrm{Gaussian}}\mathrm{\left(}{\mathit{ϵ}}_{\mathit{p}}\mathrm{\right|}{\mathit{\sigma }}_{\mathit{p}}\mathit{,}{\mathit{N}}_{\mathit{p}}\mathrm{)}\mathrm{d}{\mathit{N}}_{\mathit{p}}\mathit{.}$(11)The combined noise distribution is a convolution of the Poisson noise and Gaussian noise distributions, so it is convenient to derive the moments of our noise model via the moment-generating function of the probability distribution. The moment-generating function for the combined noise model probability density $\mathit{P}{}^{\mathrm{(}}\mathit{i}_{\mathit{p}}\mathrm{|}{\mathrm{\Lambda }}_{\mathit{p}}\mathit{,}{\mathit{\sigma }}_{\mathit{p}}{}^{\mathrm{)}}$ is proportional to the product of the moment-generating functions M_Poisson and M_Gaussian of the individual distributions: $\mathit{M}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}\mathrm{=}{\mathit{M}}_{\mathrm{Poisson}}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}{\mathit{M}}_{\mathrm{Gaussian}}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}\mathrm{=}\exp \left[\frac{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{\nu }}^{\mathrm{2}}}{\mathrm{2}}\right]\exp \mathrm{[}\mathrm{-}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathit{\nu }}\mathrm{]}\mathit{.}$(12)The first four moments of the noise model as required to derive the results given in this paper are: $\begin{array}{ccc}\mathrm{⟨}{\mathit{i}}_{\mathit{p}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{d}}{\mathrm{d}\mathit{\nu }}\mathit{M}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}{|}_{\mathit{\nu }\mathrm{=}\mathrm{0}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{p}}\mathit{,}\\ \mathrm{⟨}{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\mathit{\nu }}^{\mathrm{2}}}\mathit{M}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}{|}_{\mathit{\nu }\mathrm{=}\mathrm{0}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathit{,}\\ \mathrm{⟨}{\mathit{i}}_{\mathit{p}}^{\mathrm{3}}\mathrm{⟩}\mathrm{=}\frac{{\mathrm{d}}^{\mathrm{3}}}{\mathrm{d}{\mathit{\nu }}^{\mathrm{3}}}\mathit{M}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}{|}_{\mathit{\nu }\mathrm{=}\mathrm{0}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{3}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathit{,}\\ \mathrm{⟨}{\mathit{i}}_{\mathit{p}}^{\mathrm{4}}\mathrm{⟩}\mathrm{=}\frac{{\mathrm{d}}^{\mathrm{4}}}{\mathrm{d}{\mathit{\nu }}^{\mathrm{4}}}\mathit{M}\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}{|}_{\mathit{\nu }\mathrm{=}\mathrm{0}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{4}}\mathrm{+}\mathrm{6}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{3}}\mathrm{+}\mathrm{7}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}\mathrm{6}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\\ & \mathrm{+}& \mathrm{6}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}\mathrm{3}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{4}}\mathit{,}\end{array}$where  ⟨ ... ⟩ denotes taking the expectation over the probability distribution of the detection noise. For simplicity, only the statistics of the measured data for fixed  {Λ_p} is considered at first; averaging over multiple interferograms will be considered later. As such the angle brackets do not denote taking an average over the distribution of possible values of  {Λ_p}, due to, for example, atmospheric distortions to the interferograms.

An important assumption of the following analysis is that the detection noise on different pixels is statistically independent so that for example  ⟨ i_p1i_p2 ⟩  =  ⟨ i_p1 ⟩  ⟨ i_p2 ⟩  for p1 ≠ p2. It should be noted that this assumption will not be true in all cases of practical interest. For example when temporally-sampling an interferogram with a near-IR detector, it is common not to reset the charge accumulation capacitors at each pixel between successive time samples, so that a flux measurement i_p is derived from the difference between the voltage readings on two successive samples; in this case the final reading for one sample represents the initial reading for the next sample and so the noise on the two samples is correlated. Treating such readout schemes is beyond the scope of the current paper.

4. Bias-free power spectrum estimator

We seek an estimator ${\mathit{S}}_{\mathrm{0}}^{\mathit{ij}}$ for the noise-free power spectrum S^ij that does not suffer from statistical bias introduced by the noise on the interferogram. A bias-free estimator was calculated by Tatulli & Duvert (2007) in the presence of photon noise and read noise, but we present an alternative derivation here to demonstrate our analysis method. For the purposes of evaluating the power spectrum, we drop all telescope-pair-dependent superscripts, for example writing ${\mathit{H}}_{\mathit{p}}^{\mathit{ij}}$ as H_p.

We start by evaluating the expectation value of the biased estimator |c|². Combining Eqs. (6) and (7) we get: $\mathrm{⟨}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}⟨\sum _{\mathit{p}\mathrm{1}}\sum _{\mathit{p}\mathrm{2}}{\mathit{i}}_{\mathit{p}\mathrm{1}}{\mathit{i}}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathrm{\ast }}⟩\mathit{.}$(17)We split the sum in Eq. (17) for the cases where p1 = p2(≡ p) and for p1 ≠ p2 which gives: $\mathrm{⟨}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{p}}\mathrm{⟨}{{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathrm{+}\underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{1}}\mathrm{⟩}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{2}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathrm{\ast }}\mathit{,}$(18)noting that  ⟨ i_p1i_p2 ⟩  =  ⟨ i_p1 ⟩  ⟨ i_p2 ⟩  for p1 ≠ p2. Substitution of the moments of the noise model given in Eqs. (13) and (14) leads to: $\mathrm{⟨}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{p}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathrm{+}\underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathrm{\ast }}\mathit{.}$(19)The noise-free power spectrum, S, can be expressed in terms of the noise-free interferogram, Λ_p, as: $\mathit{S}\mathrm{=}{\left|\mathit{C}\right|}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{p}\mathrm{1}}\sum _{\mathit{p}\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathrm{\ast }}\mathit{,}$(20)this can be split into terms equivalent to those in Eq. (19): $\mathit{S}\mathrm{=}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathrm{+}\underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathrm{\ast }}\mathit{,}$(21)which can be substituted into Eq. (19) to give: $\mathrm{⟨}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathit{S}\mathrm{+}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathit{,}$(22)and written in terms of the noisy interferograms we see: $\mathit{S}\mathrm{=}\mathrm{⟨}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}⟨\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}⟩\mathit{.}$(23)The bias-free power spectrum estimator, S₀, is then: ${\mathit{S}}_{\mathrm{0}}\mathrm{=}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{-}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathit{.}$(24)It needs to be remembered that S₀ is intended to be evaluated on a per-interferogram, i.e. frame-by-frame basis as atmospheric turbulence changes the shape of the interferogram between frames. It is straightforward to show that when this estimator is averaged over a number of frames that the average is itself unbiased. This is because taking the average over frames commutes with taking the expectation over the detection statistics, so that: $⟨\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{f}}}\sum _{\mathit{f}}{\mathit{S}}_{\mathrm{0}}⟩\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{f}}}\sum _{\mathit{f}}⟨{\mathit{S}}_{\mathrm{0}}⟩\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{f}}}\sum _{\mathit{f}}\mathit{S}$(25)where N_f denotes the number of frames over which the average is taken. This result depends on the assumption that the detection noise is dependent only on the value of the flux Λ_p in each pixel and is statistically independent of the particular shape of the interferogram and hence of the atmospheric distortions affecting a particular interferogram.

For the case where read noise is negligible (σ_p = 0e⁻) and H is the DFT (|H_p|² = |e^2πiα_pf|² = 1), Eq. (24) reduces to the estimator reported by Goodman & Belsher (1976a): ${\mathit{S}}_{\mathrm{1}}\mathrm{=}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{-}\mathit{N,}$(26)where N =  ∑ _pi_p, is the total number of photons in an interferogram. For high light levels (negligible photon noise and read noise), the |c|² term dominates (as this goes as N²), and we reduce our estimator further to the uncorrected estimator (correct only under zero noise conditions): ${\mathit{S}}_{\mathrm{2}}\mathrm{=}{\left|\mathit{c}\right|}^{\mathrm{2}}\mathit{.}$(27)

5. Bias-free power spectrum variance

The derivation for the power spectrum variance requires the first four moments of the noise model, Eqs. (13) to (16), and stretches to multiple pages due to the |c^ij|⁴ term. The derivation follows the methodology presented above and in the derivation of the bias-free bispectrum estimator given in Appendix A but is not presented in this paper. For readability, the formulation of the power spectrum variance, $\mathrm{var}{}^{\mathrm{(}}\mathit{S}_{\mathrm{0}}{}^{\mathrm{)}}$, is split into terms that dominate in the case where the noise is dominated by photon noise, read noise or where both noise sources are significant: $\begin{array}{ccc}\mathrm{var}\mathrm{(}{\mathit{S}}_{\mathrm{0}}\mathrm{)}\mathrm{=}⟨{\left({\left|\mathit{c}\right|}^{\mathrm{2}}\mathrm{-}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\right)}^{\mathrm{2}}⟩\mathrm{-}{\mathrm{(}{\mathit{S}}^{\mathrm{)}}}^{\mathrm{2}}& & \end{array}$Photon noise terms: $\begin{array}{ccc}& \mathrm{=}& \mathrm{2}{\left|\mathit{C}\right|}^{\mathrm{2}}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathrm{+}\mathit{CC}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\mathrm{+}{\mathit{C}}^{\mathrm{\ast }}{\mathit{C}}^{\mathrm{\ast }}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\\ & \mathrm{+}& {\left[\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\right]}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\end{array}$Read noise terms: $\begin{array}{ccc}\mathrm{+}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{4}}\mathrm{+}{\left[\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\right]}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}& & \end{array}$Coupled terms: $\begin{array}{ccc}& \mathrm{+}& \mathrm{2}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\\ & \mathrm{+}& \sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\mathrm{+}\mathrm{2}{\left|\mathit{C}\right|}^{\mathrm{2}}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}\\ & \mathrm{+}& \mathit{CC}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\mathrm{+}{\mathit{C}}^{\mathrm{\ast }}{\mathit{C}}^{\mathrm{\ast }}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}\mathrm{-}\mathrm{2}\mathit{C}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathrm{\ast }}\\ & \mathrm{-}& \mathrm{2}{\mathit{C}}^{\mathrm{\ast }}\sum _{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}{\left|{\mathit{H}}_{\mathit{p}}\right|}^{\mathrm{2}}{\mathit{H}}_{\mathit{p}}\mathit{.}\end{array}$

For comparison with previous work we rewrite Eq. (28) assuming DFT conditions, that read noise is constant across all pixels (${\sum }_{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathrm{=}{\mathit{N}}_{\mathrm{pix}}{\mathit{\sigma }}^{\mathrm{2}}$), and that atmospheric phase fluctuations are random and as such the cross-product of independent complex quantities average to zero: $\begin{array}{ccc}\mathrm{var}\mathrm{(}{\mathit{S}}_{\mathrm{DFT}}\mathrm{)}& \mathrm{=}& \begin{array}{}\mathrm{Photon}\mathrm{noise}\mathrm{terms}\\ \stackrel{\mathrm{􏽺}0\mathrm{􏽽􏽼}0\mathrm{􏽻}}{\mathrm{2}\mathit{N}{\left|\mathit{C}\right|}^{\mathrm{2}}\mathrm{+}{\mathit{N}}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathrm{4}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathrm{4}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}}\end{array}\\ & \mathrm{+}& \stackrel{\mathrm{Read}\mathrm{noise}\mathrm{terms}}{\stackrel{\mathrm{􏽺}0\mathrm{􏽽􏽼}0\mathrm{􏽻}}{{\mathit{N}}_{\mathrm{pix}}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{pix}}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{4}}}}\\ & \mathrm{+}& \begin{array}{}\mathrm{Coupled}\mathrm{terms}\\ \stackrel{\mathrm{􏽺}0\mathrm{􏽽􏽼}0\mathrm{􏽻}}{\mathrm{2}{\left|\mathit{C}\right|}^{\mathrm{2}}{\mathit{N}}_{\mathrm{pix}}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{N}}_{\mathrm{pix}}{\mathit{\sigma }}^{\mathrm{2}}}\end{array}\mathit{.}\end{array}$(28)Separate analysis of the photon and detector noise terms, as is often done (see e.g. Perrin & Ridgway 2005; Le Bouquin 2005), means that the Gaussian-Poisson coupled terms are missed. Inspection of Eq. (28) and Fig. 1 show the importance of including these coupled terms for a wide range of photon fluxes, and therefore of incorporating both the photon noise and read noise into a single noise model. In the above example, for  ~300 source photons on the detector, the variance on the powerspectrum is underestimated by a factor of over 6 if cross-terms are excluded from the variance calculation.

Fig. 1 The relative importance of the photon, read (σp = 3e−) and coupled noise terms under DFT conditions for the power spectrum variance with respect to the number of source photons detected. The area highlighted in grey is the region where the coupled terms dominate the power spectrum variance.

6. Bias-free bispectrum estimator

The bias-free bispectrum estimator is calculated in the same manner as the bias-free power spectrum estimator. The definition of the bispectrum presented in Eq. (8) expands to: ${\mathit{B}}^{\mathit{ijk}}\mathrm{=}{\mathit{C}}^{\mathit{ij}}{\mathit{C}}^{\mathit{jk}}{\mathit{C}}^{\mathit{ki}}\mathrm{=}\sum _{\mathit{p}\mathrm{1}}\sum _{\mathit{p}\mathrm{2}}\sum _{\mathit{p}\mathrm{3}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}$(29)Splitting the sum in Eq. (29) gives 5 terms rather than 2 as for the power spectrum estimator. As such the derivation is quite lengthy and is given in Appendix A. We present here the bias-free bispectrum estimator, ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$: $\begin{array}{ccc}{\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}& \mathrm{=}& {\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{ij}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{jk}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{ki}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}\\ & \mathrm{+}& \sum _{\mathit{p}}\mathrm{(}\mathrm{2}{\mathit{i}}_{\mathit{p}}\mathrm{+}\mathrm{3}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathit{.}\end{array}$(30)This equation has been verified though computational simulations to the 0.1% level. In Sects. 7 and 8 we present the results of these simulations regarding the closure phase, derived from the bias-free bispectrum estimator in Eq. (30).

If we assume DFT conditions so that ${\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathrm{=}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}$ and an all-in-one interference pattern such that f^ij + f^jk + f^ki = 0, then we have ${\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathrm{=}{\mathit{H}}_{\mathit{p}}^{\mathit{ji}}$ where we note ${\mathit{H}}_{\mathit{p}}^{\mathit{ji}}\mathrm{=}{{}^{\mathrm{\right(}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}^{\mathrm{\left)}}}^{\mathrm{\ast }}$ and ${\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathrm{=}\mathrm{1}$.

If we also assume the read noise σ_p = σ is constant across all pixels rather than varying on a pixel by pixel basis then: $\sum _{\mathit{p}}{\mathit{\sigma }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ji}}}\mathrm{=}\mathrm{0.}$(31)The fact that the third order moment of the Gaussian distribution is zero would suggest that there should be no read noise terms under DFT conditions for constant σ, however it must be remembered that there are second order moments in the bias correction terms which will contribute to the bias under these conditions.

If we further assume that read noise is negligible (σ_p = 0), then ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ reduces to the estimator presented by Wirnitzer (1985) accounting for the bias introduced by photon noise alone: $\begin{array}{ccc}{\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}& \mathrm{=}& {\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\\ & \mathrm{-}& {\left|{\mathit{c}}^{\mathit{ij}}\right|}^{\mathrm{2}}\mathrm{-}{\left|{\mathit{c}}^{\mathit{jk}}\right|}^{\mathrm{2}}\mathrm{-}{\left|{\mathit{c}}^{\mathit{ki}}\right|}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{N,}\end{array}$(32)where N =  ∑ _pi_p is the mean number of photons in the interferogram.

For high light levels (negligible photon noise and read noise), the c^ijc^jkc^ki term dominates (as this goes as N³), and we reduce our estimator further to the uncorrected estimator (correct only under zero noise conditions): ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}\mathrm{=}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathit{.}$(33)The bispectrum variance is not presented in this paper as the |c^ijc^jkc^ki|² term leads to a long and tedious derivation. The zealous reader is invited to calculate the expression for the bispectrum variance if required using the method presented for calculating the bispectrum in Appendix A.

Fig. 2 The transform matrix, W under non-DFT conditions. The DC term shows the beam profile adopted, and the cosine (solid lines) and sine (dotted lines) components of the fringe pattern for each baseline show the deviation from integer number of fringe periods across the detector (with frequencies in the ratio 1.1 ×  {1:2:3}). This is similar (but with the addition of a DC term) to the visibility to pixel matrix introduced by Le Bouquin & Tatulli (2006).

The remainder of this paper will be concerned with demonstrating the magnitude of the bias introduced when the wrong estimator is used in the presence of photon and read noise at low light levels.

7. Closure phase bias

To illustrate the magnitude of the systematic bias in a realistic observation we present the closure phase data from a set of simple simulations. We created ideal interferograms resulting from a three telescope interferometer observing a target with a known and controllable closure phase. Photon noise and read noise were added, and the bispectrum (and hence closure phase) extracted using each of the three estimators presented above. The interferograms were created under either DFT or non-DFT conditions. The fringe frequency on the detector for each baseline was set at a ratio  {1:2:3} for the DFT case, and 1.1 ×  {1:2:3} for the non-DFT conditions. To simulate non-DFT conditions we rewrite Eq. (1) as:

${\mathit{i}}_{\mathit{p}}\mathrm{=}{\mathit{b}}_{\mathit{p}}{\mathit{C}}^{\mathrm{00}}\mathrm{+}{\mathit{b}}_{\mathit{p}}\sum _{\mathit{i}\mathit{<}\mathit{j}}^{{\mathit{N}}_{\mathrm{tel}}}\mathrm{Re}\left[{\mathit{C}}^{\mathit{ij}}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathit{\pi }{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}\right]\mathrm{+}\mathrm{noise}\mathit{,}$(34)where the beam profile, b_p is either flat: b_p = 1/N_pix for DFT conditions (as in Eq. (1)), or a cropped Gaussian of the form: ${\mathit{b}}_{\mathit{p}}\mathrm{=}\frac{\exp \left[\mathrm{-}\frac{{\left(\mathit{p}\mathrm{-}\mathit{\mu }\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\varsigma }}^{\mathrm{2}}}\right]}{{\sum }_{{\mathit{p}}^{\mathrm{\prime }}}^{{\mathit{N}}_{\mathrm{pix}}}\exp \left[\mathrm{-}\frac{{\left({\mathit{p}}^{\mathrm{\prime }}\mathrm{-}\mathit{\mu }\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\varsigma }}^{\mathrm{2}}}\right]}\mathit{,}$(35)where ς = 22 pixels and μ is centred on the 64 pixel detector. From Eq. (34) we can construct a V2PM, W, which is correct under DFT and non-DFT conditions. The form of W is shown in Fig. 2 for the non-DFT conditions described. We use the Singular Value Decomposition method to find the optimum (in terms of least squares) pseudo-inverse of W, which gives us H. An instrumental visibility loss was introduced such that the maximum visibility on a baseline is |V_i| = 1/3. The results presented in Fig. 3 are calculated using the ideal interferograms and therefore show only the closure phase bias and suffer no variance contribution. These are given by: $\begin{array}{ccc}{\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}& \mathrm{=}& {\mathit{B}}^{\mathit{ijk}}\mathrm{+}\sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ij}}}\sum _{\mathit{p}}\mathrm{(}\mathrm{\Lambda }\mathrm{+}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ji}}}\\ & \mathrm{+}& \sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{jk}}}\sum _{\mathit{p}}\mathrm{(}\mathrm{\Lambda }\mathrm{+}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{kj}}}\\ & \mathrm{+}& \sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ki}}}\sum _{\mathit{p}}\mathrm{(}\mathrm{\Lambda }\mathrm{+}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}{\mathit{\alpha }}_{\mathit{p}}{\mathit{f}}^{\mathit{ik}}}\mathrm{+}\mathit{N,}\end{array}$(36)and: ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}\mathrm{=}{\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}\mathrm{-}{\mathit{S}}^{\mathit{ij}}\mathrm{-}{\mathit{S}}^{\mathit{jk}}\mathrm{-}{\mathit{S}}^{\mathit{ki}}\mathrm{-}\mathit{N}\mathrm{-}\mathrm{3}{\mathit{N}}_{\mathrm{pix}}{\mathit{\sigma }}^{\mathrm{2}}\mathit{,}$(37)whilst the bias free estimator is simply: ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}\mathrm{=}{\mathit{B}}^{\mathit{ijk}}\mathit{.}$(38)We note the order of the indices in the DFT terms in Eqs. (36) and (37) and that the ideal source bispectrum and power spectrum are independent of the measurement system (weather or not if meets the DFT conditions), whereas the bias terms implicitly assume DFT conditions for the ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ estimators. Figure 3 shows a simulation of the bias present on the closure phase for the three estimators for simulations of 120 and 300 photons per interferogram reaching the detector.

Fig. 3 The bias in closure phase for ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ (solid line), ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ (dashed line), and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ (dotted line) estimators under a number of measurement conditions. Plots a) and b) show DFT conditions, whilst plots c) and d) show non-DFT conditions (described in the text). Read noise of σp = 3e− is present in plots b) and d). The thin and think lines show 120 and 300 photons reaching the detector respectively. The shaded regions in plots c) and d) show the range of possible biases resulting from different possible combinations of phases on the three baselines contributing to the source closure phase.

7.1. Zero read noise

In the situation where photon noise dominates the noise model (σ_p = 0e⁻) and the instrument is setup such that the interferograms are formed under DFT conditions, ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ both return the correct closure phase, however ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ is affected by a bias term resulting from photon noise. At these light levels, the bias in ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ becomes significant for the majority of source closure phases.

7.2. Non-zero read noise

At low light levels the read noise on the detector becomes important and comparable to the photon noise. Figure 3b shows simulations under DFT conditions, as in Fig. 3a, but with σ_p = 3e⁻ read noise included. This level of read noise is quite low compared to many of the optical detectors in use today. We see that failing to account for read noise introduces a significant bias in closure phase estimates for both ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$. It is clear that the magnitude of the bias term is strongly dependent on the source brightness, and that ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ is preferable to ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ under certain combinations of noise and photon flux.

7.3. Non-DFT conditions

Deviation from DFT conditions changes the way statistical bias affects the bispectrum. Inspection of Eq. (30) shows that the statistical bias can affect both the real and imaginary parts of the bispectrum when H is not the DFT. If the measurement system conforms to DFT conditions, and the read noise is negligible or constant across all pixels, the bias only affects the real part of the bispectrum estimators as shown when applying Eq. (31) to (30). It is also apparent that the magnitude of the bias is no longer dependent only on the value of the closure phase, but on the value of the individual phase of each baseline. The non-uniform beam profile causes an overlap of the fringe components in frequency space, and the non-integral fringe frequency on the detector invalidates the assumption that ${\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathrm{=}{\mathit{H}}_{\mathit{p}}^{\mathit{ji}}$ and ${\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathrm{=}\mathrm{1}$ which led to the formulation of ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$. Figure 3c is provided as an example of the impact of non-DFT conditions on the closure phase bias for a photon noise dominated instrument (σ_p = 0). Note that the non-DFT conditions change the bias terms, meaning ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ is now biased when it was not under DFT conditions. Figure 3d shows the effect of non-DFT conditions on a system with photon and read noise and again there is an alteration of the bias terms which ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ corrects but the other estimators do not. The grey regions on the non-DFT plots show the range of biases associated with each estimator depending on the individual phases of the triplet producing a given closure phase. The lines within the shaded regions represent the expectation of the closure phase bias over all possible individual phases.

We see that the bias on the closure phase does not cancel when accumulating frames; specifically, that there is no advantage to observing in phase-unstable conditions over phase-stable conditions as a means of reducing bias. Therefore, bias-wise, one should always favor fringe-tracker assisted observations whereby the longer discrete integration times result in higher SNR, even in the photon-rich regime where fringe tracking is not strictly necessary.

7.4. Closure phase bias and statistical uncertainty

The bias is only significant at low light levels, and tends to zero as the light level increases. At low light levels the variance is also high, and we must determine the relative importance of these two processes. We present Fig. 4 which shows the relative importance of bias, Φ^bias, and closure phase standard deviation, ς_Φ, using different estimators.

We show the number of frames, N_f, such that: ${\mathrm{\Phi }}^{\mathrm{bias}}\mathrm{=}\frac{{\mathit{\varsigma }}_{\mathrm{\Phi }}}{\sqrt{{\mathit{N}}_{\mathit{f}}}}\mathit{,}$(39)which is the number of interferograms required to reduce the magnitude of the statistical uncertainty to that of the bias. The closure phase standard deviation per frame, ς_Φ, is given by: ${\mathit{\varsigma }}_{\mathrm{\Phi }}\mathrm{=}\frac{\mathrm{1}}{\left|{\mathit{B}}^{\mathit{ijk}}\right|}\mathrm{\times }\sqrt{\mathrm{var}\left[\mathrm{Im}\left({\mathit{B}}_{\mathit{f}}^{\mathit{ijk}}\mathrm{\times }\frac{{\mathit{B}}^{\mathit{ijk}\mathrm{\ast }}}{\left|{\mathit{B}}^{\mathit{ijk}}\right|}\right)\right]}\mathit{,}$(40)where ${\mathit{B}}_{\mathit{f}}^{\mathit{ijk}}$ is the bispectrum calculated per frame using the chosen estimator and $\mathrm{Im}\left({\mathit{B}}_{\mathit{f}}^{\mathit{ijk}}\mathrm{\times }\frac{{\mathit{B}}^{\mathit{ijk}\mathrm{\ast }}}{\left|{\mathit{B}}^{\mathit{ijk}}\right|}\right)$ gives the projection of ${\mathit{B}}_{\mathit{f}}^{\mathit{ijk}}$ on a unit vector perpendicular to B^ijk, and. This method is used because it is insensitive to wrapping errors in the complex plane when the variance is large. In this example, the magnitude of the bias and variance are measured when the source closure phase is at 90°. The values presented in Fig. 4 were derived using numerical simulations to the 0.1% level, and the bias contribution was verified analytically using Eqs. (36) and (37).

It can be seen in Fig. 4 that the bias is important for realistic numbers of frames (e.g. 10³ frames or 25 s integration for 25 ms exposures) compared to the uncertainty for a wide range of typical light levels.

Fig. 4 The number of integrated frames required to reduce the statistical uncertainty to the level of the bias for the ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ (dashed) and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ (dotted) estimators.

8. Simulation: detection of faint companions

Closure phase nulling is an observational techniques that allows the detection of faint companions (Chelli et al. 2009). The technique has been demonstrated by Duvert et al. (2010), detecting a 5 mag fainter companion using AMBER/VLTI around HD 59717.

Closure phase nulling uses the fact that the closure phase jumps from 0 to 180 degrees as the visibility function of the primary in the system passes through a null. The presence of a companion becomes evident as a deviation from an instantaneous phase jump. The profile of the phase jump gives information about the projected separation and flux ratio of the system components. We show that statistical bias in the closure phase leads directly to incorrect estimation of the system properties.

We follow Chelli et al. (2009) and simulate a binary system where the primary is described by a uniform disc of radius, R_⋆, with a point source companion at a separation, s, with flux ratio, r. For simplicity we simulate a three telescope interferometer such that the position angle of the binary system is parallel to the frequency axis of the interferometer. The interferometer samples three frequencies  {u₁₂,u₂₃,u₁₃}  such that u₁₃ = u₁₂ + u₂₃. The instrumental parameters are given in Table 1 and we assume the instrument is functioning under DFT conditions.

The complex visibility, C^ij(u), at a given frequency, u, is described by: ${\mathit{C}}^{\mathit{ij}}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\frac{{\mathit{V}}_{\mathit{\star }}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{+}\mathit{r}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathit{\pi us}}}{\mathrm{1}\mathrm{+}\mathit{r}}\mathit{,}$(41)which is the visibility function, V_⋆(u), of the primary: ${\mathit{V}}_{\mathit{\star }}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathrm{2}\frac{{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathrm{2}\mathit{\pi u}{\mathit{R}}_{\mathit{\star }}\mathrm{\right)}}{\mathrm{2}\mathit{\pi u}{\mathit{R}}_{\mathit{\star }}}\mathit{,}$(42)plus the contribution of the secondary with a phase modulation related to the projected separation, s, and flux ratio, r. J₁ is the first order Bessel function and R_⋆ is the radius of the primary.

Fig. 5 Demonstration of the effect of bias on a closure phase measurement. Plot a): the closure phase trace resulting from a binary pair consisting of a resolved primary of radius R⋆ and an unresolved secondary with flux ratio r = 0.01 for separations (s = 10 R⋆,100 R⋆,1000 R⋆). Plot b): the closure phase trace for the 100 R⋆ system recovered by the three estimators. ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ correctly returns the ideal closure phase, whilst ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ do not. Both plots show the single source closure phase signal (thin solid line) for reference. Plots c) and d) show the closure phase calculated using ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ respectively as in plot b) along with the statistical uncertainty on the measurement after 103, 104 and 105 frames.

The ideal bispectrum is extracted from Eq. (41) using Eq. (8) and the closure phase is taken as the argument of the averaged bispectrum. The closure phase signal around the first null in the visibility profile is shown in Fig. 2 of Chelli et al. (2009) and reproduced using our simulation model in Fig. 5a in this paper for a number of companion separations (s = 10R_⋆,100 R_⋆,1000 R_⋆). The closure phase is wrapped modulo 2π and the frequencies  {u₁₂,u₂₃,u₁₃}  are set to the ratio  {1:2:3} on the detector and scaled by the radius of the primary, R_⋆. For simplicity, the baseline to fringe-frequency mapping is chosen such that it preserves the ratios of the two. The specifics of the mapping is not expected to have a significant impact on the bias. The closure phase of the primary alone (single source) is shown for comparison.

We present in Fig. 5b a simulation of the closure phase from a 100 R_⋆ separation binary with flux ratio 0.01. ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ returns the correct closure phase signal in agreement with the ideal case presented by Chelli et al. (2009). ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ are strongly affected by bias terms, meaning the closure phase is significantly different. Depending which regions of the closure phase trace are sampled by the interferometer, ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ would return an incorrect value for the primary radius and companion flux ratio. Figures 5c and d show the closure phase calculated using ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ respectively and the statistical uncertainty on the measurement after 10³, 10⁴ and 10⁵ frames calculated using Eq. (40). We see that even for very short total integration times (10³ frames is under a minute using typical discrete integration times) the bias still produces a significant contribution.

9. Conclusions

Interferometry offers exceptional angular resolution at the expense of sensitivity, as the collecting area is, by necessity, smaller than that of a single dish with the same resolution. In ground-based optical interferometry these sensitivity limitations are accentuated by the need for short observations to reduce the phase perturbations introduced by the atmosphere. The result is that noise levels on individual observation frames are often high, exacerbated by the need for fast read-out of the detector. Without understanding the statistical properties of the noise, bias terms are introduced when extracting information from the interferogram using the non-linear power spectrum and bispectrum constructs. We have shown in this paper that:

• 1.

  Neglecting the cross terms that arise in the presence of both readnoise and photon noise leads to bias on the extraction ofinterferometric observables resulting from the statistics of thenoise.

• 2.

  The statistical bias on the observables results in incorrect phisical parameters being extracted from observation data, for example the flux ratio of binary pairs or the diameter of a target as shown for the case of closure phase nulling.

To remedy this we have presented bias-free estimators for both the power spectrum and bispectrum (and hence the closure phase) which can be used in the presence of any combination of photon and read noise under any fringe encoding scheme in which the measurement equation is linear. Our bias-free estimators are presented in a format allowing easy implementation into current data reduction pipelines.

Acknowledgments

We thank the anonymous referee for constructive comments and insiteful suggestions that have undoubtedly improved this paper.

References

Appendix A: Derivation of the bias-free bispectrum estimator

The ideal, noise-free bispectrum is given by: ${\mathit{B}}^{\mathit{ijk}}\mathrm{=}{\mathit{C}}^{\mathit{ij}}{\mathit{C}}^{\mathit{jk}}{\mathit{C}}^{\mathit{ki}}\mathrm{=}\sum _{\mathit{p}\mathrm{1}}\sum _{\mathit{p}\mathrm{2}}\sum _{\mathit{p}\mathrm{3}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{,}$(A.1)which can be split into: $\begin{array}{ccc}{\mathit{B}}^{\mathit{ijk}}& \mathrm{=}& \sum _{\mathit{p}}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{3}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{1}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}\end{array}$(A.2)The equivalent expression for real data embedded with noise when average over frames is: $\mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{p}\mathrm{1}}\sum _{\mathit{p}\mathrm{2}}\sum _{\mathit{p}\mathrm{3}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{1}}{\mathit{i}}_{\mathit{p}\mathrm{2}}{\mathit{i}}_{\mathit{p}\mathrm{3}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{,}$(A.3)which splits into: $\begin{array}{ccc}\mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}& \mathrm{=}& \sum _{\mathit{p}}\mathrm{⟨}{{\mathit{i}}_{\mathit{p}}^{\mathrm{3}}}^{\mathrm{⟩}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{1}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{1}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{2}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{3}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{1}}\mathrm{⟩}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{2}}\mathrm{⟩}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{3}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}\end{array}$(A.4)Substitution of the noise model in Eqs. (13) to (15) gives the noisy bispectrum in terms of the ideal interferogram and associated noise: $\begin{array}{ccc}\mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}& \mathrm{=}& \sum _{\mathit{p}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{3}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{1}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\mathrm{1}\ne \mathit{p}\mathrm{2}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}\end{array}$(A.5)Substituting the ideal bispectrum split sum from Eq. (A.2) gives: $\begin{array}{ccc}\mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}& \mathrm{=}& {\mathit{B}}^{\mathit{ijk}}\mathrm{+}\sum _{\mathit{p}}\mathrm{(}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathit{p}}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{1}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{+}& \underset{\mathit{p}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{(}{\mathrm{\Lambda }}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Lambda }}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}\end{array}$(A.6)By rearranging and substituting in for the ensemble average over the real interferogram we find the ideal bispectrum in terms of the noisy interferogram: $\begin{array}{ccc}{\mathit{B}}^{\mathit{ijk}}& \mathrm{=}& \mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}\\ & \mathrm{-}& ⟨\sum _{\mathit{p}}\mathrm{(}\mathrm{3}{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{i}}_{\mathit{p}}\mathrm{+}\mathrm{3}{\mathit{i}}_{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}⟩\\ & \mathrm{-}& \underset{\mathit{p}\ne \mathit{p}\mathrm{1}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{1}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& \underset{\mathit{p}\ne \mathit{p}\mathrm{2}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{2}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& \underset{\mathit{p}\ne \mathit{p}\mathrm{3}}{\begin{array}{}\sum \end{array}\begin{array}{}\sum \end{array}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{⟩}}\mathrm{⟨}{\mathit{i}}_{\mathit{p}\mathrm{3}}\mathrm{⟩}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}\mathit{.}\end{array}$(A.7)Expanding out the double summations using the fact that  ∑  ∑ a ≠ b =  ∑  ∑ a b −  ∑ a = b gives:

$\begin{array}{ccc}{\mathit{B}}^{\mathit{ijk}}& \mathrm{=}& \mathrm{⟨}{\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\mathrm{⟩}\\ & \mathrm{-}& ⟨\sum _{\mathit{p}}\mathrm{(}\mathrm{3}{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{i}}_{\mathit{p}}\mathrm{+}\mathrm{3}{\mathit{i}}_{\mathit{p}}{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}⟩\\ & \mathrm{+}& \mathrm{3}⟨\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}_{\mathit{p}}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}⟩\\ & \mathrm{-}& ⟨\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\sum _{\mathit{p}\mathrm{1}}{\mathit{i}}_{\mathit{p}\mathrm{1}}{\mathit{H}}_{\mathit{p}\mathrm{1}}^{\mathit{ij}}⟩\\ & \mathrm{-}& ⟨\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\sum _{\mathit{p}\mathrm{2}}{\mathit{i}}_{\mathit{p}\mathrm{2}}{\mathit{H}}_{\mathit{p}\mathrm{2}}^{\mathit{jk}}⟩\\ & \mathrm{-}& ⟨\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}\sum _{\mathit{p}\mathrm{3}}{\mathit{i}}_{\mathit{p}\mathrm{3}}{\mathit{H}}_{\mathit{p}\mathrm{3}}^{\mathit{ki}}⟩\mathit{.}\end{array}$(A.8)Cancelling terms and noting that, ${\sum }_{\mathit{p}}{\mathit{i}}_{\mathit{p}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}\mathrm{=}{\mathit{c}}^{\mathit{ij}}$, we find the bias free bispectrum estimator:

$\begin{array}{ccc}{\mathit{B}}_{\mathrm{bias}\mathrm{-}\mathrm{free}}^{\mathit{ijk}}& \mathrm{=}& {\mathit{c}}^{\mathit{ij}}{\mathit{c}}^{\mathit{jk}}{\mathit{c}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{ij}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{jk}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\\ & \mathrm{-}& {\mathit{c}}^{\mathit{ki}}\sum _{\mathit{p}}\mathrm{(}{\mathit{i}}_{\mathit{p}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}\\ & \mathrm{+}& \sum _{\mathit{p}}\mathrm{(}\mathrm{2}{\mathit{i}}_{\mathit{p}}\mathrm{+}\mathrm{3}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{H}}_{\mathit{p}}^{\mathit{ij}}{\mathit{H}}_{\mathit{p}}^{\mathit{jk}}{\mathit{H}}_{\mathit{p}}^{\mathit{ki}}\mathit{.}\end{array}$(A.9)

All Tables

All Figures

	Fig. 1 The relative importance of the photon, read (σp = 3e−) and coupled noise terms under DFT conditions for the power spectrum variance with respect to the number of source photons detected. The area highlighted in grey is the region where the coupled terms dominate the power spectrum variance.
In the text

Fig. 2 The transform matrix, W under non-DFT conditions. The DC term shows the beam profile adopted, and the cosine (solid lines) and sine (dotted lines) components of the fringe pattern for each baseline show the deviation from integer number of fringe periods across the detector (with frequencies in the ratio 1.1 ×  {1:2:3}). This is similar (but with the addition of a DC term) to the visibility to pixel matrix introduced by Le Bouquin & Tatulli (2006).

In the text

Fig. 3 The bias in closure phase for ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ (solid line), ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ (dashed line), and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ (dotted line) estimators under a number of measurement conditions. Plots a) and b) show DFT conditions, whilst plots c) and d) show non-DFT conditions (described in the text). Read noise of σp = 3e− is present in plots b) and d). The thin and think lines show 120 and 300 photons reaching the detector respectively. The shaded regions in plots c) and d) show the range of possible biases resulting from different possible combinations of phases on the three baselines contributing to the source closure phase.

In the text

	Fig. 4 The number of integrated frames required to reduce the statistical uncertainty to the level of the bias for the ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ (dashed) and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ (dotted) estimators.
In the text

Fig. 5 Demonstration of the effect of bias on a closure phase measurement. Plot a): the closure phase trace resulting from a binary pair consisting of a resolved primary of radius R⋆ and an unresolved secondary with flux ratio r = 0.01 for separations (s = 10 R⋆,100 R⋆,1000 R⋆). Plot b): the closure phase trace for the 100 R⋆ system recovered by the three estimators. ${\mathit{B}}_{\mathrm{0}}^{\mathit{ijk}}$ correctly returns the ideal closure phase, whilst ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ do not. Both plots show the single source closure phase signal (thin solid line) for reference. Plots c) and d) show the closure phase calculated using ${\mathit{B}}_{\mathrm{1}}^{\mathit{ijk}}$ and ${\mathit{B}}_{\mathrm{2}}^{\mathit{ijk}}$ respectively as in plot b) along with the statistical uncertainty on the measurement after 103, 104 and 105 frames.

In the text