Recovering saturated images for high dynamic Kernel-Phase analysis Application to the determination of dynamical masses for the system Gl 494AB

Kernel-phase observables extracted from mid- to high-Strehl images are proving to be a powerful tool to probe within a few angular resolution elements of point sources. The attainable contrast is however limited by the dynamic range of the imaging sensors. The Fourier interpretation of images with pixels exposed beyond the saturation has so far been avoided. We show that in cases where the image is dominated by the light of a point source, we can use an interpolation to reconstruct the otherwise lost pixels with an accuracy sufficient to enable the extraction of kernel-phases from the patched image. We demonstrate the usability of our method by applying it to archive images of the Gl 494AB system, taken with the Hubble Space Telescope in 1997. Using this new data point along with other resolved observations and radial velocity measurements, we produce improved constraints on the orbital parameters of the system, and consequently the masses of its components.


Introduction
The quest for the direct imaging of exoplanets has accelerated in the recent years with the commissioning of extreme adaptive optics systems with coronagraphic capabilities, such as SPHERE (Beuzit et al. 2008), GPI (Macintosh et al. 2014) and SCExAO (Jovanovic et al. 2015). These instruments were designed to achieve high contrast detections (10 4 ) at small angular separations (down to ∼ 2λ/D where λ/D defines the resolution element, λ is the wavelength, and D is the telescope diameter). The planets discovered thus far by these instruments have however been detected at larger separations: Macintosh et al. (2015) report a detection at ∼ 10λ/D (449 mas) and Chauvin et al. (2017) at ∼ 20λ/D (830 mas). The performance of these coronagraphic devices at small angular separations is currently limited by the quality of the wavefront correction of the adaptive optics system, to which they are extremely sensitive (Guyon et al. 2006).
In the small angular separation regime, the wavefront quality requirements become so stringent that light leakage induced by instrumental phase dominates over the shot noise in the coronagraphic images. In practice, the use of interferometric techniques designed to be robust to the wavefront errors, like closure phases (Jennison 1958), or kernel-phases (KP) (Martinache 2010) are a competitive tool for discovering faint companions around nearby stars. The KERNEL project aims at developing observation and image processing means that provide observables which are intrinsically robust to small wavefront errors and can therefore be used to recover information in the image.
Kernel-phase observables (Martinache 2010) provide a reliable probe for the detection of asymmetric features at the smallest separations (down to 0.5λ/D) with moderate contrast (80-200:1). The attainable contrast is constrained by the dynamic range of the camera used. The acquisition of images for their kernel-phase analysis requires a compromise: the image must be sufficiently exposed to obtain satisfactory signal to noise ratio from the faint feature, while avoiding saturation on its brightest parts.
By destroying linearity of the intensity distribution, saturation (or clipping) breaks some of the prerequisites of the Fourier transform. As a result the shift theorem and objectimage convolution relationship become unusable, preventing the interpretation of the phase and the construction of kernels. Although the limits were pushed both by Pope et al. (2016) who had to linearize the sensor's response (soft saturation) and Martinache et al. (2016) where the uv plane sampling was truncated to avoid the problematic signals and function in a degraded mode, no kernel-phase analysis has yet been published based on hard saturated images.
The work presented in this paper attempts to circumvent this limitation and extract kernel-phases from images featuring a saturated core after using a saturation recovery algorithm. As a demonstration, we apply this method to the reduction of images taken with the Near Infrared Camera and Multi-Object Spectrograph (NICMOS) on the Hubble Space Telescope (HST) of Gl 494 (DT Vir, LHS 2665, HIP 63510, BD +13 2618) acquired in 1997 which leads to a new detection and measurement of a known companion. The novel visual astrometry data point, combined with available resolved observations at other epochs together with Article number, page 1 of 9 arXiv:1901.02824v2 [astro-ph.IM] 10 Jan 2019 A&A proofs: manuscript no. Saturation_recovery radial velocity observations, lead to improved constraints on the orbital parameters of the system and therefore on the mass of its components.

Interpolation of saturated images and kernel analysis
Typical algorithms for bad pixels correction assume that problematic pixels are few and isolated, allowing recovery by interpolating the neighboring pixels in the image, or by minimizing the power associated to the highest spatial frequencies in the Fourier plane (Ireland 2013). In the case of saturation, the pixels to recover are clustered. Our approach tackles this problem in the image plane, and uses a synthetic Point Spread Function (PSF) as a reference to interpolate the value of the saturated pixels, after which Fourier-phase can then be reliably extracted from these enhanced images.

