© R. Laugier et al. 2019

1. Introduction

The quest for the direct imaging of exoplanets has accelerated in 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⁴) 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, however, have 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, interferometric techniques designed to be robust to the wavefront errors, like closure phases (Jennison 1958) or kernel-phases (KPs; Martinache 2010), are competitive tools for discovering faint companions around nearby stars. The aim of the KERNEL project is to develop observation and image processing techniques that provide observables that are intrinsically robust to small wavefront errors, and that 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 a 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) violates some of the prerequisites of the Fourier transform. As a result the shift theorem and object-image convolution relationship become unusable, preventing the interpretation of the phase and the construction of kernels. Although the limits were pushed by Pope et al. (2016), who had to linearize the sensor’s response (soft saturation), and by 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 extracts 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 radial velocity observations, lead to improved constraints on the orbital parameters of the system and therefore on the mass of its components.

2. Interpolation of saturated images and kernel analysis

Typical algorithms for the correction of bad pixels assume that problematic pixels are few and isolated, allowing recovery by interpolating the neighboring pixels in the image or by minimizing the power associated with 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 values can then be reliably extracted from these enhanced images.

2.1. 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 nonsaturated 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 χ² variable

$$\begin{array}{c}\hfill {\chi }^2(x,y,z)=||\mathbf{i}-z·{\mathbf{p}(x,y)||}^2,\end{array}$$(1)

