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 Julien H. Girard


firstpage–lastpage
An accurate method to determine the systematics due to the peculiar velocities of galaxies in measuring the Hubble constant from gravitational wave standard sirens
    Jian-hua He
    December 30, 2023
===================================================================================================================================================================





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 (3.8 µ m) data set from the ESO archive. We detect eight low-mass companions, five of which were previously unknown, and two have angular separations of ∼0.8–1.2 λ/D (i.e. ∼80–110 mas), demonstrating that kernel phase achieves a resolution below the classical diffraction limit of a telescope. While we reach a 5σ contrast limit of ∼1/100 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








§ 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 <cit.>. It has proven particularly successful in probing our understanding of the formation of gas giant planets <cit.>, being able to estimate their mass from their luminosity and age <cit.> 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, , HD 100546, , HD 169142, ), the recent example of PDS 70 <cit.> 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 <cit.>.

However, direct imaging operates at the resolution and sensitivity limit of the most powerful instruments today <cit.>, 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, , point spread function subtraction, , principal component analysis, , ). Detecting exoplanets from the ground using these techniques has only been made possible by the recent development of extreme adaptive optics systems <cit.> and is mainly limited by non-common path aberrations which are not sensed by the wavefront control system <cit.>. 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 <cit.>. Hence, directly imaging and studying the formation of gas giant planets on solar-system scales has been extremely challenging so far <cit.> because the nearest star forming regions lie ≳100 pc away <cit.> where such orbital distances correspond to angular separations of only ≲200 mas.

In this paper, we explore the capabilities of the kernel phase technique <cit.> 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 <cit.>. 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 ≲50 mas in the near-infrared <cit.>. This gives access to objects on solar-system scales in the nearest star forming regions (i.e. projected separations of ∼40 mas for a Jupiter analog in the Scorpius Centaurus OB association, ) and has proven successful in directly imaging young exoplanets/disk features <cit.>. The caveat of sparse aperture masking is that the mask blocks ≳85 % of the light <cit.> 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 <cit.>.

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 <cit.>. 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, ). The phase ϕ of each spatial frequency can be extracted and linearly combined in a way such that the resulting closure phase θ = K·ϕ 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 K encodes this special linear combination <cit.>. 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 <cit.>, it is then convenient to define a transfer matrix A 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

    ϕ = R^-1·A·φ+ϕ_obj+𝒪(φ^3),

where R is a diagonal matrix encoding the redundancy of the spatial frequencies (i.e. the baselines of the interferometric array) and ϕ_obj is the phase intrinsic to the observed object. Multiplication with the left kernel K of R^-1·A yields

    θ   = K·ϕ
       = K·R^-1·A_= 0·φ+K·ϕ_obj+𝒪(φ^3)
       = θ_obj+𝒪(φ^3)_≈ 0,

hence the kernel θ of the measured Fourier plane phase ϕ directly represents the kernel θ_obj of the phase intrinsic to the observed object ϕ_obj, at least in the high Strehl regime (where 𝒪(φ^3) 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 <cit.>. 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.



§ METHODS





 §.§ 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 <cit.>. For this purpose, we developed our own data reduction pipeline[<https://github.com/kammerje/PyConica>] which can be fed the raw data with their associated raw calibrators from the VLT/NACO archive[<http://archive.eso.org/wdb/wdb/eso/naco/form>]. 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.




  §.§.§ 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 (16400 counts for uncorrelated high well depth mode[This is the standard imaging mode in the L' band (3.8 µ m) and all data cubes which we analyze have been taken in this mode.] according to the NACO Quality Control and Data Processing[<https://www.eso.org/observing/dfo/quality/NACO/qc/detmon_qc1.html>], with a more conservative 16000 counts used in our analysis). As kernel phase is an interferometric technique for which the fringes are coded spatially on the detector, it is very important to characterize the pixel to pixel response. Moreover, many of the data cubes which we analyze in Section <ref> feature saturated point spread functions (PSFs) which we want to correct for non-linearity before reconstructing their core (cf. Section <ref>).

In order to compute the detector linearization correction we download all frames of type “FLAT, LAMP, DETCHECK” and uncorrelated high well depth mode from 2016 March 23 and 2016 September 24 (which are closest in time to the observation of the earliest and the latest data cube which we analyze) from the VLT/NACO archive. We sort them by integration time and compute the median pixel count over all frames for each individual integration time (masking out the broken stripes in the lower left quadrant of CONICA). Then, we plot the median pixel count in dependence of the integration time t, fit a linear polynomial f(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 16000 counts (saturation threshold, cf. left panel of Figure <ref>). 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 16000 counts (cf. right panel of Figure <ref>).



  §.§.§ 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 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.



  §.§.§ Saturated pixels



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



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



  §.§.§ 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 <ref>.

Our typical performance is a pixel count standard deviation of ∼36 = 4.4+(158/s·0.2 s) outside of 10 λ/D from the center of the PSF in 0.2 s exposure, where 4.4 is the detector readout noise, λ is the observing wavelength (3.8 µ m for the L' band) and D is the diameter of the primary mirror (8.2 m for the VLT).



  §.§.§ Reconstruction of saturated pixels







Our reconstruction of saturated pixels is based on an algorithm described in Section 2.5 of <cit.>. 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 (∼2.6 arcsec^2) centered on the target. Then, we correct bad and saturated pixels for each frame separately by minimizing the Fourier plane power |f_Z| outside the region of support Z permitted by the pupil geometry. Let B_Z be the matrix which maps the bad and saturated pixel values x onto the Fourier plane domain Z, then

    f_Z = B_Z·b+ϵ_Z,

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

    b = B_Z^+·f_Z = (B_Z^*·B_Z)^-1·B_Z^*·f_Z.

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

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

    
  * 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.
    
  * 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 <ref>. 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 <ref>). 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.



  §.§.§ Frame selection





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 <ref>) 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 <cit.> because it can correctly account for a sudden drop in the Strehl ratio like shown in Figure <ref>. Note that we consider the peak pixel count after performing the PSF reconstruction (cf Section <ref>) here.



 §.§ Kernel phase extraction





  §.§.§ 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 m, which is approximately half the Nyquist sampling of λ/α≈ 0.3 m, where λ = 3.8 µ m is the observing wavelength and α = 2.610 arcsec is the image size (96 pixels). Our VLT pupil model is shown in the left panel of Figure <ref> and based on an 8.2 m primary mirror with a 1.2 m central obscuration. Another advantage of kernel phase over sparse aperture masking is the dense Fourier plane coverage which is shown in the right panel of Figure <ref>.



  §.§.§ XARA



The extraction of the Fourier plane phase and the computation of the kernel phase relies on a python package called XARA[<https://github.com/fmartinache/xara>] <cit.>. 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 <ref>), and non-saturation (restored by the procedure described in Section <ref>). The discrete achromatic representation of the VLT aperture (i.e. our pupil model) is used by XARA to compute the phase transfer matrix A and the associated left kernel operator K via a singular value decomposition of A.



With the added knowledge of the detector pixel scale and the observing wavelength, the discrete model is scaled so that the Fourier plane phase at the expected (u,v) coordinates can be extracted by a discrete Fourier transform. For the small aberration hypothesis to remain valid, the data must be properly centered prior to the Fourier transform. Failure to do so will leave a residual Fourier plane phase ramp that can wrap and lead to meaningless kernel phases (cf. left panels of Figure <ref>). 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 <ref>) due to an increased level of pupil plane phase noise from the atmosphere and the bright background (if compared to space-borne data).

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

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



ORIGINAL TEXT FROM JENS FOLLOWS:

For extracting the Fourier plane phase and computing the kernel phase we use XARA: a python package for eXtreme Angular Resolution Astronomy[<https://github.com/fmartinache/xara>] <cit.>. XARA is applicable to data from various instruments and performs sub-pixel re-centering and a linear discrete Fourier transform in order to extract the Fourier plane phase ϕ from each individual frame. For computing the kernel phase θ = K·ϕ, the kernel K of the provided pupil model R^-1·A is obtained via a singular value decomposition. Careful investigations revealed that two things are particularly important for our extensive ground-based data set. Firstly, Fourier plane positions at the outer edge of the region of support permitted by the pupil geometry suffer from low power (because they are only supported by very few baselines, i.e. have small redundancy) and must be excluded from any further analysis (cf. right panel of Figure <ref>). Secondly, the originally developed re-centering algorithm must be improved in the case of resolved and bright companions. These two issues are discussed in detail in the next two paragraphs.

For the former issue, it should be noted that simply decreasing the size of the primary mirror when constructing the pupil model is not possible because this would change the kernel K of the pupil model R^-1·A and the computed kernel phase will not represent the true kernel of the pupil geometry anymore (because the primary mirror still collects light out to a diameter of 8.2 m). Instead, we construct the pupil model R^-1·A using the full pupil diameter. Then, we remove these rows of R^-1·A which correspond to Fourier plane positions of low power (i.e. these positions outside 7.0 m from the origin of the Fourier plane, cf. right panel of Figure <ref>). This can be done using the baseline filtering option which is already implemented into XARA.



For the latter issue, we develop a new re-centering algorithm for XARA. The original algorithm minimizes the Fourier plane phase slope at the origin of the Fourier plane using a weighted centroid, which is computationally fast but fails to minimize the overall Fourier plane phase in the case of resolved and bright companions (cf. left panels of Figure <ref>). Our new re-centering algorithm uses a least squares routine in order to minimize directly the Fourier plane phase ϕ at the spatial frequencies of our pupil model. This is achieved by performing sub-pixel shifts of a frame d via

    d_shifted = ℱ^-1(ℱ(d)·exp(-2π i·(dx·U+dy·V))),

where ℱ denotes the Fourier transform, dx and dy are the (sub-pixel) shifts in x and y direction and U and V are linear ramps from -1/2 to 1/2 in x and y direction (with the same shape as d). Note that we only minimize the Fourier plane phase within 6.0 m from the origin of the Fourier plane in order to avoid errors from spatial frequencies with low power. The performance of our new re-centering algorithm in the case of a resolved and bright companion is shown in the right panels of Figure <ref>. Note that this new re-centering algorithm will be implemented into XARA soon.

Finally, before extracting its Fourier plane phase ϕ and computing its kernel phase θ = K·ϕ, we window each frame by a super-Gaussian in order to minimize edge effects when Fourier transforming it.





  §.§.§ Kernel phase uncertainties





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

    B = K·Im(F)/|F· g·w·d|.

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

    Σ_θ = B·Σ_d·B^T.


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 <ref>) 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 d_i of a data cube via

    Σ̅_θ = (∑_iΣ_θ,i^-1)^-1,

and then the weighted mean θ̅ of the kernel phase (cf. Figure <ref>) via

    θ̅ = Σ̅_θ·∑_iΣ_θ, i^-1·θ_i.

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



 §.§ Kernel phase calibration



Under perfect conditions the closure phase of a point-symmetric source, such as an unresolved star, is zero <cit.>. The same holds for the kernel phase, which is a generalization of the closure phase <cit.>. Practically however, one is limited by systematic errors caused by third order phase residuals <cit.> 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 <cit.>, 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 <cit.>, but here we want to go beyond this approach.

We use principal component analysis in the framework of a Karhunen-Loève decomposition <cit.> 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 <cit.>. We start by computing the covariance matrix E_RR of the kernel phase θ_cal,i of all calibrator data cubes i via

    E_RR,(i,j) = θ_cal,i^T·θ_cal,j.

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

    Z_(n,k) = 1/√(w_k)∑_pv_k^p·θ_cal,p^n,

where v_k^p is the p-th component of the k-th eigenvector of E_RR and θ_cal,p^n is the n-th kernel phase of the p-th calibrator data cube.

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

    P = I-Z'·Z'^T,

where I is the identity matrix and Z' is obtained from the first K_klip columns of Z. K_klip 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 P is of shape (number of kernel phases, number of kernel phases), but it has K_klip zero eigenvalues by construction. In order to properly represent the dimensions we compute another eigendecomposition of P and obtain a new projection matrix P', whose columns are those eigenvectors of P which correspond to non-zero eigenvalues. The projection matrix P' 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

    θ' = P'·θ,
       Σ_θ' = P'·Σ_θ·P'^T.

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

    θ' →θ,
       Σ_θ' →Σ_θ.




 §.§ Model fitting



From Equations <ref>–<ref> 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 θ_obj = K·ϕ_obj to θ.



  §.§.§ Binary model



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

    r_bin·e^iϕ_bin = 1+c·exp(-2π i·(Δ_RA·u/λ+Δ_DEC·v/λ)),

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

    Δ_RA = -ρ·sin(ϑ-ϑ_0),
       Δ_DEC = ρ·cos(ϑ-ϑ_0),

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



  §.§.§ Uncertainties from photon noise



Using the kernel phase covariance Σ_θ estimated from photon noise according to Section <ref> we compute the best fit contrast c_fit and its uncertainty σ_c_fit for the binary model θ_bin = K·ϕ_bin on each position of a discrete 500×500 mas square grid with spacing 13.595 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×1000 mas.

In the high-contrast regime (where c ≪ 1), the phase ϕ_bin is approximately proportional to the contrast c of the binary model, so is its kernel phase θ_bin (because K is a linear operator). Hence, the χ^2 of the binary model χ_bin^2 can be approximated as

    χ_bin^2 = (Θ-c·Θ_bin,ref)^T·Σ_Θ^-1·(Θ-c·Θ_bin,ref),

where Θ and Θ_bin,ref are vertical stacks of the kernel phase θ_i and the reference binary model θ_bin,ref,i of each data cube i and Σ_Θ^-1 is a block-diagonal matrix whose diagonal elements are the inverse kernel phase covariances Σ_θ,i^-1 of each data cube i, i.e.

    Θ = [ θ_1; θ_2;   ⋮ ],     Σ_Θ^-1 = [ Σ_θ,1^-1        0        ⋯;        0 Σ_θ,2^-1        ⋯;        ⋮        ⋮        ⋱ ].

The reference binary model θ_bin,ref is the binary model θ_bin evaluated for and normalized by a reference contrast c_ref = 0.001, i.e.

    θ_bin,ref = θ_bin(c = c_ref)/c_ref.

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

    lnL = -1/2χ_bin^2.


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

    .∂/∂clnL|_c_fit   = 0,
    ⇒ c_fit   = Θ_bin,ref^T·Σ_Θ^-1·Θ/Θ_bin,ref^T·Σ_Θ^-1·Θ_bin,ref,

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

    σ_c_fit = 1/√(Θ_bin,ref^T·Σ_Θ^-1·Θ_bin,ref).

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

    SNR_ph = c_fit/σ_ph = c_fit/σ_c_fit·√(χ_r,bin,min^2),

where we scale the uncertainty of the best fit contrast σ_c_fit by the square root of the minimal reduced χ^2 of the binary model of the entire grid (χ_r,bin,min^2). Assuming that kernel phase is proportional to contrast, this is equivalent to scaling the kernel phase covariance Σ_θ so that the minimal reduced χ^2 is 1.0. This step is necessary because kernel phase is still affected by third (or higher) order pupil plane phase noise (cf. Equations <ref>–<ref>), so that the uncertainties from photon noise σ_c_fit significantly underestimate the true errors. Note that there can be various sources of higher order phase noise <cit.>, but studying those in detail is beyond the scope of this paper.

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

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

    lnL(p|x)    = -1/2χ_bin^2
       = -1/2(Θ-Θ_bin(p))^T·Σ_Θ^-1·(Θ-Θ_bin(p)),

where p 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

    H_(i,j)   = ∂^2/∂p_i∂p_jlnL(p|x)
       ≈∂Θ_bin(p)/∂p_i·Σ_Θ^-1·∂Θ_bin(p)/∂p_j
       = -(J·Σ_Θ^-1·J^T)_(i,j),

where J and H are the Jacobian and the Hessian matrix of the binary model Θ_bin. Hence, the covariance matrix of the model parameters Σ_p can be obtained via

    Σ_p = (-H)^-1 = (J·Σ_Θ^-1·J^T)^-1,

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

    σ_p_fit = √(diag(Σ_p_fit)).

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

    corr = Σ_p_fit/σ_p_fit^T·σ_p_fit,

where ·/· denotes element-wise division.



  §.§.§ 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 analyze in Section <ref> 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 <ref>). For each of the calibrators, we then compute two contrast curves:

    
  * The azimuthal average c_RMS of the root mean square (RMS) best fit contrast c_fit.
    
  * The azimuthal average c_RMS^sub of the RMS best fit contrast c_fit^sub after subtracting the best fit binary model θ_bin 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.

    f_speckle(ρ) = c_RMS(ρ)/c_RMS^sub(ρ),

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

For each of the candidate detections, we only compute the azimuthal average c_RMS^sub of the RMS best fit contrast c_fit^sub after subtracting the best fit binary model θ_bin (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 f_speckle of all calibrators, i.e.

    σ_emp(ρ) = f̅_speckle(ρ)·c_RMS^sub(ρ),

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

    SNR_emp = c_fit/σ_emp > 5.


Furthermore, we obtain empirically motivated uncertainties on the best fit parameters p_fit by multiplying the uncertainties from photon noise σ_p_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 fit binary model θ_bin).

The kernel phase analysis tools described in Sections <ref>, <ref> and <ref> are available on GitHub[<https://github.com/kammerje/PyKernel>]. We put a strong focus on applicability to other instruments and an exchangeable kernel phase fits file format.



§ RESULTS AND DISCUSSION





 §.§ 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 <ref>.





 §.§ Detected companion candidates





Before we search the targets in Table <ref> for close companion candidates, we perform a basic vetting procedure by visually inspecting the cleaned data for wide companion candidates (cf. Section <ref>). In the field of view, which is limited to ∼1 arcsec due to the windowing, we find six wide companion candidates (cf. upper section of Table <ref>). 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 <ref>).

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 <ref>). We find two companion candidates with an empirical detection significance above the 5σ threshold, i.e. SNR_emp^can > 5 (cf. lower section of Table <ref>). 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 fit 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 analyzed 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 analyzing 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 timescales, 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 <ref>).



  §.§.§ Wide companion candidates







The wide companion candidates reported in the upper section of Table <ref> 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 fit binary model θ_bin for the measured kernel phase θ. Then, we compute the empirical detection significance SNR_emp^vis for the best fit binary model θ_bin (cf. right panels of Figures <ref> and <ref>). This is achieved using a simplification of the empirical detection method (cf. Section <ref>). Since the wide companion candidates all have a sufficiently large angular separation (i.e. ≳200 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 SNR_emp^vis as the best fit contrast c_fit divided by the azimuthal average c_RMS^sub of the RMS best fit contrast c_fit^sub after subtracting the best fit 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 <ref>).

Before we search for closer and fainter objects, we subtract the signal induced by the wide companion candidates from the measured kernel phase, i.e.

    θ→θ-θ_bin,

so that the measured kernel phase of all targets is free of wide detections. The detected wide companion candidates are shown in the left panels of Figures <ref> and <ref> and are described in more detail in the following paragraphs.

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^mag, 1.318±0.004^mag 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 <cit.> at an angular separation of ∼220 mas, a position angle of ∼201 deg and a H band contrast of ∼1.6^mag. We find a L' band contrast of 1.450±0.004^mag and an angular separation (223.5±0.4 mas) and a position angle (203.29±0.05 deg) which are in agreement with <cit.>, i.e. we can confirm the bound nature of the object.


HIP 116231 B. This object was already detected in 2004 by <cit.> at an angular separation of 641±4 mas, a position angle of 240.2±0.6 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±0.8 mas and a slightly different position angle of 254.70±0.05 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 <cit.> shows a significant PSF halo and confirms our result of ∼4^mag.


HIP 116384 C. This object was first detected in 2002 by <cit.> who found HIP 116384 (GJ 900) to be a triple system with a 510±10 mas (HIP 116384 B, ΔK = 1.61±0.03^mag) and a 760±10 mas (HIP 116384 C, ΔK = 2.38±0.04^mag) component. <cit.> resolved the system again in 2004 and 2005, finding HIP 116384 B at an angular separation of 611±2 mas and 673±2 mas respectively, and HIP 116384 C at an angular separation of 733±2 mas and 722±2 mas respectively. In the cleaned data, we only find HIP 116384 C at a slightly larger angular separation of 842.90±0.07 mas, but a position angle (346.614±0.004 deg) and a L' band contrast (2.688±0.001^mag) which are in agreement with <cit.> and <cit.>, 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 ∼1200 mas being too far away in order to be visible in our cleaned data (due to the windowing).



  §.§.§ Close companion candidates





The close companion candidates reported in the lower section of Table <ref> are all detected only by the kernel phase technique. For each target in Table <ref>, we use a grid search followed by a least squares search in order to find the best fit binary model θ_bin for the measured kernel phase θ. Then, we compute the detection significance from photon noise SNR_ph (cf. Section <ref>) at the position of the best fit binary model θ_bin from the least squares search. For this step, we always use all other targets which were observed in the same OB as calibrators for the Karhunen-Loève calibration, fixing K_klip = 4[For simplicity, we fix K_klip = 4 for all targets and regardless of the number of calibrators. Various testing has shown that subtracting off the four statistically most significant eigencomponents of the kernel phase of the calibrators mostly yields the smallest amount of significant detections, i.e. calibrates the data best. A more rigorous investigation of this relationship is foreseen for a future publication.]. 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 <ref>) for the next step and the remaining targets as calibrators.

