The Near Infrared Imager and Slitless Spectrograph for JWST. V. Kernel Phase Imaging and Data Analysis

The STScI jwstdata reduction pipeline is organized into three stages. Stage 1 performs detector level calibrations such as bias correction and jump detection, stage 2 performs photometric calibrations, and stage 3 performs high-level calibrations that are specific to each of the observing modes supported by JWST (e.g., AMI, coronagraphic imaging, integral field spectroscopy). To enable the community a straightforward access to kernel phase observables from JWST images, we develop a custom stage 3 pipeline that we name the Kpi3Pipeline. ²⁵ This pipeline can be used for extracting kernel phase observables from full pupil NIRISS and NIRCam images, support for MIRI will be added in a future update. Before the JWST images can be fed into the Kpi3Pipeline, they need to be processed with the jwststage 1 and 2 pipelines. For those, similar as with AMI data, we recommend to skip the inter-pixel capacitance, the photometry, and the resample step since Fourier plane imaging techniques in general are highly sensitive to variations in the pixel-by-pixel response of the detector (e.g., Ireland 2013). Skipping these steps avoids systematic errors being introduced by the pipeline corrections. Instead, the approach with NIRISS KPI and AMI observations is to center the science and the reference targets on the exact same detector pixel and thereby minimize flat-fielding errors. The available detector positions were carefully chosen to provide both a good uv-coverage of the pupil support (which is challenging with NIRISS given its coarse sampling) and a clean detector region with only few bad pixels (Sivaramakrishnan et al. 2023). Dithering is hence not recommended for NIRISS KPI (and AMI) observations. While the backend of the Kpi3Pipelineis based on XARA ²⁶ (Martinache 2010, 2013; Martinache et al. 2020), the front end is designed similarly to the STScI jwstpipelines to provide a uniform experience for the user. Similar to the VLT NACO kernel phase data reduction pipeline developed by Kammerer et al. (2019), the Kpi3Pipelineconsists of five steps:

• 1.  
  bad pixel fixing step,
• 2.  
  recentering step,
• 3.  
  windowing step,
• 4.  
  kernel phase extraction step,
• 5.  
  empirical uncertainties step.

An intermediate data product can be saved after each step so that the individual steps can be run separately and each of the five steps can generate a diagnostic plot for validation purposes. Each step and its diagnostic plot is described in more detail below. A list of the tunable parameters of each step can be found in Appendix B (Table 4). In principle, all steps except for the fourth one can be skipped although it is highly recommended to run at least steps 1, 2, and 4 to ensure that the kernel phase extraction step is provided with properly preprocessed data.

An additional input that is required for running the Kpi3Pipeline is a discrete model of the pupil. We note that due to different pupil plane masks, the NIRISS, NIRCam, and MIRI instruments all see slightly different pupils. We provide default pupil models for the NIRISS CLEARP and the NIRCam CLEAR pupils but the user can also provide a custom pupil model if desired. For the default pupil models, we distribute subapertures on a hexagonal grid so that 19 subapertures fall inside one JWST primary mirror segment. To obtain an isotropic grid, we also distribute subapertures on top of the gaps between the individual primary mirror segments. We are using a "gray" pupil model so that the subapertures on top of the gaps have a slightly reduced transmission. Then, we apply the transmission of the instrument-specific pupil mask and discard all subapertures with a total transmission of <0.7. This results in the NIRISS CLEARP pupil model shown in Figure 1. This model consists of 349 individual subapertures spanning 948 distinct baselines. After discarding all baselines with a redundancy of less than 10, 800 distinct baselines yielding 452 individual kernel phases remain. Discarding baselines with low redundancy helps to avoid the longest (and most noisy) baselines at the edge of the primary mirror. We find that a redundancy threshold of 10 is sufficient to avoid Fourier phases >0.5 rad and entering the non-linear regime.

Zoom In Zoom Out Reset image size

2.1.1. Bad Pixel Fixing Step

Bad pixels can have a ruinous impact on Fourier plane imaging since an individual bad pixel contributes noise to all spatial frequencies of the image. However, Ireland (2013) has shown that bad pixels do not affect kernel phase observables if they are properly corrected. In the Kpi3Pipeline, bad pixels can be fixed using two different methods: "medfilt" or "fourier". In both cases, the locations of the bad pixels are obtained from the jwst stage 1 pipeline. As a default, all pixels flagged as "DO_NOT_USE" in the data quality array are considered as bad pixels but the user may specify other data quality flags ²⁷ that shall be considered as bad pixels. With the "medfilt" method, bad pixels are simply replaced with the median filtered image using a kernel size of five pixels. With the "fourier" method, bad pixels are fixed by minimizing their Fourier power outside the region of support permitted by the pupil geometry. This method was introduced by Ireland (2013) and starts by computing the matrix B _Z which maps the bad pixel values x onto the Fourier domain Z (the complement of the pupil support), so that their Fourier power $\left|{f}_{Z}\right|$ is given by

$\begin{eqnarray}&&{f}_{Z}={{\boldsymbol{B}}}_{Z}\cdot b+{\epsilon }_{Z},\end{eqnarray} \tag{ 6 }$

where b are the corrections to be made to the bad pixels and _Z is remaining noise. These corrections are obtained by computing the Moore-Penrose pseudo inverse of B _Z and solving for

$\begin{eqnarray}&&b={{\boldsymbol{B}}}_{Z}^{+}\cdot {f}_{Z}={\left({{\boldsymbol{B}}}_{Z}^{* }\cdot {{\boldsymbol{B}}}_{Z}\right)}^{-1}\cdot {{\boldsymbol{B}}}_{Z}^{* }\cdot {f}_{Z},\end{eqnarray} \tag{ 7 }$

where the star denotes the complex conjugate. This algorithm has been successfully applied to fix bad and reconstruct saturated pixels in kernel phase observations by Kammerer et al. (2019) and is described in more detail in Sivaramakrishnan et al. (2023). While the "fourier" method is particularly suited for reconstructing saturated PSFs, we do not have any saturated PSFs in the NIRISS KPI commissioning data and therefore use the "medfilt" method for simplicity. The diagnostic plot of the bad pixel fixing step is shown in Figure 2.

Zoom In Zoom Out Reset image size

2.1.2. Recentering Step

As mentioned in Section 1, the image Fourier phase needs to be in the linear regime in order for the kernel phase technique to be applicable. Hence, recentering the images is important since a PSF being offset from the image center by only half of a pixel would already result in a linear ramp in the image Fourier phase from −π/2 to +π/2 across the support of the telescope pupil (i.e., the modulation transfer function of the pupil). While the kernel phase observables are independent of first order pupil plane phase errors (and thus a linear ramp in the image Fourier phase, Martinache 2010), there are virtually always systematic phase errors that can easily cause the image Fourier phase wrapping around the ±π discontinuity if added on top of a linear phase ramp and thereby destroying the properties of the kernel phase observables. To avoid this, the Kpi3Pipeline uses the Fourier phase norm minimization (FPNM) method implemented in XARA that recenters the images using a least-squares gradient descent minimization of the norm of the image Fourier phase within the support of the telescope pupil. This method was shown to be more robust, especially in the presence of low-contrast companions at small angular separation (Kammerer et al. 2019), albeit slightly slower than the other two methods (BCEN = centroid of the brightest speckle & COGI = center of gravity of the image) that are also implemented in XARA. We note that before recentering the images, we trim them to their maximum possible square size (constrained by the location of the detector edges) around the expected PSF center (the reference pixel position where JWST aims to place the PSF during target acquisition ²⁸ ). The diagnostic plot of the recentering step is shown in Figure 3.

Zoom In Zoom Out Reset image size

2.1.3. Windowing Step

Windowing with a smooth function is an effective method to mitigate artifacts when numerically Fourier transforming images with sharp edges. Since XARA uses a Fourier transform to extract kernel phase observables from the telescope images, we window (i.e., multiply) them with a super-Gaussian windowing function w of shape

$\begin{eqnarray}&&w(x,y)=\exp \left(-\displaystyle \frac{{\left({x}^{2}+{y}^{2}\right)}^{2}}{{r}^{4}}\right)\end{eqnarray} \tag{ 8 }$

before Fourier transforming them, where x and y are the pixel coordinates with respect to the image center and r is the radius of the super-Gaussian windowing function in units of pixels. The radius r is chosen adaptively as the smaller of 40 pixels and a fourth of the trimmed square image size to make sure that the window is not larger than the image itself. If the images are large enough, a radius of 40 pixels gives access to separations of at least 5 λ/D (sometimes more depending on the observing wavelength). At these separations, classical high-contrast imaging techniques such as coronagraphy are typically superior to KPI. The NIRISS full pupil images considered for the KPI analysis here do only have a native size of 80 by 80 pixels and we are using a constant window radius of r = 15 pixels for them resulting in a field-of-view of ∼2'', which is more than sufficient for the region where NIRISS KPI is superior to NIRCam coronagraphy (Kammerer et al. 2022). The diagnostic plot of the windowing step is shown in Figure 4.

Zoom In Zoom Out Reset image size

2.1.4. Kernel Phase Extraction Step

Once the bad pixels are fixed and the images are recentered and windowed, kernel phase observables can be extracted using XARA. First, the image Fourier phase ϕ at the discrete spatial frequencies of the pupil model is obtained using a linear discrete Fourier transform. Then, the kernel phase θ is computed as

$\begin{eqnarray}&&\theta =\tilde{{\boldsymbol{K}}}\cdot \phi ,\end{eqnarray} \tag{ 9 }$

where $\tilde{{\boldsymbol{K}}}={\boldsymbol{K}}\cdot {\boldsymbol{R}}$ and K is the kernel of the pupil model A as described in Martinache (2010). We also decided to normalize the rows of the $\tilde{{\boldsymbol{K}}}$ matrix to one so that the amplitude of the kernel phase signal is no longer depending on the geometry of the pupil model. This enables a more direct comparison between kernel phase observables obtained from different telescopes/instruments. To analytically estimate the kernel phase uncertainties, we first need to linearize the relationship between the kernel phase θ and the image I, so that

$\begin{eqnarray}&&\theta \approx {\boldsymbol{B}}\cdot I=\tilde{{\boldsymbol{K}}}\cdot \left(\displaystyle \frac{{\mathfrak{Im}}({\boldsymbol{F}})}{\left|{\boldsymbol{F}}\cdot I\right|}\right)\cdot I,\end{eqnarray} \tag{ 10 }$

