© ESO, 2015

1. Introduction

Long-baseline interferometry (LBI) enables us to spatially resolve components in close orbits, providing astrometric positions at <50 milli-arcsec (mas) scale with micro-arcsecond accuracy. When combined with spectroscopic radial velocities, we can obtain model independent estimates of the stellar masses and orbital parallaxes, which are fundamental parameters that help us study stellar properties and evolution. However, so far LBI is limited to bright stars (H < 7 mag) with typical magnitude differences ΔH< 6 mag. Adaptive optics imaging with single-dish telescopes reaches better contrasts, down to ΔH ~ 12 mag (Zurlo et al. 2014), but the angular resolution is limited to 0.2″ at these detection levels. Long-baseline interferometry is therefore a complementary technique in terms of spatial scale by probing the innermost regions.

LBI can detect components down to ~1 mas in the infrared, but the main limitation is sensitivity to high-contrast companions. We roughly know the performances reachable by current interferometric combiners. So far, the faintest companion detected with LBI has a flux ratio of 0.75% and orbits a Cepheid star of magnitude H = 3.85 mag (Gallenne et al. 2014). Objects, such as brown dwarfs and hot giant planets, are still inaccessible because of a lack of sensitivity and accuracy of the instruments. Absil et al. (2011) demonstrated a possible dynamic range of 1:500 with the VLTI/PIONIER instrument (Le Bouquin et al. 2011), but this range has not yet been reached. To achieve this detection level, several hours of multitelescope observations are required to obtain as many simultaneous interferometric measurements as possible.

Deriving the sensitivity limit from imaging can be determined directly from the noise level, however, this is not the case for interferometric observations. There are some papers in the literature that discuss detection limits and methods to search for companions from interferometric data (Absil et al. 2011; Le Bouquin & Absil 2012), but these studies have some shortcomings: the searching method is not formalized, the sigma detection is not robust, and they do not take the bandwidth smearing into account. Therefore, a robust implementation to search for components does not exist thus far. This kind of method is particularly critical to detect faint companions, as they can be at the sensitivity limit of the instrument or even not be a statistically significant detection (i.e., <3σ). We therefore created CANDID (Companion Analysis and Non-Detection in Interferometric Data) to address these aspects. This is a suite of Python tools, which contains two main functions: 1) one to perform a systematic search for faint companions (Sect. 2); and 2) one to estimate the detection limit from long-baseline interferometric observations (Sect. 3). This tool is made available to the community¹.

CANDID is made for detecting high-contrast, point-source companions orbiting a spatially resolved primary star, although it also works with an unresolved primary and contrast <50%. In this paper, we used CANDID to look for companions in binary Cepheids. We present the first main function to search for companions in Sect. 2, verify and clarify the detection level of our previously detected faint companions, and we report new detections for other Cepheids. In Sect. 3, we present the second main function of CANDID and explain our robust method to set detection limits from interferometric data. In Sect. 4, we then use this formalism to our set of Cepheid observations to derive the detection limits. We finally present our conclusions in Sect. 5.

2. Searching for companions

Our grid search is traditionally performed in a three-dimensional space. Specifically, we vary the position of the component (Δα,Δδ) and the companion/star flux ratio f, and then compute the χ² for each of these positions and flux ratios (the third dimension being f). The weakness of this method is the resolution of the grid, i.e., if the grid is to too coarse, the detection can be missed.

In addition, for high-contrast binaries and/or low accuracy data, a χ² map can show fake or nonsignificant detections. Sometimes, some authors show the estimated χ² map with the most probable location of the component without any detection level mentioned, only the reduced χ² variation is given. This parameter, however, is not optimal for checking the detection level because it depends on the number of degrees of freedom (d.o.f.). The large quantity of data required to detect a faint companion using LBI results in a large dof. For instance, we can have the most probable location with the lowest reduced chi-square, ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$, equal to 1.0, and the highest in the whole map equal to 1.1, i.e., with only a variation of 0.1. With this information, however, there is no way to know if the detection is statistically significant or not. In this section, we therefore present a systematic approach to search for components and set the detection level using the χ² and the number of dof. We then applied it to the case of binary Cepheids. We verify our previous detections for the Cepheids V1334 Cyg and AX Cir first (Gallenne et al. 2013b, 2014), and then report new detections for RT Aur and AW Per.

2.1. Detection method

A more rigorous approach is to perform a grid of fit using a least-squares algorithm, with a starting grid spacing that is guaranteed to find the global minimum. The grid in question is the 2D grid of starting points for the companion position. For each starting position, a multiparameter fit is performed: the companion position and its flux ratio (possibly the stellar diameters) are adjusted. Each position of the grid leads to a local minimum. Ideally, if the starting grid is fine enough, multiple starting points lead to the same local minima: this guarantees that all the local minima are explored and that the global minimum is indeed the global minimum. Hence, the criteria to decide if the global minimum is global is a posteriori. In CANDID, we require (a posteriori) that on average, each unique minimum is reached from two starting points of the starting grid. We also provide statistics on the “traveling” distance of the fit (between the starting and end points), compared to the size of the grid. We require that the median travel distance should be less than $\sqrt{\mathrm{2}}\mathit{/}\mathrm{2}$ of the size of the square grid. The 10 and 90% percentiles are also provided to the user to assess the typical travel distance.

A systematic search using a grid of fit is an iterative process. First, a coarse starting grid is chosen and the fits are performed. CANDID estimates, based on the traveling distances of the fits, how to refine the starting grid. A second series of fits are run, using the finer grid, which might take much longer, but the global minimum can be trusted to be global (at least within the area searched).

We searched for components with a maximum distance to the main star of 50 mas. For a wider range, the loss of coherence caused by spectral smearing of the companion is the main limitation because it degrades the dynamic range. A complete discussion about the search region is presented by Absil et al. (2011) and Le Bouquin & Absil (2012). The main limitation is the spectral sampling compared to the relative position of the component, i.e., to avoid a significant smearing we need R>ρB_p/λ, with R, B_p and ρ the spectral resolution, the projected baseline, and the separation, respectively. For PIONIER, with R = 18 and a mean projected baseline B_p = 100 m (used for our observations), we have ρ< 60 mas, and for MIRC with R = 42 and B = 200 m, we have ρ< 70 mas. Zhao et al. (2007) recommended a more stringent criteria ρ<Rλ/ (5B) to assure uncorrupted data, leading to ρ ≲ 15 mas for PIONIER and MIRC. However, a companion can still be detected at separations larger than 15 mas by taking bandwidth smearing effects into account. We chose 50 mas as our limiting range to avoid making the grid search too stringent. In addition, companions located at more than 50 mas are detected more efficiently using adaptive optics on a single-dish telescope (through imaging or sparce aperture masking).

Each point in the grid is fitted with the following binary model, representing a spatially resolved primary star with a point source component: $\begin{array}{c}\mathrm{˜}\\ \mathit{V}\end{array}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\frac{{\mathit{V}}_{\mathit{\star }}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{+}\mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathit{f}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{c}}\end{array}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{f}}\mathit{,}$(1)with, $\begin{array}{ccc}& & {\mathit{V}}_{\mathit{\star }}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{x}}\mathit{,}\\ & & \begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{c}}\end{array}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\exp \mathrm{[}\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\pi }\mathrm{(}\mathit{u}\mathrm{\Delta }\mathit{\alpha }\mathrm{+}\mathit{v}\mathrm{\Delta }\mathit{\delta }\mathrm{)}\mathit{/}\mathit{\lambda }\mathrm{]}\mathit{,}\\ & & \mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathrm{=}\left|\frac{\sin \mathit{\zeta }}{\mathit{\zeta }}\right| \mathrm{with} \mathit{\zeta }\mathrm{=}\frac{\mathit{\pi }\mathrm{(}\mathit{u}\mathrm{\Delta }\mathit{\alpha }\mathrm{+}\mathit{v}\mathrm{\Delta }\mathit{\delta }\mathrm{)}}{\mathit{R\lambda }}\mathit{,}\end{array}$where J₁(x) is the first-order Bessel function, $\mathit{x}\mathrm{=}\mathit{\pi }{\mathit{\theta }}_{\mathrm{UD}}\sqrt{{\mathit{u}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}^{\mathrm{2}}}\mathit{/}\mathit{\lambda }$, (u,v) the spatial frequencies, θ_UD the uniform disk angular diameter of the primary star, λ the wavelength, f the flux ratio between the companion and the primary star, (Δα,Δδ) the relative position of the component with respect to the primary, and R = λ/ Δλ the spectral resolution. The function G(ζ) is a corrective term to overcome the effect of bandwidth smearing (Lachaume & Berger 2013). The fitted parameters are θ_UD,f,Δα, and Δδ. Although the bandwidth smearing can be fit in CANDID, we kept it fixed for this work because we noticed that it cannot be constrained by these data.

