# Probing Photon Statistics in Adaptive Optics Images with SCExAO/MEC

Sarah Steiger, Timothy D. Brandt, Olivier Guyon, Noah Swimmer, Alexander

B. Walter, Clinton Bockstiegel, Julien Lozi, Vincent Deo, Sebastien Vievard, Nour Skaf, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Sarah Steiger, Timothy D. Brandt, Olivier Guyon, Noah Swimmer, Alexander B. Walter, et al. Probing Photon Statistics in Adaptive Optics Images with SCExAO/MEC. The Astronomical Journal, 2022, 164, 10.3847/1538-3881/ac922f . insu-03874870

HAL Id: insu-03874870

https://insu.hal.science/insu-03874870

Submitted on 28 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Probing Photon Statistics in Adaptive Optics Images with SCExAO/MEC*

Sarah Steiger<fnref n="1" /> ID, Timothy D. Brandt<fnref n="1" /> ID, Olivier Guyon<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> ID, Noah Swimmer<fnref n="1" /> ID, Alexander B. Walter<fnref n="6" />, Clinton Bockstiegel<fnref n="7" />, Julien Lozi<fnref n="2" /> ID, Vincent Deo<fnref n="2" /> ID, Sebastien Vievard<fnref n="2" /> ID, Nour Skaf<fnref n="2" />,<fnref n="8" />,<fnref n="9" /> ID, Kyohoon Ahn<fnref n="2" /> ID, Nemanja Jovanović<fnref n="10" /> ID, Frantz Martinache<fnref n="11" /> ID, and Benjamin A. Mazin<fnref n="1" /> ID

## Abstract

We present an experimental study of photon statistics for high-contrast imaging with the Microwave Kinetic Inductance Detector (MKID) Exoplanet Camera (MEC) located behind the Subaru Coronagraphic Extreme Adaptive Optics System at the Subaru Telescope. We show that MEC measures the expected distributions for both on-axis companion intensity and off-axis intensity, which manifests as quasi-static speckles in the image plane and currently limits high-contrast imaging performance. These statistics can be probed by any MEC observation due to the photon-counting capabilities of MKID detectors. Photon arrival time statistics can also be used to directly distinguish companions from speckles using a post-processing technique called Stochastic Speckle Discrimination (SSD). Here, we give an overview of the SSD technique and highlight the first demonstration of SSD on an extended source—the protoplanetary disk AB Aurigae. We then present simulations that provide an in-depth exploration as to the current limitations of an extension of the SSD technique called Photon-Counting SSD (PCSSD) to provide a path forward for transitioning PCSSD from simulations to on-sky results. We end with a discussion of how to further improve the efficacy of such arrival-time-based post-processing techniques applicable to both MKIDs, as well as other high-speed astronomical cameras.

*Unified Astronomy Thesaurus concepts: Adaptive optics (2281); Exoplanet detection methods (489); Direct imaging (387); Astronomical detectors (84); Astronomical instrumentation (799); Astronomy data analysis (1858)*

## 1. Introduction

Direct imaging is an extremely technologically challenging exoplanet detection technique with necessary contrasts easily exceeding 10⁻<fnref n="6" /> for even the largest self-luminous planets with orbits wider than that of Saturn (>10 au). For this reason, nearly all of the ~10–20 directly imaged planets have separations of 10–150 au, typically ρ ~ 0''4–2'' on the sky (e.g., Marois et al. 2008b; Lagrange et al. 2009; Kuzuhara et al. 2013; Rameau et al. 2013; Currie et al. 2014; Macintosh et al. 2015; Chauvin et al. 2017). In order to directly image a new regime of planets at Jupiter-to-Saturn-like separations, new technologies and techniques need to be developed to push for higher contrasts at smaller inner working angles (IWA).

The current limitation for ground-based direct imaging is point-spread function (PSF)-sized “speckles”, which are caused by uncontrolled diffracted and scattered starlight and have a range of correlation timescales (τ) and sources. Rapidly

evolving atmospheric speckles (τ ~ 1–20 ms) result from aberrations left uncorrected by an adaptive optics (AO) system and average out over the course of long-exposure images, forming a smooth halo (e.g., Perrin et al. 2003; Soummer & Ferrari 2007). These “fast” speckles are easier to remove and can be corrected by improved AO control loops, which will mitigate the temporal bandwidth error and measurement (photon noise) error (e.g., Guyon 2005). Alternatively, quasi-static speckles result from instrument imperfections such as noncommon path errors, and telescope vibrations (Guyon 2005; Lozi et al. 2018). These speckles do not average out over long exposures, can interfere with atmospheric speckles, and be pinned to the diffraction rings, easily masquerading as point sources in an image (Soummer & Ferrari 2007). This quasi-static speckle noise is temporally well correlated (τ ~ 10–60 minutes), presenting a fundamental obstacle in exoplanet direct imaging (e.g., Marois et al. 2008a).