where F denotes the linear discrete Fourier transform at the spatial frequencies of the pupil model and the fraction denotes element-wise division (Kammerer et al. 2019). Then, the kernel phase covariance Σ can be computed from the image covariance Σ_image as

$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={\boldsymbol{B}}\cdot {{\boldsymbol{\Sigma }}}_{\mathrm{image}}\cdot {{\boldsymbol{B}}}^{T},\end{eqnarray} \tag{ 11 }$

where Σ_image is assumed to be uncorrelated with the square of the "ERR" extension of the stage 2-reduced pipeline product on the diagonal.

If the recentering step is being run in conjunction with the kernel phase extraction step, a more accurate method is used to extract the kernel phase observables where the integer pixel recentering is performed by rolling the image and the subpixel recentering is performed by adding a linear ramp to the extracted image Fourier phase. This method circumvents interpolation errors when performing subpixel shifts of the images. We note that for extracting the kernel phase uncertainties, we still use the images that have been recentered with subpixel precision. The diagnostic plot of the kernel phase extraction step is shown in Figure 5.

Zoom In Zoom Out Reset image size

2.1.5. Empirical Uncertainties Step

A single JWST data product (here called an exposure) typically consists of a large number (≳10) of individual integrations. While all previous pipeline steps can be performed on hundreds of individual images within timescales of seconds to minutes, it can be very time- and memory-consuming to perform model fitting on such a large number of individual data points, especially if accounting for correlated uncertainties. Hence, it is advisable to average individual integrations within a single exposure into one final data point. This is achieved in the Kpi3Pipeline by computing the covariance-weighted mean of the kernel phase observables according to

$\begin{eqnarray}&&{\theta }_{\mathrm{mean}}={{\boldsymbol{\Sigma }}}_{\mathrm{mean}}\cdot \sum _{i=1}^{{N}_{{\rm{f}}}}{{\boldsymbol{\Sigma }}}_{i}^{-1}\cdot {\theta }_{i},\end{eqnarray} \tag{ 12 }$

where θ_i and Σ_i are the kernel phase and its covariance of the individual integrations and

$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mathrm{mean}}={\left(\sum _{i=1}^{{N}_{{\rm{f}}}}{{\boldsymbol{\Sigma }}}_{i}^{-1}\right)}^{-1}\end{eqnarray} \tag{ 13 }$

is the mean covariance matrix (Kammerer et al. 2019).

Another advantage of having several individual integrations is that the uncertainties can be estimated empirically as the standard deviation of the kernel phase observables over these individual integrations. This can be especially helpful if the kernel phase observables are dominated by systematic errors (e.g., Kammerer et al. 2019; Wallace et al. 2020) which are not reflected in the analytical error estimation from the previous pipeline step. An empirical estimation for the kernel phase covariance can then be obtained by computing the correlation matrix C _mean of Σ_mean (which only includes image pixel uncertainties) and multiplying it with the empirically estimated kernel phase standard deviation σ_emp, that is

$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mathrm{emp},{nm}}={{\boldsymbol{C}}}_{\mathrm{mean},{nm}}{\sigma }_{\mathrm{emp},m}{\sigma }_{\mathrm{emp},n}.\end{eqnarray} \tag{ 14 }$

Here, m and n denote the indices of the matrices. This procedure ensures that the kernel phase uncertainties remain consistent with the square root of the diagonal of the kernel phase covariance matrix. The diagnostic plot of the empirical uncertainties step is shown in Figure 6.

Zoom In Zoom Out Reset image size

A common and well-defined data exchange format has proven highly valuable for sharing astronomical data and reproducing observational results. In the optical and near-infrared (long-baseline) interferometry community, the OIFITS file format (Pauls et al. 2005; Duvert et al. 2017) has established itself as a golden standard that is used by most major observatories and instruments (e.g., CHARA/MIRC, Monnier et al. 2004, VLTI/PIONIER, Le Bouquin et al. 2011, VLTI/GRAVITY, Lapeyrere et al. 2014) and data analysis tools (e.g., LITpro, Tallon-Bosc et al. 2008, CANDID, Gallenne et al. 2015). Due to the very similar nature of the data, the same OIFITS file format is also being used for AMI data (e.g., in ImPlaneIA, Greenbaum et al. 2015, AMICAL, Soulain et al. 2020, SAMpip ²⁹ ). Although KPI and AMI are closely related techniques (see, e.g., Ireland 2016), there are a few subtle differences that strongly motivate the need for a distinct file format for KPI data. First, as shown in Section 1, the kernel phase θ is obtained from a special linear combination $\tilde{{\boldsymbol{K}}}$ of the image Fourier phase ϕ. Thus, the kernel matrix $\tilde{{\boldsymbol{K}}}$ is of fundamental importance in KPI and is repeatedly used in the model fitting process to project image Fourier phase onto kernel phase. Fast and easy access to this kernel matrix $\tilde{{\boldsymbol{K}}}$ is therefore vital. While it is in principle possible (albeit time-consuming) to recompute this matrix from the geometry of the pupil array that can be saved in the OI_ARRAY extension of an OIFITS file, we strongly prefer to save this matrix directly in one of the FITS file extensions, not least because handling hundreds of individual subapertures becomes impractical using the OI_ARRAY extension. Second, recent work from Martinache et al. (2020) has shown that the use of "gray" pupil models featuring subapertures with a continuous transmission 0 < T ≤ 1 is highly desirable in KPI as they help to reduce the systematic kernel phase signal of unresolved and point-like PSF reference targets. However, gray pupil models are not supported by the OIFITS file format due to the third dimension required to store transmission. Finally, the Karhunen-Loève calibration method applied in Kammerer et al. (2019); Wallace et al. (2020), and Kammerer et al. (2021) can be used implicitly if the kernel phase θ and the kernel matrix $\tilde{{\boldsymbol{K}}}$ are saved in matrix format. This is because the Karhunen-Loève calibration (Soummer et al. 2012) is a linear projection represented by a matrix P so that replacing the kernel phase and the kernel matrix with P · θ and ${\boldsymbol{P}}\cdot \tilde{{\boldsymbol{K}}}$, respectively, will ensure that any model fitting code will automatically run correctly with a calibrated data set.