Saturation recovery algorithm
The restoration of saturated pixels assumes that the local signal is dominated by the flux of a bright point source. This makes it possible to fit a parametric theoretical (or empirical) PSF to the non-saturated parts of the image and to replace the saturated pixels by interpolation of this adjusted PSF. The best solution of position (x, y) and amplitude z minimizes the following χ 2 variable: where i = Σ − 1 2 ·s img and p(x, y) = Σ − 1 2 ·s P SF (x, y) are, respectively, error-normalized (whitened) vectors of the image signal s img , and the (x, y) shifted PSF signal s P SF (x, y). The synthetic PSF is obtained with the TinyTim (Krist et al. 2011) simulator, and both vectors exclude saturated pixels, which ensures the fitting is not biased by their arbitrary value. A centroid algorithm provides a starting point within a few pixels of the final result, scale at which the problem is convex. Minimizing χ 2 (x, y, z) should give the best performance assuming a good estimator of the covariance matrix Σ is used.
This covariance has contributions from the image sensor, and from the PSF model. For the image, the pixels are independent (and the contribution diagonal), and a good estimator can be built from the data in the reduced image file, with Eq. 6 for the terms on the diagonal. For the PSF however, the deviations in the image plane are strongly correlated, and their estimations would rely on strong hypotheses on the spectral distribution of the aberration modes of the wavefront errors. In the case of a space telescope, thanks to the good stability of the PSF, one would expect to be able to neglect this error term; effectively using the sensor noise as the only source of deviation between the image and the PSF. However it was found that a uniform error estimator made the interpolation less sensitive to the evolution of the small instrumental phase over long timescales, and therefore produced smaller biases in the observables. In practice, since the images are large, an exponential windowing function (sometimes called super-Gaussian) was used to exclude the pixels that have low SNR: where r is the distance from the pixel to the approximated center of the star, and r 0 is a radius parameter chosen depending on the SNR of the image. From the error estimation point of view this can be seen as using the inverse of the values of the mask as the corresponding σ 2 ii terms of a diagonal matrix Σ. A constant value of r 0 = 40 pixels gave satisfactory results with the NICMOS images and was used for the whole dataset.
An analytical expression ofẑ, the optimal value of z can be obtained from Eq. 1: leading, through the substitution in Eq. 1 to the definition of a new criterion, function of x and y only: Minimizing ε(x, y) is therefore equivalent to minimizing χ 2 (x, y, z) but is computationally more efficient.

Simulation of realistic
NICMOS images The simulation of images is necessary for two purposes: the evaluation of the fidelity of the algorithm in Sect. 2.4, and the bootstrapping of the kernel-phase covariance matrix in Sect. 3.3.
Images of single and binary stars were simulated based on a PSF of HST obtained with the TinyTim software. Binary stars were constructed by adding a shifted and contrasted copy of the original PSF. The noise behavior of the non-destructive reads was also emulated in order to reproduce the typical behavior of the sensor in MULTIAC-CUM/STEP128 mode. First the number of correct readouts and exposure time is estimated for each pixel through an iterative process using the full well capacity and the flux. Then the readout noise map is estimated as: where σ ro,single is the readout noise value for a single read and N s is the number of successful samples. This is useful when simulating images in order to build sensible approximations for the metadata maps included in the cal.fits files, but isn't necessary when bootstrapping the errors for existing data as in Sect. 3.3. For the application of photon (shot) noise, we follow the directions (and notations) provided in the NICMOS data handbook (Thatte & Dahlen 2009) and consider a total error map of: where SCI is the flux in count per second, G is the ADC inverse gain (in electron/ADU), ERR is the noise map, and T IM E is the map of the total exposure time, all available in the FITS file. Figure 1 shows an example of a simulated binary star image. The images are then clipped to a maximum value, interpolated using the algorithm described in Sect. 2.1, and compared to ideal (non saturated) images to judge of the fidelity of the algorithm. Example of a simulated HST NICMOS2 128 seconds image through the F222M filter (2.22µm) including shot, readout noise and saturation. The plate scale is 76.5 mas per pixel and the image also features a companion at 10 3 contrast located 450mas left and 450mas down from the primary. The clipped value of nine pixels at the core of the PSF (displayed in white) would usually prevent the use of Fourier analysis.

The Kernel-Phase method and pipeline
Kernel analysis relies on the discrete model of the pupil of the telescope which was built using a dedicated python package called XARA 1 , developed in the context of the KER-NEL project. A discrete pupil model of 140 subapertures was constructed following a Cartesian grid covering the aperture of the telescope with an outer diameter of 1.95 m, with an inner obstruction diameter of 0.71 m and spider 0.07 m thickness spider arms. Figure 2 shows the sampling of the model both in the pupil plane and in the UV plane. According to the Nyquist-Shannon sampling theorem, the grid's pitch of b min = 0.13 m allows for detection up to ρ max = λ/2b min ≈ 1.76 arcsec at 2.22µm without aliasing.
The extraction of kernel-phases requires five steps: 1. Integer pixel recentering of the image. This procedure is applied to each image of a given dataset, and produces each time a vector κ of n k = 262 kernelphases.

Fidelity of the saturation recovery procedure
Since the aim of this study is to enable the extraction of kernel-phases, the main concern is with the fidelity of the 1 XARA is available at github.com/fmartinache/xara  A binary star image with a 10 3 contrast companion is simulated as described in 2.2 with a central peak intensity corresponding to four times the dynamic range of the sensor, resulting in nine clipped pixels. Figure 3 shows the Fourier phase compared between the saturated image where the outer edge of the uv plane is unusable, the ideal nonsaturated image, and the recovered image for which the phase information was restored.
Our recovery algorithm assumes that the image is dominated by the PSF of the primary star. The end result is Fig. 4. Illustration of both the deviation from the single star PSF that influence the fitting, and the bias of recovered pixel due to omitted feature signal. This configuration at a very low contrast of 3 and ∼ 1λ/D separation exaggerates the situation. therefore likely to be affected by the three following potential biases: 1. The signal of the companion acquired by the saturated pixels is lost and omitted by the recovery process as illustrated in orange by Fig. 4. 2. The signal of the companion still present in the image, but not in the ideal fitted PSF, therefore biasing the interpolated value; showed in blue in Fig. 4. 3. The instrumental PSF is distorted by a small amount of instrumental phase, therefore biasing the interpolated value. The impact of this bias was mitigated through the use of a uniform deviation estimation in section 2.1.
At high contrast, the effect of bias 1 is negligible compared to shot noise. For bias 2, it is much more difficult to evaluate as the effect is non-linear and highly dependent on the geometry of the PSF.
We characterized the combined effect of biases 1 and 2 on kernel-phases through a controlled noiseless simulation in a case with 4 saturated pixels, representative of the GL 494 case highlighted in Section 4. We measured the amplitude of the relative error e = ||κe−κr || ||κe|| , where κ e is the expected kernel-phase signal, and κ r is the recovered kernel-phase signal. In the high contrast regime, this metric sets an upper bound to the fitted contrast, but also generally encompasses fitting residuals in the whole kernel-phase subspace. As shown in Fig. 5, this metric mostly depends on the separation, and it drops below 2% at separations ≥ 200 mas.
The influence of biases 2 and 3 are both affected by the choice of weighting scheme mentioned in Sect. 2.1. Although the weighting function used produced satisfying results, further optimization work could lead to more improvements in the future.

Kernel analysis of NICMOS archival images
The dataset that we propose to examine comes from the Hubble Space Telescope proposal # 7420. The aim of the study was to observe nearby stars within 10pc of the Sun in order to find very low mass (M < 0.2M ) stars. Observations were conducted using the HST Near Infrared Camera and Multi-Object Spectrograph camera 2 through the F110W, F180M, F207M and F222M filters during from 1997 to 1998.
Analysis of the data was published by Krist et al. (1998) and Dieterich et al. (2012) where PSF subtraction was performed through manual adjustments. They determined detection limits by trying to identify artificially inserted companion stars over the subtraction residual. Contrast detection limits are given for the F180M filter from 0.2" to 4". In comparison, our model is sensitive between 0.12" to 1.76" due to its extent and resolution. In this work, we focus on the images taken in the F222M filter consisting of 159 images of 80 targets.