While focal-plane wave-front control methods offer a pathway to suppress these speckles on-sky (e.g., Give'on et al. 2007; Martinache et al. 2016), they are most commonly removed in post-processing. Unfortunately, the most common post-processing techniques that utilize advanced PSF subtraction methods (e.g., Lafrenière et al. 2007; Soummer et al. 2012) become less effective at small IWAs where direct detections are most challenging.

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹Department of Physics, University of California, Santa Barbara, Santa Barbara, CA, USA; steiger@physics.ucsb.edu</footnote>

<footnote n="2">²Subaru Telescope, National Astronomical Observatory of Japan, 650 North A'ohokū Place, Hilo, HI 96720, USA</footnote>

<footnote n="3">³Steward Observatory, The University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="4">⁴College of Optical Sciences, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="5">⁵Astrobiology Center, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan</footnote>

<footnote n="6">⁶Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷IQM, Keilaranta 19, FI-02150 Espoo, Finland</footnote>

<footnote n="8">⁸LESIA, Observatoire de Paris, Univ. PSL, CNRS, Sorbonne Univ., Univ. de Paris, 5 pl. Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="9">⁹Department of Physics and Astronomy, University College London, London, UK</footnote>

<footnote n="10">¹⁰Department of Astronomy, California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="11">¹¹Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France</footnote>

<footnote>Received 2022 August 11; revised 2022 September 11; accepted 2022 September 13; published 2022 October 12</footnote>

<footnote>* Based in part on data collected by Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.</footnote>
Angular differential imaging (ADI; Marois et al. 2006) uses the rotation of the Earth or, analogously, the field-of-view rotation of an altitude–azimuth telescope with the field derotator turned off, to separate speckles (which will remain stationary) from planetary companions (which will rotate with the sky). The arc traced by an object across the image, however, scales proportionally with angular separation for a given unit time, resulting in less rotation at smaller IWAs. Additionally, the rotation in λ/D units is smaller within a few diffraction beamwidths, resulting in severe self-subtraction and partial subtraction of a planet signal at these small separations (Lafrenière et al. 2007; Mawet et al. 2012). Spectral differential imaging (SDI; Marois et al. 2000) utilizes differences in the wavelength-dependent scaling of the diffracted light speckles and companions in polychromatic images, but requires broad spectral coverage close to the primary; otherwise it also suffers from self-subtraction effects. Reference star differential Imaging (RDI/RSDI; Soummer et al. 2012) does not inherently suffer at small IWAs, but requires careful matching between the target of interest and the reference star for optimal performance (Ruane et al. 2019). A method to suppress quasi-static speckles that is free of the limitations of ADI, SDI, and RDI would significantly improve our ability to detect Jovian planets at Jupiter-to-Saturn-like separations.

Post-processing techniques that take advantage of photon arrival time statistics, such as stochastic speckle discrimination (SSD; Gladysz et al. 2008), currently show much promise as they use differences in the intensity distributions of speckles and companions at millisecond timescales to distinguish the two populations. As this technique uses only temporal information in the form of the “instantaneous” intensity, it requires no PSF reference subtraction and is free of the spectral, spatial, and reference matching issues of SDI, ADI, and RDI, respectively. Effectively utilizing this information however can be challenging as you need many photons to arrive within a single speckle decorrelation time ($\tau \sim 0.1$ s). Additionally, cameras with a frame rate faster than this timescale and very low read noise at these high speeds are needed.

Microwave kinetic inductance detector (MKID) instruments are therefore very appealing for these techniques as they are photon-counting instruments with µs precision and no read noise or dark current (see Mazin et al. 2012; Zypryt et al. 2017 for more details on MKIDs). Specifically, the MKID Exoplanet Camera (MEC; Walter et al. 2019) is a y–J band MKID IFU located behind the Subaru Coronagraphic Extreme Adaptive Optics System (SCExAO; Jovanovic et al. 2015b) at the Subaru Telescope on Maunakea that has been recently commissioned and has shown its ability to use the temporal resolution afforded by its MKID detector to enable the discovery of low-mass companions using SSD (Steiger et al. 2021).

In this work, we demonstrate the ability of MKIDs to accurately probe the distinct intensity distributions of quasi-static speckles and companions using data from SCExAO/MEC. We also highlight more generally how others can use MEC to access photon arrival time information and use the techniques discussed here, or to develop new time-based analysis techniques in the future. We then give an overview of the status of current photon statistical post-processing techniques for high-contrast imaging with on-sky MEC SSD results shown for both point sources and diffuse sources. Next we describe a new post-processing technique first shown in Walter

et al. (2019) called Photon-Counting SSD (PCSSD). Using both simulations and on-sky data, we identify key features that have prevented this technique from achieving desired performance on-sky and the steps that can be taken to improve performance in the future. Finally we end by discussing key paths forward toward developing new and more effective photon arrival-time-based post-processing techniques for exoplanet direct imaging.

## 2. Photon Statistics in Millisecond AO Images

Before discussing how to utilize photon arrival time information with MEC and photon-arrival-time-based post-processing techniques, we present an overview of the statistics that govern the different parts of an image downstream of an AO system at millisecond frame rates.

### 2.1. Speckles

Originally derived by Goodman (1975) and experimentally verified by Cagigal & Canales (2001) and Fitzgerald & Graham (2006), the underlying probability density function (PDF) that estimates the intensity distribution of off-axis stellar speckles in the image plane is given by a modified Rician (MR)

$$p_{\text{MR}}(I) = \frac{1}{I_s} \exp\left(-\frac{I + I_c}{I_s}\right) I_0\left(\frac{2\sqrt{I_c}}{I_s}\right), \quad (1)$$

where $I_0(x)$ denotes the zero-order modified Bessel function of the first kind, $I_c$ describes the coherent intensity component attributed to the unaberrated PSF of the primary, and $I_s$ is the time variable component of the total intensity that describes the speckle field (see also Marois et al. 2008a). It is important to note that the shape of this distribution is always positively skewed (i.e., the distribution tail falls to the right-hand side of the mean) and that the skewness of the MR can be parameterized by the ratio of $I_c/I_s$ (the larger the $I_c/I_s$ ratio, the less skewed the distribution).

For a sequence of exposures shorter than the decorrelation time of atmospheric speckles (~10–100 ms), a histogram of the image plane intensity follows an MR, and $I_c$ and $I_s$ can be determined for each pixel in an image (Fitzgerald & Graham 2006). A key feature of MEC is that it stores the arrival time information of every photon, and so all time binning can be done in post-processing. This is important as the bin size that ideally samples the MR is difficult to determine a priori and may also vary spatially across the image.

The optimal bin size should be shorter than the decorrelation timescale of the speckles in the image but long enough that each bin contains many photons. If too large of a bin size is chosen, many realizations of the speckle intensity will be averaged over. Conversely, if too small of a bin size is selected, then not enough photons will arrive per bin and the distribution will tend toward Poisson statistics.

### 2.2. Companions

On-axis sources (nondiffracted light, i.e., astronomical objects in the image) behind an extreme AO (ExAO) system do not follow MR statistics because of the ExAO system itself, which can be thought of as a high-pass spatial filter acting on the phase (Sivaramakrishnan et al. 2001). At the center of an image, the complex amplitude is simply the integral of the pupil complex amplitude as given by Equation (2)
(Soummer & Ferrari 2007).

$$ \Psi(0) = \int P(x) e^{i\phi(x)} dx. \qquad (2) $$

Here $P(x)$ denotes the pupil function, and $\phi(x)$ is the phase of the wave front in the pupil plane.

With no (or low) AO correction, $\phi(x)$ is large and so the phase vectors, $e^{i\phi(x)}$, can take any orientation in the complex plane. Summing a large number of these vectors will produce a random walk resulting in circular Gaussian statistics by the central limit theorem. This adherence to circular Gaussian statistics in phase will result in an MR distribution in intensity and be indistinguishable from off-axis speckles.

At high-correction levels (i.e., behind an ExAO system), $\phi(x)$ is small, and so the vectors are not oriented randomly in the complex plane. Summing them will not produce a random walk, the corresponding distribution is not a circular Gaussian, and the resulting on-axis intensity distribution will not follow an MR. The spatial extent where the transition between the on-axis and off-axis intensity distributions occurs is at $\ll 1\lambda/D$ (Soummer & Ferrari 2007), which can be seen qualitatively in Figure 1—see also Section 2.3. For reference, at Subaru $1\lambda/D$ is equal to 27.6 mas (2.75 MEC pixels) at $1.1\mu$m.

Instead of following the MR, at these high-correction levels the Strehl ratio (SR) distribution (which is proportional to the intensity) instead follows the PDF described by Equation (3), which was derived independently by Soummer & Ferrari (2007) and Gladysz et al. (2008):

$$ p_{SR}(sr) = \frac{p_{\delta^2}(-\ln sr)}{sr}, \qquad (3) $$

where sr is the instantaneous SR, and $p_{\delta^2}$ is given by

$$ p(x; k, \theta, \mu) = \frac{\left(\frac{x-\mu}{\theta}\right)^{k-1} \exp\left(-\frac{x-\mu}{\theta}\right)}{\Gamma(k)\theta}. \qquad (4) $$

Here $k > 0$ is the shape parameter, $\theta > 0$ is the scale parameter, and $\mu$ is the location parameter, which shifts the PDF left and right. $\Gamma(k)$ is a Gamma function. For ease of discussion we will refer to the entirety of Equation (3) as the “Gamma” distribution for the remainder of this work.

In contrast to the MR distribution, this Gamma distribution is negatively skewed (i.e., the distribution tail falls to the left-hand side of the mean). Figure 2 illustrates these differences in skewness for various on and off-axis sources with a SCExAO/MEC observation of $\Theta$ 1 Orionis (the Trapezium Cluster).

## 2.3. Satellite Spots

Astrometric and spectrophotometric calibrations are very difficult when the target star (typically the only reference in the field of view) is obscured by a coronagraph. In order to perform these calibrations, faint copies of the obscured stellar PSF called satellite spots may be intentionally placed into the image plane. For SCExAO, these are generated by placing a “waffle” pattern (two orthogonal sine waves) on the AO system’s deformable mirror (DM).

In SCExAO, the satellite spots are additionally modulated at a rate of 2 kHz by flipping the sign of the two sine waves (equivalent to phase shifting them by $\pi$). This is done to make the speckles incoherent with the underlying background and improve photometric performance (Jovanovic et al. 2015a). If there is a coherent speckle underneath the spots then these two

polarities will not have the same brightness as one will interact constructively and the other destructively. In the regime of these millisecond images, this phase switching is unlikely to affect the statistics at the spot locations as 10–100 s of these transitions are being averaged over in a single time bin. It is then expected that the satellite spots should follow the same on-axis statistics of the primary.

Figure 1 shows intensity histograms moving from the speckle field toward the satellite spot center for a single 25 s MEC observation of HIP 109427. The statistics shift from positively to negatively skewed, showing that the satellite spot statistics are not only distinct from that of the speckle field, but also qualitatively follow the same distribution expected for the primary. Additionally, this transition between the on- and off-axis distributions occurs at the expected location of $\ll 1\lambda/D$ ($<2.75$ MEC pixels).

The satellite spot statistics themselves are important because if the distribution of the satellite spots matches that of the on-axis source, then they can be used to help fit the free parameters of the on-axis Gamma distribution (Equation (3)) when observing with a coronagraph. As the brightness of the satellite spots can be controlled by changing the amplitude of the sine waves on the DM, this would be the simplest way to determine the shape of the companion distribution using on-sky data at high signal-to-noise ratio (SNR). Once the shape of this distribution is known, it could then enable the use of post-processing techniques that try to explicitly separate the companion and speckle PDFs such as PDF deconvolution (see Gladysz & Yaitskova 2010).

## 3. Utilizing Photon Arrival Time Information with MEC

Arrival time information can be easily accessed using the MKID Pipeline<fnref n="12" />, which is open-source and can be downloaded via GitHub (Steiger et al. 2022). As MKID detectors record the arrival time and energy of each incident photon, the format of raw MKID data is a time- and energy-tagged photon list that can be queried using the Photontable class on pixel location, time range, photon wavelength, or any combination thereof. A result of this is that all spectral and temporal binning is done in post-processing, and MEC has no set “exposure time” for its observations. This is especially beneficial for post-processing techniques that leverage differences in arrival time statistics, like those that are discussed in the following section, as many different timescales can be probed from a single observation.

MKIDs temporal resolution limit is determined by the readout speed of the detector ($\sim 1$ MHz) and the firmware-imposed dead time. This dead time is set by material properties of the MKID array and for MEC has a value of 10 $\mu$s. During this time, no additional photons are able to be recorded for that pixel to allow it time to return to its idle state (Fruitwala et al. 2020).

Immediately after this dead time, a pile-up of photon events has been observed, which is likely insignificant for total integrated observations but can cause unintended effects when using photon arrival time information directly. Empirically this effect decays rapidly after 40 $\mu$s (Figure 3), and any work done with MEC that uses arrival time information should take this into account so as to not contaminate results.

<footnote n="12">¹² https://github.com/MazinLab/MKIDPipeline</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Histograms of the binned satellite spot intensity moving from three pixels away from the center of the spot (purple, top) to the center (yellow, bottom). Histograms shift from positive to negative skewness as they approach the center of the satellite spot as is shown by the best-fit MR $I_c/I_s$ ratio (dashed lines). Here a higher value indicates less skew. As these data were taken in the $y+J$ band, the satellite spots are elongated in the image and the plotted intensities were found by adding the intensities for six pixels at the same specified distance from the spot center (three on each side). Insets: satellite spot image where the pixels used to generate the histograms are color coded by their distance from the spot center. Dead pixels were purposefully avoided.</figcaption></figure>

# 4. Photon Arrival-time-based Post-processing Techniques

## 4.1. Stochastic Speckle Discrimination

Stochastic speckle discrimination (SSD) is a post-processing technique first demonstrated by Gladysz & Christou (2008) that relies solely on photon arrival time statistics to distinguish between speckles and faint companions in coronagraphic images.

It is achieved with MEC by fitting an MR to every pixel in an image with a user-specified temporal bin size. While the components of the MR distribution do not inherently describe

the signal from a faint companion, the ratio of the coherent component to the time variable component $-I_c/I_s$—may reveal faint companions from a comparably bright speckle field (Gladysz & Christou 2009; Meeker et al. 2018; Steiger et al. 2021). This is because the addition of light from a companion (whose statistics follow a negatively skewed Gamma distribution) will make the best-fit MR more negatively skewed at that location. This is analogous to increasing the best-fit $I_c$ and results in a larger $I_c/I_s$ ratio at the location of the companion compared to other pixels at the same angular separation from the primary.

SSD was recently shown to be effective on-sky with MEC and facilitated the discovery of a low-mass (~0.25 $M_\odot$) stellar companion to the nearby A star HIP 109427 (tet Peg); see Figure 4 and Steiger et al. (2021) for more details.

### 4.1.1. SSD on AB Aur

Some post-processing techniques, such as ADI, struggle to reveal structures with azimuthal symmetry. SSD does not suffer the same limitations for these sources because the light from the disk will still be on-axis (and thus follow a Gamma distribution) even if it is spread over a region instead of contained within a single PSF.

Here we tested the performance of SSD on extended sources using a SCExAO/MEC y-band observation of AB Aurigae (AB Aur) taken on 2021 October 14. AB Aur is a well-studied system with a known protoplanetary disk that also harbors a recently discovered protoplanet (Boccaletti et al. 2020; Currie et al. 2022; Zhou et al. 2022). AB Aur was observed using SCExAO/MEC for 16 minutes in exceptional seeing conditions of ~0.3. The results from the SSD reduction of the protoplanetary disk surrounding AB Aur can be seen in Figure 5 where inner disk structures are revealed in the $I_c/I_s$ map (right) not seen in the y-band total intensity image (left). This bears a strong resemblance to images taken of this system in polarized intensity by the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument at the VLT (SPHERE; Beuzit et al. 2019) as described in Boccaletti et al. (2020)—see their Figure 4. In contrast with the SPHERE observations, here only the millisecond intensity distributions were used to generate these images with no polarization information or PSF subtraction techniques employed. It is important to note that this result was in a large part facilitated by the exceptional seeing conditions as better seeing leads to a less intense speckle halo and allows the disk's Gamma-distributed intensity to significantly modulate the underlying MR distribution of the speckles.

## 4.2. Photon-counting SSD (PCSSD)

PCSSD is an extension of the SSD formalism where contributions from an incoherent source of constant intensity ($I_p$) are accounted for in addition to $I_c$ and $I_s$, which define the shape of the MR. Given a list of photon interarrival times, the maximum likelihood values of $I_c$, $I_s$, and $I_p$ are determined. As all interphoton arrival times are used in this technique, no binning is done, and it has been shown to perform twice as well as perfect PSF subtraction on simulated data where the companions were modeled as constant, incoherent sources (Walter et al. 2019).

One of the main motivations for expanding the SSD formalism is that the SSD $I_c/I_s$ maps—while helpful for
<figure><img src="image_6.png" /></figure>

**Figure 2.** Single-pixel arrival time histograms taken on-sky with SCExAO/MEC of the unocculted star Θ 1 Orionis (purple, bottom left), a faint secondary companion (teal, upper right), and at a random location in the field with comparable brightness to the faint secondary companion (yellow, bottom right). The best-fit MR was plotted for each distribution and is shown with the dashed line. Note that due to the brightness of the primary, 2 ms time bins were used to generate that intensity histogram instead of the 10 ms time bins used for the secondary and field locations. It is clear to see that the MR adequately describes the field location, but breaks down at the location of the primary. The companion is seen to be a convolution of the primary and field PDFs.

<figure><img src="image_7.png" /></figure>

**Figure 3.** Interarrival time histogram for a single pixel. An excess of interarrival times between 10 and 40 µs can be clearly seen, which are inconsistent with Poisson statistics.

extracting companion astrometry and disk morphology—only quantify the skewness of the fit MR to the data. This output is therefore not easily converted to physically meaningful units, and the images are unhelpful for performing spectroscopy or photometry as would typically be desired to determine key companion properties such as temperature, composition, and mass. PCSSD attempts to solve for this by calculating a likelihood for each interphoton arrival time so that the output units can be reported in counts or counts s⁻¹. Additionally, PCSSD is able to leverage the photon-counting nature of MKID detectors to beat the long-exposure noise limit by not temporally binning.

### 4.2.1. PCSSD on HIP 109427 B

In the form described by Walter et al. (2019), PCSSD makes the following assumptions:

1. The speckle halo intensity is entirely described by the MR distribution.

2. $I_C$, $I_S$, and $I_P$ (the intensity of a companion) remain constant over an observation.

3. Chromaticity is ignored.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Left: Total intensity image of HIP 109427 B taken with SCExAO/MEC at the *y* and *J* bands where the location of the companion has been circled in red. Right: SSD *I*<sub>C</sub>/*I*<sub>S</sub> image of HIP 109427 B. Here the companion is plainly visible as well as dark regions at the edge of the coronagraph showing the removal of pinned speckles from the total intensity image. This figure is reproduced from Steiger et al. (2021).</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Left: Total *y*-band intensity dither combined image of AB Aurigae taken with SCExAO/MEC. Disk features here are largely obscured. The satellite spots can be seen as radially smeared bright patches on the edge of this image due to the corrected sky rotation. Right: *I*<sub>C</sub>/*I*<sub>S</sub> image of the disk around AB Aurigae clearly showing some of the inner disk features. These include the two main spirals that have roots to the North and South as well as the extended point-like source to the South at ρ = 0.16 as found in Boccaletti et al. (2020). A Gaussian filter has been applied over this image to smooth over small scale interpixel variations.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

The assumptions that are perhaps the most problematic are that the off-axis intensity is entirely described by the MR distribution and that the companion intensity (*I*<sub>P</sub>) remains constant over the course of an observation. MEC has a known infrared (IR) background that can cause the intensity at any pixel to not be fully described by the MR PDF. Additionally, for realistic observing conditions the on-axis companion intensity (which is proportional to the SR) varies quite considerably and we know is described by the Gamma PDF in Equation (3).

The fact that the on-axis light follows this negatively skewed distribution is in fact one of the main bases that allow traditional SSD to work. Companions stand out in *I*<sub>C</sub>/*I*<sub>S</sub> images due to the addition of the negatively skewed companion PDF,

which shifts the whole distribution or, analogously, increases the best-fit MR *I*<sub>C</sub> (Gladysz & Christou 2008; Meeker et al. 2018; Steiger et al. 2021).

For this reason, using the PCSSD technique on real data does not accurately separate the companion light into *I*<sub>P</sub>, but instead attributes that light largely to *I*<sub>C</sub>. See Figure 6 for an example of PCSSD run on the same HIP 109427 B data set as Figure 4.

4.2.2. PCSSD on More Realistic Simulated Data

To verify that a limitation for performing PCSSD on on-sky data is the assumption that the companion intensity (*I*<sub>P</sub>) is constant, we generated new mock photon lists following the procedure as described in Walter et al. (2019) but with the notable exception that the companion intensities were sampled
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Output of PCSSD on HIP 109427 B (see also Figure 4). The red arrows in each image point to the same location of the low-mass stellar companion. Here the companion light clearly shows up in $I_C$ with corresponding minima in the intensity in the $I_P$ and $I_S$ images.</figcaption></figure>

from the Gamma PDF. The exact distribution from which the companion intensities were sampled can be seen in Figure 7—here sensible values of $k$, $\Theta$, and $\mu$ were chosen to roughly match observed values. These companion photons were also assumed to be correlated in time with a decorrelation timescale ($\tau$) of 0.1 s—the same as that of the MR. A suite of companion separations and contrasts was tested ranging from 3.5 to 12.5 $\lambda/D$ and $5 \times 10^{-5}$ to $4 \times 10^{-4}$, respectively. A suite of brighter companions was also tested (with contrasts ranging from $7.5 \times 10^{-4}$ to $6 \times 10^{-3}$) to more closely replicate the on-sky PCSSD results for binary stellar companions like HIP 109427 B, which has a J-band contrast of $1.27 \times 10^{-3}$.

The results are summarized in Figure 8. The addition of the Gamma-distributed companion flux causes more companion light to be misattributed to $I_C$ over $I_P$ in the case of bright companions ($<1 \times 10^{-3}$ contrast—see right most column) matching observations.

Interestingly, in the case of higher-contrast sources (left column), the shape of the companion distribution does not seem to make much of a difference in the performance of the PCSSD. We believe that in this regime there are not enough companion photons to shift the skewness of the underlying MR distribution at the location of companion pixels and so the light does not incorrectly end up in $I_C$. In other words, at these low count rates it becomes hard to distinguish between the Poisson and Gamma distributions.

At these higher contrasts, however, other factors such as background noise sources are likely to become more significant. The effect of background counts on the performance of the PCSSD is summarized in Figure 9. Here the PCSSD was run on mock photon lists using the suite of higher-contrast companions (see Figure 8, left panels) with an added constant and uniform 50 photons $s^{-1}$ pixel$^{-1}$ background with uncorrelated Poisson distributed arrival times. This count rate approximately matches the current IR background count rate for MEC and significantly degrades PCSSD performance (center column).

We also wanted to explore the effect of performing a wavelength cut on MEC data before running PCSSD. As MEC stores the arrival time and wavelength of each incident photon, a wavelength range can be specified and only the arrival times of photons with wavelengths within that range are returned. This is typically done in part to remove background counts as MEC’s thermal background lies outside of our science bands

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Gamma distribution (Equation (3)) from which the companion intensities were sampled in Section 4.2.2 and Figures 8 and 9. Here $k = 5$, $\Theta = 30$, and $\mu$ corresponds to the median strehl ratio used (0.7).</figcaption></figure>

($y-J$). While these out-of-band photons themselves can be removed from analysis, their effect on temporally proximate photons still remains. Even if it is not used, an out-of-band photon incident on the detector still activates the 10 $\mu$s dead time of that pixel meaning that it is no longer photosensitive over that time range in a way that is unaccounted for by the PCSSD code. For the typical background count rates observed with MEC, however, this effect does not appear significant enough to impact PCSSD performance. This can be seen in the right column of Figure 9 as the background removed $I_P$ images are consistent with the case of no thermal background present (left column).

For all of these tests, the resulting $I_P$ SNRs are summarized in Table 1. The SNR was calculated using $10^3$ independently generated 30 s photon lists and is given by $\text{SNR} = (\langle I_P \rangle - \langle \text{Background} \rangle)/(\text{std. dev.} \langle I_P \rangle)$. Here the “Background” is determined using photon lists without any injected companions and “std.dev. $\langle I_P \rangle$” is the standard deviation of the mean companion intensity.

In summary, a constant 50 photons $s^{-1}$ pixel$^{-1}$ background significantly degrades the $I_P$ SNR of the faint companions at all
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Simulated *I*<sub>C</sub>, *I*<sub>S</sub>, and *I*<sub>P</sub> images generated by running PCSSD on constant and Gamma-distributed companion intensities for two different suites of companion contrasts. Left: 4 × 10<sup>-4</sup>, 2 × 10<sup>-4</sup>, 1 × 10<sup>-4</sup>, 5 × 10<sup>-5</sup>. Right: 6 × 10<sup>-3</sup>, 3 × 10<sup>-3</sup>, 1.5 × 10<sup>-3</sup>, 7.5 × 10<sup>-4</sup>. The “True” columns are the input *I*<sub>C</sub>, *I*<sub>S</sub>, and *I*<sub>P</sub> images.</figcaption></figure>

separations when not removed. The perfect removal of these
photons in post-processing recovers the results as if there had
never been a background present, but this type of perfect
subtraction can be challenging with real data due to the current
energy resolution of MEC (R ~ 5). An exploration of how
imperfect background subtraction effects PCSSD performance
will be left for future work.

5. Discussion

To use photon arrival time statistics as a metric for
differentiating photons from the bright speckle halo from those
that originate from a faint companion, understanding and
incorporating information about the distributions from which
they originate is essential. This is especially true in regimes
where the companion intensity is comparable to that of the
underlying speckle field and deconvolving those two distribu-
tions becomes more important.

Additionally, identifying and removing any possible sources
of light outside of these distributions, or incorporating them
explicitly into the model being used, is important in very-high-
contrast regimes where background photons can easily out-
number photons coming from the faint source of interest.

Thankfully hardware fixes can remove out-of-band back-
ground light. Future work will identify and mitigate MEC's
known IR background at wavelengths longer than its science
band (y--J band). This will significantly reduce the number of
background counts and should aid the performance of even the
current PCSSD algorithm on faint sources.

PCSSD has been derived only for a linear combination of an
MR intensity distribution and a constant intensity component
(Walter et al. 2019). Section 4.2.2 shows that a companion
whose intensity instead follows a Gamma distribution can
significantly weaken the sensitivity of PCSSD. PCSSD could
be expanded within its current analytical framework to account
for a Gamma-distributed intensity component. Doing so could

require fitting an additional three parameters, rather than one
for a constant-intensity companion, unless some of the Gamma
distribution's parameters are independently known (e.g., from
AO telemetry). This could make performing the fit untenable
for even typical MKID data sets that can easily run in excess of
10⁹ photons.

Nonanalytic approaches can be computationally less expen-
sive than PCSSD and assume no knowledge about the
underlying intensity distributions except for the fact that they
are measurably different. This work, and the works cited
herein, have already demonstrated that intensity distributions
differ for speckles and incoherent on-axis sources. The most
prevalent nonanalytic SSD techniques overlap with machine-
learning-based approaches for detecting patterns in time-series
data such as long short-term memory recurrent neural networks
(Hochreiter & Schmidhuber 1997). To facilitate further
development of machine-learning-based SSD, the MKID
Exoplanet Direct Imaging Simulator (MEDIS) can be used to
generate more realistic simulated MKID data sets (Dodkins
et al. 2020).

In high-contrast imaging there is unique information
available at millisecond and microsecond timescales that is
lost when taking long-exposure images. New instruments
deploying fast, noiseless detectors (like MEC) are now able to
access this information and use it to start suppressing/
differentiating quasi-static speckles and companions. Though
not discussed here, work using millisecond images (sometimes
combined with wave-front sensor telemetry) to directly
measure and remove quasi-static speckles on-sky is another
promising path forward. However, these techniques have
additional challenges needing to run in conjunction with the
AO loop during observing instead of in post-processing
(Martinache et al. 2014; Gerard et al. 2018; Rodack et al.
2021). Moving forward, the advancement of both real-time and
post-processing techniques that leverage the information
available at these fast timescales will be necessary for reaching
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Simulated *I*<sub>C</sub>, *I*<sub>S</sub>, and *I*<sub>P</sub> images generated by running PCSSD on mock photon lists that sampled a Gamma PDF for the companion intensities. Left: No background included (same as the second column of Figure 8, left). Middle: Adding an additional constant and uncorrelated 50 photons s<sup>-1</sup> pixel<sup>-1</sup> background. Right: The 50 photons s<sup>-1</sup> pixel<sup>-1</sup> background is inserted but then removed before running PCSSD. Any photons that were excluded from the analysis due to their proximity to a background photon (i.e., falling within the 10 µs dead time) are not accounted for, which simulates performing a wavelength cut on MEC data in post-processing.</figcaption></figure>

<table><caption>**Table 1**
Figure 9 Companion Signal-to-Noise Ratio</caption><thead><tr><th rowspan="2">Contrast</th><th colspan="3">Separation = 3.5 λ/D</th><th colspan="3">Separation = 6.5 λ/D</th><th colspan="3">Separation = 9.5 λ/D</th><th colspan="3">Separation = 12.5 λ/D</th></tr><tr><th>Control</th><th>50 cps</th><th>Removed</th><th>Control</th><th>50 cps</th><th>Removed</th><th>Control</th><th>50 cps</th><th>Removed</th><th>Control</th><th>50 cps</th><th>Removed</th></tr></thead><tbody><tr><td>4·10<sup>-4</sup></td><td>4.5</td><td>3.4</td><td>4.2</td><td>6.3</td><td>4.4</td><td>6.8</td><td>7.4</td><td>3.0</td><td>7.0</td><td>5.8</td><td>1.3</td><td>6.4</td></tr><tr><td>2·10<sup>-4</sup></td><td>2.6</td><td>1.9</td><td>2.6</td><td>3.5</td><td>2.6</td><td>3.6</td><td>4.7</td><td>2.1</td><td>4.8</td><td>5.5</td><td>1.1</td><td>5.4</td></tr><tr><td>1·10<sup>-4</sup></td><td>1.6</td><td>1.0</td><td>1.7</td><td>2.2</td><td>1.4</td><td>2.2</td><td>2.6</td><td>1.5</td><td>2.7</td><td>3.4</td><td>0.8</td><td>3.6</td></tr><tr><td>5·10<sup>-5</sup></td><td>0.8</td><td>0.6</td><td>0.8</td><td>1.2</td><td>0.7</td><td>1.4</td><td>1.6</td><td>0.8</td><td>1.6</td><td>1.9</td><td>0.3</td><td>1.9</td></tr></tbody></table>

the best achievable contrasts with current and next-generation telescopes.

S.S. is supported by a grant from the Heising-Simons Foundation. T.B. gratefully acknowledges support from the Heising-Simons foundation and from NASA under grant #80NSSC18K0439. N.Sk. and V.D. acknowledge support from NASA grant #80NSSC19K0336. K.A. acknowledges funding from the Heising-Simons foundation.

The development of SCExAO was supported by the Japan Society for the Promotion of Science (grant-in-Aid for

Research #23340051, #26220704, #23103002, #19H00703 & #19H00695), the Astrobiology Center of the National Institutes of Natural Sciences, Japan, the Mt Cuba Foundation and the directors contingency fund at Subaru Telescope. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community and are most fortunate to have the opportunity to conduct observations from this mountain.

*Facilities:* Subaru (SCExAO Jovanovic et al. 2015b + MEC Walter et al. 2020).
*Software:* The MKID Pipeline (Steiger et al. 2022), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020).

## ORCID iDs

Sarah Steiger https://orcid.org/0000-0002-4787-3285

Timothy D. Brandt https://orcid.org/0000-0003-2630-8073

Olivier Guyon https://orcid.org/0000-0002-1097-9908

Noah Swimmer https://orcid.org/0000-0001-5721-8973

Julien Lozi https://orcid.org/0000-0002-3047-1845

Vincent Deo https://orcid.org/0000-0003-4514-7906

Sebastien Vievard https://orcid.org/0000-0003-4018-2569

Nour Skaf https://orcid.org/0000-0002-9372-5056

Kyohoon Ahn https://orcid.org/0000-0002-1094-852X

Nemanja Jovanovic https://orcid.org/0000-0001-5213-6207

Frantz Martinache https://orcid.org/0000-0003-1180-4138

Benjamin A. Mazin https://orcid.org/0000-0003-0526-1114

## References

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Boccaletti, A., Di Folco, E., Pantin, E., et al. 2020, A&A, 637, L5

Cagigal, M. P., & Canales, V. F. 2001, OptEn, 40, 2690

Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, L9

Currie, T., Daemgen, S., Debes, J., et al. 2014, ApJL, 780, L30

Currie, T., Lawson, K., Schneider, G., et al. 2022, NatAs, 6, 751

Dodkins, R. H., Davis, K. K., Lewis, B., et al. 2020, PASP, 132, 104503

Fitzgerald, M. P., & Graham, J. R. 2006, ApJ, 637, 541

Fruitwala, N., Strader, P., Cancelo, G., et al. 2020, RScI, 91, 124705

Gerard, B. L., Marois, C., & Galicher, R. 2018, AJ, 156, 106

Give'on, A., Kern, B., Shaklan, S., Moody, D. C., & Pueyo, L. 2007, Proc. SPIE, 6691, 66910A

Gladysz, S., & Christou, J. C. 2008, ApJ, 684, 1486

Gladysz, S., & Christou, J. C. 2009, ApJ, 698, 28

Gladysz, S., Christou, J. C., Bradford, L. W., & Roberts, L. C. 2008, PASP, 120, 1132

Gladysz, S., & Yaitskova, N. 2010, Proc. SPIE, 7387, 73870Z

Goodman, J. W. 1975, in Laser Speckle and Related Phenomena, ed. J. C. Dainty, Vol. 9 (Berlin: Springer), 9

Guyon, O. 2005, ApJ, 629, 592

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hochreiter, S., & Schmidhuber, J. 1997, Neural Comput., 9, 1735

Hunter, J. D. 2007, CSE, 9, 90

Jovanovic, N., Guyon, O., Martinache, F., et al. 2015a, ApJL, 813, L24

Jovanovic, N., Martinache, F., Guyon, O., et al. 2015b, PASP, 127, 890

Kuzuhara, M., Tamura, M., Kudo, T., et al. 2013, ApJ, 774, 11

Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770

Lagrange, A. M., Desort, M., Galland, F., Udry, S., & Mayor, M. 2009, A&A, 495, 335

Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, JATIS, 4, 049001

Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Sci, 350, 64

Marois, C., Doyon, R., Racine, R., & Nadeau, D. 2000, PASP, 112, 91

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2008a, ApJ, 673, 647

Marois, C., Macintosh, B., Barman, T., et al. 2008b, Sci, 322, 1348

Martinache, F., Guyon, O., Jovanovic, N., et al. 2014, PASP, 126, 565

Martinache, F., Jovanovic, N., & Guyon, O. 2016, A&A, 593, A33

Mawet, D., Pueyo, L., Lawson, P., et al. 2012, Proc. SPIE, 8442, 844204

Mazin, B. A., Bumble, B., Meeker, S. R., et al. 2012, OExpr, 20, 1503

Meeker, S. R., Mazin, B. A., Walter, A. B., et al. 2018, PASP, 130, 065001

Perrin, M. D., Sivaramakrishnan, A., Makidon, R. B., Oppenheimer, B. R., & Graham, J. R. 2003, ApJ, 596, 702

Rameau, J., Chauvin, G., Lagrange, A. M., et al. 2013, ApJL, 779, L26

Rodack, A. T., Frazin, R. A., Males, J. R., & Guyon, O. 2021, JOSAA, 38, 1541

Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118

Sivaramakrishnan, A., Koresko, C. D., Makidon, R. B., Berkefeld, T., & Kuchner, M. J. 2001, ApJ, 552, 397

Soummer, R., & Ferrari, A. 2007, ApJL, 663, L49

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Steiger, S., Bailey, J. I., Zobrist, N., et al. 2022, AJ, 163, 193

Steiger, S., Currie, T., Brandt, T. D., et al. 2021, AJ, 162, 44

Szypryt, P., Meeker, S. R., Coiffard, G., et al. 2017, OExpr, 25, 25894

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walter, A. B., Bockstiegel, C., Brandt, T. D., & Mazin, B. A. 2019, PASP, 131, 114506

Walter, A. B., Fruitwala, N., Steiger, S., et al. 2020, PASP, 132, 125005

Zhou, Y., Sanghi, A., Bowler, B. P., et al. 2022, ApJL, 934, L13