Based on the aforementioned issues with the OIFITS file format, the kernel phase community has expressed the need for a dedicated data exchange format in the framework of the Masking/Kernel Phase Hackathon ³⁰ held virtually in mid 2021 to ensure compatibility among a variety of data reduction, calibration, and model fitting tools. The kernel phase FITS (KPFITS) file format proposed here is based on the XARA file format developed by Martinache (2010) and Martinache (2013). However, based on recent developments from Kammerer et al. (2019) and Martinache et al. (2020) several modifications were made to the original XARA file format in consultation with kernel phase experts from the global community. For completeness, we also added an optional extension for saving kernel amplitude observables as introduced in Pope (2016). The structure of the KPFITS file format is described in Table 1 and shall be used as a standard for saving and exchanging kernel phase data in the future. The exact order of the FITS file extensions is in principle irrelevant as extensions shall be referred to in the code with their extension names, the numbering in Table 1 can hence be regarded a suggestion. We note that in principle, AMI data can also be exchanged using the KPFITS file format in a more efficient way.

Table 1. Structure of the KPFITS File Format

Ind	Opt	Name	Type	Dimensions	Description
0	(N)	PRIMARY	Image	Nf × Nλ × Npix × Npix	Telescope images (can be empty array)
1	N	APERTURE	Table	Nsap × 3	Description of pupil model
	N	XXC	Column		Subaperture x-coordinate [m]
	N	YYC	Column		Subaperture y-coordinate [m]
	N	TRM	Column		Subaperture transmission (0 < T ≤ 1)
2	N	UV-PLANE	Table	Nuv × 3	Fourier plane coverage of pupil model
	N	UUC	Column		Fourier u-coordinate [m]
	N	VVC	Column		Fourier v-coordinate [m]
	N	RED	Column		Redundancy of uv-position (integer)
3	N	KER-MAT	Image	${N}_{\ker }\times {N}_{\mathrm{uv}}$	Matrix $\tilde{{\boldsymbol{K}}}$ mapping image Fourier phase onto kernel phase
4	N	BLM-MAT	Image	Nuv × Nsap	Matrix A mapping pupil plane phase onto image Fourier phase
5	N	KP-DATA	Image	${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }$	Kernel phase data [rad]
6	N	KP-SIGM	Image	${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }$	Kernel phase uncertainties [rad]
7	N	CWAVEL	Table	Nλ × 2	Description of bandpass
	N	CWAVEL	Column		Central wavelength of bandpass [m]
	N	BWIDTH	Column		Best available estimate of effective half-power bandwidth [m]
8	N	DETPA	Image	Nf	Detector position angle E of N [deg]
9	N	CVIS-DATA	Image	2 × Nf × Nλ × Nuv	Complex visibility data (dim 1 = real part, dim 2 = imag. part)
>9	Y	KA-DATA	Image	${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }$	Kernel amplitude data
>9	Y	KA-SIGM	Image	${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }$	Kernel amplitude uncertainties
>9	Y	CAL-MAT	Image	$({N}_{\ker }-{K}_{\mathrm{klip}})\times {N}_{\ker }$	Karhunen-Loève projection matrix ${\boldsymbol{P}}^{\prime}$ as in K19
>9	Y	KP-COV	Image	Flexible, up to ${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }\times {N}_{\ker }$	Kernel phase covariance [rad2]
>9	Y	KA-COV	Image	Flexible, up to ${N}_{{\rm{f}}}\times {N}_{\lambda }\times {N}_{\ker }\times {N}_{\ker }$	Kernel amplitude covariance
>9	Y	FULL-COV	Image	Flexible, up to ${N}_{{\rm{f}}}\times {N}_{\lambda }\times 2{N}_{\ker }\times 2{N}_{\ker }$	Kernel phase [rad 2] and kernel amplitude covariance
>9	Y	IMSHIFT	Table	Nf × 2	Shift to recenter images
	Y	XSHIFT	Column		Shift along x-axis (1-axis in Python) [pix]
	Y	YSHIFT	Column		Shift along y-axis (0-axis in Python) [pix]
>9	Y	WINMASK	Image	Npix × Npix	Super-Gaussian windowing mask

Note. "Ind" denotes the index of the FITS file extension, "Opt" specifies whether an extension is optional or not, N_f denotes the number of frames, N_λ denotes the number of spectral channels, N_pix denotes the number of pixels per image axis, N_sap denotes the number of subapertures, N_uv denotes the number of uv-points and ${N}_{\ker }$ denotes the number of kernels. We note that the telescope images are optional and the PRIMARY extension can in principle be an empty array. K19 = Kammerer et al. (2019). All data should be of type float.