Saturated and bad pixels
The first step of the reduction process is the identification of problematic pixels. When considering the raw output of a camera, this might seem trivial: identify the pixels that have reached the maximum ADC value of the camera, often in the form 2 n bits in scientific sensors. However, this will prove erroneous when the gain of the camera is chosen for maximal dynamic range and the value is actually limited by the full well capacity of the pixel. This value can vary from pixel to pixel and is often more of a soft maximum showing a non-linear region (that is often linearized by advanced reduction pipelines). Furthermore, dark subtraction and flat-fielding will also play a role in making the clipped pixels harder to identify. For all these reasons, it is more appropriate to have the image reduction pipeline identify the saturated pixels and provide the information along with the meta-data. Here we relied on the metadata provided along with the image in the cal.fits output of the CALNICA re-duction pipeline. This is explained in more detail for our case in chapter 2.2.1 of the NICMOS Data Handbook by Thatte & Dahlen (2009).
All the images of the dataset were obtained using MUL-TIACCUM mode with a STEP128 sampling scheme for 128 seconds exposures. This is a mode that uses multiple nondestructive readouts to compute the flux on each pixel. If saturation is detected in some of the last samples, the first non-saturated samples can still be used to compute the flux. In this case, a pixel is only considered saturated if the first 0.303 seconds of exposure are sufficient to saturate the pixel. This explains the remarkable dynamic range of single images.
Bad pixels are also identified during this step. Bad pixels near the core of the PSF would be treated as saturated pixels and recovered, but pixels further out are simply interpolated with their nearest neighbors.

