Kernel phase imaging with VLT/NACO

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

Jens Kammerer, Michael J. Ireland, Frantz Martinache and Julien H. Girard
Research School of Astronomy & Astrophysics, Australian National University, ACT 2611, Australia
Laboratoire Lagrange, Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Parc Valrose, Bât. H. FIZEAU, 06108 Nice, France
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
Accepted XXX. Received YYY; in original form ZZZ

Directly imaging exoplanets is challenging because quasi-static phase aberrations in the pupil plane (speckles) can mimic the signal of a companion at small angular separations. Kernel phase, which is a generalization of closure phase (known from sparse aperture masking), is independent of pupil plane phase noise to second order and allows for a robust calibration of full pupil, extreme adaptive optics observations. We applied kernel phase combined with a principal component based calibration process to a suitable but not optimal, high cadence, pupil stabilized L’ band () data set from the ESO archive. We detect eight low-mass companions, five of which were previously unknown, and two have angular separations of (i.e. ), demonstrating that kernel phase achieves a resolution below the classical diffraction limit of a telescope. While we reach a contrast limit of at such angular separations, we demonstrate that an optimized observing strategy with more diversity of PSF references (e.g. star-hopping sequences) would have led to a better calibration and even better performance. As such, kernel phase is a promising technique for achieving the best possible resolution with future space-based telescopes (e.g. JWST), which are limited by the mirror size rather than atmospheric turbulence, and with a dedicated calibration process also for extreme adaptive optics facilities from the ground.

planets and satellites: detection – planets and satellites: formation – techniques: high angular resolution – techniques: image processing – techniques: interferometric – binaries: close
pubyear: 2018pagerange: Kernel phase imaging with VLT/NACO: high-contrast detection of new candidate low-mass stellar companions at the diffraction limitLABEL:lastpage

1 Introduction

Direct imaging is vital for studying the outer regions of extrasolar systems which are inaccessible to transit observations and can only be revealed by decades-long, time consuming radial velocity surveys (e.g. fischer2014). It has proven particularly successful in probing our understanding of the formation of gas giant planets (e.g. d'angelo2010), being able to estimate their mass from their luminosity and age (e.g. spiegel2012) and resolve their orbit. Although the majority of detected companion candidates are arguably consistent with being emission or scattering from disk material (e.g. LkCa 15, kraus2012, HD 100546, quanz2013, HD 169142, biller2014), the recent example of PDS 70 (keppler2018) 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. bowler2018).

However, direct imaging operates at the resolution and sensitivity limit of the most powerful instruments today (e.g. pepe2014), placing demanding requirements on the observing and the post-processing techniques which are used to uncover faint companions at high contrasts (e.g. angular differential imaging, marois2006, point spread function subtraction, lafreniere2007a, principal component analysis, amara2012, soummer2012). Detecting exoplanets from the ground using these techniques has only been made possible by the recent development of extreme adaptive optics systems (e.g. milli2016) and is mainly limited by non-common path aberrations which are not sensed by the wavefront control system (e.g. sauvage2007). These aberrations manifest themselves as quasi-static speckles on the detector images which can mimic the signal of a companion and place a strong constraint on the achievable contrast at small angular separations (e.g. fitzgerald2006). Hence, directly imaging and studying the formation of gas giant planets on solar-system scales has been extremely challenging so far (e.g. bowler2016) because the nearest star forming regions lie away (e.g. loinard2007) where such orbital distances correspond to angular separations of only .

In this paper, we explore the capabilities of the kernel phase technique (martinache2010) 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 (tuthill2000). 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 in the near-infrared (i.e. the L’ band, cheetham2016). This gives access to objects on solar-system scales in the nearest star forming regions (i.e. projected separations of for a Jupiter analog in the Scorpius Centaurus OB association, preibisch2008) and has proven successful in directly imaging young exoplanets/disk features (e.g. kraus2012). The caveat of sparse aperture masking is that the mask blocks of the light (for VLT/NACO, tuthill2010) 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. pope2016; sallum2019).

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. readhead1988). In the Fourier transform of the detector image (hereafter referred to as Fourier plane), these sub-apertures map onto their auto-correlation (i.e. their spatial frequencies, ireland2016). The phase of each spatial frequency can be extracted and linearly combined in a way such that the resulting closure phase 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 encodes this special linear combination (e.g. ireland2016). For observations from the ground, the pupil plane phase is affected by noise from atmospheric seeing and non-common path aberrations which ultimately cause quasi-static speckles. Being more robust with respect to these systematic effects, sparse aperture masking achieves a superior angular resolution.