Download table as:  ASCIITypeset image

To ensure that the information necessary to consistently reprocess the kernel phase data is present, the KPFITS file format also requires a PRIMARY FITS header with the header keywords specified in Table 2. While several of these keywords are optional, others are strictly required to rerun the Kpi3Pipeline and extract kernel phase observables or to aid the model fitting procedures described in Section 2.3.

There are several publicly available model fitting toolkits that understand OIFITS files and can be used to analyze long-baseline interferometry and AMI data (e.g., LITpro, Tallon-Bosc et al. 2008, CANDID, Gallenne et al. 2015, PMOIRED, Mérand 2022). However, none of these toolkits accounts for correlated uncertainties in the fits. Recently, Lachaume et al. (2019) and Kammerer et al. (2020) developed methods to extract and model correlations in long-baseline interferometry data and Kammerer et al. (2020) showed that accounting for them in the fits can improve VLTI/GRAVITY companion detection limits by a factor of up to ∼2. More importantly, they also found that widely used detection criteria based on χ²-statistics are only valid when accounting for correlations and yield an excessive number of false positive detections otherwise.

In light of these findings, we develop the fouriever toolkit which provides a single solution for analyzing KPI, AMI, and long-baseline interferometry data while accounting for correlated uncertainties in the fits. The toolkit also enables modeling correlations in long-baseline interferometry and AMI data as described in Kammerer et al. (2020) and calibrating science data using a Karhunen-Loève projection based on calibrator data as introduced in Kammerer et al. (2019). The fouriever toolkit is written in Python and can be obtained from GitHub. ³¹ Many functionalities were inspired by CANDID but have been modified to improve the performance given that fits accounting for correlations involve significantly more complex matrix multiplications. In the following, we briefly outline the search for companions and the estimation of detection limits with fouriever and highlight similarities and improvements with respect to CANDID. We also note that fouriever comes with tutorials and test data for JWST NIRISS AMI, VLT NACO and Keck NIRC2 KPI, and VLTI PIONIER and VLTI GRAVITY long-baseline interferometry applications.

2.3.1. Companion Search

The first step in the analysis of high-contrast imaging data usually consists of the search for one or multiple companions. For this purpose, fouriever can compute a χ²-detection map based on binary model fits similar to the fitMap in CANDID but accounting for correlated uncertainties in the fits. This χ²-detection map is obtained from a set of least squares gradient descent minimizations initialized on a grid around the science target. The optimizer aims to minimize

$\begin{eqnarray}&&{\chi }^{2}={R}^{T}\cdot {{\boldsymbol{\Sigma }}}^{-1}\cdot R,\end{eqnarray} \tag{ 15 }$

where R = D − M are the residuals between data and model and Σ is the data covariance matrix. For simplicity, other toolkits assume that there are no correlations and that the matrix Σ is diagonal. The model observables are obtained from the complex visibility of a binary source

$\begin{eqnarray}&&{V}_{\mathrm{bin}}=\displaystyle \frac{{V}_{1}+{V}_{2}f\exp \left(-2\pi i\left(\tfrac{{{\rm{\Delta }}}_{\mathrm{RA}}u}{\lambda }+\tfrac{{{\rm{\Delta }}}_{\mathrm{DEC}}v}{\lambda }\right)\right)}{1+f},\end{eqnarray} \tag{ 16 }$

where V₁ and V₂ are the complex visibilities of the primary and the secondary, f is the relative flux of the secondary with respect to the primary, Δ_RA and Δ_DEC are the R.A. and decl. offset between the primary and the secondary, u and v are the Fourier plane coordinates for which the complex visibility shall be obtained, and λ is the observing wavelength (e.g., Berger 2003). For the case of an unresolved companion considered here, we set

$\begin{eqnarray}&&{V}_{1}=2\displaystyle \frac{{J}_{1}(\pi \vartheta b)}{\pi \vartheta b},\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}&&{V}_{2}=1,\end{eqnarray} \tag{ 18 }$

where J₁ denotes the first-order Bessel function of the first kind and $b=\sqrt{{u}^{2}+{v}^{2}}$ denotes the length of the baseline for which the complex visibility shall be obtained, so that the primary is described by a uniform disk with angular diameter ϑ and the secondary is described by an unresolved point-source (e.g., Berger 2003). The model observables M in the form of squared visibility amplitudes (v2), closure phases (cp), or kernel phases (kp) are then obtained via

$\begin{eqnarray}&&{\rm{v}}2={\left|V\right|}^{2},\end{eqnarray} \tag{ 19 }$

$\begin{eqnarray}&&\mathrm{cp}={\boldsymbol{C}}\cdot \angle V,\end{eqnarray} \tag{ 20 }$

$\begin{eqnarray}&&\mathrm{kp}=\tilde{{\boldsymbol{K}}}\cdot \angle V,\end{eqnarray} \tag{ 21 }$

where ∠ denotes the phase of a complex number and C and $\tilde{{\boldsymbol{K}}}$ denote the closure and kernel matrices, respectively, whose rows contain the linear combinations of Fourier phases forming closure and kernel phases. The detection significance is obtained using χ²-statistics and computed via

