Angular differential kernel phases

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. Our aim is to experimentally demonstrate and validate the efficacy of an angular differential kernel-phase approach, a new method for self-calibrating interferometric observables that operates similarly to angular differential imaging, while retaining their statistical properties. We used linear algebra to construct new observables that evolve outside of the subspace spanned by static biases. On-sky 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. 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. 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.


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. Non-redundant masking, in particular, has been established as an observing mode in most of the high-resolution instruments available (Tuthill et al. 2010). Kernel-phase 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 Karhunen-Loè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 quasi-static 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 alt-azimuth 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

Analytical formulation
In this section, we describe a mathematical framework for generating self-calibrated ADKs building up from the traditional kernel-phase formalism.
Subsection §2.1 provides the context for kernel-phase observables and discusses the nature of the biases that arise. Construction of the new self-calibrating observables is handled in two steps: construction of the projection and whitening.

Classical calibration
The fundamental approximation of kernel-phase techniques is based on the finding that at high Strehl ratios (behind ExAO or in space), the phase error in the telescope's Pupil-plane ϕ will, to first-order, propagate linearly into the Fourier plane phase (Martinache 2010). This linear application is described by the matrix A, which leads, through the van Cittert-Zernike 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 model-fitting 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 third-order residual effects of wavefront error leaking through the kernel matrix.
Here, both are treated as Gaussian ) 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 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 Equation (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.

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 cutting-edge 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.

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 self-calibration. 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 kernel-phase 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 kernel-phase observables. The model observables κ 0,s and the residuals ε s are treated in the same way, leading to the new model-fitting 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 kernel-phase 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 model-fitting equation can then be formulated as follows: where any evolution of the quasi-static 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 non-zero 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.

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 Equation (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 block-diagonal: where Σ i is the covariance matrix of the i th kernel-phase 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 Equation (5), leading to the model-fitting 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 post-processing techniques. An alternative approach consists in deriving the same reasoning on a whitened version of Equation (12). We show in appendix C that this approach, although it is more complicated mathematically, leads to a projection matrix that is equivalent.
Having established the mathematical principles for a selfcalibration procedure, we now propose a practical implementation and on-sky 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 full-pupil kernel-phase observations, which requires an imaging instrument that provides a good wavefront correction.

The SCExAO instrument
The Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) (Jovanovic et al. 2015;Lozi et al. 2018) instrument is a planethunting, high-contrast imaging instrument on the Subaru telescope located on Maunakea, Hawaii. It is equipped with a 50x50 actuators deformable mirror controlled in closed-loop 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 on-sky 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 near-infrared camera without any coronagraphic mask. This camera is a C-RED2 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 Nyquist-Shannon sampling requirement in the H band.
The non-common-path aberrations (NCPA) on SCExAO are generally corrected for the CHARIS module (Groff et al. 2015), which is the prime near-infrared scientific instrument for highcontrast observations. As a consequence, the data was affected by large amounts of low-order static phase error dominated by astigmatism and trefoil for a total of 2 radians peak-to-valley in the Fourier phase, which is very challenging for kernel-phase 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 35cm. This model, represented in Fig. 1, is composed of 329 sub-apertures 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.

The data reduction pipeline
The fast frame rate of the camera resulted in millions of images. All were dark-subtracted and flat-fielded using dome flats. Images were co-added to increase the signal-to-noise 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 co-adds and frame selection rates were selected 1 https://github.com/fmartinache/xara  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 re-centered on the brightest speckle to integer pixel. A fourth-order 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 kernel-phase 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 kernel-phase 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 ma-trix. The mean parallactic angle is also evaluated for each bin to be used in the model fitting.
As described in §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 Equations (C.3), (C.4), and (C.7) and to be used in the ADK model-fitting Equation (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 matrixW s . This provided what will be subsequently referred to as the classical calibration analysis: Equations (19), (18), and (20) are used in parallel for comparison of the three approaches.

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 main-sequence 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 2019-03-21T13:40:00.
Here, we acquire 308µs exposures at 3,193 Hz for 40 minutes 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.
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 2019-06-25T14: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 2ms at 500 Hz, and co-added 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.

Basic analysis
Lucky imaging using a stack of the images around parallactic angle -0.4 degrees is given in Fig. 3 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 model-fitting and colinearity maps of the residual of this fit shown in Fig. 4 do not reveal any other strong binary signal.

Re-evaluation of the covariance
The amplitude of the residuals of this model-fitting 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 kernel-phase 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 model-fitting result.

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.
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 on-sky observations.

Performance comparison
The distribution of the residuals shown in Fig. 5 demonstrates a comparable distribution between ADK and classical calibration residuals. 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 under-evaluated covariance in the case of the classical calibration.

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 kernel-phase 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 minutes timescale. However, 90% of consistency is maintained over 20 minutes and around 80% over 30 minutes. Although the sample is too small to draw conclusions at longer timescales, the stability seems to deteriorate after that point.
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 . 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.
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 kernel-phase 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 3.2, for all the possible pairs of bins -as per the pairwise solution proposed with Equation (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 y l..m = W L,l..m · L l..m · κ s,l..m + κ s,th,l..m and x l..m = W L,l..m · L l..m · κ s,th,l..m .
Article number, page 7 of 11 A&A proofs: manuscript no. main   The result, presented in Fig. 7, shows that the signal increases steadily during the acquisition and only starts to taper-off after 30 minutes when the bias decorrelates to a significant extent.

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 quasi-static errors. These effects impact both the classical and ADK calibration procedures.
In the case of both observing nights, the static non-common path aberrations are expected to degrade the rejection performance of the kernel-phase 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 quasi-static and quickly varying errors. This has an substantial effect on the sensitivity of the kernel-phase analysis in general.
Because of the effect of these errors, our test does not provide a definitive reference for the performance of kernel-phase observations with SCExAO. However, the comparison between raw observables, classically calibrated observables, and ADK observables remains a topic of interest.

Non-stationarity 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 use-cases 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 point-like 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 non-stationarity 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.

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 non-redundant masking and in the case of kernel phases for full-pupil 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 kernel-phase 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 minutes 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 kernel-phase models, as well as the use of dedicated NCPA corrections for the internal infrared camera.
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 ) that provides both weighting and decorrelation in a single step. In the context of the concatenated observables, this translates to a block-diagonal W s matrix that can be computed block by block as follows: where Σ i is the covariance matrix of the i th kernel-phase observable vector. Working with this new observable, the steps are similar, although they are more complicated.
W s κ s = W s κ 0,s + W s U f · κ bias + W s ε sc , (C.2) 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 block-diagonal) because of the structure of both matrices. The result is a n k by n k matrix B. L 2 W s κ s = L 2 W s κ 0,s + L 2 W s ε s .

(C.5)
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 L 2 that has full column-rank. 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 non-zero eigenvalues of L 2 (here by using Λ which is a correspondingly cropped version of Λ): This time, the covariance of the error parameter L 2 W s ε s is: With this new projection matrix, Equation (C.5) becomes: L 2 W s κ s = L 2 W s κ 0,s + L 2 W s ε s . (C.9) Once again, we have successfully built a new projected observable in a smaller subspace that is robust to static errors (L 2 W s U f κ bias = 0) and has identity covariance (Cov L 2 W s ε s = I). 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 L 2 that transforms Equation (C.2) into Equation (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 drop-in replacement of usual whitened observables.
The matrices L and L 2 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 Moore-Penrose pseudo-inverse. 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 : χ 2 L 2 W s κ s = χ 2 MW L L κ s , (C.11) which develops: