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

Abstract

1 INTRODUCTION

Direct imaging is vital for studying the outer regions of extrasolar systems that are inaccessible to transit observations and can only be revealed by decades-long, time-consuming radial velocity surveys (e.g. Fischer et al. 2014). It has proven particularly successful in probing our understanding of the formation of gas giant planets (e.g. D’Angelo, Durisen & Lissauer 2010), being able to estimate their mass from their luminosity and age (e.g. Spiegel & Burrows 2012) and resolve their orbit. Although the majority of detected companion candidates are arguably consistent with being emission or scattering from disc material (e.g. LkCa 15, Kraus & Ireland 2012; HD 100546, Quanz et al. 2013; HD 169142, Biller et al. 2014), the recent example of PDS 70 (Keppler et al. 2018) demonstrates that direct imaging of wide-separation but still Solar-system scale planets is possible at relatively moderate contrasts in the vicinity of young stars. This is spurring an ongoing discussion about the nature of planet formation and the commonness of gas giant planets with large orbital distances (e.g. Bowler & Nielsen 2018).

However, direct imaging operates at the resolution and sensitivity limit of the most powerful instruments today (e.g. Pepe, Ehrenreich & Meyer 2014), placing demanding requirements on the observing and the post-processing techniques that are used to uncover faint companions at high contrasts (e.g. angular differential imaging, Marois et al. 2006; point spread function (PSF) subtraction, Lafrenière et al. 2007a; principal component analysis, Amara & Quanz 2012; Soummer, Pueyo & Larkin 2012). Detecting exoplanets from the ground using these techniques has only been made possible by the recent development of extreme adaptive optics systems (e.g. Milli et al. 2016) and is mainly limited by non-common path aberrations that are not sensed by the wavefront control system (e.g. Sauvage et al. 2007). These aberrations manifest themselves as quasi-static speckles on the detector images that can mimic the signal of a companion and place a strong constraint on the achievable contrast at small angular separations (e.g. Fitzgerald & Graham 2006). Hence, directly imaging and studying the formation of gas giant planets on Solar-system scales has been extremely challenging so far (e.g. Bowler 2016) because the nearest star-forming regions lie |${\gtrsim }100\, \text{pc}$| away (e.g. Loinard et al. 2007) where such orbital distances correspond to angular separations of only |${\lesssim }200\, \text{mas}$|⁠.

In this paper, we explore the capabilities of the kernel phase technique (Martinache 2010) for high-contrast imaging at the diffraction limit from the ground. This post-processing technique can be seen as refinement of sparse aperture masking and the closure phase technique (Tuthill et al. 2000). By probing only certain linear combinations of the phase of the Fourier transformed detector images, kernel phase, and sparse aperture masking allow for a robust calibration of the time-varying optical transfer function of the system and a significant mitigation of quasi-static speckles and achieve an angular resolution of |${\lesssim }50\, \text{mas}$| in the near-infrared (i.e. the L’ band, Cheetham et al. 2016). This gives access to objects on Solar-system scales in the nearest star-forming regions (i.e. projected separations of |${\sim }40\, \text{mas}$| for a Jupiter analogue in the Scorpius Centaurus OB association, Preibisch & Mamajek 2008) and has proven successful in directly imaging young exoplanets/disc features (e.g. Kraus & Ireland 2012). The caveat of sparse aperture masking is that the mask blocks |${\gtrsim }85~{{\ \rm per\ cent}}$| of the light (for VLT/NACO; Tuthill et al. 2010) and therefore significantly decreases the sensitivity and hence the contrast limit of the observations for relatively faint targets. However, kernel phase uses the light collected by the entire pupil and should perform better in the high Strehl regime and the bright limit (e.g. Pope et al. 2016; Sallum & Skemer 2019).

For sparse aperture masking, a mask is placed at the Lyot stop of an instrument in order to split the primary mirror into a discrete interferometric array of real sub-apertures (e.g. Readhead et al. 1988). In the Fourier transform of the detector image (hereafter referred to as Fourier plane), these sub-apertures map on to their autocorrelation (i.e. their spatial frequencies, Ireland 2016). The phase ϕ of each spatial frequency can be extracted and linearly combined in a way such that the resulting closure phase |$\boldsymbol{\theta} = \boldsymbol{\mathsf{K}}\cdot\boldsymbol{\phi}$| is independent of the pupil plane (or instrumental) phase φ to second order (i.e. terms of first and second order in φ are vanishing), where the matrix |$\boldsymbol{\mathsf{K}}$| encodes this special linear combination (e.g. Ireland 2016). For observations from the ground, the pupil plane phase φ is affected by noise from atmospheric seeing and non-common path aberrations that ultimately cause quasi-static speckles. Being more robust with respect to these systematic effects, sparse aperture masking achieves a superior angular resolution.