The interferometric observables, i.e., the squared visibility V², the closure phase CP, and the bispectrum amplitude B_amp, are then estimated from the squared modulus and the bispectrum in closed triangles as follows: ${\mathit{V}}^{\mathrm{2}}\mathrm{=}\mathrm{|}\begin{array}{c}\mathrm{˜}\\ \mathit{V}\end{array}{\mathrm{|}}^{\mathrm{2}} \mathrm{and} \begin{array}{c}\mathrm{˜}\\ \mathit{B}\end{array}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{12}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{23}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{31}}^{\mathrm{\ast }}\end{array}\mathit{,}$(5)which provides the bispectrum amplitude and the closure phase from the definition $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{B}\end{array}\mathrm{=}{\mathit{B}}_{\mathrm{amp}}{\mathrm{e}}^{\mathrm{-}\mathrm{iCP}}$: ${\mathit{B}}_{\mathrm{amp}}\mathrm{=}\mathrm{\left|}\begin{array}{c}\mathrm{˜}\\ \mathit{B}\end{array}\mathrm{\right|} \mathrm{and} \mathrm{CP}\mathrm{=}\arg \mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ \mathit{B}\end{array}\mathrm{\right)}\mathit{.}$(6)The data are then fitted simultaneously using a Levenberg-Marquardt least-squares minimization algorithm with $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \begin{array}{c}\sum \end{array}\mathrm{(}{\mathrm{CP}}_{\mathrm{o}}\mathrm{-}{\mathrm{CP}}_{\mathrm{m}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{\mathrm{CP}}^{\mathrm{2}}\mathrm{+}\begin{array}{c}\sum \end{array}{\mathrm{(}{\mathit{B}}_{\mathrm{amp}\mathit{,}\mathrm{o}}\mathrm{-}{\mathit{B}}_{\mathrm{amp}\mathit{,}\mathrm{m}}\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{{\mathrm{B}}_{\mathrm{amp}}}^{\mathrm{2}}\\ & & \mathrm{+}\begin{array}{c}\sum \end{array}{\mathrm{(}{\mathit{V}}_{\mathrm{o}}^{\mathrm{2}}\mathrm{-}{{\mathit{V}}_{\mathrm{m}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{{\mathit{V}}^{\mathrm{2}}}^{\mathrm{2}}\mathit{,}\end{array}$(7)where the indexes o and m denote the data and the model, respectively. CANDID can fit all observables or just one, depending on the data available and the user. We then divided by the number of degrees of freedom to obtain a map of the chi-square minima, and interpolated these minima on a regular grid.

The grid resolution is critical as it depends on the distance explored from the initial and final positions. This is why we implemented an estimate of the optimum grid resolution in our code. It is worth mentioning that a denser starting grid is not necessarily better as the fit would not be improved, and we would lose in computation time (CANDID performs a 30 × 30 grid with 1 mas steps in 50 s using six cores and a data set with 879 degrees of freedom). To make the grid search faster, CANDID was developed for parallel processing on multicore machines.

Once the map of the minima is computed, we can check the variation in the whole map, however this does not tell us about the statistical significance: specifically, whether the most probable location is detected at 1σ or more. Although it might not be important for low-contrast companions for which the variation is large enough, this is critical for components with a flux ratio <5% with a variation of the minima of a few percent. The number of dof is an important parameter in that context. Assuming that the data follow Gaussian statistics, we implemented in CANDID an estimate of the number of sigma for each point in the grid in order to obtain a nσ detection map. The formalism we used is based on the probability P (or confidence interval) with ν degrees of freedom, as already employed by Absil et al. (2011). The number of sigma demonstrates how our binary model is significant compared to a uniform disk model (i.e., a single star). We used the following formula for the probability: $\mathit{P}\mathrm{\left(}\mathrm{\Delta }\mathit{\alpha ,}\mathrm{\Delta }\mathit{\delta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathrm{CDF}}_{\mathit{\nu }}\left(\frac{\mathit{\nu }{\mathit{\chi }}_{\mathrm{UD}}^{\mathrm{2}}}{{\mathit{\chi }}_{\mathrm{r}\mathit{,}\mathrm{bin}}^{\mathrm{2}}\mathrm{\left(}\mathrm{\Delta }\mathit{\alpha ,}\mathrm{\Delta }\mathit{\delta }\mathrm{\right)}}\right)\mathit{,}$(8)where ${\mathit{\chi }}_{\mathrm{bin}}^{\mathrm{2}}$ and ${\mathit{\chi }}_{\mathrm{UD}}^{\mathrm{2}}$ are the minimum chi-square for the binary and the uniform disk model (i.e., fitting Eqs. (1)and (2), respectively), and CDF denotes the χ² cumulative probability distribution function with ν degrees of freedom. We then convert the probability into the number of sigmas, nσ (e.g., 99.73% = 3σ, 99.99% = 4σ,...). To avoid big float numbers, we limited the maximum value to 50σ.

This formalism therefore provides a ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$ map to find the most probable location of a companion, if any, and a nσ map giving the detection level at each point in the grid.

Fig. 1 ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$ map of the local minima (left) and detection level map (right) of V1334 Cyg for the observations on 2012-10-01. The yellow lines represent the convergence from the starting points to the final fitted position. The maps were reinterpolated in a regular grid for clarity. The axis limit was chosen according to the location of the companion.

Fig. 2 Same as Fig. 1, except for AX Cir for the observations on 2013-07-14.

Fig. 3 Same as Fig. 1, except for RT Aur using only the closure phase.

2.2. Previously published detections

To validate our method, we computed the maps of the χ² minima and estimated the detection level for our previously detected companions around the Cepheids V1334 Cyg and AX Cir (Gallenne et al. 2013b, 2014, observed with the instruments CHARA/MIRC and VLTI/PIONIER, respectively). The journal of these previous observations are reported in Table 1. The maps are presented in Figs. 1 and 2, for which all observables were fitted, except for PIONIER where only CP and V² were used because there is no good estimator of B_amp so far. The central part has been hidden to improve the clarity of the intensity map level, which can be biased by the primary star. The companion orbiting V1334 Cyg is detected at more than 25σ, and at more than 13σ for AX Cir. We summarized the detection levels for these two stars in Table 2, including fitting only the closure phase signal. The CP is more sensitive to faint off-axis companions and is also less affected by instrumental and atmospheric perturbations than the other observables (i.e., V² and B_amp). Fitting all of the observables can improve the detection level because we add more information, but it can also affect the results, depending on the magnitude of the biases. We notice from Table 2 that including V² and B_amp degrades the detection level for V1334 Cyg, although it is still significant, and including these observables improves the detection for AX Cir. A possible explanation that addresses why the detection level decreases sometimes when we add the V² is that we add correlated noise, which is not consistent with our hypothesis of uncorrelated noise. When only the CP is used, the angular diameter of the primary is first determined by fitting only a uniform disk model to the square visibilities, and then kept fixed during the grid search.

The resulting fitted parameters (i.e., the astrometric position, flux ratio, and the angular diameter) are in good agreement with the values determined in Gallenne et al. (2013b, 2014).

2.3. New detections of the companions around the Cepheids RT Aur and AW Per

Our interferometric program on Galactic binary Cepheids, which started two years ago, is promising to directly measure the dynamical masses. We have obtained several observing nights in both hemispheres with the multitelescope combiners CHARA/MIRC and VLTI/PIONIER to detect the close companions of a few Cepheids (Gallenne et al. 2014, 2013b,a). Here we report new detections for the Cepheids AW Per and RT Aur.

The observations were performed in 2012 using the Michigan InfraRed Combiner (MIRC) installed at the CHARA Array (ten Brummelaar et al. 2005), located on Mount Wilson, California. The CHARA Array consists of six 1 m aperture telescopes in a Y-shaped configuration (two telescopes on each branch), oriented to the east (E1, E2), west (W1,W2) and south (S1, S2), providing good coverage of the (u,v) plane. The baselines range from 34 m to 331 m, providing high angular resolution down to 0.5 mas at H. The MIRC instrument (Monnier et al. 2004, 2010) is an image-plane combiner, which enables us to combine the light coming from all six telescopes in H or K. MIRC also offers three spectral resolutions (R = 42,150 and 400), which provide 15 visibility and 20 closure phase measurements across a range of spectral channels.

Fig. 4 Same as Fig. 1, except for AW Per.

We observed the short-period Cepheids RT Aur (HD 45412, P_puls = 3.73 d) and AW Per (HD 30282, P_puls = 6.46 d) with all six telescopes. We used the H-band filter with the lowest spectral resolution in which the light is split into eight spectral channels. Table 1 lists the journal of our observations. We followed a standard observing procedure in which we monitored the interferometric transfer function by observing a calibrator before and/or after the Cepheids. The calibrators, listed in Table 1, were selected using the SearchCal² software (Bonneau et al. 2006, 2011) provided by the Jean-Marie Mariotti Center³.

We reduced the data using the standard MIRC pipeline (Monnier et al. 2007), which consists of computing the squared visibilities and triple products for each baseline and spectral channel, and correcting for photon and readout noise. Squared visibilities are estimated using Fourier transforms, while the triple products are evaluated from the amplitudes and phases between three baselines forming a closed triangle.

We used CANDID to search for a component within ±50 mas. For RT Aur, we might have detected a companion at 3.8σ using only the closure phases, while we only have a 2.3σ detection using all observables. The possible companion is detected at ρ = 2.1 mas and PA = −136° with a flux ratio f = 0.22%. The component orbiting AW Per is detected at >15σ at ρ = 32 mas and PA = 67° with a flux ratio f = 1.22%. The grids of fit and detection level maps are shown in Fig 3 and 4, and the final fitted parameters are listed in Table 3. For the uncertainties, we used the conservative formalism of Boffin et al. (2014) for all the fitted parameters, i.e., ${\mathit{\sigma }}_{\mathrm{X}}^{\mathrm{2}}\mathrm{=}{\mathit{N}}_{\mathrm{sp}}{\mathit{\sigma }}_{\mathrm{stat}}^{\mathrm{2}}\mathrm{+}\mathrm{0.0001}{\mathit{X}}^{\mathrm{2}}\mathit{,}$(9)where N_sp is the number of spectral channels and X denotes the fitted parameters (i.e., Δx,Δy,f and θ_UD). The first term takes into account that the spectral channels are almost perfectly correlated, and the second term comes from the fact that the wavelength calibration is only precise at a 1% level. The parameter σ_stat is the statistical error from the bootstrapping technique (bootstrap on the calibrated data with replacement) using 10 000 bootstrap samples (also included in CANDID). We then took from the distributions the median and the maximum value between the 16% and 84% percentiles as the uncertainty (although the distributions were roughly symmetrical).

There are additional significant peaks in the RT Aur maps, with a detection level >3σ, however, these are spuriously produced by the (u,v) coverage and the presence of the companion. It is worth mentioning that RT Aur was only observed for one hour (one sequence), and we need more data to confirm the presence of the companion. A more complete discussion about the detected companions is presented in Sect. 4.

Fig. 5 Same as Fig. 3, except we analytically removed the companion.

2.4. Undetected companions

We also observed with MIRC the Cepheids SU Cas (HD 17463, P = 1.95 day) and T Vul (HD 198726, P = 4.44 day). We used five and six telescope configurations, using the same instrument setup and calibration procedures as explained in the previous section. Table 1 lists the journal of these observations.

We did not detect any significant companions around these stars, either using all of the observables or only the CP, i.e., no more than 1.9σ and 2.9σ for SU Cas and T Vul, respectively. However, we were able to measure their angular diameters, which are listed in Table 3. A more detailed discussion of these binary systems is presented in Sect. 4.

3. Detection limit of high-contrast binaries

It is important to be able to check the dynamic range that can be reached with a given set of data and any interferometric combiner. Absil et al. (2011) already presented a method to set detection limits for the VLTI/PIONIER instrument, but it has some shortcomings, that is why we propose a more robust formalism.

Absil’s method:

their method is based on comparing a uniform disk model with a binary model for each position (Δα,Δδ) in the grid. They then checked whether the probability of the binary model is consistent with the data using, $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathrm{\Delta }\mathit{\alpha ,}\mathrm{\Delta }\mathit{\delta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathrm{CDF}}_{\mathit{\nu }}\left(\frac{\mathit{\nu }{\mathit{\chi }}_{\mathit{r,}\mathrm{bin}}^{\mathrm{2}}\mathrm{\left(}\mathrm{\Delta }\mathit{\alpha ,}\mathrm{\Delta }\mathit{\delta }\mathrm{\right)}}{{\mathit{\chi }}_{\mathrm{UD}}^{\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(10)We notice that this equation has a different ratio in the CDF than Eq. (8). This is because we assumed that the binary model is the true model, while Absil et al. (2011) assumed the uniform disk as the true model. In theory, both equations should lead to the same results, however, as we see in subsequent sections, their method is more sensitive to biased data and can sometimes lead to under- or overestimated detection limits.

Our method:

we suggest an alternative method, which is based on the injection of a companion into the data at each astrometric position with different flux ratios. As we inject a companion, we therefore know that the binary model should be the true model, and we can use Eq. (8)to obtain the probability of the binary model to be the true model. We introduced this method because we think it is more robust, as we demonstrate in the next section.

In this section, we first introduce approximate formulae for high-contrast companions. We then explain how we inject an additional component into the data and derive the detection limits.

3.1. High-contrast approximation

The complex visibility for a binary system composed of a resolved primary star and an unresolved component is given by Eq. (1). The squared visibility is given by, $\begin{array}{ccc}{\mathit{V}}^{\mathrm{2}}& \mathrm{=}& \mathrm{\left|}\begin{array}{c}\mathrm{˜}\\ \mathit{V}\end{array}\mathit{.}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}^{\mathrm{\ast }}\end{array}\mathrm{\right|}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{f}}^{\mathrm{2}}}\mathrm{\left|}\mathrm{\right(}{\mathit{V}}_{\mathit{\star }}\mathrm{+}\mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathit{f}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }}\mathrm{\left)}\mathrm{\right(}{\mathit{V}}_{\mathit{\star }}\mathrm{+}\mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathit{f}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\varphi }}\mathrm{\left)}\mathrm{\right|}\\ & \mathrm{=}& \left|\frac{{\mathit{V}}_{\mathit{\star }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathit{f}{\mathit{V}}_{\mathit{\star }}\cos \mathit{\varphi }\mathrm{+}\mathit{G}\mathrm{(}\mathit{\zeta }{\mathrm{)}}^{\mathrm{2}}{\mathit{f}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{2}}}\right|\mathit{,}\end{array}$(11)where we kept the f² term to avoid having an (small) offset (see next section). For a high-contrast companion, i.e., for f ≪ 1, we can approximate the bispectrum at the first order in f as, $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathit{B}\end{array}& \mathrm{=}& \frac{\mathrm{(}{\mathit{V}}_{\mathrm{1}}\mathrm{+}\mathit{G}\mathrm{(}{\mathit{\zeta }}_{\mathrm{1}}\mathrm{\left)}\mathit{f}{\mathrm{e}}^{\mathrm{i}{\mathit{\varphi }}_{\mathrm{1}}}\mathrm{\right)}\mathrm{(}{\mathit{V}}_{\mathrm{2}}\mathrm{+}\mathit{G}\mathrm{(}{\mathit{\zeta }}_{\mathrm{2}}\mathrm{\left)}\mathit{f}{\mathrm{e}}^{\mathrm{i}{\mathit{\varphi }}_{\mathrm{2}}}\mathrm{\right)}\mathrm{(}{\mathit{V}}_{\mathrm{3}}\mathrm{+}\mathit{G}\mathrm{(}{\mathit{\zeta }}_{\mathrm{3}}\mathrm{\left)}\mathit{f}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\varphi }}_{\mathrm{3}}}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{3}}}\\ & \mathrm{~}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{B}}_{\mathit{\star }}\end{array}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{3}}}\left[\mathrm{1}\mathrm{+}\mathit{f}\left(\frac{\mathit{G}\mathrm{\left(}{\mathit{\zeta }}_{\mathrm{1}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\varphi }}_{\mathrm{1}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{1}\mathit{\star }}\end{array}}\mathrm{+}\frac{\mathit{G}\mathrm{\left(}{\mathit{\zeta }}_{\mathrm{2}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\varphi }}_{\mathrm{2}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{2}\mathit{\star }}\end{array}}\mathrm{+}\frac{\mathit{G}\mathrm{\left(}{\mathit{\zeta }}_{\mathrm{3}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}{\mathit{\varphi }}_{\mathrm{3}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{3}\mathit{\star }}^{\mathrm{\ast }}\end{array}}\right)\right]\\ & \mathrm{~}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{B}}_{\mathit{\star }}\end{array}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{3}}}\mathit{Z,}\end{array}$(12)where $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{B}}_{\mathit{\star }}\end{array}\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{V}}_{\mathrm{1}\mathit{\star }}\end{array}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{V}}_{\mathrm{2}\mathit{\star }}\end{array}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{V}}_{\mathrm{3}\mathit{\star }}^{\mathrm{\ast }}\end{array}$.