where $\mathbf{i}={\mathbf{\Sigma }}^{-\frac{1}{2}}·{\mathbf{s}}_{\mathrm{img}}$ and $\mathbf{p}(x,y)={\mathbf{\Sigma }}^{-\frac{1}{2}}·{\mathbf{s}}_{\mathrm{PSF}}(x,y)$ are, respectively, error-normalized (whitened) vectors of the image signal s_img and the (x, y) shifted PSF signal s_PSF(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, a scale at which the problem is convex. Minimizing χ²(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, we 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 S/N,

$$\begin{array}{c}\hfill g={\mathrm{e}}^{-{(\frac{r}{r_0})}^4},\end{array}$$(2)

where r is the distance from the pixel to the approximated center of the star and r₀ is a radius parameter chosen depending on the S/N 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 ${\sigma }_{ii}^2$ terms of a diagonal matrix Σ. A constant value of r₀ = 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),

$$\begin{array}{c}\hfill \widehat{z}(x,y)=\frac{\mathbf{p}{(x,y)}^t·\mathbf{i}}{\mathbf{p}{(x,y)}^t·\mathbf{p}(x,y)},\end{array}$$(3)

leading, through the substitution in Eq. (1) to the definition of a new criterion, a function of x and y only:

$$\begin{array}{c}\hfill \epsilon (x,y)=-\frac{{(\mathbf{p}{(x,y)}^t·\mathbf{i})}^2}{\mathbf{p}{(x,y)}^t·\mathbf{p}(x,y)}·\end{array}$$(4)

Minimizing ε(x, y) is therefore equivalent to minimizing χ²(x, y, z), but is computationally more efficient.

2.2. Simulation of realistic NICMOS images

The simulation of images is necessary for two purposes: the evaluation of the fidelity of the algorithm (Sect. 2.4) and the bootstrapping of the kernel-phase covariance matrix (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 nondestructive reads was also emulated in order to reproduce the typical behavior of the sensor in MULTIACCUM/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

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{ro}}=\frac{{\sigma }_{\mathrm{ro},\mathrm{single}}}{\sqrt{N_{\mathrm{s}}}},\end{array}$$(5)

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 is not 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

$$\begin{array}{c}\hfill {\mathrm{ERR}}_{\mathrm{total}}=\sqrt[]{\frac{\mathrm{SCI}}{G\times \mathrm{TIME}}+{\mathrm{ERR}}^2},\end{array}$$(6)

where SCI is the flux in counts per second, G is the ADC inverse gain (in electrons/ADU), ERR is the noise map, and TIME 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.

Fig. 1. Example of a simulated HST NICMOS2 128 s 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 103 contrast located 450 mas to the left and 450 mas below 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 images are then clipped to a maximum value, interpolated using the algorithm described in Sect. 2.1, and compared to ideal (nonsaturated) images to judge the fidelity of the algorithm.

2.3. Kernel-phase method and pipeline

Kernel analysis relies on the discrete model of the pupil of the telescope that was built using a dedicated python package called XARA¹, developed in the context of the KERNEL 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 detection up to ρ_max = λ/2b_min ≈ 1.76 arcsec at 2.22 μm without aliasing.

Fig. 2. Discrete representation of the HST for kernel analysis of the NICMOS images. Left panel: representation of the pupil model defined by the cold mask of the instrument. The model takes into account central obstruction and spider vanes, but neglects finer structures such as support pads. The grid pitch is 0.13 m, the outer diameter is 1.95 m, the inner diameter is 0.71 m, and the thickness of the spider vanes is 0.07 m. Right panel: corresponding uv-plane coverage.

The extraction of kernel-phases requires five steps:

1. Integer pixel recentering of the image;

2. Application of an apodization mask to reduce the effect of readout noise and avoid Fourier domain ringing;

3. Subpixel recentering of the image by application of a wedge phasor in the Fourier domain;

4. Extraction of the complex visibility vector;

5. Multiplication of the phase of the visibility by the Kernel matrix.

This procedure is applied to each image of a given dataset, and produces each time a vector κ of n_k = 262 kernel-phases.

2.4. 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 phase of the Fourier transform of the image. The discretized pupil model (see Sect. 2.3) provides a comprehensive sampling scheme of the Fourier plane.

A binary star image with a 10³ contrast companion is simulated, as described in Sect. 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.

Fig. 3. Representation of the Fourier plane phase (in radians) sampled for an example of saturated image of a binary star (shown in Fig. 1) for the ideal nonsaturated image, for the recovered image, and for the residual. For the original image (upper left) the saturation results in large amplitude aberrations in the outer ring. After recovery (upper right), the phase is very close to what it would be if the sensor had remained linear (lower left). The difference between the two (lower right) shows small residuals concentrated in the outermost ring.

Our recovery algorithm assumes that the image is dominated by the PSF of the primary star. The end result is 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 (in orange in 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 (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 (Sect. 2.1).

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.

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 nonlinear 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 four saturated pixels, representative of the GL 494 case highlighted in Sect. 4. We measured the amplitude of the relative error $e=\frac{||{\mathit{\kappa }}_{\mathrm{e}}-{\mathit{\kappa }}_{\mathrm{r}}||}{||{\mathit{\kappa }}_{\mathrm{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.

Fig. 5. Relative norm of the kernel-phase bias due to saturation-recovery process (normalized by the amplitude of the companion signal) obtained by noiseless simulations. The colored region indicates the range of values for contrasts between 10 and 250, while the solid line indicates their mean. As we can see, the bias quickly drops to below 2% past the ∼1λ/D = 235 mas separation.

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.

3. 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 10 pc of the Sun in order to find very low-mass stars (M <  0.2M). Observations were conducted from 1997 to 1998 using the HST NICMOS camera 2 through the F110W, F180M, F207M, and F222M filters.

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″ and 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.

3.1. 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 maximum 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 nonlinear region (which 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 metadata. Here we relied on the metadata provided along with the image in the cal.fits output of the CALNICA reduction 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 in the dataset were obtained using MULTIACCUM mode with a STEP128 sampling scheme for 128 s 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 nonsaturated samples can still be used to compute the flux. In this case, a pixel is only considered saturated if the first 0.303 s 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 farther out are simply interpolated with their nearest neighbors.

3.2. 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. These correlations must be accounted for before model-fitting or hypothesis testing. Decorrelation can be performed using a whitening transformation described by the matrix

$$\begin{array}{c}\hfill \mathbf{W}={\mathbf{\Sigma }}_{\mathrm{K}}^{-\frac{1}{2}},\end{array}$$(7)

where Σ_K is an estimate of the covariance of κ, the kernel-phase signal vector, therefore ensuring that

$$\begin{array}{c}\hfill \text{cov}(\mathbf{W}·\mathit{\kappa })\approx \mathbf{I}.\end{array}$$(8)

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 ${\chi }_{\nu }^2$ writes as

$$\begin{array}{c}\hfill {\chi }_{\nu }^2=\frac{||\mathbf{W}·{\mathit{\kappa }}_{\mathrm{o}}-\mathbf{W}·{\mathit{\kappa }}_{\mathrm{c}}-\mathbf{W}·{\mathit{\kappa }}_{\mathrm{m}}{||}^2}{n_{\mathrm{k}}-n_{\mathrm{p}}},\end{array}$$(9)

where κ_o, κ_c, and κ_m respectively represent the kernel-phase for the object of interest, the calibrator (explained in 3.5) and the model. The value of κ_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 ${\chi }_{\nu }^2$ with a Levenberg–Marquardt algorithm allows us to identify the best fit model.

3.3. 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 target-calibrator pair as

$$\begin{array}{c}\hfill {\mathbf{\Sigma }}_{\mathrm{K}}=\text{cov}({\mathit{\kappa }}_{\mathrm{o}})+\text{cov}({\mathit{\kappa }}_{\mathrm{c}}).\end{array}$$(10)

Although these covariances can 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 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. These realizations are then pushed through the same pipeline as the science image, as described in Sect. 2.1 for the recovery and in 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; the same covariance estimator is used for both images of the same target because they share the same exposure conditions, sensor location, and aberrations.

3.4. Colinearity maps

In addition to the three parameters ${\chi }_{\nu }^2$ used for model fitting, colinearity maps can be built that constitute a matched filter defined by Scharf (1990) as

$$\begin{array}{c}\hfill M(\alpha ,\delta )=\frac{{(\mathbf{W}·{\mathit{\kappa }}_{\mathrm{o}}-\mathbf{W}·{\mathit{\kappa }}_{\mathrm{c}})}^T·\mathbf{W}·{\mathit{\kappa }}_{\mathrm{m}}(\alpha ,\delta ,c_0)}{||\mathbf{W}·{\mathit{\kappa }}_{\mathrm{m}}(\alpha ,\delta ,c_0)||},\end{array}$$(11)

where c₀ is a fixed high-contrast value (10² is used in our case).

In the high-contrast regime, the uv phase signal ${\mathbf{\Phi }}_0=\arg (1+\frac{1}{c}{\mathrm{e}}^{-i\frac{2\pi }{\lambda }(u{\alpha }_0+v{\delta }_0)})$ and therefore of W ⋅ κ are inversely proportional to c. As a consequence, M(α, δ) will peak at coordinates [α, δ]=[α₀, δ₀] 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.

3.5. 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, one for the detection of binary signals among all targets of the sample and one 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 ten remaining kernel-phase signals constitutes the generic calibration signal, subtracted from all raw extracted kernel-phases, that allowed the identification of the companion in the saturated images of Gl 494, as shown in Fig. 6.

Fig. 6. Correlation plot showing the good agreement between the Gl494 data and the fitted binary model. Here the data is whitened and calibrated with the mean of list of ten selected calibrators. The covariance matrix used for whitening takes into account the covariance of the target, the covariances of the calibrators, and the variance between calibrators.

Fig. 7. Examples of colinearity maps described in Sect. 3.4 for rejected calibrators (left panels) where the signal from the calibrator is prominent, and selected calibrators (right panels) where the signal of interest is dominant and visible on the right. In the high-contrast approximation, the value of each pixel peaks for the calibrated signal vector colinear to the signal of a companion at the pixel’s location, regardless of its contrast. The map is antisymmetric, reflecting that the kernel-phase is – like the closure phases – a measure of the asymmetries of the target.

Evaluating the impact of this calibration-induced bias on the extracted binary parameters requires as many calibrators as possible. Suitable ones were selected using visual examination of colinearity maps for Gl 494 after subtraction of each of the 79 other possible calibrators (see Fig. 7). Candidates for which the map was not dominated by the binary signal identified at the first stage were discarded, leaving the 32 targets listed in Table 1. The distribution of the parameters fitted after subtraction of each calibrator signal was used to estimate the confidence intervals for the final measurement.

4. Gl 494

The binary system Gl 494AB was first identified by Heintz (1994) as an astrometric binary with a 14.5 year period. The system was later resolved with the PUEO adaptive optics (Véran et al. 1999) at the Canada-France-Hawaii Telescope (CFHT) with a brief summary of the properties of the system (Beuzit et al. 2004). The dynamical masses derived remained inconclusive because they relied on the Heintz orbital parameters, which are given without standard errors, and show some inconsistency with the position angle. The system was later resolved again with VLT NaCo (Lenzen et al. 2003; Rousset et al. 2003) in 2005 and 2006 as published by Ward-Duong et al. (2015), but the orbital parameters were only recently updated by Mann et al. (2019).

The Hubble Space Telescope archive images of the system used here went overlooked in the early studies, probably for two combined reasons: the system is close to the resolution limit of the telescope (≈1.3λ/D) with a high contrast (≈55), and the central peak of the images is saturated. Contrast sensitivities for this target in the same dataset in the F180M filter were later reexamined by Dieterich et al. (2012) to be Δ_H = 0.1 at 0.2″ and Δ_H = 2.6 at 0.4″ on this target, which is insufficient for this system. Figure 8 shows a crop of the saturated image used here, and the corresponding unusable uv phase signal. Overcoming the saturation problem allows us to use the more sensitive kernel analysis and adds a valuable high-precision visual astrometric point. With the recent observations reported by Mann et al. (2019), the visual observations of this binary system now span 18 years, covering more than one orbital period. Together with the radial velocity data, these observations give strong constraints on the orbital parameters of the system, and therefore on the masses of its components.

Fig. 8. Left panel: cropped image of Gl 494 taken in F222M filter in log scale before recovery. The four saturated pixels appear in white. The companion is embedded within the first airy ring, to the right of the primary, and not discernible by eye. Right panel: Fourier phase of Gl 494 image with saturation. The spurious phase signal visible in the outer ring would usually prevent interpretation.

4.1. Extraction of the visual astrometry

The method was applied to images of the system Gl 494 taken in August 1997 with the HST NICMOS camera 2 through the F222M filter. This configuration provides 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 were followed to obtain a detection, then a distribution of the parameters of the binary 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.3deg in position angle, which is significantly lower than the value obtained for the same target with larger ground-based telescopes. We also measured the contrast at Δ_K = 4.36 ± 0.15, which is consistent with the initial PUEO measurement by Beuzit et al. (2004) and with later measurements by Mann et al. (2019) 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.

4.2. Determination of orbital parameters

In order to determine orbital parameters for the system, we used the radial velocity data obtained with the ELODIE spectrograph (Baranne et al. 1996) at the 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 2018, 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 that 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-squares fit for a 10° of freedom problem, resulting in the parameters presented in Table 3 and Fig. 9.

Fig. 9. Orbital fitting for the Gl 494AB system. Upper panel: astrometric measurements and corresponding fitted model results. The 1997 point is the novel measurement provided by this work. Lower panel: radial velocity measurements obtained with CORAVEL (1983–1995) and ELODIE (2000–2006), and the corresponding fitted model results.

The apparent magnitude of the system from 2MASS m_K = 5.578 ± 0.016 together with the distance and contrast provide us with an 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.

Fig. 10. Comparison of the mass and absolute magnitudes obtained from this work with the mass-luminosity relationship for the low-mass stars in the solar neighborhood, showing reasonable agreement.

5. Conclusion

We demonstrated via 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. 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 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 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 obtaining 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 large numbers of binary systems that will be discovered by the Gaia mission. The possibility of using saturated images could lead us to reevaluate how to optimize the observation strategy for this goal, and for the discovery of new stellar binary systems, especially in the context of the James Webb Space Telescope.

Acknowledgments

The authors wish to thank Xavier Delfosse, Stephane Udry, Damien Ségransan, and the CORAVEL team at the Observatoire Astronomique de l’Université de Genève for allowing us to use the data from the CORAVEL instrument; and Sébastien Peretti for helping us with the data. 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). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

References

All Tables

All Figures

Fig. 1. Example of a simulated HST NICMOS2 128 s 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 103 contrast located 450 mas to the left and 450 mas below 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.

In the text

Fig. 2. Discrete representation of the HST for kernel analysis of the NICMOS images. Left panel: representation of the pupil model defined by the cold mask of the instrument. The model takes into account central obstruction and spider vanes, but neglects finer structures such as support pads. The grid pitch is 0.13 m, the outer diameter is 1.95 m, the inner diameter is 0.71 m, and the thickness of the spider vanes is 0.07 m. Right panel: corresponding uv-plane coverage.

In the text

Fig. 3. Representation of the Fourier plane phase (in radians) sampled for an example of saturated image of a binary star (shown in Fig. 1) for the ideal nonsaturated image, for the recovered image, and for the residual. For the original image (upper left) the saturation results in large amplitude aberrations in the outer ring. After recovery (upper right), the phase is very close to what it would be if the sensor had remained linear (lower left). The difference between the two (lower right) shows small residuals concentrated in the outermost ring.

In the text

	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.
In the text

	Fig. 5. Relative norm of the kernel-phase bias due to saturation-recovery process (normalized by the amplitude of the companion signal) obtained by noiseless simulations. The colored region indicates the range of values for contrasts between 10 and 250, while the solid line indicates their mean. As we can see, the bias quickly drops to below 2% past the ∼1λ/D = 235 mas separation.
In the text

	Fig. 6. Correlation plot showing the good agreement between the Gl494 data and the fitted binary model. Here the data is whitened and calibrated with the mean of list of ten selected calibrators. The covariance matrix used for whitening takes into account the covariance of the target, the covariances of the calibrators, and the variance between calibrators.
In the text

Fig. 7. Examples of colinearity maps described in Sect. 3.4 for rejected calibrators (left panels) where the signal from the calibrator is prominent, and selected calibrators (right panels) where the signal of interest is dominant and visible on the right. In the high-contrast approximation, the value of each pixel peaks for the calibrated signal vector colinear to the signal of a companion at the pixel’s location, regardless of its contrast. The map is antisymmetric, reflecting that the kernel-phase is – like the closure phases – a measure of the asymmetries of the target.

In the text

	Fig. 8. Left panel: cropped image of Gl 494 taken in F222M filter in log scale before recovery. The four saturated pixels appear in white. The companion is embedded within the first airy ring, to the right of the primary, and not discernible by eye. Right panel: Fourier phase of Gl 494 image with saturation. The spurious phase signal visible in the outer ring would usually prevent interpretation.
In the text

	Fig. 9. Orbital fitting for the Gl 494AB system. Upper panel: astrometric measurements and corresponding fitted model results. The 1997 point is the novel measurement provided by this work. Lower panel: radial velocity measurements obtained with CORAVEL (1983–1995) and ELODIE (2000–2006), and the corresponding fitted model results.
In the text

	Fig. 10. Comparison of the mass and absolute magnitudes obtained from this work with the mass-luminosity relationship for the low-mass stars in the solar neighborhood, showing reasonable agreement.
In the text