For the next step, we compute the empirical detection significance SNR_emp^can (cf. Section <ref>) for each of the candidate detections from the previous step. For this step, we always use all remaining targets which 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. SNR_emp^can > 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 Figure <ref> and are described in more detail in the following paragraphs. Please note that we report the correlation of the best fit parameters in Appendix <ref> and present model-data correlation plots in Appendix <ref>.

HIP 50156 B. This object was already detected in 2011 by <cit.> at an angular separation of ∼90 mas and a K band contrast of ∼1.8^mag. Just nine month later, <cit.> cannot resolve this companion and report an upper limit of ∼20 mas for its angular separation. We find HIP 50156 B at an angular separation of 77.3±0.8 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 which was not reported before. It has a 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 ∼23.1 arcsec companion of almost equal mass (HIP 37923, M ≈ 0.95 M_⊙, ). Together with our companion candidate, this would make the system triple.



 §.§ Detection limits





In Section <ref>, we present our empirical approach to find meaningful detection limits for the data analyzed 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 fit 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 <ref>). 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 panel of Figure <ref>.

At the small angular separations which are inaccessible by classical high-contrast imaging techniques (i.e. within ∼200 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 (∼10 Myr) gas giants <cit.>. 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 <cit.>. 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 <ref>).



 §.§ Windowing correction