The bispectrum amplitude and the closure phase are then estimated as, $\begin{array}{ccc}{\mathit{B}}_{\mathrm{amp}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{3}}}\mathrm{\left|}\begin{array}{c}\mathrm{˜}\\ {\mathit{B}}_{\mathit{\star }}\end{array}\mathrm{\right|}\mathrm{\left|}\mathit{Z}\mathrm{\right|}\\ \mathrm{CP}& \mathrm{=}& \arg \mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {\mathit{B}}_{\mathit{\star }}\end{array}\mathrm{\right)}\mathrm{-}\arg \mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathit{.}\end{array}$

3.2. Adding/removing a component

If no companion is detected in the data, we can assume that the measured values only represent the primary star, i.e., a uniform disk, and we can substitute the index “⋆” in the previous equations by the index “obs” (note that the function Z also depends on the primary star visibility). It is now simple, from Eqs. (11), (13)and (14), to inject a companion to the observed data, which corresponds to the following equations for the following new observables: $\begin{array}{ccc}& & {\mathit{V}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{obs}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{G}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\mathit{f}{\mathit{V}}_{\mathrm{obs}}\cos \mathit{\varphi }\mathrm{+}\mathit{G}\mathrm{(}\mathit{\zeta }{\mathrm{)}}^{\mathrm{2}}{\mathit{f}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{B}}_{\mathrm{amp}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}{\mathrm{)}}^{\mathrm{3}}}{\mathit{B}}_{\mathrm{amp}\mathit{,}\mathrm{obs}}\mathrm{\left|}\mathit{Z}\mathrm{\right|}\mathit{,}\\ & & \end{array}$The original “oifits” files provide all the necessary spatial and spectral information to reconstruct the individual phases and visibilities. We also take the coherence loss effect caused by the spectral smearing of the companion into account.

Inversely, if a component is detected, we can use this formula to analytically remove the companion (with a negative flux ratio), and then check for another possible fainter component or estimate the detection limit to rule out any other companions. This step is critical to obtain unbiased detection limits. As an example in Fig. 5, we analytically removed the detected component orbiting RT Aur, and notice that there is no significant detection of a third component and that the other χ² minima are not additional companions.

To check our approximation, we created a model of a single star with a uniform disk of 1 mas (3h of PIONIER observations with a point every 30 min), for which we added a companion at a position Δα = Δδ = 50 mas with a flux ratio f = 5%. We then compare this model with a true binary model using the same parameters (Eq. (1)). The difference between the approximation and the true model is shown in Fig. 6. We found for the amplitude of the bispectrum a standard deviation of the relative error <0.3%, and a difference <0.005° for the closure phase, which is lower than the achievable interferometric accuracy. Closer companions with the same or higher contrast (i.e., <5%) lead to smaller errors.

Fig. 6 Difference between the true binary model and the approximation model. For the bispectrum, the abscissa is the maximum of the three spatial frequencies.

We kept the f² term in the squared visibility, otherwise it gives an offset in the relative error of about 0.3% for f = 5%, and decreases with decreasing flux ratios. Although this value is negligible compared to the current possible data accuracy (ΔV²/V² = ΔB_amp/B_amp ~ 2%, ΔCP ~ 0.5°), the offset is larger for lower contrasts (~1.2% with f = 10%). As no approximation is made for V², the formula is valid for any flux ratio. We therefore also checked the bispectrum approximation for brighter components, and found that we can use it up to a flux ratio of 50% only using CP and V². We listed in Table 4 the relative error for different flux ratios and various positions (1 <ρ< 50 mas and 0 < PA < 2π). We see that CANDID can also be used for low-contrast companions.

3.3. Estimating the detection limit

We implemented in CANDID a tool to estimate the detection limit for a given set of interferometric data. The method is to compare the χ² obtained for a model without a companion and the χ² obtained for a model with an injected companion. The method does not work, a priori, with an already detected companion, first because the previous approximation might not be valid and also because the resulting detection limit would be biased by a systematic offset related to the flux ratio between the components. This means that any detected component has to be analytically removed first.

The method is based on an N × N grid with a range of ±50 mas, for which the minimum N depends on the optimum resolution estimated from the χ² map (see Sect. 2). At each point in the grid, we inject a companion with various flux ratios and we compute the ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$. As we know that the true model is the binary model (because we injected a companion), we used our previous equation, Eq. (8), to estimate the number of sigmas for each flux ratio. We then interpolated the flux ratio values at 3σ, which we set as the significance level. This means that lower flux ratios are not detected significantly. Doing this for all points in the grid, we then have a 3σ detection limit map for the flux ratio. To have a quantitative estimate of the sensitivity limit with respect to the separation, we estimated a radial profile, f_3σ(r), using the 90% completeness level (i.e., 90% of all possible positions) from the cumulated histogram in rings for all azimuths. This tool also includes parallel processing to make the calculation faster.