$\begin{eqnarray}&&{P}_{\mathrm{bin}}=1-{\mathrm{CDF}}_{\nu }\left(\displaystyle \frac{\nu {\chi }_{{\rm{r}},\mathrm{ud}}^{2}}{{\chi }_{{\rm{r}},\mathrm{bin}}^{2}}\right),\end{eqnarray} \tag{ 22 }$

where CDF_ν denotes the cumulative distribution function of a χ²-distribution with ν degrees of freedom, ${\chi }_{{\rm{r}},\mathrm{ud}}^{2}$ is the reduced χ² of the uniform disk only model (i.e., without a companion), and ${\chi }_{{\rm{r}},\mathrm{bin}}^{2}$ is the reduced χ² of the binary model. As in CANDID, fouriever also supports numerical bandwidth smearing which helps to prevent underestimating the companion flux. Moreover, an on-sky search region can be specified where fouriever is looking for companions. This feature is particularly useful if an already known companion resides outside the diffraction field-of-view (FoV) of $0.5{\lambda }_{\min }/{B}_{\min }$ and creates aliasing artifacts that could be mistaken for a true companion, where ${\lambda }_{\min }$ is the shortest observing wavelength and ${B}_{\min }$ is the smallest baseline of the pupil or interferometric array. Searching for additional companions is possible by repeating the computation of the χ²-detection map after analytically subtracting the best fit companion model from the data.

The uncertainties derived from least squares gradient descent minimizations are usually unreliable if the uncertainties in the underlying data have been wrongly estimated. While CANDID employs the bootstrapping (with replacement) method (Efron & Tibshirani 1986) to extract the uncertainties and correlations of the model parameters, fouriever employs the MCMC method (using emcee, Foreman-Mackey et al. 2013) with a temperature ${T}_{{\chi }^{2}}$ by default equaling the reduced χ² of the best fit binary model obtained from the χ²-detection map to account for potential over- or underestimation of the uncertainties in the underlying data (e.g., Andrae 2010). The log-likelihood function that is being sampled by the MCMC method is thus given by

$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=-\displaystyle \frac{1}{2}\displaystyle \frac{{\chi }^{2}}{{T}_{{\chi }^{2}}}.\end{eqnarray} \tag{ 23 }$

Both the bootstrapping and the MCMC method yield comparable results and an advantage of the latter is that it can also be used to fit more complex models such as multiple companions simultaneously or geometric disk and ring models to the data (see, e.g., Kammerer et al. 2021; Blakely et al.2022).

Figure 13 shows a companion search in VLTI/PIONIER data of AX Cir from Gallenne et al. (2015) with fouriever and CANDID. For a more direct comparison, we did assume uncorrelated uncertainties and numerical bandwidth smearing evaluated at three uniformly spaced wavelength nodes across the observing bandpass in both cases. The recovered model parameters (host star uniform disk diameter ϑ/diam*, relative companion flux f, R.A. offset of the companion Δ_RA/x, and decl. offset of the companion Δ_DEC/y) from the MCMC fit with fouriever and the bootstrapping with CANDID agree well within their uncertainties. Both approaches also reveal a correlation in the model parameters between the host star uniform disk diameter and the relative companion flux. We observe a similar computation time for fouriever and CANDID, although we note that fouriever natively uses a higher resolution grid than CANDID when computing the χ²-detection map and also achieves a similar performance with correlated uncertainties.

2.3.2. Detection Limits

Estimating companion detection limits (also known as contrast curves) is one of the most fundamental pathways to assess the sensitivity of high-contrast imaging observations and to quantitatively interpret non-detections. For this purpose, fouriever offers the same two methods to estimate these limits as CANDID, namely the "Absil" and the "Injection" method.

"Absil" method: this method has been introduced by Absil et al. (2011) and refined by Gallenne et al. (2015) and computes the detection limit as the relative companion flux f at which the binary model deviates by a certain probability (e.g., 99.73% for 3σ) from the uniform disk only model according to

$\begin{eqnarray}&&{P}_{\det }=1-{\mathrm{CDF}}_{\nu }\left(\displaystyle \frac{\nu {\chi }_{{\rm{r}},\mathrm{bin}}^{2}}{{\chi }_{{\rm{r}},\mathrm{ud}}^{2}}\right),\end{eqnarray} \tag{ 24 }$

that is assuming that the binary model is the true model (see Gallenne et al. 2015 for a more detailed discussion of this choice). This condition is evaluated on a grid around the science target and azimuthally averaged to obtain a contrast curve.

"Injection" method: this method has been introduced by Gallenne et al. (2015) and has been reported to be more robust than the "Absil" method if the data is biased or affected by correlations. It consists of analytically injecting a fake companion into the data and then computing the significance of the best fit uniform disk only model over the best fit binary model fitted to the synthetic data with the injected companion. In this case, P_bin (Equation (22)) sets the detection limit in the form of relative companion flux f. This method is computationally more expensive than the "Absil" method but it resembles more closely the procedure that is undertaken if a real companion is detected.