For observations from space, which do not suffer from atmospheric seeing, kernel phase has proven to be successful in resolving close companions at the diffraction limit (Martinache 2010; Pope, Martinache & Tuthill 2013). It is our goal to determine if, under good observing conditions, kernel phase also is a competitive alternative to sparse aperture masking from the ground.

2 METHODS

2.1 Data reduction

A basic direct imaging data reduction (such as dark, flat, background subtraction and bad pixel correction) is also essential for the kernel phase technique (e.g. Sallum & Eisner 2017). For this purpose, we developed our own data reduction pipeline¹ that can be fed the raw data with their associated raw calibrators from the VLT/NACO archive.² Our data reduction pipeline performs the following steps which are described in more detail in the following sections:

• Linearize the raw frames

• Compute master darks and their bad pixel maps

• Compute master flats and their bad pixel maps

• Flag saturated pixels

• Apply dark, flat, background and bad pixel corrections

• Perform a dither subtraction

• Reconstruct saturated pixels

• Select frames with sufficient Strehl ratio

2.1.1 Detector linearization correction

Like most photon counting devices, NACO’s infrared detector CONICA suffers from a non-linear response when the pixel counts approach the saturation threshold (16 400 counts for uncorrelated high well depth mode³ according to the NACO Quality Control and Data Processing,⁴ with a more conservative 16 000 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 analyse in Section 3 feature saturated PSFs which we want to correct for non-linearity before re-constructing 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 analyse) 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(t) 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(t) to all data points with less than 16 000 counts (saturation threshold, cf. the left-hand panel of Fig. 1). We linearize the detector using a continuously differentiable piecewise polynomial approach h to the correction curve f(g) with a linear function up to 8500 counts and a cubic polynomial between 8500 and 16 000 counts (cf. the right-hand panel of Fig. 1).

Figure 1.

Left-hand panel: median pixel count in dependence of the integration time t for uncorrelated high well depth mode from the detector monitoring (blue curve) and the linear (orange curve) and cubic (green curve) polynomials f(t) and g(t) that we fit to it. Right-hand panel: correction curve (blue) f(g) and the cubic polynomial h (orange curve) that we fit to it and use for linearizing all pixels with measured counts between 8500 and 16 000. In both panels, the solid red lines mark the end of the linear regime and the saturation threshold. Note that very low (i.e. negative) pixel counts occur due to the use of a narrow-band filter (⁠|$\Delta \lambda = 0.018\, {\mu \rm m}$|⁠) for the detector monitoring, whereas the L’ science frames are taken with a wide-band filter (⁠|$\Delta \lambda = 0.62\, {\mu \rm m}$|⁠).

Open in new tabDownload slide

2.1.2 Master darks and master flats

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:

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

• 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 that 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.

2.1.3 Saturated pixels

The data cubes that we analyse 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 that are all pixels with more than h(16 000) counts.

2.1.4 Dark, flat, background, and bad pixel correction

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-hand panel of Fig. 2, where residual systematic noise (mainly from the detector) is still clearly visible.

Figure 2.

Left-hand panel: median frame of a data cube of HIP 47425 after dark, flat and a simple background subtraction. The pixel counts are scaled by an arcsinh stretch so that both the PSF and the background are visible in the image. Right-hand panel: same median frame after performing the dither subtraction described in Section 2.1.5. This second step is essential to remove residual systematic noise from the detector that can be seen as grid-like structure in the left-hand panel. Note that the two panels have the same colour scale.

Open in new tabDownload slide

2.1.5 Dither subtraction

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-hand panel of Fig. 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 λ/D from the centre of the PSF in |$0.2\, \text{s}$| exposure, where 4.4 is the detector readout noise, λ is the observing wavelength (⁠|$3.8\, {\mu \rm m}$| for the L’ band) and D is the diameter of the primary mirror (⁠|$8.2\, \text{m}$| for the VLT).

2.1.6 Reconstruction of saturated pixels

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

• computing the Fourier transform of the corrected frame from the previous iteration;

• windowing this frame by the Fourier domain Z;

• computing the inverse Fourier transform of this frame; and