For full pupil kernel phase imaging, there is no mask and the entire primary mirror is discretized into an interferometric array of virtual sub-apertures. According to martinache2010, it is then convenient to define a transfer matrix which maps the baselines between each pair of virtual sub-apertures onto their corresponding spatial frequency. In the high Strehl regime, where the pupil plane phase can be linearized, we obtain the relationship


where is a diagonal matrix encoding the redundancy of the spatial frequencies (i.e. the baselines of the interferometric array) and is the phase intrinsic to the observed object. Multiplication with the left kernel of yields


hence the kernel of the measured Fourier plane phase directly represents the kernel of the phase intrinsic to the observed object , at least in the high Strehl regime (where is negligible). This is why frame selection based on the Strehl ratio is essential. Note that the kernel phase is a generalization of the closure phase to the case of redundant apertures.

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 (martinache2010; pope2013). 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. sallum2017). For this purpose, we developed our own data reduction pipeline111 which can be fed the raw data with their associated raw calibrators from the VLT/NACO archive222 Our data reduction pipeline performs the following steps which are described in more detail in the following sections:

  1. Linearize the raw frames.

  2. Compute master darks and their bad pixel maps.

  3. Compute master flats and their bad pixel maps.

  4. Flag saturated pixels.

  5. Apply dark, flat, background and bad pixel corrections.

  6. Perform a dither subtraction.

  7. Reconstruct saturated pixels.

  8. Select frames with sufficient Strehl ratio.

2.1.1 Detector linearization correction

Figure 1: Left panel: median pixel count in dependence of the integration time for uncorrelated high well depth mode from the detector monitoring (blue curve) and the linear (orange curve) and cubic (green curve) polynomials and which we fit to it. Right panel: correction curve (blue) and the cubic polynomial (orange curve) which we fit to it and use for linearizing all pixels with measured counts between 8500 and 16000. 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 () for the detector monitoring, whereas the L’ science frames are taken with a wide-band filter ().

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 mode333This is the standard imaging mode in the L’ band () and all data cubes which we analyze have been taken in this mode. according to the NACO Quality Control and Data Processing444, 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 , fit a linear polynomial to all data points with less than 8500 counts (end of the linear regime for uncorrelated high well depth mode) and a cubic polynomial 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 to the correction curve with a linear function up to 8500 counts and a cubic polynomial between 8500 and 16000 counts (cf. right panel of Figure 1).

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:

  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.

2.1.3 Saturated pixels

The data cubes which we analyze in Section 3 consist of 100 frames of 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 counts.

2.1.4 Dark, flat, background and bad pixel correction

Figure 2: Left 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 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 which can be seen as grid-like structure in the left panel. Note that the two panels have the same color scale.

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.

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 panel of Figure 2.

Our typical performance is a pixel count standard deviation of outside of from the center of the PSF in exposure, where is the detector readout noise, is the observing wavelength ( for the L’ band) and is the diameter of the primary mirror ( for the VLT).

2.1.6 Reconstruction of saturated pixels

Figure 3: Left panel: mean over a horizontal and a vertical cross-section through the center of the median frame shown in the right panel of Figure 2. Right 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 panel in order to illustrate the reconstruction of the peak in the PSF core.
Figure 4: Left panels: Fourier plane phase of the median frame shown in the right panel of Figure 2 (top). The phase is flat in the center, but the cutoff 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 panels: same as in the left 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.

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 () centered on the target. Then, we correct bad and saturated pixels for each frame separately by minimizing the Fourier plane power outside the region of support permitted by the pupil geometry. Let be the matrix which maps the bad and saturated pixel values onto the Fourier plane domain , then


where are the corrections to the bad and saturated pixel values (i.e. the corrected pixel values are ) and is remaining Fourier plane noise. We solve for using the Moore-Penrose pseudo-inverse of , 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 of 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 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 .

  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.

2.1.7 Frame selection

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.

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.

2.2 Kernel phase extraction