In theory, for uncorrelated data with Gaussian statistics, our method should lead to the same results as Absil et al. (2011) in terms of the detection limit. However, real data are often biased by different sources (atmospheric turbulence, mechanical vibrations, ...). We performed two tests to compare both methods, one including uncorrelated Gaussian noise and another with a noise model, i.e., correlated non-Gaussian noise. We used all observables and three data sets for each test. The Aspro2 software⁴ were used to create these synthetic data sets, however, it does not have implemented the bandwidth smearing effect. We therefore did not take it into account in CANDID (this would not change the conclusion of our test).

First test:

we created the first “ideal” synthetic data set (i.e., without noise) representing a uniform disk of 1 mas (3 h observation with three spectral channels with the PIONIER instrument). We then added uncorrelated Gaussian noise and estimated the detection limit using the Absil’s method and our formalism. As expected, we see in Fig. 7 that both methods provide the same results. We then created a second data set by adding a component at Δα = Δδ = 5 mas with f = 2% as a source of bias to the previous model. The companion creates a departure from the ideal measurements (like noise) but in such a way that the observations at different spatial frequencies have a correlated departure (hence, correlated noise). We then estimated the detection limits for this second data set, and although the trend between the two methods are similar, there is an offset of ~2% because of the presence of the faint companion. Finally, we analytically removed the companion to get our third data set. Comparing the three data sets, i.e., estimating the total variation of the detection limit for each radius, we noticed that both methods vary in a similar way. This was expected as all of the noise sources inserted have Gaussian statistics. However, it is worth mentioning that the scatter tends to be a bit larger for the Absil’s method for increasing flux ratios and wider separations of the companion (due to the correlated noise introduced by the component).

Fig. 7 Comparison between the CANDID detection limit tool and the formalism of Absil et al. (2011) for a uniform disk model with Gaussian noise.

Fig. 8 Comparison between the CANDID detection limit tool and the formalism of Absil et al. (2011) for a uniform disk model with a noise model, biased by the presence of a faint companion.

Second test:

we created two new data sets similar to the first test, i.e., one for a uniform disk of 1 mas and another for a binary system with Δα = Δδ = 5 mas and f = 2% (with average atmospheric conditions, which is an option for the noise model). The only difference is that the noise is no longer Gaussian; it is represented by a more complicated model that includes the instrument response, atmospheric turbulence, photon, and detector noise⁵. As shown in Fig. 8, the absil’s method is more sensitive to the presence of the companion, which in this case gives a higher sensitivity limit. The third data set is created by removing the component using our approximation. We then compared the variation of the detection limit between all of the data sets, as described previously, i.e., the minimum and maximum value between the three synthetic data sets at each separation. Figure 9 shows that the scatter from the Absil’s method is larger. We performed additional tests with lower flux ratios and at different positions, and the scatter from the Absil’s method is larger than our formalism most of the time. We therefore conclude that our method is more robust to biased data when estimating detection limits.

In our tests, the limits derived from the injection method are lower than when we use Absil’s method, but this is not always the case for real data. The injection method occasionally results in higher sensitivity limits, depending on the magnitude of the biases. We find that Absil’s method may under- or overestimate the detection limits depending on the data set. However, in some cases when the data are not affected by significant biases, both methods give similar results.

Fig. 9 Total variation (minimum and maximum values) between the three synthetic data sets with a noise model for both methods: a single star, a single star + an unresolved component as source of bias, and the previous model with the component analytically removed. The variation has been normalized to unity.

Fig. 10 Flux ratio detection limit at 3σ for a second companion around V1334 Cyg (for 2012 Oct. 01).

4. Detection limit of our binary Cepheids

In this section, we set the detection limit for the sample of Cepheids previously presented. For V1334 Cyg, AX Cir, RT Aur, and AW Per, where a companion is detected, we removed analytically the companion first. Although both methods are estimated by CANDID, all of the given detection limits are derived from our injection method. We present two detection limits for each star, one using all of the observables, and one using only the closure phases. We listed three different values, the average for r< 50 mas, r< 25 mas and r> 25 mas, which can be relevant when the limit increases with r. All of the final detection limits are listed in Table 5, where the values are conservative as they correspond to the mean plus the standard deviation for the given radius range.

From an evolutionary timescale point of view, most of the companions should be stars close to the main sequence. We therefore set upper limits for the spectral type of the companion assuming it is on the main sequence, based on their H-band luminosities.

V1334 Cyg:

this binary system contains a visual and a spectroscopic component. While the visual companion is >150 mas, the close component was spatially resolved using interferometry by Gallenne et al. (2013b), and has a flux ratio of ~3.1% in H (more details on this companion are presented in Gallenne et al. 2013b). There is no evidence of a third component so far. After removing the close companion analytically, we estimated the dynamic range as explained in the previous section. Between the two epochs, we reached a maximum average sensitivity limit of f_3σ = 1.54%. The limit is lower for the first epoch because the atmospheric conditions were better than the second epoch. The contrast upper limit at 3σ is shown in Fig. 10, with the average values for the two epochs in Table 5. This converts to magnitude difference of Δm_H> −2.5log f_3σ = 4.5 mag. Using the distance d = 683 pc from the K-band P–L relation for first overtone (FO) pulsators (Bono et al. 2002) and the average magnitudes of the Cepheid H = 4.66 mag (Cutri et al. 2003), we can exclude the presence of additional companions with a spectral type earlier than B7V stars.

AX Cir:

the spectroscopic companion was first detected by Lloyd Evans (1982), and was spatially resolved for the first time with LBI by Gallenne et al. (2014, with a more detailed discussion about this companion). The average detection limits are listed in Table 5. We reached a dynamic range of f_3σ = 0.36%, corresponding to Δm_H> 6.1 mag if another companion is present. With d = 500 pc (from the K-band P–L relation of Storm et al. 2011) and H = 4.66 mag for the Cepheid (Cutri et al. 2003), we can rule out any other component with a spectral type earlier than A5V star.

RT Aur:

the binary nature of this short-period Cepheid is still uncertain, however we might have detected it for the first time. Some authors suggested the presence of an early-type companion (Janot-Pacheco 1976; Balona 1977), while others did not see evidence of an additional component (Harris 1981; Gieren 1985). Leonard & Turner (1986) summarized various studies from that time about the possible companion orbiting RT Aur. Based on spectra from the International Ultraviolet Explorer satellite (IUE), Evans (1992a) did not report any detection and showed that any main-sequence secondary has to be cooler than an A4 star. Radial velocity measurements do not show any variations from orbital motion. Recently, Turner et al. (2007) studied long-term photometric light curves and reported a sinusoidal trend consistent with a light travel time effect in the binary system. The companion we might have resolved has a flux ratio of 0.21 ± 0.12%, i.e., Δm_H = 6.7 ± 0.6 mag, leading to m_H = 10.6 ± 1.0 mag (using for the Cepheid at our given phase m_H = 3.94 ± 0.01 mag from Monson & Pierce 2011). Using the K-band P–L relation to get d = 428 pc (Storm et al. 2011), we estimate its spectral type to be later than an F1 star. This is compatible with the A4V star upper limit determined by Evans (1992a) from the International Ultraviolet Explorer (IUE) spectra. However, we need more observations to confirm the existence of the companion as we are at the sensitivity limit of this data set.

Removing this possible companion from our interferometric data, we estimated the detection limit. We reached a minimum dynamic range of 0.47%, corresponding to Δm_H> 5.8 mag. This is a 90% completeness azimuthal value, which is a conservative value, and the component at a given position with a slightly higher contrast might be detected (as in this case). This limit allows us to exclude the presence of an additional companion with a spectral type earlier than F0V.

AW Per:

this Cepheid is a spectroscopic binary with an orbital period of ~40 yr. First discovered by Miller & Preston (1964), it took several years to derive the first orbit from radial velocity data (Evans 1989). It is likely that this companion is itself a binary, as argued by Evans (1989), Evans et al. (2000), because the magnitude difference between the Cepheid and its companion is not consistent with equal masses and predictions from evolutionary tracks. Unfortunately, we do not have enough angular resolution and sensitivity to detect this third companion. The properties of the brightest companion were studied based on IUE spectra by Evans (1994, 1995), who found its spectral type to be B8.3V. This is in agreement with our detection with a flux ratio f = 1.22 ± 0.30%, i.e., a spectral type in the range B6-B9V (using d = 853 pc from a P–L relation and K = 4.63 mag for the Cepheid). We estimated for the companion m_H = 9.6 ± 0.3 mag (using for the Cepheid at our given phase m_H = 4.84 ± 0.01 mag from Monson & Pierce 2011). Massa & Evans (2008) also determined the angular separation of the component for another epoch, which allows us later, with more astrometric points from interferometry, to estimate all the orbital elements, including the inclination and the semi-major axis (see also Gallenne et al. 2013a).

We derived a maximum sensitivity limit at 3σ of 0.62%, i.e., Δm_H> 5.5 mag. We can therefore exclude any additional components with a spectral type earlier than B9V.

SU Cas:

the binary nature of this Cepheid is still ambiguous. A component was first detected by Evans (1985) based on studying the CAII H and K lines, and then a spectral type of B9.5V was determined from IUE spectra (Evans 1991). Although the location of SU Cas in a two color diagram is consistent with the presence of a companion, the radial velocity measurements do not show convincing evidence. Szabados (1991) found four possible orbital periods from the observations available at that time. Later, (Gorynya et al. 1996) compiled more data and suggested SU Cas as a possible spectroscopic binary. They derived an orbital period of 408 days with an eccentricity of e = 0.43, but Groenewegen (2008) could not confirm the eccentricity with a larger data set and found a period of 407 days, fixing e = 0. Recently, Evans et al. (2013) reanalyzed special dates of Gorynya data, where the velocity difference was supposed to be the largest from the derived orbit, but they concluded that the orbital motion from radial velocity data do not show any convincing detection.

