Kernel phase imaging with VLT/NACO: high-contrast detection of new candidate low-mass stellar companions at the diffraction limit

Like most photon counting devices, NACO’s infrared detector CONICA suffers from a non-linear response when the pixel counts approach the saturation threshold (16400 counts for uncorrelated high well depth mode³³3This is the standard imaging mode in the L’ band ($3.8\text{textmu m}$) and all data cubes which we analyze have been taken in this mode. according to the NACO Quality Control and Data Processing⁴⁴4https://www.eso.org/observing/dfo/quality/NACO/qc/detmon_qc1.html, with a more conservative 16000 counts used in our analysis). As kernel phase is an interferometric technique for which the fringes are coded spatially on the detector, it is very important to characterize the pixel to pixel response. Moreover, many of the data cubes which we analyze in Section 3 feature saturated point spread functions (PSFs) which we want to correct for non-linearity before reconstructing their core (cf. Section 2.1.6).

In order to compute the detector linearization correction we download all frames of type “FLAT, LAMP, DETCHECK” and uncorrelated high well depth mode from 2016 March 23 and 2016 September 24 (which are closest in time to the observation of the earliest and the latest data cube which we analyze) from the VLT/NACO archive. We sort them by integration time and compute the median pixel count over all frames for each individual integration time (masking out the broken stripes in the lower left quadrant of CONICA). Then, we plot the median pixel count in dependence of the integration time $t$, fit a linear polynomial $f\left(t\right)$ to all data points with less than 8500 counts (end of the linear regime for uncorrelated high well depth mode) and a cubic polynomial $g\left(t\right)$ to all data points with less than 16000 counts (saturation threshold, cf. left panel of Figure 1). We linearize the detector using a continuously differentiable piecewise polynomial approach $h$ to the correction curve $f\left(g\right)$ with a linear function up to 8500 counts and a cubic polynomial between 8500 and 16000 counts (cf. right panel of Figure 1).

For each observation block (OB) we compute master darks from the associated dark frames as the median of all dark frames with a unique set of size and exposure time. Then we compute a bad pixel map for each master dark based on the frame by frame median and variance of each pixel’s count. Therefore, we first compute two frames:

1. The absolute difference between the master dark and the median filtered master dark.

2. The absolute of the median subtracted variance dark.

Then, we identify bad pixels in each of these frames based on their difference to the median of these frames. For frame (i) we classify pixels which are above 10 times the median as bad, for frame (ii) pixels which are above 75 times the median. Note that these thresholds were identified empirically. From the median subtracted dark frames, we estimate the readout noise as the mean over each frame’s pixel count standard deviation.

We proceed similar for the flat frames, but also group them by filter as well as size and exposure time, subtract a master dark with matching properties (i.e. similar size and exposure time) from each master flat and normalize it by its median pixel count.

The data cubes which we analyze in Section 3 consist of 100 frames of $0.2\text{s}$ exposure. For each data cube, we reject the first frame (which we find to consistently suffer from a bias), so that there are 99 frames left. Note that NACO appends the median of all 100 frames at the end of each data cube which is also rejected here. Before proceeding, we also flag the saturated pixels in each frame which are all pixels with more than $h\left(16000\right)$ counts.

We clean each frame of a data cube individually by subtracting a master dark with matching properties (i.e. similar size and exposure time), dividing it by a master flat with matching properties (i.e. similar size, exposure time and filter), correcting bad pixels (which are bad pixels from the master dark or the master flat) with a median filter of size five pixels and performing a simple background subtraction by subtracting the median pixel count of the frame from each pixel. A typical result is shown in the left panel of Figure 2, where residual systematic noise (mainly from the detector) is still clearly visible.

In order to mitigate the residual systematic noise from the detector and the sky background we perform a dither subtraction. After cleaning all data cubes within one OB, we find for each data cube (which we will here call data cube A) the data cube B with the target furthest away (on the detector) and subtract its median frame from each frame of data cube A. The new bad and saturated pixel maps are then the logical sums of those from both involved data cubes. After this step the residual noise appears like Gaussian random noise as is shown in the right panel of Figure 2.

Our typical performance is a pixel count standard deviation of $\sim 36=4.4+(158/\text{s}\cdot 0.2\text{s})$ outside of $10\lambda /D$ from the center of the PSF in $0.2\text{s}$ exposure, where $4.4$ is the detector readout noise, $\lambda$ is the observing wavelength ($3.8\text{textmu m}$ for the L’ band) and $D$ is the diameter of the primary mirror ($8.2\text{m}$ for the VLT).

Our reconstruction of saturated pixels is based on an algorithm described in Section 2.5 of ireland2013. This technique also identifies and corrects residual bad pixels, with no more than 10 additional bad pixels corrected in a typical frame. First, we crop all frames to a size of 96 by 96 pixels ($\sim 2.6{\text{arcsec}}^2$) centered on the target. Then, we correct bad and saturated pixels for each frame separately by minimizing the Fourier plane power $\left|f_Z\right|$ outside the region of support $Z$ permitted by the pupil geometry. Let $B_Z$ be the matrix which maps the bad and saturated pixel values $x$ onto the Fourier plane domain $Z$, then

where $b$ are the corrections to the bad and saturated pixel values $x$ (i.e. the corrected pixel values are $x-b$) and ${ϵ}_Z$ is remaining Fourier plane noise. We solve for $b$ using the Moore-Penrose pseudo-inverse of $B_Z$, i.e.

Since a broad-band filter was used for the observations, but we use a monochromatic central filter wavelength in our analysis and also blur the edge of the pupil through the use of a windowing function, we use a sightly larger pupil diameter to define this region $Z$ of $10\text{m}$ here. In fact, the only important thing for recovering the Fourier plane phase is that the Fourier plane power outside the region of support permitted by the pupil geometry is minimized, so using a larger pupil diameter just assures this in case of low quality data and is a conservative choice, especially in the case of our data which is far from the Nyquist sampling criterion.

Sometimes, the remaining Fourier plane noise ${ϵ}_Z$ can be significant, which is why we repeat the entire correction process up to 15 times for each frame. After each iteration, we look for remaining bad pixels by:

1. Computing the Fourier transform of the corrected frame from the previous iteration.

2. Windowing this frame by the Fourier domain $Z$.

3. Computing the inverse Fourier transform of this frame.

4. Identifying remaining bad pixels in this frame based on their difference to the median filtered frame.

If no remaining bad pixels are identified, we terminate the iteration.

A cross-section of a saturated PSF before and after performing the reconstruction is shown in Figure 3. Obviously, this reconstruction cannot reveal any structure or companions hidden behind saturated pixels, but it allows us to perform our kernel phase analysis on saturated data cubes which would otherwise suffer from high Fourier plane phase noise (cf. Figure 4). Please note that a method from the class of least squares spectral analysis techniques (i.e. image plane fringe fitting) may be more robust in dealing with bad pixels, but would require the simultaneous fitting of all Fourier plane phases and amplitudes and is therefore beyond the scope of this paper, although it is a promising approach for future work.

As explained in the Introduction, a high Strehl ratio is essential for the kernel phase technique in order for the mathematical framework (i.e. the linearization of the Fourier plane phase, cf. Equation 1) to be valid. Therefore, we select frames with sufficient Strehl ratio based on their peak pixel count. For each data cube, we first compute the median peak count of the 10% best frames. Then, we reject all frames with a peak count below 75% of this value. Using this dynamic threshold is better than simply rejecting a fixed fraction of the frames (e.g. law2006) because it can correctly account for a sudden drop in the Strehl ratio like shown in Figure 5. Note that we consider the peak pixel count after performing the PSF reconstruction (cf Section 2.1.6) here.

In order to extract the kernel phase from VLT/NACO data we first need to construct a model for the VLT pupil (i.e. split the primary mirror into an interferometric array of virtual sub-apertures). We sample 140 virtual sub-apertures on a square grid with a pupil plane spacing of $0.6\text{m}$, which is approximately half the Nyquist sampling of $\lambda /\alpha \approx 0.3\text{m}$, where $\lambda =3.8\text{textmu m}$ is the observing wavelength and $\alpha =2.610\text{arcsec}$ is the image size (96 pixels). Our VLT pupil model is shown in the left panel of Figure 6 and based on an $8.2\text{m}$ primary mirror with a $1.2\text{m}$ central obscuration. Another advantage of kernel phase over sparse aperture masking is the dense Fourier plane coverage which is shown in the right panel of Figure 6.

The extraction of the Fourier plane phase and the computation of the kernel phase relies on a python package called XARA⁵⁵5https://github.com/fmartinache/xara (eXtreme Angular Resolution Astronomy, martinache2010; martinache2013). XARA has been designed to process data produced by multiple instruments assuming that the images comply to the kernel phase analysis requirements of proper sampling, high-Strehl (boosted by our frame selection procedure described in Section 2.1.7), and non-saturation (restored by the procedure described in Section 2.1.6). The discrete achromatic representation of the VLT aperture (i.e. our pupil model) is used by XARA to compute the phase transfer matrix $A$ and the associated left kernel operator $K$ via a singular value decomposition of $A$.

With the added knowledge of the detector pixel scale and the observing wavelength, the discrete model is scaled so that the Fourier plane phase at the expected $(u,v)$ coordinates can be extracted by a discrete Fourier transform. For the small aberration hypothesis to remain valid, the data must be properly centered prior to the Fourier transform. Failure to do so will leave a residual Fourier plane phase ramp that can wrap and lead to meaningless kernel phases (cf. left panels of Figure 7). XARA offers several centering algorithms. It is crucial to carefully choose from the available options depending on the requirements coming from the data. For our extensive ground-based data set for example, we find that minimizing directly the Fourier plane phase which is extracted by XARA is most robust and the offered sub-pixel re-centering is very valuable (cf. right panels of Figure 7) due to an increased level of pupil plane phase noise from the atmosphere and the bright background (if compared to space-borne data).

Moreover, virtual baselines near the outer edge of the Fourier coverage suffer from low power as they are only supported by very few baselines, i.e. have small redundancy. The phase measured for these baselines is systematically noisier and needs to be excluded from the model to prevent the noise to propagate into the estimation of all kernel phases. This can be achieved using the baseline filtering option implemented in XARA. In our case, baselines of length greater than $7.0\text{m}$ and the corresponding rows of $A$ are eliminated prior to the computation of $K$. Some of the theoretically available kernel phases are lost but the remaining kernel phases can nevertheless be used just like for the complete model.

Finally, to limit the impact of readout noise in regions of the image where little signal is present, frames are windowed by a super-Gaussian ($g\left(r\right)=\exp -(r/r_0{)}^4$) with a radius $r_0=25$ pixels, effectively limiting our field of view to $\sim 1000\text{mas}$. Note that Section LABEL:sec:windowing_correction will further comment on the effect of this window and how it can affect contrast estimates for detections at large separations.

For estimating the uncertainties, we compute the kernel phase covariance ${\Sigma }_{\theta }$ for each frame $d$ from its photon count variance ${\Sigma }_d=g\cdot w^2\cdot (d+b)$ in units of (photo-electrons)${}^2$, where $g$ is the detector gain ($g=9.8$ for uncorrelated high well depth mode), $w$ is the super-Gaussian window, $d$ is the cleaned and re-centered frame and $b$ is its background (from the simple background subtraction, cf. Section 2.1.4). Therefore, we first need to find a linear operator $B$ which maps each frame $g\cdot w\cdot d$ in units of photo-electrons to its kernel phase $\theta$. The linear discrete Fourier transform $F$ and the kernel $K$ of the pupil model $R^{-1}\cdot A$ are already linear operators, and the Fourier plane phase $\varphi \left(z\right)$ (of a complex number $z$) can be approximated as $Im\left(z\right)/|z|$ for small angles. Hence, we compute

Note that $B\cdot g\cdot w\cdot d$ would be a small-angle approximation for the kernel phase. Then, we obtain an estimate for the kernel phase covariance by propagating the photon count variance according to

Note that this kernel phase covariance model includes the contribution of shot noise only. Any residual calibration errors not taken into account in the following section are therefore expected to increase the reduced ${\chi }^2$ of any model fitting, potentially to much more than 1.0 in the case of high signal-to-noise data with highly imperfect calibration.

Using the kernel phase covariance ${\Sigma }_{\theta }$ estimated from photon noise according to Section 2.2.3 we compute the best fit contrast $c_{\text{fit}}$ and its uncertainty ${\sigma }_{c_{\text{fit}}}$ for the binary model ${\theta }_{\text{bin}}=K\cdot {\varphi }_{\text{bin}}$ on each position of a discrete $500\times 500\text{mas}$ square grid with spacing $13.595\text{mas}$ (which is half the detector pixel scale of CONICA). In some cases, where we suspect a companion candidate at a larger angular separation, we also extend the grid to $1000\times 1000\text{mas}$.

In the high-contrast regime (where $c\ll 1$), the phase ${\varphi }_{\text{bin}}$ is approximately proportional to the contrast $c$ of the binary model, so is its kernel phase ${\theta }_{\text{bin}}$ (because $K$ is a linear operator). Hence, the ${\chi }^2$ of the binary model ${\chi }_{\text{bin}}^2$ can be approximated as

where $\Theta$ and ${\Theta }_{\text{bin,ref}}$ are vertical stacks of the kernel phase ${\theta }_i$ and the reference binary model ${\theta }_{\text{bin,ref},i}$ of each data cube $i$ and ${\Sigma }_{\Theta }^{-1}$ is a block-diagonal matrix whose diagonal elements are the inverse kernel phase covariances ${\Sigma }_{\theta ,i}^{-1}$ of each data cube $i$, i.e.

The reference binary model ${\theta }_{\text{bin,ref}}$ is the binary model ${\theta }_{\text{bin}}$ evaluated for and normalized by a reference contrast $c_{\text{ref}}=0.001$, i.e.

Finally, we obtain the log-likelihood $\ln L$ for the binary model ${\theta }_{\text{bin}}$ as

The final parameters $p_{\text{fit}}$ for the best fit binary model are then obtained from a least squares search which maximizes the log-likelihood $\ln L$ of the binary model under varying angular separation, position angle and contrast simultaneously. For the least squares search, we use the grid position with the maximal log-likelihood as prior and restrict the search box for the angular separation $\rho$ to $50\text{mas}\le \rho \le 1000\text{mas}$.

Using only the uncertainties from photon noise, it is still difficult to distinguish between residual speckle noise (i.e. third order phase noise in the pupil plane) and real detections at small angular separations. This is the case because the data set which we analyze in Section 3 is very limited in terms of diversity of calibrator PSFs. For this reason, we use an empirical approach as the primary method to determine whether a detection is real or not.

First, we split our targets into candidate detections and calibrators based on their detection significance from photon noise ${\text{SNR}}_{\text{ph}}$ (cf. Section LABEL:sec:close_companion_candidates). For each of the calibrators, we then compute two contrast curves:

1. The azimuthal average $c_{\text{RMS}}$ of the root mean square (RMS) best fit contrast $c_{\text{fit}}$.

2. The azimuthal average $c_{\text{RMS}}^{\text{sub}}$ of the RMS best fit contrast $c_{\text{fit}}^{\text{sub}}$ after subtracting the best fit binary model ${\theta }_{\text{bin}}$ from the measured kernel phase $\theta$.

Here, the assumption is that the calibrators are single stars, so that the ratio of the two RMS contrast curves computed above, i.e.

is a correction factor for the relative contrast of the residual speckle noise. This is illustrated in the left panel of Figure 9.

For each of the candidate detections, we only compute the azimuthal average $c_{\text{RMS}}^{\text{sub}}$ of the RMS best fit contrast $c_{\text{fit}}^{\text{sub}}$ after subtracting the best fit binary model ${\theta }_{\text{bin}}$ (which might or might not be a real detection) from the measured kernel phase $\theta$. Then, we multiply this RMS contrast curve with the mean of the relative speckle contrast $f_{\text{speckle}}$ of all calibrators, i.e.

where the bar denotes the mean, in order to obtain an empirical contrast uncertainty ${\sigma }_{\text{emp}}$ as a function of the angular separation $\rho$ for each candidate detection (cf. right panel of Figure 9). We classify a candidate detection as real if its empirical detection significance ${\text{SNR}}_{\text{emp}}$ is above the $5\sigma$ threshold, i.e.

Furthermore, we obtain empirically motivated uncertainties on the best fit parameters $p_{\text{fit}}$ by multiplying the uncertainties from photon noise ${\sigma }_{p_{\text{fit}}}$ with the ratio $f_{\text{err}}$ of the empirical contrast uncertainty ${\sigma }_{\text{emp}}$ to the contrast uncertainty from photon noise ${\sigma }_{\text{ph}}$ (at the position of the best fit binary model ${\theta }_{\text{bin}}$).

The kernel phase analysis tools described in Sections 2.2, 2.3 and 2.4 are available on GitHub⁶⁶6https://github.com/kammerje/PyKernel. We put a strong focus on applicability to other instruments and an exchangeable kernel phase fits file format.