2.2.1 VLT pupil model

Figure 6: Left panel: our VLT pupil model consisting of 140 virtual sub-apertures sampled on a square grid with a pupil plane spacing of . The cyan circles show the size of the primary mirror and the central obscuration. Right 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 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 from the origin (i.e. these which do not suffer from low power, cf. Section 2.2.2) are shown.

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 , which is approximately half the Nyquist sampling of , where is the observing wavelength and is the image size (96 pixels). Our VLT pupil model is shown in the left panel of Figure 6 and based on an primary mirror with a 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.

2.2.2 Xara

The extraction of the Fourier plane phase and the computation of the kernel phase relies on a python package called XARA555 (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 and the associated left kernel operator via a singular value decomposition of .

Figure 7: Left panels: Fourier plane phase of the median frame of a data cube of TYC 6849 1795 1 (resolved and bright binary) after imperfect re-centering of the frames (top). The phase is flat in the center, 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 panels: same as in the left panels, but after proper re-centering 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 Figure 4.

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 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 and the corresponding rows of are eliminated prior to the computation of . 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 () with a radius pixels, effectively limiting our field of view to . 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.

2.2.3 Kernel phase uncertainties

Figure 8: Measured mean kernel phase over all data cubes of HIP 47425 (typical calibrator) and TYC 6849 1795 1 (resolved and bright binary) as well as of its best fit binary model (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 and that the raw binary parameters reported here are not corrected for the windowing.

For estimating the uncertainties, we compute the kernel phase covariance for each frame from its photon count variance in units of (photo-electrons), where is the detector gain ( for uncorrelated high well depth mode), is the super-Gaussian window, is the cleaned and re-centered frame and is its background (from the simple background subtraction, cf. Section 2.1.4). Therefore, we first need to find a linear operator which maps each frame in units of photo-electrons to its kernel phase . The linear discrete Fourier transform and the kernel of the pupil model are already linear operators, and the Fourier plane phase (of a complex number ) can be approximated as for small angles. Hence, we compute


Note that 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


Now, we have a kernel phase and a kernel phase covariance for each frame. In order to save computation time for the model fitting (cf. Section 2.4) we compute a weighted mean of the kernel phase for each data cube. Therefore, we first compute the average kernel phase covariance over all frames of a data cube via


and then the weighted mean of the kernel phase (cf. Figure 8) via


For the rest of this paper, we omit the bar for better readability, i.e.


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. monnier2007). The same holds for the kernel phase, which is a generalization of the closure phase (e.g. ireland2016). Practically however, one is limited by systematic errors caused by third order phase residuals (e.g. ireland2013) and even point-symmetric sources have non-zero kernel phase.

For this reason, calibration is of fundamental importance when analyzing interferometric measurables (like closure or kernel phase). The systematic errors are expected to be quasi-static (e.g. ireland2013), 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 martinache2010, but here we want to go beyond this approach.

We use principal component analysis in the framework of a Karhunen-Loève decomposition (soummer2012; pueyo2016) in order to subtract the statistically most significant phase residuals of the calibrator kernel phase from that of the science target. Note that the following technique is new, but very similar to the POISE observables in ireland2013. We start by computing the covariance matrix of the kernel phase of all calibrator data cubes via


Then, we compute an eigendecomposition of this matrix in order to obtain its sorted (in descending order) eigenvalues and eigenvectors . Finally, we compute the Karhunen-Loève transform of shape (number of kernel phases, number of calibrator data cubes) via


where is the p-th component of the k-th eigenvector of and is the n-th kernel phase of the p-th calibrator data cube.

From the Karhunen-Loève transform we obtain a projection matrix via


where is the identity matrix and is obtained from the first columns of . is an integer representing the order of the correction, i.e. how many eigencomponents of the calibrator kernel phase should be corrected for. The projection matrix is of shape (number of kernel phases, number of kernel phases), but it has zero eigenvalues by construction. In order to properly represent the dimensions we compute another eigendecomposition of and obtain a new projection matrix , whose columns are those eigenvectors of which correspond to non-zero eigenvalues. The projection matrix is of shape (number of “good” kernel phases, number of kernel phases), where “good” means statistically independent of systematic errors, and can be used to project the measured kernel phase and its covariance into a sub-space of dimension (number of “good” kernel phases), which is more robust with respect to quasi-static errors, via


For the rest of this paper, we omit the prime for better readability, i.e.


2.4 Model fitting

From Equations 24 it becomes clear that the measured kernel phase directly represents the kernel phase intrinsic to the observed object . Hence, we can infer information about the spatial structure of the observed object by fitting models for to .

2.4.1 Binary model

In order to search for companion candidates we use the binary model


where is the contrast ratio between secondary and primary, and are the coordinates of the sampled Fourier plane positions (i.e. the spatial frequencies of the pupil model), is the observing wavelength and


where is the angular separation between primary and secondary, is the position angle of the secondary with respect to the primary and is the detector position angle during the observation. Figure 8 shows the best fit binary model for the measured kernel phase of TYC 6849 1795 1 (resolved and bright binary).

2.4.2 Uncertainties from photon noise

Using the kernel phase covariance estimated from photon noise according to Section 2.2.3 we compute the best fit contrast and its uncertainty for the binary model on each position of a discrete square grid with spacing (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 .

In the high-contrast regime (where ), the phase is approximately proportional to the contrast of the binary model, so is its kernel phase (because is a linear operator). Hence, the of the binary model can be approximated as


where and are vertical stacks of the kernel phase and the reference binary model of each data cube  and is a block-diagonal matrix whose diagonal elements are the inverse kernel phase covariances of each data cube , i.e.


The reference binary model is the binary model evaluated for and normalized by a reference contrast , i.e.


Finally, we obtain the log-likelihood for the binary model as


The best fit contrast for the binary model is then obtained by maximizing for each grid position, i.e.


and its uncertainty is the square root of its variance, i.e.


Finally, the detection significance based on photon noise is computed for each grid position as


where we scale the uncertainty of the best fit contrast by the square root of the minimal reduced of the binary model of the entire grid (). Assuming that kernel phase is proportional to contrast, this is equivalent to scaling the kernel phase covariance so that the minimal reduced is 1.0. This step is necessary because kernel phase is still affected by third (or higher) order pupil plane phase noise (cf. Equations 24), so that the uncertainties from photon noise significantly underestimate the true errors. Note that there can be various sources of higher order phase noise (e.g. ireland2013), but studying those in detail is beyond the scope of this paper.

The final parameters for the best fit binary model are then obtained from a least squares search which maximizes the log-likelihood 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 .

The uncertainties of the best fit parameters follow from the likelihood function for Gaussian errors (which are applicable to high confidence detections)


where represents the three-dimensional parameter space of angular separation, position angle and contrast. Differentiating twice and neglecting terms containing second order derivatives of a single parameter yields


where and are the Jacobian and the Hessian matrix of the binary model . Hence, the covariance matrix of the model parameters can be obtained via


and the uncertainties of the model parameters for the best fit binary model are


We also compute the correlation of the best fit model parameters as


where denotes element-wise division.

2.4.3 Empirical uncertainties

Figure 9: Left panel: mean of the azimuthal average of the RMS best fit contrast of all non-detections of OB 2 (cf. column “Det” of Table 3.1) before (solid blue curve) and after (dashed blue curve) subtracting the best fit binary model from the measured kernel phase. The dotted black curve represents the correction factor for the relative contrast of the residual speckle noise . Right panel: same as in the left panel, but for HIP 50156 (close binary). The empirical detection limit which we use for our analysis (solid green curve) is obtained by multiplying the azimuthal average of the RMS best fit contrast after subtracting the best fit binary model from the measured kernel phase (dashed orange curve) with the correction factor .

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

  1. The azimuthal average of the root mean square (RMS) best fit contrast .

  2. The azimuthal average of the RMS best fit contrast after subtracting the best fit binary model from the measured kernel phase .

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 of the RMS best fit contrast after subtracting the best fit binary model (which might or might not be a real detection) from the measured kernel phase . Then, we multiply this RMS contrast curve with the mean of the relative speckle contrast of all calibrators, i.e.


where the bar denotes the mean, in order to obtain an empirical contrast uncertainty as a function of the angular separation for each candidate detection (cf. right panel of Figure 9). We classify a candidate detection as real if its empirical detection significance is above the threshold, i.e.


Furthermore, we obtain empirically motivated uncertainties on the best fit parameters by multiplying the uncertainties from photon noise with the ratio of the empirical contrast uncertainty to the contrast uncertainty from photon noise (at the position of the best fit binary model ).

The kernel phase analysis tools described in Sections 2.22.3 and 2.4 are available on GitHub666 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 analyze program 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 3.1.

OB Target SpT [pc] K [mag] [s] Vis Can Det
1 HIP 68994 F3/5V 71.7 6.715 395.8 N 46.6 Y 4.7 N
HIP 63734 F7/8V 54.1 6.436 389.2 N 44.3 N N
HIP 55052 K7V 23.7 6.808 389.2 N 45.1 N N
2 HIP 44722 K7V 14.5 5.757 395.6 N 22.1 N N
HD 108767 B K0V 26.7 6.235 310.2 N 21.4 N N
HIP 47425 M3V 9.6 6.056 388.8 N 32.4 Y 1.0 N
HIP 50156 M0.7V 23.4 6.261 395.2 N 292.7 Y 33.2 Y
HD 102982 G3V 53.2 6.605 316.6 N 23.9 N N
HIP 58029 G7V 42.2 6.78 395.8 N 32.9 Y 1.4 N
HIP 61804 G3V 59.2 6.869 395.8 N 27.1 N N
HD 110058 A0V 130.0 7.583 383.4 N 30.3 N N
HIP 72053 G3V 59.7 6.994 382.4 N 29.4 N N
3 HIP 58241 G4V 35.5 6.24 256.0 N 16.7 N N
TYC 8312 0298 1 K0II 804.5 6.475 162.0 N 18.0 N N
HIP 78747 F5V 41.1 4.859 280.8 N 22.8 Y 2.0 N
4 HIP 37918 K0IV-V 34.4 6.275 389.2 N 336.5 Y 20.3 Y
HIP 36985 M1.0V 14.1 5.934 334.4 Y 182.2 31.1 N N
TYC 7401 2446 1 K0V 42.2 6.778 117.4 Y 195.0 14.4 N N
5 TYC 6849 1795 1 K5V 27.6 6.911 305.4 Y 250.1 13.3 N N
HIP 92403 M3.5V 3.0 5.370 750.8 N 39.1 Y 2.5 N
HIP 94020 B K5V 29.1 6.999 657.0 N 23.9 N N
6 BDp19 3532 K0 240.2 5.842 1361.2 N N N
HIP 108085 B8IV-V 64.7 3.45 401.8 N N N
7 HIP 116231 B9.5III 53.4 4.611 285.2 Y 4.3 N N
HIP 116258 K2V 34.0 6.685 367.0 N N N
8 HIP 11484 B9III 60.4 4.392 279.6 N N N
9 HIP 3203 B K5V 26.5 6.834 181.6 N N N
10 TYC 5835 0469 1 G8V 60.9 6.997 465.0 Y 95.8 N N
TYC 9339 2158 1 K3V 30.3 6.712 461.2 N N N
11 HIP 7554 M0V 22.2 6.621 637.4 N N N
HIP 13754 K2V 38.6 6.883 503.8 N N N
12 HIP 116384 K7V 20.8 6.044 739.4 Y 99.7 26.2 N N
HIP 12925 F8 57.1 6.52 595.4 N 24.0 N N
HIP 13008 F5V 39.1 5.442 617.8 N 127.4 Y N/A N
13 HIP 14555 M1V 19.6 6.367 609.6 N 35.1 N N
HIP 20737 G9.5V 35.6 6.742 626.6 N 31.1 N N
HIP 22506 G9V 50.8 6.876 620.0 N 35.8 Y 4.3 N
HIP 23362 B9V 60.5 4.974 311.8 N 28.7 N N
Table 1: Target list grouped by OB for the VLT/NACO program 097.C-0972(A), PI J. Girard. For each target, we report the spectral type (SpT), the distance (), the apparent K band magnitude (K) and the total integration time after frame selection (). Whether we find any wide (visual) companion candidates, close (kernel phase) candidate detections and real detections is highlighted in columns “Vis”, “Can” and “Det”. We further report the empirical detection significance for the wide (visual) companion candidates (), the detection significance from photon noise for all targets () and the empirical detection significance for the close (kernel phase) candidate detections ().
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description