Using the distance d = 392 pc (from the K-band P–L relation of first overtone pulsators of Bono et al. 2002) and a Cepheid magnitude H = 4.27 mag (Cutri et al. 2003), a B9.5V companion should give a flux ratio ~1.8% in the H band, which should be detectable with MIRC. Taking the velocity amplitude of K₁ = 1 km s^-1 from Groenewegen (2008, while Gorynya et al. 1996 derived a velocity amplitude of 3 km s^-1, we estimated a sini> 0.1 mas. Our nondetection could be explained if the companion was located at < 0.5 mas because it would not have been spatially resolved, but this type of close component would have an effect on the radial velocities (unless the orbit is face-on, but the probability of this kind of an orbit is low). The other possibility would be that this companion has a wider orbit. We therefore searched within a 100 mas radius range, but we did not find a significant detection. We then estimated the sensitivity limit for r< 100 mas, and found a maximum 3σ flux ratio of 1.65%. This means that if the B9.5V companion was within 100 mas, we would have detected it. Therefore, if this companion exists, it should have a wider orbit.

The average sensitivity limits for a given radius range are tabulated in Table 5. We reached a mean contrast of 1.37%, which means that Δm_H> 4.6 mag. This converts to an upper limit of the spectral type of an A0V star.

T Vul:

as with SU Cus, the radial velocities of this Cepheid do not show any signature of an orbiting companion, while a hot A0.8V component was detected by Evans (1992b). In the literature, we found contradictory estimates of the orbital period from radial velocity measurements; for instance, Kovacs et al. (1990) found a long period modulation of 738 d, Szabados (1991) estimated P = 1745 d from a larger data set, while Bersier et al. (1994) did not find any orbital motion larger 0.55 km s^-1 using additional more accurate data. They showed that the long period of Szabados (1991) is not compatible and argued that the 738 d period might be an artifact of the time sampling because the observations were only made in autumn. Kiss & Vinkó (2000) reached the same conclusion with additional measurements showing no signature of orbital motion in the radial velocity curve.

From our interferometric observations, we did not detect any companion within 50 mas. The A0.8V component detected from IUE spectra would correspond to a flux ratio in H of ~0.7%. According to our estimated interferometric detection limit, listed in Table 5, this kind of a component is below our detection level. We reached an average sensitivity limit of 1.04%, and we can therefore exclude any other possible companion with a spectral type earlier than a B9V star.

5. Conclusion

We presented an overview of CANDID, a new tool to search for point-source companions and estimate the sensitivity level from interferometric observations using the squared visibilities, closure phases, and amplitude of the bispectrum, when available. CANDID allows us to:

• efficiently detect companions using a grid of fit and determine thedetection level by giving the number of sigmas;

• set the detection limit for a companion in data where a companion has not been detected; and

• set the detection limit for a tertiary companion, in the case where a companion has been detected.

We used CANDID to investigate a sample of binary Cepheids. We first determined the detection level for our previous detections (Gallenne et al. 2013b, 2014) and showed that the components were detected at >13σ for AX Cir and >25σ for V1334 Cyg. We also reported a new detection for AW Per, with a detection level >15σ; the companion is located at ρ = 32 mas and PA = 67°, with a flux ratio of f = 1.22%. The companion orbiting RT Aur might have been detected at 3.8σ, using only the closure phase signal, however, more observations are needed to confirm the presence of this component. Any additional companions were not detected signficantly (i.e., with a detection level >3σ) around these stars. Likewise, no companions were detected around SU Cas and T Vul. From these interferometric data, we were able to set upper limits for the spectral types; we found no components with a spectral type earlier than B7V, A5V, F0V, B9V, A0V, and B9V for V1334 Cyg, AX Cir, RT Aur, AW Per, SU Cas, and T Vul, respectively.

The fitting procedure of CANDID also allowed us to measure the uniform disk angular diameters of the new Cepheids observed. We found for RT Aur θ_UD = 0.699 ± 0.011 mas (at the pulsation phase φ = 0.32), for AW Per θ_UD = 0.627 ± 0.018 mas (at φ = 0.52), for SU Cas θ_UD = 0.609 ± 0.043 mas (at φ = 0.77), and θ_UD = 0.608 ± 0.013 mas and θ_UD = 0.635 ± 0.018 mas for T Vul (at φ = 0.27 and φ = 0.12).

We demonstrated that the approximation we used to analytically inject a companion and estimate the detection limits is valid (i.e., error <0.5%) for contrasts % if we use all of the observables, and up to % using only the squared visibilities and the closure phases. This makes CANDID a useful tool for analyzing long-baseline interferometric observations of binary star systems.

Finally, this work demonstrates the capabilities of the MIRC and PIONIER instruments, which can reach a dynamic range of 1:200, depending on the angular distance of the companion and the (u,v) plane coverage.

Acknowledgments

The authors would like to thank the CHARA Array and Mount Wilson Observatory staff for their support. Research conducted at the CHARA Array is funded by the National Science Foundation through NSF grant AST-0908253, by Georgia State University, the W. M. Keck Foundation, the Packard Foundation, and the NASA Exoplanet Science Institute. The authors also thank all the people involved in the VLTI project. A.G. acknowledges support from FONDECYT grant 3130361. J.D.M. acknowledges funding from the NSF grants AST-0707927, AST-0807577, and AST 1108963. We also acknowledge the STSCI grant HST-GO13841.006-A and HST-GO-13454.07-A. W.G. and G.P. gratefully acknowledge financial support for this work from the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) PFB-06/2007. Support from the Polish National Science Centre grant MAESTRO DEC-2012/06/A/ST9/00269 and the Polish Ministry of Science grant Ideas Plus (awarded to G. P.) is also acknowledge. W.G. also acknowledges financial support from the Millenium Institute of Astrophysics (MAS) of the Iniciativa Cientifica Milenio del Ministerio de Economía, Fomento y Turismo de Chile, project IC120009. A.G. also acknowledge computational resources from FONDECYT grant 1130521 and 1130721. We also acknowledge support from the ECOS/Conicyt grant C13U01. PIONIER is funded by the Université Joseph Fourier (UJF), the Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), the Agence Nationale pour la Recherche (ANR-06-BLAN-0421, ANR-10-BLAN-0505, ANR-10-LABX56), and the Institut National des Science de l’Univers (INSU PNP and PNPS). The integrated optics beam combiner is the result of a collaboration between IPAG and CEA-LETI based on CNES R&T funding. This research received the support of PHASE, the high angular resolution partnership between ONERA, Observatoire de Paris, CNRS, and University Denis Diderot Paris 7. This work made use of the SIMBAD and VIZIER astrophysical database from CDS, Strasbourg, France and the bibliographic informations from the NASA Astrophysics Data System. This research has made use of the Jean-Marie Mariotti Center SearchCal and ASPRO services, co-developed by FIZEAU and LAOG/IPAG, and of CDS Astronomical Databases SIMBAD and VIZIER.

References

All Tables

All Figures

	Fig. 1 ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$ map of the local minima (left) and detection level map (right) of V1334 Cyg for the observations on 2012-10-01. The yellow lines represent the convergence from the starting points to the final fitted position. The maps were reinterpolated in a regular grid for clarity. The axis limit was chosen according to the location of the companion.
In the text

	Fig. 2 Same as Fig. 1, except for AX Cir for the observations on 2013-07-14.
In the text

	Fig. 3 Same as Fig. 1, except for RT Aur using only the closure phase.
In the text

	Fig. 4 Same as Fig. 1, except for AW Per.
In the text

	Fig. 5 Same as Fig. 3, except we analytically removed the companion.
In the text

	Fig. 6 Difference between the true binary model and the approximation model. For the bispectrum, the abscissa is the maximum of the three spatial frequencies.
In the text

	Fig. 7 Comparison between the CANDID detection limit tool and the formalism of Absil et al. (2011) for a uniform disk model with Gaussian noise.
In the text

	Fig. 8 Comparison between the CANDID detection limit tool and the formalism of Absil et al. (2011) for a uniform disk model with a noise model, biased by the presence of a faint companion.
In the text

	Fig. 9 Total variation (minimum and maximum values) between the three synthetic data sets with a noise model for both methods: a single star, a single star + an unresolved component as source of bias, and the previous model with the component analytically removed. The variation has been normalized to unity.
In the text

	Fig. 10 Flux ratio detection limit at 3σ for a second companion around V1334 Cyg (for 2012 Oct. 01).
In the text