The fouriever toolkit can compute these detection limits accounting for correlated uncertainties in the fits and uses a more aggressive smoothing when azimuthally averaging the detection maps because the sparse uv-coverage especially in long-baseline interferometry observations often results in detection maps showing strong changes in sensitivity at high spatial frequencies (aliasing). This can be seen in the top panels of Figure 14 which also compares the detection limits of the VLTI/PIONIER data of AX Cir from Gallenne et al. (2015) obtained with fouriever and CANDID. For better comparability, we apply the same azimuthal averaging that is being used in fouriever when computing the CANDID detection limits shown in Figure 14. To illustrate how accounting for correlated uncertainties in the fits improves the detection limits, we also model the correlations among the closure phase observables of the AX Cir data considering that two telescope triplets of closing triangles at the VLTI always share one of their three baselines and are therefore not mathematically independent (e.g., Monnier 2007). This results in the correlation models from Section 2.2 of Kammerer et al. (2020) with x = y = 0. Then, we repeat the computation of the detection limits with fouriever using these correlation models in the fits. For the AX Cir data with only four closing triangles and three spectral channels, the correlations are small and the detection limits improve by only ∼5% (Figure 14).

When inspecting the jwst pipeline-processed commissioning data by eye it is already possible to identify a low-contrast close-in companion candidate around CPD-66 562. This observation is supported by the raw kernel phase observables extracted with the Kpi3Pipeline which show a much larger scatter (±0.152 rad) for CPD-66 562 than for the other three targets (±0.007 rad, Figure 8). We note that the theoretically expected kernel phase signal of a centro-symmetric point-source is zero so that larger scatter in the kernel phase observables is always indicative of a detection (e.g., Martinache 2010).

Zoom In Zoom Out Reset image size

A companion search with fouriever reveals the precise parameters of the detected companion candidate. Before, though, we calibrate the raw kernel phase observables of CPD-66 562 using the other three targets as point-source references. Here, we consider the data from both dither positions of the first run on 2022 May 23. Using the Karhunen-Loève calibration class (klcal) in fouriever, the Karhunen-Loève basis of the kernel phase observables of the three reference targets is computed and clipped to K_klip = 3 principal components (Soummer et al. 2012; Kammerer et al. 2019). Then, the kernel phase observables of CPD-66 562 are projected into an ${N}_{\ker }-3$-dimensional subspace that is orthogonal to (and thus independent of) these first three principal components of the reference target kernel phase observables. Then, the location of the companion candidate is derived by computing a χ²-detection map from the calibrated kernel phase observables of CPD-66 562 as described in Section 2.3.1. Finally, the companion parameters including their uncertainties are estimated using an MCMC fit initialized around the best fit position from the χ²-detection map (Figure 9).

Zoom In Zoom Out Reset image size

The χ²-detection map shows a highly significant detection reaching the numerically set limit of N_σ = 8.0. The companion candidate is fairly bright with a relative flux of ∼20% of the primary and separated by ∼151 mas which corresponds to only ∼1 λ/D at λ = 4.8 μm. The small relative uncertainties in the flux of <1% and in the position of <1 permille obtained from the MCMC fit demonstrate the high precision at which the kernel phase technique can resolve companions at the diffraction limit.

To evaluate the performance of JWST NIRISS full pupil KPI, we compute 5σ companion detection limits using the Absil method in fouriever. For this purpose, we only consider the three point-source reference targets that were reobserved on 2022 June 5 with successful target acquisition achieving a precision of <0.1 pixels (Rigby et al. 2022). This also means that we exclude the low-contrast binary CPD-66 562 reported in Section 4.1.

First, we investigate the raw kernel phase signal of the three point-source reference targets from the second run averaged over both dither positions (Figure 10). The average observed scatter is ±0.008 rad and thus comparable, albeit slightly larger, to the first run. This small difference can stem from unflagged (and thus uncorrected) bad pixels in the data from the second run (e.g., from cosmic rays) or the temporally varying wave front quality of JWST (due to e.g., wing tilt events, Rigby et al. 2022). Since kernel phase observables are often dominated by systematic errors from higher order phase aberrations in the system (e.g., Kammerer et al. 2019; Martinache et al. 2020), it is common practice to use reference targets for calibrating out these systematic errors. The stability of the system can then be assessed by considering the difference between kernel phase observables from different reference targets. For this purpose, the right panel of Figure 10 shows the kernel phase residuals between each pair of targets observed in the second run (again averaged over both dither positions) and reveals that the signal measured on 2MASS J062802.01-663738.0 seems to be an outlier.

Zoom In Zoom Out Reset image size

Since two slightly different dither positions were used for the two exposures on each target, we repeat the same analysis for each dither position separately. Figure 15 shows 2MASSJ062802.01-663738.0 as an outlier in both dither positions. A companion search with fouriever using the other two targets as calibrators reveals that this outlier is roughly consistent with a source at ∼240 mas separation and ∼0.6% flux ratio (Figure 17). At a contrast of ∼5.6 mag, this companion candidate is above the 5σ detection threshold measured for the other two targets. To exclude the possibility that this detection is caused by wave front drift between the observations of the different targets, we also analyze the data from the first run for an additional companion candidate around 2MASS J062802.01-663738.0. Figure 16 compares the kernel phase residuals between each pair of targets from the first run and also reveals 2MASS J062802.01-663738.0 as an outlier in both dither positions. A companion search using only the data from the first run also yields a detection toward the South-East of 2MASS J062802.01-663738.0, albeit at significantly different separation and flux ratio.

