Issue 
A&A
Volume 636, April 2020



Article Number  A21  
Number of page(s)  11  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201937121  
Published online  08 April 2020 
Angular differential kernel phases
^{1}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Lagrange, France
email: romain.laugier@oca.eu
^{2}
Research School of Astronomy & Astrophysics, Australian National University, Canberra, ACT 2611, Australia
^{3}
Subaru Telescope, National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 650 North Aohoku Place, Hilo, HI 96720, USA
^{4}
Steward Observatory, University of Arizona, Tucson, AZ 85721, USA
^{5}
College of Optical Sciences, University of Arizona, Tucson, AZ 85721, USA
^{6}
Astrobiology Center of NINS, 2211, Osawa, Mitaka, Tokyo 1818588, Japan
^{7}
European Southern Observatory, KarlSchwarzschildStr 2, 85748 Garching, Germany
Received:
15
November
2019
Accepted:
24
February
2020
Context. To reach its optimal performance, Fizeau interferometry requires that we work to resolve instrumental biases through calibration. One common technique used in high contrast imaging is angular differential imaging, which calibrates the point spread function and flux leakage using a rotation in the focal plane.
Aims. Our aim is to experimentally demonstrate and validate the efficacy of an angular differential kernelphase approach, a new method for selfcalibrating interferometric observables that operates similarly to angular differential imaging, while retaining their statistical properties.
Methods. We used linear algebra to construct new observables that evolve outside of the subspace spanned by static biases. Onsky observations of a binary star with the SCExAO instrument at the Subaru telescope were used to demonstrate the practicality of this technique. We used a classical approach on the same data to compare the effectiveness of this method.
Results. The proposed method shows smaller and more Gaussian residuals compared to classical calibration methods, while retaining compatibility with the statistical tools available. We also provide a measurement of the stability of the SCExAO instrument that is relevant to the application of the technique.
Conclusions. Angular differential kernel phases provide a reliable method for calibrating biased observables. Although the sensitivity at small separations is reduced for small field rotations, the calibration is effectively improved and the number of subjective choices is reduced.
Key words: instrumentation: high angular resolution / techniques: image processing / techniques: interferometric / methods: numerical
© Romain Laugier et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Since the advent of speckle interferometry (Labeyrie 1970), Fizeau interferometry techniques that use the aperture of a single telescope have proven to be a reliable way to obtain measurements at and beyond the classically defined resolution limit of telescopes. By working with the Fourier transform of images, they exploit observables that were originally developed for long baseline interferometry, such as closure phases (Jennison 1958; Baldwin et al. 1986), to provide observables that are robust for instrumental phase errors. Nonredundant masking, in particular, has been established as an observing mode in most of the highresolution instruments available (Tuthill et al. 2010). Kernelphase observables (Martinache 2010) rely on a generalization of the notion of closure phase for redundant apertures that is applicable in the high Strehl regime.
These interferometric techniques rely on simplifications, such as the monochromatic approximation, the short exposure approximation, or the absence of scintillation and instrumental amplitude errors. As described by Ireland (2013), deviations from these model conditions induce biases in the data that may sometimes prevent direct interpretations. Calibration observations, where calibrator targets are observed in the same conditions as the science target, are routinely used to remove these instrumental biases.
Techniques exist for improving the performance of the calibration, such as optimal calibration weighting (Ireland 2013), or KarhunenLoève projection (Kammerer et al. 2019). In both of those approaches, a number of different calibrators are used either to better interpolate the calibration signal or to project the target signal outside of the subspace containing the most common calibration signals. Both of these methods attempt to circumvent the problem of hidden variables that affect the quality of the calibration, but they tend to be limited by the small number of calibration sources available and the diversity they represent. Sampling this diversity requires the acquisition of a lot of data on calibration sources which is very inefficient. Among the hidden variables are spectral properties of the targets, inaccuracies of the model, and quasistatic instrumental biases. An ideal solution to the problem would be, therefore, to have a method to calibrate the science target with itself during the same observing session, thus preserving the consistency of the aforementioned parameters.
Systematic biases also affect coronagraphic observations, and similar approaches are used to characterize the flux leakage through coronagraphs. Angular differential imaging (ADI), proposed by Marois et al. (2006), uses the opportunity provided by the fact that the telescope uses an altazimuth mount and the relay train provides a fixed pupil relative to the detector along with a slowly rotating field of view. In consequence, the technique is capable of distinguishing instrumental errors that are fixed with the pupil from features of the target that are rotating with the field. The technique has been used extensively as part of different analysis processes such as TLOCI (Marois et al. 2014), KLIP (Soummer et al. 2012), or PYNPOINT (Amara & Quanz 2012).
The application of such a technique to closure phases or kernel phases has, to our knowledge, not yet been reported. We therefore introduce a mathematical expression for angular differential kernels (ADK), a selfcalibration technique that modifies the observables to completely remove the static bias while preserving their properties, both in terms of instrumental phase rejection and in terms of versatility, therefore building upon the current tools that have been developed for kernelphase analyses.
2. Analytical formulation
In this section, we describe a mathematical framework for generating selfcalibrated ADKs building up from the traditional kernelphase formalism.
Section 2.1 provides the context for kernelphase observables and discusses the nature of the biases that arise. Construction of the new selfcalibrating observables is handled in two steps: construction of the projection and whitening.
2.1. Classical calibration
The fundamental approximation of kernelphase techniques is based on the finding that at high Strehl ratios (behind ExAO or in space), the phase error in the telescope’s Pupilplane φ will, to firstorder, propagate linearly into the Fourier plane phase (Martinache 2010). This linear application is described by the matrix A, which leads, through the van CittertZernike theorem, to a comparison of the observations and the model as follows:
where Φ and Φ_{0} are the vectors of observed and theoretical Fourier phases and ε′ are deviations between model and observation. Their dimension is the number of considered baselines defined by a discrete representation of the pupil. Identifying the left null space of A (Martinache 2010) provides the matrix K:
therefore leading to the modelfitting equation:
The origins of different error sources that comprise K ⋅ ε′ have previously been described by Ireland (2013) and are also described further below. To simplify our notations, we introduce κ = K ⋅ Φ as the observables, κ_{0} = K ⋅ Φ_{0} as the model. For the purpose of this work, we decompose the errors into two types based on their temporal properties, namely, K ⋅ ε′ = κ_{bias} + ε, where κ_{bias} regroups the errors that remain stable over the timescale of the observations (the bias) and ε regroups all the errors varying during the observation. Therefore, the modelfitting equation is written as follows:
The varying error can encompass both correlated effects of sensor noise (readout noise and shot noise) or by thirdorder residual effects of wavefront error leaking through the kernel matrix. Here, both are treated as Gaussian (Ceau et al. 2019) and described by a covariance matrix Σ.
This covariance is usually estimated through either analytic propagation of the shot noise (Kammerer et al. 2019), imageplane bootstrapping of the shot noise (Ceau et al. 2019; Laugier et al. 2019), or empirically, when a sufficient number of short exposure images is available as in the case of the example presented here. The covariance Σ is used to construct a square whitening matrix
corresponding to a subspace where the considered errors are uncorrelated. This transforms Eq. (4) into:
which decorrelates uncertainties (Cov(Wε) = I), where W is the whitening matrix.
This approach, which lies at the heart of our statistical treatment, allows for optimal performance if the evaluation of the covariance is representative of the actual experimental errors.
A calibration signal κ_{calibrator} acquired on an unresolved source (κ_{0} = 0) is subtracted from the target signal as follows:
In this manner, the biases κ_{bias}, which are present in both observations, are canceled, whereas ε″ now contains the combined errors from the target and the calibrator, and W must therefore be calculated based on the covariance Σ that is written as:
with Σ_{target} and Σ_{calibrator} as the covariances of the target and the calibrator, respectively.
2.2. Discussion of the classical approach
It should be noted here that ε″ also contains residual calibration errors as described by Ireland (2013) that arise from a number of factors including some intrinsic to the choice of calibrators. Furthermore, we make the assumption that the calibrator is an unresolved point source, so that κ_{0} only represents the signal from the target. In practice, verifying this hypothesis is difficult, especially when the observing campaign has cuttingedge sensitivity and attempts to make new discoveries. Thus, any resolved signature in the calibrator propagates through the entire pipeline and either imprints artifacts or false detections on the final results, or it suppresses a real companion signal.
While a valid approach is to use a second calibrator to lift the ambiguity, requiring two calibrators for each target is inefficient as a smaller fraction of observation time is spent on the target.
Experience shows that observations of different calibration sources differ by amounts larger than can be explained by the statistics observed for each (Kraus et al. 2008; Martinache et al. 2009; Laugier et al. 2019); this is either because of the differences in spectral properties between the calibrators or because of slow variations in observing conditions. For this reason, a larger number of calibration sources are often used to try to improve the performance of the calibration step as described by Ireland (2013) and Kammerer et al. (2019).
Unfortunately, the distribution of the calibration errors are very difficult to evaluate both because of the small number of realizations (i.e., of calibrators) and because they combine variations of the spectral properties of the target, sky location leading to different air mass, adaptive optics performance, and telescope flexures. Our approach mitigates this effect while preserving the possibility to use the whitening transform to manage the fast varying errors.
2.3. Angular differential kernels
An alternative approach would be to completely remove the errors introduced by the intrinsic differences between target and calibrator by using only observations of the target itself in the logic of selfcalibration. As is the case in ADI, the signal of interest would then be identified through the diversity brought on by the field rotation.
For this purpose, we consider an arbitrary n_{f} number of frames acquired during an observing sequence, each taken at a known field position angle (parallactic angle for ground based telescopes, or roll angle for space telescopes, or other).
We consider a series of n_{f} kernelphase observables vectors κ_{i} extracted from the series of images, each corresponding to a different field rotation angle. We concatenate the whole dataset into a single long observable vector κ_{s} of length n_{k} × n_{f}, where n_{k} is the number of kernelphase observables. The model observables κ_{0, s} and the residuals ε_{s} are treated in the same way, leading to the new modelfitting equation, now written with the concatenated observables:
where κ_{bias, s} is the contribution of the static bias on all the observable vectors, and ε_{s} is the residual in this concatenated form. Since it is the same for all the frames, this contribution can be expressed as a function of an arbitrary single kernelphase signal of dimension n_{k}:
where κ_{bias} is the static calibrator signal of length n_{k}, and U_{f} is an unfolding matrix that maps the constant calibrator signal into a series of repeated signals (a concatenation of n_{f} identity matrices I, each of size n_{k} by n_{k})
The modelfitting equation can then be formulated as follows:
where any evolution of the quasistatic component has to be considered as part of ε_{s}. In a manner reminiscent of the ADI approach, we build a matrix L that cancels out the static bias in the observables. This matrix subtracts from each of the observable vectors κ_{i}, the mean of the observable vectors of each frame. We call this matrix L and express it as:
Details behind the construction of this matrix and its properties can be found in Appendix A. Since, by design, LU_{f} = 0, multiplying by L brings to zero any purely static contribution to the signal. Although the signal of interest κ_{0, s} may contain some static part that will be lost, the antisymmetric nature of phase signatures relative to the image plane guarantees that any feature with a nonzero phase signature will never be perfectly static with field rotation. As a consequence, some sensitivity loss is expected at small separations and small field rotations, the same as ADI.
2.4. Statistical whitening
The fact that the application described by L is not a surjection poses problems to subsequent statistical treatment as it implies a singular covariance for the resulting observables. In order to make direct use of statistical tools, we use a modified version L′ of the matrix L, which spans the same subspace but has been made into a surjection by removing its last n_{k} lines.
This turns Eq. (12) into:
This is better described in Appendix B; also, it does not affect the information gathered by the observables.
In the particular case n_{f} = 2, the L′ matrix can be written as two blocks:
which reveals that it operates a subtraction of the observables in the case of two frames. This algebraic approach provides flexibility of the number of frames that are acquired during a session, as well as the moment they are taken.
Assuming the statistical errors affecting the measurements are independent from one frame to the other, the covariance matrix of the concatenated observables is blockdiagonal:
where Σ_{i} is the covariance matrix of the ith kernelphase observable vector. Therefore, the covariance of the ADK observables built with the matrix L′ becomes:
In this case, the covariance matrix is invertible, and a whitening matrix W_{ADK} can be computed with Eq. (5), leading to the modelfitting equation to be used in the data reduction:
In line with the goal of this work, this equation can be used in all the typical applications of kernel phases, such as modelfitting and hypothesis testing. Observable quantities that are thus formed are insensitive to static biases and make up a generalization of simpler intuitive approaches. While similar in principle to what is carried out with ADI in the image plane, their covariance can be inverted for more rigorous postprocessing techniques.
An alternative approach consists in deriving the same reasoning on a whitened version of Eq. (12). We show in Appendix C that this approach, although it is more complicated mathematically, leads to a projection matrix that is equivalent.
3. Onsky validation
Having established the mathematical principles for a selfcalibration procedure, we now propose a practical implementation and onsky demonstration of its usability to show that an observing sequence of a rotating field can be exploited without the need to observe a calibration source. A classical calibration approach will be followed in parallel to provide a qualitative comparison.
Although the approach can also be used with other kinds of observables, we use fullpupil kernelphase observations, which requires an imaging instrument that provides a good wavefront correction.
3.1. The SCExAO instrument
The Subaru Coronagraphic Extreme Adaptive Optics (SCExAO; Jovanovic et al. 2015; Lozi et al. 2018) instrument is a planethunting, highcontrast imaging instrument on the Subaru telescope located on Maunakea, Hawaii. It is equipped with a 50 × 50 actuators deformable mirror controlled in closedloop with a pyramid wavefront sensor at several kHz. It operates as a second stage to the facility adaptive optics system AO188 which corrects 188 modes with a larger stroke using a curvature wavefront sensor.
Our observations took place within some of the dedicated onsky engineering time of the instrument, in parallel with other characterization and commissioning tasks, which was decisive in the choice of targets and observing time. We used the internal nearinfrared camera without any coronagraphic mask. This camera is a CRED2 camera providing low readout noise (30 electrons) at high frame rates (up to several kHz in cropped mode). It provides a plate scale of 16.2 mas per pixel, which satisfies the NyquistShannon sampling requirement in the H band.
The noncommonpath aberrations (NCPA) on SCExAO are generally corrected for the CHARIS module (Groff et al. 2015), which is the prime nearinfrared scientific instrument for highcontrast observations. As a consequence, the data was affected by large amounts of loworder static phase error dominated by astigmatism and trefoil for a total of 2 radians peaktovalley in the Fourier phase, which is very challenging for kernelphase observables.
The camera was also used to capture the pupil of the instrument using the internal calibration source. The pupil is then discretized using the tools provided in the xara^{1} package with a step of 35 cm. This model, represented in Fig. 1, is composed of 329 subapertures of equal transmission that normally generate 781 baselines to be sampled in the uv plane. However, in order to avoid the problem of π radians degeneracy of the phase that appears as a combination of the large NCPA and the phase noise, the longest baselines were discarded. The values of these parameters, and others described in this section, can be found in Table 1.
Fig. 1. Discrete representation of the pupil of the SCExAO instrument (left) and corresponding uv plane samples (right) used for the analysis of 3 Ser and the calibrator 31 Boo. Some more of the longest baselines were discarded for HD 211976, leading to a slightly smaller sampled uv plane. 
Value of relevant parameters for the extraction of the kernel phases.
3.2. The data reduction pipeline
The fast frame rate of the camera resulted in millions of images. All were darksubtracted and flatfielded using dome flats. Images were coadded to increase the signaltonoise ratio (S/N), then the frame selection was applied based on a criterion of empirical Strehl ratio, while the image timestamps were tracked. The number of coadds and frame selection rates were selected to prevent the π radians degeneracy mentioned before, yet to keep a large enough number of images for the statistical treatment to be valid.
Each frame was then recentered on the brightest speckle to integer pixel. A fourthorder exponential mask was numerically applied to the images both to reduce the contribution of the read noise in the unused parts of the image and to avoid the Fourier space aliasing (Laugier et al. 2019). Its radius r_{0} was chosen to correspond to a radius of 0.5λ/b, where b is the smallest baseline.
The complex visibilities were extracted for the pupil model using the xara package, which computes a discrete Fourier transform at the exact uv coordinates of the model, directly providing a vector of complex visibilities. A gradient descent algorithm was used to find a phase wedge that minimizes the square of the phase in the uv plane, thus providing a subpixel centering. The phase was then multiplied by the kernel matrix to obtain kernelphase observables.
The parallactic angle was determined for each frame and the kernel phases were binned together with a resolution of a few degrees of field rotation into n_{f} individual chunks. The value must be selected to be smaller than the rotation angle at which the kernelphase signal decorrelates, which depends on the extent of the region of interest. Larger separations therefore require smaller angular bins. This allows us to consider the signal as stationary for the statistical analysis, as well as to keep the calculation within reasonable scales. For the bins containing a number of realizations at least three times larger than n_{k}, the mean κ_{i} and its covariance Σ_{i} are empirically determined. The other bins are discarded so that they do not produce a biased covariance matrix. The mean parallactic angle is also evaluated for each bin to be used in the model fitting.
As described in Sect. 2.3 yielding κ_{s} and W_{s} to be used in what will be subsequently referred to as the raw analysis:
which neglects any bias, where Cov(W_{s}ε_{s}) = I, and where κ_{0, s} is a concatenated signal model evaluated for each of the mean field rotation parallactic angles.
The projection matrix is built based on Eqs. (C.3), (C.4), and (C.7) and to be used in the ADK modelfitting Eq. (3).
For the standard calibration strategy, the mean observables of the calibrator κ_{calibrator} were subtracted from each of the target observables and the covariance adjusted accordingly as the sum of the covariance matrices of the two signals, leading to the corresponding whitening matrix . This provided what will be subsequently referred to as the classical calibration analysis:
Equations (18)–(20) are used in parallel for comparison of the three approaches.
3.3. Targets
We observed targets of opportunity that were easily accessible over the engineering time of the instrument. 3 Ser (HDS 2143, HIP 74649) is a K0III star with a mainsequence stellar companion expected at a separation of ≈280 mas and a contrast in the infrared of ΔH ≈ 5 according to the orbit determined by Horch et al. (2015). The observations took place during the engineering time of the SCExAO instrument, starting 20190321T13:40:00. Here, we acquire 308 μs exposures at 3193 Hz for 40 min during its transit. The total rotation during this period was 35 degrees.
In order to compare the approach with a standard calibration procedure, a calibrator star of the same spectral type and the same apparent magnitude was selected in order to observe it with the same exposure times and similar AO performance. Proximity on sky is also necessary in order to replicate the airmass and instrumental flexures and for it to be observed just before or after the target. These constraints and the subjective decision they entail are a large part of the problems caused by classical calibration. In our case, we observed 31 Boo, a G7III star with the same frame rate and exposure times for a duration of two minutes between 13:09 and 13:40 UTC on the same night.
During this observing session, the wind was low, resulting in a high prominence of low wind effect (Sauvage et al. 2016) with devastating consequences for the PSF at high spatial frequencies and evolutions in the timescales of seconds. This was dealt with during the processing through frame selection by selecting only the 5% best frames. Figure 2 shows a few examples of images of the calibrator 31 Boo which were kept and rejected at this stage. In our case, for a companion located at the expected separation of 280 mas, the kernel phases were binned into steps of 1.8 degrees. Some of those bins that contained too few images had to be discarded, leaving 11 bins covering a total field rotation of 29 degrees.
Fig. 2. First two rows: images of the calibrator 31 Boo rejected at frame selection due to the lowwind effect. Last two rows: images kept for the analysis. Despite the selection, the appropriate images still display some variability on top of a static aberration revealed by the aspect of the first diffraction ring. 
For the second observing session, the aim was to study the stability of the instrument by observing a single star for a long period during its transit in a more typical set of observing conditions. We observed HD 211976 (HIP 110341) during an engineering time starting at 20190625T14:21:00 for a duration of 38.5 minutes, resulting in a total field rotation of 40 degrees. Because the star was fainter, the integration time was pushed to 2 ms at 500 Hz, and coadded four by four before a selection of the 50% best images based on empirical Strehl ratio. The field rotation binning was kept at 1.8 degrees, providing 23 bins with a total field rotation of 40 degrees. Parameters for both targets are available in Table 1.
3.4. Basic analysis
Lucky imaging using a stack of the images around parallactic angle −0.4 degrees is given in Fig. 3, showing the instrumental PSF and the region of interest. Careful examination of the upper right quadrant reveals the companion hidden in a field of speckles.
Fig. 3. Lucky imaging integration extracted from the parallactic angle bin with PA ≈ −0.4deg. The white square represents the region of interest of the analysis and the axis are labelled in mas. 
Colinearity maps are plotted as in Laugier et al. (2019) to highlight the possible location of the companion. The maximum value of this map is then used as a starting point for the extraction of the parameters of the binary by maximizing the likelihood with a Levenberg–Marquardt algorithm.
Figure 4 shows the colinearity maps where a companion is clearly visible, compatible with HDS 2143B. The companion’s position and contrast are evaluated by modelfitting and colinearity maps of the residual of this fit shown in Fig. 4 do not reveal any other strong binary signal.
Fig. 4. Left: colinearity maps of the region of interest around the target 3 Ser. The orientation corresponds to the orientation on the detector at a parallactic angle of 0 deg. Middle: correlograms corresponding to the final fit with adjusted covariance matrix. The corresponding χ^{2} is indicated for each plot. Right: colinearity map of the residual, after subtraction of the fitted binary signal. No additional companion is visible. The raw data has a significant amount of residual visible in the background and residual of the colinearity maps, as well as on the residual map. The level of bias is also visible with the poor fitting quality of the fit. ADK and classical calibration show comparable amount of residual background in the colinearity maps, but the χ^{2} value closer to 1 seem to indicate an even lower level of bias with ADK. 
3.5. Reevaluation of the covariance
The amplitude of the residuals of this modelfitting are much higher than expected. They do not show a strong binary signature, as revealed by colinearity maps of the residuals. We introduce the hypothesis that the extracted companion signal is the only significant kernelphase signal on 3 Ser. Based on this assumption, we evaluate the residual variance of the observables along the course of the observation. This variance is then added to the diagonal of the covariance matrix of each slice in order to better represent the errors. In the case of the classical calibration, the signal of the calibrator was appended to the series in order to account for the evolution of the signal in the relevant timescale.
Throughout the entire process, the different calibration strategies were applied independently from each other so that the residuals would not be affected by small differences in the modelfitting result.
4. Results
4.1. Revised analysis
Model fitting was performed anew with the corrected covariance estimation and the result is shown in Table 2. The results show consistent parameters for the three reduction pipelines. The estimated standard deviation for the position angle is significantly smaller in the case of ADK.
Extracted binary parameters compared for the different calibration approaches.
In the case of classical calibration, the higher value of χ^{2} is a sign of under estimation of the errors. This may indicate that a source of error was not taken into account. In this case, evolution of the bias between the calibrator and the target is a probable cause. This uncertainty is a typical source of problems in calibrating onsky observations.
4.2. Performance comparison
The distribution of the residuals shown in Fig. 5 demonstrates a comparable distribution between ADK and classical calibration residuals.
Fig. 5. Comparison of the distribution of residuals after the final model fitting. The residual from classical calibration appears skewed to the left, whereas the ADK residual looks very compatible with a centered Gaussian distribution. 
Based on the revised covariance estimation, sensitivity maps were drawn following the energy detector T_{E} test proposed by Ceau et al. (2019). They were determined for a detection probability of 0.95 and false alarm rate of 0.01 for the region of interest, which is plotted in Fig. 6. They show a sensitivity lower in ADK than in classical calibration, especially at small separations. Although sensitivity loss at small separation is expected in ADK, a significant part of the overall lower performance shown here can be attributed to the underevaluated covariance in the case of the classical calibration.
Fig. 6. Classical calibration and ADK sensitivity maps compared. The sensitivity is expressed in logcontrast on a common scale and the dashed line marks a separation of 1λ/D. The performance in classical calibration are expected to be optimistic, as hinted by a larger value of χ^{2} = 1.68. 
4.3. Stability analysis
Stability of the measurement is key to any calibration process. In classical calibration, the bias must remain static over the timescale necessary to switch from the calibration reference to the target, whereas in ADK, it must remain static during the necessary time to allow sufficient field rotation. We made use of the measurements on the single star HD 211976 to evaluate the evolution timescale of the bias signal. Since it is a single star, the only signal present is the bias signal. We considered each possible pair of kernelphase signals κ_{l} and κ_{m}, and built a common whitening matrix based on the mean of the two covariance matrices Σ_{l} and Σ_{m}:
This defines the preferred basis in which we can compare the two vectors. Then we computed the correlation c of the two vectors in this base:
where y_{l} = W_{l, m} ⋅ Σ_{l} and y_{m} = W_{l, m} ⋅ Σ_{m}. This correlation value is plotted in mean and standard deviation in Fig. 7. It shows that the signal loses the first 10% of consistency on the 5 min timescale. However, 90% of consistency is maintained over 20 min and around 80% over 30 min. Although the sample is too small to draw conclusions at longer timescales, the stability seems to deteriorate after that point.
Fig. 7. Representations of the evolution of the biases and signals with field rotation and time based on the observations of the single star HD 211976. Top: correlation between pairs of the kernelphase vectors evolving with their temporal separation. Bottom: series ADK signal for different target separations. The colored regions represent the standard deviation of the realizations in our sample, therefore, it is not relevant for the largest field rotation samples that only have one realization. The contrast for the binary signals are adjusted to provide signals of similar amplitudes. 
In the case of ADK, the bias must remain static in the timescale necessary for the field rotation to generate a significant signal in the differential kernel phases. As a consequence, the ADK sensitivity becomes the result of a race between the decorrelation of the signal and the decorrelation of the bias. This result is mostly influenced by the separation and contrast of the target and by the peak elevation of the target as seen from a given observatory (as it conditions the rate of field rotation). The amplitude of the bias signal is best evaluated in comparison to the signal of interest. Here, we add to the kernelphase value of each bin, the corresponding theoretical signal κ_{s, th} of a companion. Then we produce a reduction of the signal similar to what is described in Sect. 3.2, for all the possible pairs of bins – as per the pairwise solution proposed with Eq. (15) – and for all the possible series of consecutive bins (as per what is proposed in this work). We then project this signal onto the injected signal, in the same way as in the colinearity maps presented previously:
where
and
The result, presented in Fig. 7, shows that the signal increases steadily during the acquisition and only starts to taperoff after 30 min when the bias decorrelates to a significant extent.
5. Discussion
5.1. Adverse conditions
In the comparative test, the performance of both methods are affected by two adverse effects that are contingent upon the observing conditions. The low wind effect introduces strong correlated wavefront errors and the use of frame selection introduces selection bias in the residual wavefront errors, which can dominate the quasistatic errors. These effects impact both the classical and ADK calibration procedures.
In the case of both observing nights, the static noncommon path aberrations are expected to degrade the rejection performance of the kernelphase analysis as they impose a deviation from the linear regime. We expect this effect to increase the errors on the observables and, therefore, increase the covariance due to both quasistatic and quickly varying errors. This has an substantial effect on the sensitivity of the kernelphase analysis in general.
Because of the effect of these errors, our test does not provide a definitive reference for the performance of kernelphase observations with SCExAO. However, the comparison between raw observables, classically calibrated observables, and ADK observables remains a topic of interest.
5.2. Nonstationarity of the problem
In all of the test cases, despite the improvements brought by the ADK approach, our efforts at calibration were not sufficient to bring the residuals down to a level that could be explained by the covariances derived from each temporal bins. As in most of the published usecases of robust observables, ad hoc adjustment was employed for the final analysis. This could be explained either by the presence of astrophysical signal different from a pointlike companion (no response in the binary colinearity maps) or by the evolution of the bias at intermediate timescales.
In the general case (not only in ADK), the rotation of the field of view introduces an added difficulty to the evaluation and treatment of this error because of the nonstationarity of the problem. The observables we use evolve in time, which imposes some constraints to the analysis. The temporal binning of the parallactic angle in particular must be selected with care as a compromise. Although its step has to be small enough that the signal of interest can be considered static in its timescale, a longer timescale may improve the quality of the evaluation of the covariance. The signature of binary stars decorrelates faster at larger separations, therefore requiring a finer temporal binning scheme than for smaller separations. Since the density of the pupil model also has to increase in order to reach larger separations, this may render explorations farther than 10λ/D in long continuous sequences which are computationally demanding in ADK, where the projection matrix grows as (n_{k} ⋅ n_{f})^{2}.
As is the case with ADI, the sensitivity is slightly reduced at very small separations, especially in the case where small field rotation is available.
5.3. Conclusion
ADK is a new approach proposed for the calibration of robust observables that is aimed at removing most of the subjective choices that may often affect classical calibration techniques. It can be used both in the case of closure phases for nonredundant masking and in the case of kernel phases for fullpupil imaging. In all our test cases, it was the approach that provided the smallest and most Gaussian residuals.
Although our proposition for the matrix L′ may seem mathematically simplistic compared to some of the more advanced forms of ADI, it was designed for the purposes of a flawless integration with our kernelphase pipeline, especially with regard to the statistical whitening and hypothesis testing proposed by Ceau et al. (2019). Furthermore, the mathematical framework provided in Appendix A can be used even with more elaborate definitions of the bias signal by replacing the matrix U_{f} with the appropriate matrices.
The performance of the method currently appears to be limited by the opposition between the decorrelation of the signal of interest by field rotation and the decorrelation of the bias signal through the stability of the instrument. In the case of SCExAO, we showed that a steady accumulation of signal of interest was possible over more than 30 min and more than 33 degrees of field rotation, which makes the technique competitive even down to 1λ/D.
Further improvements are expected with the coming improvements of the kernelphase models, as well as the use of dedicated NCPA corrections for the internal infrared camera.
Acknowledgments
KERNEL has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG – 683029). Based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. The development of SCExAO was supported by the Japan Society for the Promotion of Science (GrantinAid for Research #23340051, #26220704, #23103002, #19H00703 & #19H00695), the Astrobiology Center of the National Institutes of Natural Sciences, Japan, the Mt Cuba Foundation and the director’s contingency fund at Subaru Telescope. The authors wish to recognise and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.
References
 Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [NASA ADS] [CrossRef] [Google Scholar]
 Baldwin, J. E., Haniff, C. A., Mackay, C. D., & Warner, P. J. 1986, Nature, 320, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Ceau, A., Mary, D., Greenbaum, A., et al. 2019, A&A, 120, 1 [Google Scholar]
 Groff, T. D., Kasdin, N. J., Limbach, M. A., et al. 2015, Techniques and Instrumentation for Detection of Exoplanets VII, 9605, 96051C [NASA ADS] [CrossRef] [Google Scholar]
 Horch, E. P., Van Belle, G. T., Davidson, J. W., et al. 2015, AJ, 150, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, M. J. 2013, MNRAS, 433, 1718 [NASA ADS] [CrossRef] [Google Scholar]
 Jennison, R. 1958, MNRAS, 118, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890 [NASA ADS] [CrossRef] [Google Scholar]
 Kammerer, J., Ireland, M. J., Martinache, F., & Girard, J. H. 2019, MNRAS, 486, 639 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kraus, A. L., Ireland, M. J., Martinache, F., & Lloyd, J. P. 2008, ApJ, 679, 762 [NASA ADS] [CrossRef] [Google Scholar]
 Labeyrie, A. O. 1970, A&A, 6, 85 [NASA ADS] [Google Scholar]
 Laugier, R., Martinache, F., Ceau, A., et al. 2019, A&A, 164, 1 [Google Scholar]
 Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, Proc., SPIE, 10703, 1070359 [Google Scholar]
 Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
 Marois, C., Correia, C., Galicher, R., et al. 2014, Proc. SPIE, 9148, 13 [Google Scholar]
 Martinache, F. 2010, ApJ, 724, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Martinache, F., RojasAyala, B., Ireland, M. J., Lloyd, J. P., & Tuthill, P. G. 2009, ApJ, 695, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Sauvage, J.F., Fusco, T., Lamb, M., et al. 2016, Adaptive Optics Systems V, 9909, 990916 [CrossRef] [Google Scholar]
 Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, 1 [Google Scholar]
 Tuthill, P. G., Lacour, S., Amico, P., et al. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, 7735, 77351O [CrossRef] [Google Scholar]
Appendix A: The matrix L via projection
There are several ways to obtain a matrix L that has the required properties. We will detail one and suggest a few approaches that eventually lead to mathematically equivalent results. We start with Eq. (12) describing the application of a static signal in the data:
where U_{f} is the unfolding matrix described in Sect. 2.3.
Considering that U_{f} is full column rank, the matrix is invertible, and a matrix P_{U} can be built that projects the signal into the space spanned by U_{f}.
And conversely, a projection into the orthogonal complement to this subspace can be computed as:
This matrix, therefore, projects the concatenated observables outside of the subspace reached by static bias as LU_{f} = 0. This can easily be demonstrated as:
where .
In this particular case,
which leads us to write
which is used in Eq. (13).
It is fair to note a resemblance to an approach that would evaluate an estimate of κ_{bias}, then subtract it from the observables of the series. Indeed, a rigorous least squares approach written in linear algebra, leads to the construction of the same matrix.
This reasoning would lead one to consider following this approach with a whitened observable W_{s} ⋅ κ_{s} in order to improve the estimate by benefiting from the weighting and decorrelation provided the whitening transform. This approach is provided in Appendix C for reference, but it does not improve the result, since the subspace is exactly the same, as it is perfectly defined by the matrix U_{f}.
Appendix B: Managing the correlations
Typical statistical analyses such as model fitting or detection tests rely on the computation of a likelihood functions which, for multivariate distributions, require the inversion of their covariance. In the typical kernelphase analyses, this is commonly done by the construction of the whitening matrix. In both cases, this requires that the covariance matrix of the final observables be invertible. In the case of the basic L matrix, the covariance of the observables is written as:
where, in the bestcase scenario, Cov(ε_{s}) = I, leads to:
because L is symmetric (L^{T} = L) and, as a projection, idempotent (L^{2} = L). Here L, as the matrix of a projection, is not invertible. In the general case, this comes from the fact that L is nonsurjective as it does not span the entire subspace of dimension n_{f} × n_{k}. As a consequence, the observables cannot be directly used through the rest of the usual kernelphase pipeline.
Since this is caused by the nonsurjective nature of L, an obvious solution to this problem is to reduce the dimensions of the output vectors to match the span of the matrix in order to make it into a surjection, and obtain an invertible covariance matrix. In the general case, an eigenvalue decomposition can be used to identify the subspace that is not reached by the transformation as will be defined by the eigenvectors that correspond to zero eigenvalues.
Here, we remark that the matrix is constituted of blocks corresponding to each frames, where each row of blocks can be expressed as a linear combination of all the other rows (minus their sum). We remove the last row and call this new matrix L′. This reduces the number of rows to the rank of the matrix which is (n_{f} − 1) × n_{k}, and therefore making it a surjection matrix.
Appendix C: Computing the projection on whitened observables
Since our construction of the matrix L is based around the expression of an estimator of κ_{bias}, one might expect to obtain better result using a weighted estimation. In our case the error ε_{s} is correlated. A good way to deal with this is to use a whitening transformation W (Ceau et al. 2019) that provides both weighting and decorrelation in a single step. In the context of the concatenated observables, this translates to a blockdiagonal W_{s} matrix that can be computed block by block as follows:
where Σ_{i} is the covariance matrix of the ith kernelphase observable vector. Working with this new observable, the steps are similar, although they are more complicated.
where U_{f} is still the unfolding matrix that maps the constant calibrator signal into a series of repeated signals; and κ_{bias} is the fixed calibrator signal.
As was done in Appendix A, we can compute the matrix of the projection into the orthogonal complement of U_{f}. This time, the method requires inverting (W_{s}U_{f})^{T}W_{s}U_{f} which is a bit more complicated, but can still be done very efficiently in the case where the errors between frames are independent (W_{s} is blockdiagonal) because of the structure of both matrices. The result is a n_{k} by n_{k} matrix B.
By defining the matrix as
we can simplify the equation into the following:
This approach uses the hypothesis on the calibration signal to project the observables in a subspace that is beyond its reach. Now since the error term L_{2}W_{s}ε_{s} is correlated, there is an additional step to perform as its covariance L_{2} is not invertible. However, since L_{2} is a projection, we can build a different projection matrix that has full columnrank. As a real symmetric matrix, L_{2} is diagonalized as per the form:
where Λ is a diagonal matrix containing 0s and 1s (because L_{2} is a projection). In the subspace defined by V^{T}, the operation corresponds to cropping the observable vector by n_{k}. We can therefore define a new projection matrix with full column rank with the rows of V^{T} corresponding to the nonzero eigenvalues of L_{2} (here by using Λ′ which is a correspondingly cropped version of Λ):
This time, the covariance of the error parameter is:
With this new projection matrix, Eq. (C.5) becomes:
Once again, we have successfully built a new projected observable in a smaller subspace that is robust to static errors () and has identity covariance (). This allows for the construction of a likelihood function for κ_{0, s} and the application of all the statistical tools defined in Ceau et al. (2019). In other terms, it means that based simply on the covariance Σ_{i} of the observables in the different frames, we can easily build the matrix that transforms Eq. (C.2) into Eq. (C.9), effectively removing the bias signal; and because of the precautions taken here, it can be used in any model fitting or statistical test as a dropin replacement of usual whitened observables.
The matrices L′ and are both designed to explore the same subspace that is purely defined as the orthogonal complement of the subspace of U_{f}. Although the observables are described in a different basis, the likelihood function that they provide is expected to be the same.
A transfer matrix M can be computed to pass from one set of observables to the other:
where the + sign indicates a MoorePenrose pseudoinverse. Although we do not provide a proof for it, we have found in our cases that M is a real unitary matrix which both confirms that they span the same subspace and that they provide the same sensitivity since they produce the same χ^{2}:
which develops:
and since M^{T}M = I, we obtain:
All Tables
All Figures
Fig. 1. Discrete representation of the pupil of the SCExAO instrument (left) and corresponding uv plane samples (right) used for the analysis of 3 Ser and the calibrator 31 Boo. Some more of the longest baselines were discarded for HD 211976, leading to a slightly smaller sampled uv plane. 

In the text 
Fig. 2. First two rows: images of the calibrator 31 Boo rejected at frame selection due to the lowwind effect. Last two rows: images kept for the analysis. Despite the selection, the appropriate images still display some variability on top of a static aberration revealed by the aspect of the first diffraction ring. 

In the text 
Fig. 3. Lucky imaging integration extracted from the parallactic angle bin with PA ≈ −0.4deg. The white square represents the region of interest of the analysis and the axis are labelled in mas. 

In the text 
Fig. 4. Left: colinearity maps of the region of interest around the target 3 Ser. The orientation corresponds to the orientation on the detector at a parallactic angle of 0 deg. Middle: correlograms corresponding to the final fit with adjusted covariance matrix. The corresponding χ^{2} is indicated for each plot. Right: colinearity map of the residual, after subtraction of the fitted binary signal. No additional companion is visible. The raw data has a significant amount of residual visible in the background and residual of the colinearity maps, as well as on the residual map. The level of bias is also visible with the poor fitting quality of the fit. ADK and classical calibration show comparable amount of residual background in the colinearity maps, but the χ^{2} value closer to 1 seem to indicate an even lower level of bias with ADK. 

In the text 
Fig. 5. Comparison of the distribution of residuals after the final model fitting. The residual from classical calibration appears skewed to the left, whereas the ADK residual looks very compatible with a centered Gaussian distribution. 

In the text 
Fig. 6. Classical calibration and ADK sensitivity maps compared. The sensitivity is expressed in logcontrast on a common scale and the dashed line marks a separation of 1λ/D. The performance in classical calibration are expected to be optimistic, as hinted by a larger value of χ^{2} = 1.68. 

In the text 
Fig. 7. Representations of the evolution of the biases and signals with field rotation and time based on the observations of the single star HD 211976. Top: correlation between pairs of the kernelphase vectors evolving with their temporal separation. Bottom: series ADK signal for different target separations. The colored regions represent the standard deviation of the realizations in our sample, therefore, it is not relevant for the largest field rotation samples that only have one realization. The contrast for the binary signals are adjusted to provide signals of similar amplitudes. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.