As mentioned in Section <ref>, we window all frames by a super-Gaussian (with a FWHM of 1240 mas) in order to minimize edge effects when computing their Fourier transform. Due to this windowing, the brightness of companions at angular separations ≳215 mas deviates by more than 1% 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 fit 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 over-simplification of the problem. Nevertheless, this method agrees fairly well with the contrasts which we measure in the cleaned fits files and we use it to correct the contrast of all wide companion candidates (cf. right panel of Figure <ref>). We add an additional contrast correction error in quadrature based on injection-recovery tests to companions wider than 500 mas to account for limitations in this technique.



§ 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-centering algorithm in the case of resolved and bright companions. Furthermore, we apply a principal component analysis based calibration to the data <cit.> and develop a suite of analytic model fitting algorithms in order to search for point source companions with the kernel phase technique[<https://github.com/kammerje/PyKernel>].

For the archival data set which we analyze in Section <ref>, 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 mas) companion candidates which 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 which we analyze, 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 prinipal 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 JWST, as it gives access to an exciting parameter space which 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 <cit.>. 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 which aims to push the kernel phase technique towards higher contrasts is the VIKiNG instrument <cit.>, 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 <cit.>, which aims to detect dozens of Earth-like exoplanets in the solar neighborhood.



§ 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 program (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.







mnras



@urlcharsothermakeother $&#_%
@doi@urlcharsother ifnextchar [ @doi@
  @doi@[]
@doi@[#1]#2tempa#1tempaempty http://dx.doi.org/#2 doi:#2http://dx.doi.org/#2 #1
@eprint#1#2@eprint@#1:#2::nil
@eprint@arXiv#1http://arxiv.org/abs/#1 arXiv:#1
@eprint@dblp#1http://dblp.uni-trier.de/rec/bibtex/#1.xml
  dblp:#1
@eprint@#1:#2:#3:#4niltempa #1tempb #2tempc
  #3tempc empty tempc tempb tempb tempa tempb empty tempb arXivifundefined
  mn@eprint@tempbtempb:tempcmn@eprint@tempbtempc

[Amara & QuanzAmara &
  Quanz2012]amara2012
Amara A.,  Quanz S. P.,  2012, @doi []
  10.1111/j.1365-2966.2012.21918.x, http://adsabs.harvard.edu/abs/2012MNRAS.427..948A 427, 948

[Biller et al.,Biller
  et al.2014]biller2014
Biller B. A.,  et al., 2014, @doi [] 10.1088/2041-8205/792/1/L22,
  http://adsabs.harvard.edu/abs/2014ApJ...792L..22B 792, L22

[BowlerBowler2016]bowler2016
Bowler B. P.,  2016, @doi [] 10.1088/1538-3873/128/968/102001,
  http://adsabs.harvard.edu/abs/2016PASP..128j2001B 128, 102001

[Bowler & NielsenBowler &
  Nielsen2018]bowler2018
Bowler B. P.,  Nielsen E. L.,  2018, preprint, http://adsabs.harvard.edu/abs/2018arXiv180210132B  (@eprint arXiv
  1802.10132)

[Bowler, Liu, Shkolnik  &
  TamuraBowler et al.2015]bowler2015
Bowler B. P.,  Liu M. C.,  Shkolnik E. L.,   Tamura M.,  2015, @doi
  [] 10.1088/0067-0049/216/1/7, http://adsabs.harvard.edu/abs/2015ApJS..216....7B 216, 7

[Brandt et al.,Brandt
  et al.2014]brandt2014
Brandt T. D.,  et al., 2014, @doi [] 10.1088/0004-637X/786/1/1,
  http://adsabs.harvard.edu/abs/2014ApJ...786....1B 786, 1

[Cantalloube et al.,Cantalloube
  et al.2015]cantalloube2015
Cantalloube F.,  et al., 2015, @doi [] 10.1051/0004-6361/201425571,
  http://adsabs.harvard.edu/abs/2015A

[Cheetham, Girard, Lacour, Schworer,
  Haubois  & BeuzitCheetham et al.2016]cheetham2016
Cheetham A. C.,  Girard J.,  Lacour S.,  Schworer G.,  Haubois X.,
  Beuzit J.-L.,  2016, in Optical and Infrared Interferometry and Imaging V.
  p. 99072T, @doi10.1117/12.2231983

[D'Angelo, Durisen  &
  LissauerD'Angelo et al.2010]d'angelo2010
D'Angelo G.,  Durisen R. H.,   Lissauer J. J.,  2010, Giant Planet
  Formation.
pp 319–346

[Desidera, Gratton, Lucatello  &
  ClaudiDesidera et al.2006]desidera2006
Desidera S.,  Gratton R. G.,  Lucatello S.,   Claudi R. U.,  2006,
  @doi [] 10.1051/0004-6361:20064896, http://adsabs.harvard.edu/abs/2006A

[Fischer, Howard, Laughlin, Macintosh,
  Mahadevan, Sahlmann  & YeeFischer et al.2014]fischer2014
Fischer D. A.,  Howard A. W.,  Laughlin G. P.,  Macintosh B.,
  Mahadevan S.,  Sahlmann J.,   Yee J. C.,  2014, @doi [Protostars and
  Planets VI] 10.2458/azu_uapress_9780816531240-ch031, http://adsabs.harvard.edu/abs/2014prpl.conf..715F pp 715–737

[Fitzgerald & GrahamFitzgerald &
  Graham2006]fitzgerald2006
Fitzgerald M. P.,  Graham J. R.,  2006, @doi [] 10.1086/498339,
  http://adsabs.harvard.edu/abs/2006ApJ...637..541F 637, 541

[Foreman-MackeyForeman-Mackey2016]corner
Foreman-Mackey D.,  2016, @doi [The Journal of Open Source Software]
  10.21105/joss.00024, 24

[Foreman-Mackey, Hogg, Lang  &
  GoodmanForeman-Mackey et al.2013]foreman-mackey2013
Foreman-Mackey D.,  Hogg D. W.,  Lang D.,   Goodman J.,  2013, @doi
  [] 10.1086/670067, http://adsabs.harvard.edu/abs/2013PASP..125..306F 125, 306

[Galicher et al.,Galicher
  et al.2016]galicher2016
Galicher R.,  et al., 2016, VizieR Online Data Catalog, http://adsabs.harvard.edu/abs/2016yCat..35940063G 359

[IrelandIreland2013]ireland2013
Ireland M. J.,  2013, @doi [] 10.1093/mnras/stt859, http://adsabs.harvard.edu/abs/2013MNRAS.433.1718I 433, 1718

[IrelandIreland2016]ireland2016
Ireland M. J.,  2016, in Boffin H. M. J.,  Hussain G.,  Berger J.-P.,
  Schmidtobreick L.,  eds,  Astrophysics and Space Science Library Vol. 439,
  Astronomy at High Angular Resolution. p. 43,
  @doi10.1007/978-3-319-39739-9_3

[Kammerer & QuanzKammerer &
  Quanz2018]kammerer2018
Kammerer J.,  Quanz S. P.,  2018, @doi []
  10.1051/0004-6361/201731254, http://adsabs.harvard.edu/abs/2018A

[Keppler et al.,Keppler
  et al.2018]keppler2018
Keppler M.,  et al., 2018, @doi [] 10.1051/0004-6361/201832957,
  http://adsabs.harvard.edu/abs/2018A

[Kraus & IrelandKraus &
  Ireland2012]kraus2012
Kraus A. L.,  Ireland M. J.,  2012, @doi []
  10.1088/0004-637X/745/1/5, http://adsabs.harvard.edu/abs/2012ApJ...745....5K 745, 5

[Lafrenière, Marois, Doyon, Nadeau
  & ArtigauLafrenière et al.2007a]lafreniere2007a
Lafrenière D.,  Marois C.,  Doyon R.,  Nadeau D.,   Artigau
  É.,  2007a, @doi [] 10.1086/513180, http://adsabs.harvard.edu/abs/2007ApJ...660..770L 660, 770

[Lafrenière et al.,Lafrenière
  et al.2007b]lafreniere2007b
Lafrenière D.,  et al., 2007b, @doi [] 10.1086/522826, http://adsabs.harvard.edu/abs/2007ApJ...670.1367L 670, 1367

[Law, Mackay  & BaldwinLaw
  et al.2006]law2006
Law N. M.,  Mackay C. D.,   Baldwin J. E.,  2006, @doi []
  10.1051/0004-6361:20053695, http://adsabs.harvard.edu/abs/2006A

[Loinard, Torres, Mioduszewski,
  Rodríguez, González-Lópezlira, Lachaume, Vázquez
  & GonzálezLoinard et al.2007]loinard2007
Loinard L.,  Torres R. M.,  Mioduszewski A. J.,  Rodríguez
  L. F.,  González-Lópezlira R. A.,  Lachaume R.,  Vázquez
  V.,   González E.,  2007, @doi [] 10.1086/522493, http://adsabs.harvard.edu/abs/2007ApJ...671..546L 671, 546

[Marois, Lafrenière, Doyon,
  Macintosh  & NadeauMarois et al.2006]marois2006
Marois C.,  Lafrenière D.,  Doyon R.,  Macintosh B.,   Nadeau D.,
   2006, @doi [] 10.1086/500401, http://adsabs.harvard.edu/abs/2006ApJ...641..556M 641, 556

[MartínMartín2003]martin2003
Martín E. L.,  2003, @doi [] 10.1086/376742, http://adsabs.harvard.edu/abs/2003AJ....126..918M 126, 918

[MartinacheMartinache2010]martinache2010
Martinache F.,  2010, @doi [] 10.1088/0004-637X/724/1/464, http://adsabs.harvard.edu/abs/2010ApJ...724..464M 724, 464

[MartinacheMartinache2013]martinache2013
Martinache F.,  2013, @doi [] 10.1086/670670, http://adsabs.harvard.edu/abs/2013PASP..125..422M 125, 422

[Martinache & IrelandMartinache &
  Ireland2018]martinache2018
Martinache F.,  Ireland M. J.,  2018, @doi []
  10.1051/0004-6361/201832847, http://adsabs.harvard.edu/abs/2018A

[Milli, Mawet, Mouillet, Kasper  &
  GirardMilli et al.2016]milli2016
Milli J.,  Mawet D.,  Mouillet D.,  Kasper M.,   Girard J. H.,  2016,
  in Boffin H. M. J.,  Hussain G.,  Berger J.-P.,   Schmidtobreick L.,
  eds,  Astrophysics and Space Science Library Vol. 439, Astronomy at High
  Angular Resolution. p. 17 (@eprint arXiv 1701.00836),
  @doi10.1007/978-3-319-39739-9_2

[MonnierMonnier2007]monnier2007
Monnier J. D.,  2007, @doi [] 10.1016/j.newar.2007.06.006, http://adsabs.harvard.edu/abs/2007NewAR..51..604M 51, 604

[Pepe, Ehrenreich  & MeyerPepe
  et al.2014]pepe2014
Pepe F.,  Ehrenreich D.,   Meyer M. R.,  2014, @doi []
  10.1038/nature13784, http://adsabs.harvard.edu/abs/2014Natur.513..358P 513, 358

[Pope, Martinache  & TuthillPope
  et al.2013]pope2013
Pope B.,  Martinache F.,   Tuthill P.,  2013, @doi []
  10.1088/0004-637X/767/2/110, http://adsabs.harvard.edu/abs/2013ApJ...767..110P 767, 110

[Pope, Tuthill, Hinkley, Ireland,
  Greenbaum, Latyshev, Monnier  & MartinachePope
  et al.2016]pope2016
Pope B.,  Tuthill P.,  Hinkley S.,  Ireland M. J.,  Greenbaum A.,
  Latyshev A.,  Monnier J. D.,   Martinache F.,  2016, @doi []
  10.1093/mnras/stv2442, http://adsabs.harvard.edu/abs/2016MNRAS.455.1647P 455, 1647

[Preibisch & MamajekPreibisch &
  Mamajek2008]preibisch2008
Preibisch T.,  Mamajek E.,  2008, The Nearest OB Association:
  Scorpius-Centaurus (Sco OB2).
p. 235

[PueyoPueyo2016]pueyo2016
Pueyo L.,  2016, @doi [] 10.3847/0004-637X/824/2/117, http://adsabs.harvard.edu/abs/2016ApJ...824..117P 824, 117

[Quanz, Amara, Meyer, Kenworthy,
  Kasper  & GirardQuanz et al.2013]quanz2013
Quanz S. P.,  Amara A.,  Meyer M. R.,  Kenworthy M. A.,  Kasper M.,
  Girard J. H.,  2013, @doi [] 10.1088/2041-8205/766/1/L1, http://adsabs.harvard.edu/abs/2013ApJ...766L...1Q 766, L1

[Quanz, Kammerer, Defrère, Absil,
  Glauser  & KitzmannQuanz et al.2018]quanz2018
Quanz S. P.,  Kammerer J.,  Defrère D.,  Absil O.,  Glauser
  A. M.,   Kitzmann D.,  2018, in Optical and Infrared Interferometry and
  Imaging VI. p. 107011I (@eprint arXiv 1807.06088),
  @doi10.1117/12.2312051

[Readhead, Nakajima, Pearson,
  Neugebauer, Oke  & SargentReadhead et al.1988]readhead1988
Readhead A. C. S.,  Nakajima T. S.,  Pearson T. J.,  Neugebauer G.,
  Oke J. B.,   Sargent W. L. W.,  1988, @doi [] 10.1086/114724,
  http://adsabs.harvard.edu/abs/1988AJ.....95.1278R 95, 1278

[Sallum & EisnerSallum &
  Eisner2017]sallum2017
Sallum S.,  Eisner J.,  2017, @doi [] 10.3847/1538-4365/aa90bb,
  http://adsabs.harvard.edu/abs/2017ApJS..233....9S 233, 9

[Sallum & SkemerSallum &
  Skemer2019]sallum2019
Sallum S.,  Skemer A.,  2019, arXiv e-prints, http://adsabs.harvard.edu/abs/2019arXiv190101266S 

[Sauvage, Fusco, Rousset  &
  PetitSauvage et al.2007]sauvage2007
Sauvage J.-F.,  Fusco T.,  Rousset G.,   Petit C.,  2007, @doi
  [Journal of the Optical Society of America A] 10.1364/JOSAA.24.002334,
  http://adsabs.harvard.edu/abs/2007JOSAA..24.2334S 24, 2334

[Schöller, Correia, Hubrig  &
  AgeorgesSchöller et al.2010]schoeller2010
Schöller M.,  Correia S.,  Hubrig S.,   Ageorges N.,  2010, @doi
  [] 10.1051/0004-6361/201014246, http://adsabs.harvard.edu/abs/2010A

[Soummer, Pueyo  & LarkinSoummer
  et al.2012]soummer2012
Soummer R.,  Pueyo L.,   Larkin J.,  2012, @doi []
  10.1088/2041-8205/755/2/L28, http://adsabs.harvard.edu/abs/2012ApJ...755L..28S 755, L28

[Spiegel & BurrowsSpiegel &
  Burrows2012]spiegel2012
Spiegel D. S.,  Burrows A.,  2012, @doi []
  10.1088/0004-637X/745/2/174, http://adsabs.harvard.edu/abs/2012ApJ...745..174S 745, 174

[Tuthill, Monnier, Danchi, Wishnow  &
  HaniffTuthill et al.2000]tuthill2000
Tuthill P. G.,  Monnier J. D.,  Danchi W. C.,  Wishnow E. H.,
  Haniff C. A.,  2000, @doi [] 10.1086/316550, http://adsabs.harvard.edu/abs/2000PASP..112..555T 112, 555

[Tuthill et al.,Tuthill
  et al.2010]tuthill2010
Tuthill P.,  et al., 2010, in Ground-based and Airborne Instrumentation for
  Astronomy III. p. 77351O (@eprint arXiv 1006.2586),
  @doi10.1117/12.856806

[Wenger et al.,Wenger
  et al.2000]wenger2000
Wenger M.,  et al., 2000, @doi [] 10.1051/aas:2000332, http://adsabs.harvard.edu/abs/2000A




















§ PARAMETER CORRELATION







§ CORRELATION PLOTS






