Whitening for the image plane errors
Usually in Fourier phase analysis, the image noise (mainly shot and readout noise) translates into correlated phase noise and, in turn, correlated kernel-phase noise. Such correlations must be accounted for before model-fitting or hypothesis testing. De-correlation can be performed using a whitening transformation described by the matrix where Σ K is an estimate of the covariance of κ, the kernelphase signal vector; therefore ensuring that: This is very similar to the approach of Ireland (2013) who uses the finite-dimensional spectral theorem to compute a unitary matrix that diagonalizes the covariance, then normalizes the observables by the corresponding standard deviation, which are the square roots of the terms of the diagonal matrix. In our case the same goal is reached through a non-unitary matrix W that also applies normalization, as shown by Eq. 8. Using this whitening matrix, the reduced goodness-of-fit parameter χ 2 ν writes as: where κ o , κ c and κ m respectively represent the kernelphase for the object of interest, the calibrator (explained in 3.5 and the model. κ m is computed using a parametric binary object model of n P = 3 parameters: angular separation ρ, position angle θ and contrast c, which leaves ν = n k − n p degrees of freedom at the denominator. The minimization of this χ 2 ν with a Levenberg-Marquardt algorithm allows us to identify the best fit model.

Estimation of the covariance of the Kernel-Phase
The evaluation of the covariance matrix Σ K of the calibrated kernel-phase vector κ = κ o − κ c is necessary for the whitening step (3.2). It is obtained for each targetcalibrator pair as: Although those covariances could be approximated analytically for normal images, the saturation recovery process used here is bound to introduce covariance even within the image plane, making the problem more complicated. Building a reliable covariance matrix estimator for the 262 observables vectors requires the acquisition of a large number of realizations. Since the data available only consists of two snapshots, we use a bootstrapping Monte Carlo approach. A number of realizations are simulated by adding noise to the real science image as described in Sect. 2.2. Those realizations are then pushed through the same pipeline as the science image as described in Sect. 2.1 for the recovery and Sect. 2.3 for the extraction of the observables. The only exception is the saturation mask which is kept identical. Since this process is rather computationally intensive because of the necessity to interpolate each image, we use 5000 realizations for each target, and the same covariance estimator is used for both images of the same target because they share the same exposure conditions, sensor location, and aberrations.

Colinearity maps
In addition to the three parameters χ 2 ν used for model fitting, colinearity maps can be built that constitute a matched filter as defined by Scharf (1990): where c 0 is a fixed high contrast value (10 2 is used in our case).
In the high contrast regime, the uv phase signal Φ 0 = arg(1+ 1 c e −i 2π λ (uα0+vδ0) ), and therefore of W·κ are inversely proportional to c. As a consequence, M (α, δ) will peak at coordinates [α, δ] = [α 0 , δ 0 ], regardless of their contrast. For this reason, this map is useful for 2D visualization of binary signals present in the kernel signature. It can be computed very quickly over a grid, and plotted as in Fig.  7.

Choice of calibrators
The calibration of interferometric observations requires the selection of point sources observed under similar conditions and with spectral and photometric properties similar to those of the target. The dataset associated with proposal # 7420 consists of images acquired on a wide variety of targets with respect to intrinsic brightness and spectral types; with observations separated by several weeks. We used two different procedures for: the detection of binary signals among all targets of the sample and; for the binary model fitting and parameter estimation.
For the detection, we built a generic calibrator through the elimination of images producing outlier kernel signals (in the sense of Euclidean distance to the median of the remaining subset) one at a time until 10 out of the original 159 remained. This conservative approach excludes occurrences of spurious signals caused by potential failures of the interpolation process, sensor defects and the mask alignment inconsistency mentioned by Krist et al. (1998), as well as targets with resolved features. The mean of the 10 remaining kernel-phase signals constitutes the generic Fig. 7. Examples of colinearity maps described in 3.4 for rejected calibrators (left) where signal from the calibrator is prominent, and selected calibrators (right) where the signal of interest is dominant and visible on the right hand side of the center. In the high contrast approximation, the value of each pixel peaks for calibrated signal vector colinear to the signal of a companion at the pixel's location, regardless of its contrast. Note that the map is antisymmetric, reflecting that kernel-phase is -like closure phases -a measure of the asymmetries of the target.
3.1 pixels per resolution element. The kernel-phase vectors of both images were then averaged to reduce noise. The steps described in Sect. 3.5 are followed to obtain a detection, then a distribution of the parameters of the binary which are provided in the first row of table 2 for use in the orbital determination. The kernel-phase analysis provides a standard deviation 5 mas in separation and 2.3 deg in position angle which is significantly lower than what was obtained for the same target with larger ground based telescopes. We also measure the contrast at ∆ K = 4.36 ± 0.15 which is consistent with the initial PUEO measurement by Beuzit et al. (2004), and later measurements by Mann et al. (2018) showing ∆ Ks = 4.269 ± 0.017, but not with the NACO measurements presented by Ward-Duong et al. (2015). An instrumental bias in the NACO measurements, or the variability of the active primary could explain these differences.

Determination of orbital parameters
In order to determine orbital parameters for the system, we use the radial velocity data obtained with the ELODIE spectrograph (Baranne et al. 1996) at Observatoire de Haute-Provence between 2000 and 2006 obtained through the ELODIE archive, and the CORAVEL spectrograph between 1983 and 1995. The standard deviation figure used for ELODIE data was 75 m/s. Two data points acquired with the SOPHIE spectrograph are also available but were not included in the model because it was not worth calibrating for one more instrument. The parallax of the system is given by the Gaia DR2 catalog at π = 86.86 ± 0.1515 mas (Gaia Collaboration et al. 2018; Gaia Collaboration 2016) which gives a distance of d = 11.513 ± 0.02 pc. The technique employed to fit the orbital parameters with both the visual orbit and the radial velocity data is similar to the one described by Martinache et al. (2007). The lmfit package (Newville et al. 2014), a python implementation of the Levenberg-Marquardt algorithm is used to obtain a least-square fit for a 10 degrees of freedom problem, resulting in the parameters presented in table 3 and Fig. 9.
The apparent magnitude of the system from 2MASS m K = 5.578±0.016 together with the distance and contrast provide us with absolute magnitude of M K = 5.29 ± 0.017 for Gl 494A and M K = 9.65 ± 0.15 for Gl 494B. Figure 10 shows how this compares to the expected mass-luminosity relationship provided by Benedict et al. (2016) for the M dwarfs in the solar neighborhood. The two components are found to be in reasonable agreement with the model.

Conclusion
We have demonstrated using both simulation and on-sky data that it is possible to use our knowledge of the PSF to recover the saturated pixels of images. Unlike the original images, the mended images can then be used successfully in kernel-phase analysis. The proposed method works particularly well in the case of images dominated by the PSF of a single star, like in the case of high contrast binary stars or exoplanetary systems. Its primary purpose is to expand the field of applicability of kernel-phase methods, especially to the larger body of archival space telescope images.
As an example, we have used the technique on images from the HST of the Gl 494AB system with a saturated core, and produced a new visual astrometry data point at an interesting older epoch and small uncertainties. Combining this information with radial velocity data, we have improved the orbital model of the system and determined the masses of its two components.
As remarked by Torres (1999), providing high precision visual measurement of astrometric and spectroscopic binary systems is key to obtain accurate orbital parameters and masses, and this work shows once more how kernel analysis can play a major role in the follow-up of the binary systems --0.553 ± 0.007 that will be discovered in large quantities by the Gaia mission. The possibility of using saturated images could lead us to reevaluate how to optimize the observation strategy for this goal, as well as for the discovery of new stellar binary systems, especially in the context of the James Webb Space Telescope.  (1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) and ELODIE (2000ELODIE ( -2006, and the corresponding fitted model results.