As a final check, we use classical PSF subtraction methods to search for the companion candidate in the image plane using only the data from the second run. We build a PSF library with the images of the other two targets (TYC 8906-1660-1 and CPD-67 607) and use the publicly available pyKLIP package (Wang et al. 2015) to subtract the first 20 Karhunen-Loève modes of the PSF library from the science target (2MASS J062802.01-663738.0). The bottom panel of Figure 17 shows the PSF-subtracted image and clearly reveals a point-source at a position that is consistent with the best fit from the kernel phase reduction. Hence, the detection around 2MASS J062802.01-663738.0 appears to be real and could be either a companion or a background source. The agreement between the kernel phase and the classical PSF subtraction techniques is a great confirmation of our methodology and a more detailed comparison between Fourier plane and classical imaging techniques is planned for a future publication.

Finally, comparing the four dither positions tried during instrument commissioning suggests that the first dither position from the second run should be avoided since it yields a ∼3 times larger scatter in the raw kernel phase and a ∼2 times larger scatter in the kernel phase residuals. We note that the rms of the kernel phase uncertainty estimated from the standard deviation over the data cubes is on the order of ∼0.0034 rad for all targets and thus slightly larger than the calibration residuals between the best two targets, meaning that with 2e8 collected photons the observations are not yet limited by systematic calibration errors.

Next, we compute companion detection limits for each of the three point-source reference targets from the second run. For each of the three targets, we use the other two targets as references and apply the Karhunen-Loève calibration in fouriever (Soummer et al. 2012; Kammerer et al. 2019) with K_klip = 2 to calibrate out systematic errors. Then, we use the Absil method in fouriever to compute companion detection limits from the calibrated kernel phase observables of each target. For 2MASS J062802.01-663738.0, we analytically remove the best fit companion before computing the detection limits. Figure 11 shows the detection limits obtained from this procedure. The KPI detection limits show a prominent bump just inside of 300 mas separation which is caused by increased photon noise from the first Airy ring of the PSF. This is highlighted by the gray shaded background showing the azimuthal average of a NIRISS CLEARP F480M PSF computed with WebbPSF ³⁴ (Perrin et al. 2012) in a logarithmic color stretch. For the best target, the detection limits reach ∼6.5 mag at ∼200 mas and ∼7 mag at ∼400 mas. These limits agree well with the expectation from the kernel phase calibration residuals of σ_KP ∼ 0.002 rad between TYC 8906-1660-1 and CPD-67 607 as shown in Figure 10 and translating into a contrast limit of ∼6.75 mag. Compared to the theoretical KPI detection limits predicted by Ceau et al. (2019) using their binary test T_B (which is similar to our detection criterion), our on-sky limits are ∼1 mag worse. This is consistent with the ∼1 mag difference between the theoretical AMI detection limits from Ireland (2013) and the limits measured on-sky (see Figure 12). We expect that some fraction of this difference is caused by charge migration which has not been accounted for in the simulations and theoretical predictions. However, a more detailed analysis of the systematic errors is still ongoing.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For comparison, Figure 12 shows the detection limits achieved with the NIRISS AMI F480M commissioning observations of AB Dor together with our KPI detection limits. Both the KPI and the AMI observations were designed to collect 2e8 photons on the detector and use the same F480M filter. The AMI data (observations 12 and 13 of PID 1093) were reduced with the jwst stage 1 and 2 pipelines and closure phases and visibility amplitudes were extracted from the cleaned interferograms using a custom version of ImPlaneIA ³⁵ (Greenbaum et al. 2015). Next, the observables of AB Dor were calibrated against those of the two observed reference targets (HD 37093, observations 15 and 16, and HD 36805, observations 18 and 19) using median subtraction. Then, the best fit parameters of the known close-in companion AB Dor C (separation ∼326 mas) were obtained with fouriever as discussed in Section 2.3.1 using both closure phases and visibility amplitudes in the fit. Finally, the best fit companion was analytically subtracted from the data before computing the detection limits using the Absil method in fouriever. The NIRISS AMI commissioning data reduction and analysis is described in more detail in Sivaramakrishnan et al. (2023).

Figure 12 shows that AMI (light blue curve) reaches deeper contrasts than KPI, especially at small angular separations ≲400 mas. This is the case since the NRM only collects photons on non-redundant baselines whereas the vast majority of the photons collected in KPI with the full pupil are affected by redundancy noise, so that the AMI observations achieve better detection limits than the KPI observations with the same number of collected photons. For a fair comparison, it needs to be considered that the NRM needs 5.6 times more time to collect the same number of photons than the full pupil since the NRM has an optical throughput of only 15% while the CLEARP full pupil mask has an optical throughput of 84%. Assuming that the contrast scales with the square root of the number of collected photons (e.g., Ireland 2013) and multiplying the AMI detection limits by $\sqrt{0.84/0.15}$ results in more comparable limits between AMI (dark blue curve) and KPI. While the scaled AMI contrast curve still achieves better limits at ∼250–325 mas where KPI is suffering from increased photon noise from the first Airy ring of the PSF, KPI now achieves better limits beyond ∼325 mas where AMI is limited by the reduced throughput and uv-sampling. Finally, we note that the unscaled AMI detection limits (light blue curve) are still ∼1 mag above the fundamental noise floor according to Ireland (2013) and future efforts will aim at further improving the performance of NIRISS AMI and KPI (Sivaramakrishnan et al. 2023).