• 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 Fig. 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 that would otherwise suffer from high Fourier plane phase noise (cf. Fig. 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.

Figure 3.

Left-hand panel: mean over a horizontal and a vertical cross-section through the centre of the median frame shown in the right-hand panel of Fig. 2. Right-hand panel: same cross-section, but after reconstructing bad and saturated pixels as described in Section 2.1.6. The dashed black line marks the maximum of the cross-section in the left-hand panel in order to illustrate the reconstruction of the peak in the PSF core.

Open in new tabDownload slide

Figure 4.

Left-hand panels: Fourier plane phase of the median frame shown in the right-hand panel of Fig. 2 (top). The phase is flat in the centre, but the cut-off spatial frequency is smaller than the region of support permitted by the pupil geometry (magenta circle). Median Fourier plane phase at the spatial frequencies of our pupil model (bottom). Right-hand panels: same as in the left-hand panels, but after reconstructing bad and saturated pixels as described in Section 2.1.6. In both upper panels, the magenta line traces out the spatial frequencies of our pupil model (from left to right) in order to illustrate how the patterns observed in the lower panels are obtained.

Open in new tabDownload slide

2.1.7 Frame selection

As explained in the Introduction section, 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 per cent best frames. Then, we reject all frames with a peak count below 75 per cent of this value. Using this dynamic threshold is better than simply rejecting a fixed fraction of the frames (e.g. Law, Mackay & Baldwin 2006) because it can correctly account for a sudden drop in the Strehl ratio like shown in Fig. 5. Note that we consider the peak pixel count after performing the PSF reconstruction (cf. Section 2.1.6) here.

Figure 5.

Peak count for all 99 frames of a data cube of HIP 116258. The horizontal red line marks the rejection threshold computed according to Section 2.1.7. Around frame 70, the observing conditions suddenly become worse and a clear drop in the peak count can be observed.

Open in new tabDownload slide

2.2 Kernel phase extraction

2.2.1 VLT pupil model

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\, {\mu \rm 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-hand panel of Fig. 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-hand panel of Fig. 6.

Figure 6.

Left-hand panel: our VLT pupil model consisting of 140 virtual sub-apertures sampled on a square grid with a pupil plane spacing of |$0.6\, \text{m}$|⁠. The cyan circles show the size of the primary mirror and the central obscuration. Right-hand panel: Fourier plane coverage of the same pupil model. The magenta circle shows the region of support permitted by the pupil geometry in the left-hand panel. Since the Fourier transform is symmetric, we only use the phase measured in one half-plane. Note that only these Fourier plane positions within |$7.0\, \text{m}$| from the origin (i.e. these which do not suffer from low power, cf. Section 2.2.2) are shown.

Open in new tabDownload slide

2.2.2 XARA

The extraction of the Fourier plane phase and the computation of the kernel phase rely on a python package called XARA⁵ (eXtreme Angular Resolution Astronomy; Martinache 2010, 2013). 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 |$\boldsymbol{\mathsf{A}}$| and the associated left kernel operator |$\boldsymbol{\mathsf{K}}$| via a singular value decomposition of |$\boldsymbol{\mathsf{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 centred 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. the left-hand panels of Fig. 7). XARA offers several centring 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-centring is very valuable (cf. the right-hand panels of Fig. 7) due to an increased level of pupil plane phase noise from the atmosphere and the bright background (if compared to space-borne data).

Figure 7.

Left-hand panels: Fourier plane phase of the median frame of a data cube of TYC 6849 1795 1 (resolved and bright binary) after imperfect re-centring of the frames (top). The phase is flat in the centre, but there is an overall phase ramp from bottom to top caused by the resolved and bright companion. Median Fourier plane phase at the spatial frequencies of our pupil model (bottom). Right-hand panels: same as in the left-hand panels, but after proper re-centring of the frames. The residual Fourier plane phase is of considerably reduced amplitude and can be properly assembled to form meaningful kernel phases. The magenta circles and lines represent the same as in Fig. 4.

Open in new tabDownload slide

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 |$\boldsymbol{\mathsf{A}}$| are eliminated prior to the computation of |$\boldsymbol{\mathsf{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(r) = exp  − (r/r₀)⁴) with a radius r₀ = 25 pixels, effectively limiting our field of view to |${\sim }1000\, \text{mas}$|⁠. Note that Section 3.4 will further comment on the effect of this window and how it can affect contrast estimates for detections at large separations.

2.2.3 Kernel phase uncertainties

Figure 8.

Measured mean kernel phase |$\bar{\boldsymbol{\theta}}$| over all data cubes of HIP 47425 (typical calibrator) and TYC 6849 1795 1 (resolved and bright binary) as well as of its best-fitting binary model |$\boldsymbol{\theta}_\text{bin} = \boldsymbol{\mathsf{K}}\cdot\boldsymbol{\phi}_\text{bin}$| (cf. Section 2.4.1). Data and model agree very well, so that the green curve overlaps with the orange curve. Note that we normalize each kernel phase by the norm of its corresponding row of |$\boldsymbol{\mathsf{K}}$| and that the raw binary parameters reported here are not corrected for the windowing.

Open in new tabDownload slide

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 χ² of any model fitting, potentially to much more than 1.0 in the case of high signal-to-noise data with highly imperfect calibration.

2.3 Kernel phase calibration

Under perfect conditions the closure phase of a point-symmetric source, such as an unresolved star, is zero (e.g. Monnier 2007). The same holds for the kernel phase, which is a generalization of the closure phase (e.g. Ireland 2016). Practically, however, one is limited by systematic errors caused by third-order phase residuals (e.g. Ireland 2013) and even point-symmetric sources have non-zero kernel phase.

For this reason, calibration is of fundamental importance when analysing interferometric measurables (like closure or kernel phase). The systematic errors are expected to be quasi-static (e.g. Ireland 2013), i.e. slowly varying with time, so that the kernel phase of a well-known point source measured close in time to that of the science target can serve as a calibrator. The simplest calibration technique would be to subtract the kernel phase of a well-known point source from that of the science target. This technique was for example used successfully in Martinache (2010), but here we want to go beyond this approach.

2.4 Model fitting

From equations (2)–(4), it becomes clear that the measured kernel phase θ directly represents the kernel phase intrinsic to the observed object θ_obj. Hence, we can infer information about the spatial structure of the observed object by fitting models for |$\boldsymbol{\theta}_\text{obj} = \boldsymbol{\mathsf{K}}\cdot\boldsymbol{\phi}_\text{obj}$| to θ.

2.4.1 Binary model

2.4.2 Uncertainties from photon noise

Using the kernel phase covariance |$\boldsymbol{\mathsf{\Sigma}}_\theta$| estimated from photon noise according to Section 2.2.3, we compute the best-fitting contrast c_fit and its uncertainty |$\sigma _{c_\text{fit}}$| for the binary model |$\boldsymbol{\theta}_\text{bin} = \boldsymbol{\mathsf{K}}\cdot\boldsymbol{\phi}_\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}$|⁠.

The final parameters p_fit for the best-fitting binary model are then obtained from a least squares search that 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 ρ to |$50\, \text{mas} \le \rho \le 1000\, \text{mas}$|⁠.

2.4.3 Empirical uncertainties

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 analyse 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 SNR_ph (cf. Section 3.2.2). For each of the calibrators, we then compute two contrast curves:

• The azimuthal average c_rms of the root mean square (rms) best-fitting contrast |$c_\text{fit}$|

• The azimuthal average |$c_\text{RMS}^\text{sub}$| of the rms best-fitting contrast |$c_\text{fit}^\text{sub}$| after subtracting the best-fitting binary model θ_bin from the measured kernel phase θ

Figure 9.

Left-hand panel: mean of the azimuthal average c_rms of the rms best-fitting contrast |$c_\text{fit}$| of all non-detections of OB 2 (cf. column ‘Det’ of Table 1) before (solid blue curve) and after (dashed blue curve) subtracting the best-fitting binary model from the measured kernel phase. The dotted black curve represents the correction factor for the relative contrast of the residual speckle noise f_speckle. Right-hand panel: same as in the left-hand panel, but for HIP 50156 (close binary). The empirical 1σ detection limit that we use for our analysis (solid green curve) is obtained by multiplying the azimuthal average |$c_\text{RMS}^\text{sub}$| of the rms best-fitting contrast |$c_\text{fit}^\text{sub}$| after subtracting the best-fitting binary model from the measured kernel phase (dashed orange curve) with the correction factor f_speckle.

Open in new tabDownload slide

Furthermore, we obtain empirically motivated uncertainties on the best-fitting parameters p_fit by multiplying the uncertainties from photon noise |$\sigma _{p_\text{fit}}$| with the ratio f_err of the empirical contrast uncertainty σ_emp to the contrast uncertainty from photon noise σ_ph (at the position of the best-fitting binary model θ_bin).

The kernel phase analysis tools described in Sections 2.2–2.4 are available on GitHub.⁶ We put a strong focus on applicability to other instruments and an exchangeable kernel phase fits file format.

3 RESULTS AND DISCUSSION

3.1 Target list

We test our methods on an archival data set because the kernel phase technique is optimized for detecting companions at much smaller angular separations to their host star than conventional high-contrast imaging techniques (such as ADI and reference star differential imaging, i.e. RDI). Hence, the parameter space that we are looking at is still unexplored. We search the VLT/NACO archive for L’-band RDI surveys and decide to analyse programme 097.C-0972(A) (PI: J. Girard) due to a large number of observed targets and therefore potential calibrators. A target list together with our detections is reported in Table 1.

3.2 Detected companion candidates

Before we search the targets in Table 1 for close companion candidates, we perform a basic vetting procedure by visually inspecting the cleaned data for wide companion candidates (cf. Section 3.2.1). In the field of view, which is limited to |${\sim }1\, \text{arcsec}$| due to the windowing, we find six wide companion candidates (cf. the upper section of Table 2). Three of them are already known and we classify our detections as confirmed, whereas the other three have not been reported before and therefore are new detections. Note that we correct the contrast of the wide companion candidates for the windowing (cf. Section 3.4).

After detecting and subtracting off the signal induced by the wide companion candidates, we use the kernel phase technique in order to search for closer and fainter objects (cf. Section 3.2.2). We find two companion candidates with an empirical detection significance above the 5σ threshold, i.e. |$\text{SNR}_\text{emp}^\text{can} \gt 5$| (cf. lower section of Table 2). One of them is already known and we classify our detection as confirmed, whereas the other one has not been reported before and therefore is a new detection. For HIP 13008, we note that the empirical detection significance is 9.4σ when using only HIP 116384 as calibrator, but only 1.9σ when using HIP 12925 due to high residuals and a very large f_err correction. Therefore, HIP 12925 seems to be a bad calibrator and we do not report any best-fitting parameters for HIP 13008 due to a lack of credibility. Follow-up observations are required to confirm the true nature of this object. Also note that OBs 6–11 contain only one or two targets and are not analysed with the kernel phase technique because the diversity of kernel phase amongst calibrators is essential for our empirical detection method. As there are systematic differences between the individual nights in the measured kernel phase, for this paper we are only analysing OBs which contain at least two PSF calibrators (observed in the same night). Although this choice was made for simplicity and it might be possible to calibrate targets over longer time-scales, this adds significant additional complexity which is beyond the scope of this paper.

From the targets for which we detect neither a wide nor a close companion candidate, we compute a contrast curve (i.e. the detection limit as a function of the angular separation) for the kernel phase technique (cf. Section 3.3).

3.2.1 Wide companion candidates

The wide companion candidates reported in the upper section of Table 2 are all detected by visually inspecting the cleaned data. When we find a companion candidate, we use a grid search followed by a least squares search in order to find its best-fitting binary model θ_bin for the measured kernel phase θ. Then, we compute the empirical detection significance |$\text{SNR}_\text{emp}^\text{vis}$| for the best-fitting binary model θ_bin (cf. right-hand panels of Figs 10 and 11). This is achieved using a simplification of the empirical detection method (cf. Section 2.4.3). Since the wide companion candidates all have a sufficiently large angular separation (i.e. |${\gtrsim }200\, \text{mas}$|⁠) and are sufficiently bright (otherwise we could not detect them by eye), we can skip the use of any calibrators and compute the empirical detection significance |$\text{SNR}_\text{emp}^\text{vis}$| as the best-fitting contrast |$c_\text{fit}$| divided by the azimuthal average |$c_\text{RMS}^\text{sub}$| of the rms best-fitting contrast |$c_\text{fit}^\text{sub}$| after subtracting the best-fitting binary model θ_bin from the measured kernel phase θ. Note that we do not use any Karhunen–Loève calibration for this step either, i.e. θ′ = θ (cf. Section 2.3).

Figure 10.

Left-hand panels: median frame of a cleaned data cube of the targets for which we detect a wide companion candidate. The magenta star indicates the position of the host star and the magenta circle indicates the position of the companion candidate, obtained from a least squares fit of the binary model θ_bin to the measured kernel phase θ. Note that the colour scale is logarithmic, reaching from 1e+1.5 to 1e+3.5 pixel counts. Right-hand panels: map of the empirical detection significance |$\text{SNR}_\text{emp}^\text{vis}$| (cf. Section 3.2.1) for the same targets as in the left-hand panels. The number in the lower left corner of each panel reports the empirical detection significance at the position of the best-fitting binary model θ_bin (note that this is not necessarily the position with the highest detection significance) and the dashed cyan circle indicates the 99 per cent threshold of the super-Gaussian window (i.e. the brightness of objects outside this circle is decreased by more than 1 per cent by the windowing).

Open in new tabDownload slide

Figure 11.

Left-hand panels: median frame of a cleaned data cube of the targets for which we detect a wide companion candidate. The magenta star indicates the position of the host star and the magenta circle indicates the position of the companion candidate, obtained from a least squares fit of the binary model θ_bin to the measured kernel phase θ. Note that the colour scale is logarithmic, reaching from 1e+1.5 to 1e+3.5 pixel counts. Right-hand panels: map of the empirical detection significance |$\text{SNR}_\text{emp}^\text{vis}$| (cf. Section 3.2.1) for the same targets as in the left-hand panels. The number in the lower left corner of each panel reports the empirical detection significance at the position of the best-fitting binary model θbin (note that this is not necessarily the position with the highest detection significance) and the dashed cyan circle indicates the 99 per cent threshold of the super-Gaussian window (i.e. the brightness of objects outside this circle is decreased by more than 1 per cent by the windowing).

Open in new tabDownload slide

HIP 36985 B, TYC 7401 2446 1 B, TYC 5835 0469 1 B. These objects are new companion candidates which were not reported before. They have L’-band contrasts of 2.619 ± 0.005, 1.318 ± 0.004, and 2.399 ± 0.003 mag, respectively, and therefore are candidates for stellar mass companions.

TYC 6849 1795 1 B. This object was already detected in 2005 by Galicher et al. (2016) at an angular separation of |${\sim }220\, \text{mas}$|⁠, a position angle of |${\sim }201\, \text{deg}$| and a H-band contrast of ∼1.6 mag. We find an L’-band contrast of 1.450 ± 0.004 mag and an angular separation (⁠|$223.5\pm 0.4\, \text{mas}$|⁠) and a position angle (⁠|$203.29\pm 0.05\, \text{deg}$|⁠) that are in agreement with Galicher et al. (2016), i.e. we can confirm the bound nature of the object.

HIP 116231 B. This object was already detected in 2004 by Schöller et al. (2010) at an angular separation of |$641\pm 4\, \text{mas}$|⁠, a position angle of |$240.2\pm 0.6\, \text{deg}$|⁠, and a K-band contrast of 2.75 ± 0.01 mag. We find a L’-band contrast of 4.43 ± 0.02 mag, a slightly larger angular separation of |$874.6\pm 0.8\, \text{mas}$|⁠, and a slightly different position angle of |$254.70\pm 0.05\, \text{deg}$|⁠, but (allowing for orbital motion) we can confirm the bound nature of the object. Note that there is a huge disagreement in the contrast, but a brief look at the raw data from Schöller et al. (2010) shows a significant PSF halo and confirms our result of ∼4 mag.

HIP 116384 C. This object was first detected in 2002 by Martín (2003) who found HIP 116384 (GJ 900) to be a triple system with a |$510\pm 10\, \text{mas}$| (HIP 116384 B, ΔK = 1.61 ± 0.03 mag) and a |$760\pm 10\, \text{mas}$| (HIP 116384 C, ΔK = 2.38 ± 0.04 mag) component. Lafrenière et al. (2007b) resolved the system again in 2004 and 2005, finding HIP 116384 B at an angular separation of |$611\pm 2$| and |$673\pm 2\, \text{mas}$|⁠, respectively, and HIP 116384 C at an angular separation of |$733\pm 2$| and |$722\pm 2\, \text{mas}$|⁠, respectively. In the cleaned data, we only find HIP 116384 C at a slightly larger angular separation of |$842.90\pm 0.07\, \text{mas}$|⁠, but a position angle (⁠|$346.614\pm 0.004\, \text{deg}$|⁠) and a L’-band contrast (2.688 ± 0.001 mag) that are in agreement with Martín (2003) and Lafrenière et al. (2007b), so that we can confirm the bound nature of the object. Looking at the raw data, we also find HIP 116384 B (which is the brighter of the two companions), noticing that it has moved to an angular separation of |${\sim }1200\, \text{mas}$| being too far away in order to be visible in our cleaned data (due to the windowing).

3.2.2 Close companion candidates

The close companion candidates reported in the lower section of Table 2 are all detected only by the kernel phase technique. For each target in Table 1, we use a grid search followed by a least squares search in order to find the best-fitting binary model θ_bin for the measured kernel phase θ. Then, we compute the detection significance from photon noise SNR_ph (cf. Section 2.4.2) at the position of the best-fitting binary model θ_bin from the least squares search. For this step, we always use all other targets that were observed in the same OB as calibrators for the Karhunen–Loève calibration, fixing K_klip = 4.⁷ Knowing that the majority of VLT/NACO targets do not have any close companions, we then classify the ∼1/3 of the targets with the highest SNR_ph in each OB as candidate detections (cf. column ‘Can’ of Table 1) for the next step and the remaining targets as calibrators.

For the next step, we compute the empirical detection significance |$\text{SNR}_\text{emp}^\text{can}$| (cf. Section 2.4.3) for each of the candidate detections from the previous step. For this step, we always use all remaining targets that were classified as calibrators in the previous step for the Karhunen–Loève calibration, again fixing K_klip = 4. If the empirical detection significance is above the 5σ threshold, i.e. |$\text{SNR}_\text{emp}^\text{can} \gt 5$|⁠, we classify the candidate detection as real. If not, we add the candidate detection to the list of calibrators and redo the computation of the empirical detection significance (this time with one calibrator more than before). We repeat this process until all candidate detections are real. The detected close companion candidates are shown in Fig. 12 and are described in more detail in the following paragraphs. Please note that we report the correlation of the best-fitting parameters in Appendix  A and present model-data correlation plots in Appendix  B.

Figure 12.

Map of the empirical detection significance |$\text{SNR}_\text{emp}^\text{can}$| for the targets for which we detect a close companion candidate. The cyan star indicates the position of the host star and the solid cyan circle indicates the position of the companion candidate, obtained from a least squares fit of the binary model θ_bin to the measured kernel phase θ. The number in the lower left corner of each panel reports the empirical detection significance at the position of the best-fitting binary model θ_bin and the dashed cyan circle indicates the 99 per cent threshold of the super-Gaussian window (like in Fig. 10).

Open in new tabDownload slide

HIP 50156 B. This object was already detected in 2011 by Bowler et al. (2015) at an angular separation of |${\sim }90\, \text{mas}$| and a K-band contrast of ∼1.8 mag. Just nine month later, Brandt et al. (2014) cannot resolve this companion and report an upper limit of |${\sim }20\, \text{mas}$| for its angular separation. We find HIP 50156 B at an angular separation of |$77.3\pm 0.8\, \text{mas}$| and an L’-band contrast of ∼1.91 ± 0.03 mag, confirming the detection and notable orbital motion.

HIP 37918 B. This object is a new companion candidate that was not reported before. It has an L’-band contrast of ∼3.29 ± 0.05 mag, and therefore is a candidate for a stellar mass companion. Furthermore, HIP 37918 (M ≈ 0.98 M_⊙) is known to have a |${\sim }23.1\, \text{arcsec}$| companion of almost equal mass (HIP 37923, M ≈ 0.95 M_⊙, Desidera et al. 2006). Together with our companion candidate, this would make the system triple.

3.3 Detection limits

In Section 2.4.3, we present our empirical approach to find meaningful detection limits for the data analysed in this paper. Based on this approach, we compute the contrast limit of the kernel phase technique as a function of the angular separation as the azimuthal mean of the rms best-fitting contrast c_rms of all targets for which we do not detect any companions with the kernel phase technique (i.e. all non-detections, cf. column ‘Det’ of Table 1). Note that we already subtracted off the signal induced by the wide companion candidates. The mean, the best and the worst contrast limit are shown in the left-hand panel of Fig. 13.

Figure 13.

Left-hand panel: 5σ empirical contrast limit, i.e. rms contrast curve c_rms(ρ) multiplied by 5, for all non-detections (cf. column ‘Det’ of Table 1). Shown are the mean, the best and the worst contrast limit. Right-hand panel: value of the super-Gaussian windowing function depending on the angular separation. The brightness of companions outside of |${\sim }200\, \text{mas}$| is decreased significantly. We use this curve to recover the true contrast of the detected wide (visual) companions (cf. the upper section of Table 2). For reference, their position on this curve is indicated by the circles.

Open in new tabDownload slide

At the small angular separations that are inaccessible by classical high-contrast imaging techniques (i.e. within |${\sim }200\, \text{mas}$| in the L’ band), the kernel phase technique achieves contrast limits of ∼1e−2. This is not yet deep enough to detect companions in the planetary-mass regime, which would start between 1e−3 and 1e−4 for young (⁠|${\sim }10\, \text{Myr}$|⁠) gas giants (e.g. Bowler 2016). However, our closest detections prove that the resolution which is required to resolve Solar-system scales in the nearest star-forming regions can be achieved with the kernel phase technique. At larger angular separations, our best contrast limit is comparable with the limits achieved by RDI (e.g. Cantalloube et al. 2015). The large spread in the contrast limit comes from the fact that the amplitude of the background noise is nearly the same for all data cubes, whereas the peak value of the PSF varies heavily due to the PSF reconstruction (cf. Section 2.1.6).

3.4 Windowing correction

As mentioned in Section 2.2.2, we window all frames by a super-Gaussian (with an FWHM of |$1240\, \text{mas}$|⁠) in order to minimize edge effects when computing their Fourier transform. Due to this windowing, the brightness of companions at angular separations |${\gtrsim }215\, \text{mas}$| deviates by more than 1 per cent from the true value. In order to correct for this effect, we again assume that kernel phase is proportional to contrast in the high-contrast regime, so that we can obtain the true contrast of a companion by dividing its measured contrast (i.e. the best-fitting contrast from the binary model) by the value of the super-Gaussian windowing function. We are aware that this method has its limits, as each PSF has a spatial extent on the detector and assuming that the entire PSF is multiplied by the same value is an oversimplification of the problem. Nevertheless, this method agrees fairly well with the contrasts that we measure in the cleaned fits files and we use it to correct the contrast of all wide companion candidates (cf. the right-hand panel of Fig. 13). We add an additional contrast correction error in quadrature based on injection-recovery tests to companions wider than |$500\, \text{mas}$| to account for limitations in this technique.

4 CONCLUSIONS

We use the kernel phase technique in order to search for close companions at the diffraction limit in an archival VLT/NACO RDI L’-band data set. Therefore, we develop our own data reduction pipeline for VLT/NACO data, which performs a basic dark, flat, bad pixel and background (i.e. dither) subtraction, but also reconstructs saturated PSFs in order to reduce their Fourier plane noise. Furthermore, we select frames with sufficiently high Strehl ratio, which is essential for the kernel phase technique as it relies on a linearization of the Fourier plane phase. Then, we use XARA for extracting the kernel phase and improve its re-centring algorithm in the case of resolved and bright companions. Furthermore, we apply a principal component analysis based calibration to the data (i.e. Karhunen–Loève decomposition; Soummer et al. 2012) and develop a suite of analytic model fitting algorithms in order to search for point source companions with the kernel phase technique.⁸

For the archival data set that we analyse in Section 3, we find that our kernel phase covariance model (which only takes into account shot noise) is not sufficient and significantly underestimates the true errors. This is still the case after calibrating the data, because the diversity of calibrator PSFs is not sufficient. Hence, we develop an empirical method for estimating the relative contrast of the residual speckle noise and finding meaningful detection limits for the data. With this empirical approach, we detect six wide companion candidates by visually inspecting the cleaned data and two close (∼80–|$110\, \text{mas}$|⁠) companion candidates that are detected only by the kernel phase technique. All eight companion candidates lie in the stellar-mass regime and five of them were previously unknown.

In order to reach the planetary-mass regime, a better library of calibrator PSFs is required. Therefore, it is extremely important that the targets and their calibrators are observed as close in time as possible. This becomes very clear from the archival data set that we analyse, where there are in fact multiple calibrators observed in one night, but not close enough in time, so that the kernel phase calibration does not reduce the quasi-static errors satisfyingly. In order to make better use of our principal component analysis based calibration, we propose star-hopping sequences of ∼6 targets, and to revisit each target at least twice. Star-hopping is an observing strategy for which the instrument (and in particular the AO system) acquisition is only performed once at the beginning of each sequence. Then, one slews (‘hops’) from target to target without interrupting the AO system. Furthermore, we aim to examine more extensive Keck data sets where we are hopeful that the significant investment of telescope resources gives adequate calibrator diversity to characterize the systematic errors and possibly use Bayesian Monte Carlo techniques.

In this paper, we have shown that kernel phase is able to achieve a resolution below the classical diffraction limit of a telescope under good observing conditions (i.e. sufficiently high Strehl ratio). This is of particular interest for future space-based observatories, such as the James Webb Space Telescope, as it gives access to an exciting parameter space that could otherwise not be explored due to the limited mirror size (and therefore resolution). Space-based telescopes do not suffer from atmospheric turbulence, what makes the calibration much less challenging than for the ground-based VLT/NACO data (e.g. Martinache 2010). Nevertheless, with an optimized observing strategy, kernel phase is also a competitive high-contrast imaging technique from the ground.

The application of kernel phase is of course not limited to imaging telescopes. One concept that aims to push the kernel phase technique towards higher contrasts is the VIKiNG instrument (Martinache & Ireland 2018), which proposes kernel phase nulling interferometry with the VLTI. By combining kernel phase with a high-contrast booster (i.e. a nulling interferometer), it would allow for self-calibrating the observables and achieving a better robustness with respect to residual wavefront errors. This would in turn also be an option to reduce the demanding stability requirements on space-based nulling interferometers, such as the LIFE concept (Kammerer & Quanz 2018; Quanz et al. 2018), which aims to detect dozens of Earth-like exoplanets in the solar neighbourhood.

ACKNOWLEDGEMENTS

MJI was supported by the Australian Research Council Future Fellowship (FT130100235). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement CoG #683029). JHG gratefully acknowledges support from the Director’s Research Funds at the Space Telescope Science Institute. The manuscript was also substantially improved following helpful comments from an anonymous referee.

Footnotes

REFERENCES

APPENDIX A: PARAMETER CORRELATION

See Fig. A1.

Figure A1.

For the targets for which we detect a close companion candidate (i.e. HIP 50156, top and HIP 37918, bottom) we report the correlation of the best-fitting parameters using a corner plot from Foreman-Mackey (2016). Here, we use an MCMC technique (emcee, Foreman-Mackey et al. 2013) with six random walkers initialized at the best-fitting position and a temperature of |$f_\text{err}^2$| in order to find the best-fitting parameters including their correlated uncertainties by maximizing the log-likelihood ln L of the binary model.

Open in new tabDownload slide

APPENDIX B: CORRELATION PLOTS

See Fig. B1.

Figure B1.

Correlation of the measured kernel phase and the best-fitting binary model kernel phase for the targets for which we detect a close companion candidate in blue. The presented errorbars are computed based on photon noise (cf. Section 2.2.3) and scaled up by f_err according to our empirical uncertainties (cf. Section 2.4.3). The orange line indicates the identity which would represent perfect agreement between measured and model kernel phase. Similar to Fig. 8, we normalize each kernel phase by the norm of its corresponding row of |$\boldsymbol{\mathsf{P}}^{\prime}\cdot\boldsymbol{\mathsf{K}}$| since we are dealing with calibrated kernel phase here.

Open in new tabDownload slide