Stochastic Speckle Discrimination with Time-Tagged Photon Lists: Digging Below the Speckle Noise Floor

TL;DR

This paper introduces a photon timing-based algorithm that distinguishes speckles from incoherent sources in high contrast images, enabling detection of exoplanets closer to stars by overcoming speckle noise limitations.

Contribution

The novel algorithm leverages photon arrival time distributions to separate speckles from incoherent sources, improving high contrast imaging capabilities at small angular separations.

Findings

Can overcome speckle noise limit in simulated data

Effective when a coronagraph reduces diffraction patterns

Allows subsequent application of differential imaging techniques

Abstract

We present an algorithm that uses the distribution of photon arrival times to distinguish speckles from incoherent sources, like planets and disks, in high contrast images. Using simulated data, we show that our approach can overcome the noise limit from fluctuating speckle intensity. The algorithm is likely to be most advantageous when a coronagraph limits the coherent diffraction pattern in the image plane but the intensity is still strongly modulated by fast-timescale uncorrected stellar light, for example from atmospheric turbulence. These conditions are common at small inner working angles and will allow probing of exoplanet populations at smaller angular separations. The technique requires a fast science camera that can temporally resolve the speckle fluctuations, and the detection of many photons per speckle decorrelation time. Since the algorithm directly extracts the incoherent light, standard differential imaging post-processing techniques can be performed afterwards to further boost the signal.

Tables

Table 1. Table 1 Companion SSD Signal-to-Noise Ratio

		$\mathrm{Separation}=3.5\lambda /D$				$\mathrm{Separation}=6.5\lambda /D$				$\mathrm{Separation}=9.5\lambda /D$				$\mathrm{Separation}=12.5\lambda /D$
Contrast		MLE	MAP	Limit$\dagger$$\dagger$Detection limit for a 30 s long exposure with perfect PSF subtraction of the stellar light.		MLE	MAP	Limit$\dagger$$\dagger$Detection limit for a 30 s long exposure with perfect PSF subtraction of the stellar light.		MLE	MAP	Limit$\dagger$$\dagger$Detection limit for a 30 s long exposure with perfect PSF subtraction of the stellar light.		MLE	MAP	Limit$\dagger$$\dagger$Detection limit for a 30 s long exposure with perfect PSF subtraction of the stellar light.
$4\cdot {10}^{-4}$		4.8	5.3	2.3		6.8	7.0	3.8		7.6	8.3	5.7		7.0	9.4	8.8
$2\cdot {10}^{-4}$		2.8	3.0	1.2		4.0	4.0	1.9		4.9	4.9	2.9		5.6	5.9	4.5
$1\cdot {10}^{-4}$		1.6	1.7	0.6		2.3	2.3	1.0		2.8	2.8	1.4		3.5	3.5	2.3
$5\cdot {10}^{-5}$		0.9	0.9	0.3		1.3	1.2	0.5		1.6	1.5	0.7		2.0	2.0	1.1

Equations

$${\rho }_{\mathrm{MR}}\left[I|I_c,I_s\right]=\frac{1}{I_s}\exp \left[-\frac{I+I_c}{I_s}\right]I_0\left[\frac{2\sqrt{I{I}_c}}{I_s}\right],$$

$${\sigma }_{I,\mathrm{tot}}^2\approx \frac{2{\tau }_s\left(I_s^2+2I_c{I}_s\right)+I_c+I_s+I_p}{T_{\mathrm{tot}}}.$$

$$p_{\star }[n|I_c,I_s]$$

$$=\frac{1}{I_s+1}{\left(1+\frac{1}{I_s}\right)}^{-n}\exp \left[-\frac{I_c}{I_s}\right]$$

$$\times L_n\left[\frac{-I_c}{I_s^2+I_s}\right]\exp \left[\frac{I_c}{I_s^2+I_s}\right],$$

$${\mathcal{L}}_i=\sum _{m=0}^{k_i}{p}_p[m|I_p]p_{\star }[k_i-m|I_c,I_s],$$

$${\mathcal{L}}_i=$$

$$\times \sum _{m=0}^{k_i}\frac{I_p^m}{m!}{\left(1+\frac{1}{I_s}\right)}^{-(k_i-m)}{L}_{k_i-m}\left[\frac{-I_c}{I_s^2+I_s}\right].$$

$$p[\delta t|I]=I{e}^{-I\delta t}.$$

$$p[\delta t|I_c,I_s]={\int }_0^{\infty }p[\delta t|I]{\rho }_{\mathrm{MR}}[I|I_c,I_s]dI,$$

$$p[\delta t|I_c,I_s]\propto {\int }_0^{\infty }p[\delta t|I]{\rho }_{\mathrm{MR}}[I|I_c,I_s]IdI.$$

$$p[\delta t|I_c,I_s,I_p]\propto {\int }_0^{\infty }p[\delta t|I+I_p]{\rho }_{\mathrm{MR}}[I|I_c,I_s](I+I_p)dI,$$

$${\int }_{{\tau }_0}^{\infty }p[\tau |I_c,I_s,I_p]d\tau =1.$$

$$u_i=\frac{1}{1+I_s\delta t_i}\text{leavevmode},\text{leavevmode}\text{leavevmode}{u}_{max}=\frac{1}{1+I_s{\tau }_0}.$$

$$\log (\mathcal{L})=$$

$$=$$

$$\displaystyle+\sum_{i=1}^{N}\log\big{[}I_{c}^{2}u_{i}^{5}+4I_{c}I_{s}u_{i}^{4}+(2I_{s}^{2}+2I_{p}I_{c})u_{i}^{3}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ +2I_{p}I_{s}u_{i}^{2}+I_{p}^{2}u_{i}\big{]}$$

$$-N\frac{(u_{max}-1)(I_p+I_c{u}_{max})}{I_s{u}_{max}}$$

$$-N\log \left[I_c{u}_{max}^3+I_s{u}_{max}^2+I_p{u}_{max}\right].$$

$$\log (\mathcal{L})\to \log (\mathcal{L})-\frac{1}{2}{\left(\left(I_c-{\tilde{I}}_c\right)/\sigma [{\tilde{I}}_c]\right)}^2$$

$$⟨I⟩=I_c+I_s$$

$${\sigma }_I^2=⟨I^2⟩-⟨I{⟩}^2=I_s^2+2I_c{I}_s$$

$$⟨(I_i-⟨I⟩)(I_j-⟨I⟩)⟩={\sigma }_I^2\exp \left[-\frac{|t_i-t_j|}{{\tau }_s}\right].$$

$$\overline{I}=\frac{1}{N}\sum _i{I}_i.$$

$$⟨{\overline{I}}^2⟩-⟨\overline{I}{⟩}^2$$

$$\displaystyle=\frac{1}{N^{2}}\bigg{\langle}\sum_{i}\sum_{j}\left(I_{i}-\langle I\rangle\right)\left(I_{j}-\langle I\rangle\right)\bigg{\rangle}$$

$$=\frac{{\sigma }_I^2}{N^2}\sum _i\sum _j\exp \left[-\frac{|t_i-t_j|}{{\tau }_s}\right].$$

$$\frac{1}{N^2}\sum _i\sum _j\exp \left[-\frac{|t_i-t_j|}{{\tau }_s}\right]$$

$$=\frac{1}{N^2}\sum _{i=0}^{N-1}\left(\sum _{j=0}^{N-i-1}\exp \left[-\frac{j\delta t}{{\tau }_s}\right]+\sum _{j=0}^i\exp \left[-\frac{j\delta t}{{\tau }_s}\right]-1\right)$$

$$=\frac{1}{N^2}\sum _{i=0}^{N-1}\left(\frac{{\tau }_s}{\delta t}\left(2-\exp \left[-\frac{(N-i)\delta t}{{\tau }_s}\right]-\exp \left[-\frac{(i+1)\delta t}{{\tau }_s}\right]\right)-1\right)$$

$$=\frac{2{\tau }_s}{N\delta t}-\frac{2{\tau }_s}{N\delta t}\frac{1}{N}\left(\sum _{i=1}^N\exp \left[-\frac{i\delta t}{{\tau }_s}\right]\right)-\frac{1}{N}$$

$$=\frac{2{\tau }_s}{N\delta t}-2{\left(\frac{{\tau }_s}{N\delta t}\right)}^2\left(1-\exp \left[-\frac{(N+1)\delta t}{{\tau }_s}\right]\right)+\frac{2{\tau }_s}{N^2\delta t}-\frac{1}{N}.$$

$$⟨{\overline{I}}^2⟩-⟨\overline{I}{⟩}^2$$

$$=\left(\frac{2{\tau }_s}{T_{\mathrm{tot}}}-2{\left(\frac{{\tau }_s}{T_{\mathrm{tot}}}\right)}^2\left(1-\exp \left[-\frac{T_{\mathrm{tot}}}{{\tau }_s}\right]\right)\right)\left(I_s^2+2I_c{I}_s\right).$$

$$⟨{\overline{I}}^2⟩-⟨\overline{I}{⟩}^2$$

$${\sigma }_{I,\mathrm{tot}}^2\approx \frac{2{\tau }_s\left(I_s^2+2I_c{I}_s\right)+I_c+I_s+I_p}{T_{\mathrm{tot}}}.$$

Full text

1. Introduction

Direct imaging is a challenging exoplanet discovery and characterization technique due to the extreme contrast ($<$10*-4* for ground based targets) and small angular separations ($\lesssim$1*′′*) between the planetary companion and its stellar host. Despite this, adaptive optics (AO) and coronagraphy have enabled the discovery of planets up to $\sim$106 times fainter than their host stars (Marois et al., 2008; Lagrange et al., 2010; Kuzuhara et al., 2013; Macintosh et al., 2015; Keppler et al., 2018). Imaging an exoplanet requires subtracting the light of its host star in the form of the point-spread function (PSF). If this background were static and could be subtracted perfectly, exoplanet imaging would be limited only by the photon shot noise of the bright host star. Instead, high-contrast imaging is limited by uncontrolled scattered and diffracted light, which produces a coherent speckle halo in the image plane (Guyon, 2005).

Fast atmospheric speckles average down over an observation, while slower, quasistatic speckles must be removed using post-processing techniques. Angular differential imaging (ADI, Marois et al., 2006) exploits the rotation of the Earth, and hence the field-of-view of an altitude-azimuth telescope, to distinguish diffraction speckles from astrophysical sources. Spectral differential imaging (SDI, Racine et al., 1999; Marois et al., 2000; Sparks & Ford, 2002) uses the scaling of diffraction speckles with wavelength. Since the initial development of ADI and SDI, a variety of post-processing algorithms have refined their approaches to dig deeper into the stellar PSF (e.g. Lafrenière et al., 2007; Soummer et al., 2012; Marois et al., 2014).

The time variability and chromaticity of quasi-static speckles limit the performance of ADI and SDI (Gerard et al., 2019). Both techniques also suffer at small separations where exoplanets are more likely to hide. The speckle spectral dispersion used by SDI is proportional to the separation: close to the star, it becomes smaller than the planet’s PSF. For ADI, the arclength traced by the companion’s sky rotation is proportional to the separation. Furthermore, the precision of the background estimate for PSF subtraction is limited by low counting statistics at small separations (Mawet et al., 2014). Even without these issues, the variability induced by speckle fluctuations can dominate the photon noise and be well above the shot noise expected from the total number of photons.

Stochastic Speckle Discrimination (SSD, Gladysz & Christou, 2008) is a post-processing technique designed to reduce the additional noise caused by speckle fluctuations by temporally resolving them. It relies on the difference in photon arrival time statistics between a planet/extended source and the off-axis residual stellar speckles (Canales & Cagigal, 1999; Cagigal & Canales, 2001; Aime & Soummer, 2004; Fitzgerald & Graham, 2006; Soummer et al., 2007). With millisecond imaging cameras, SSD-like approaches have been shown to reduce speckle noise from even fast atmospheric speckles and improve the contrast limit for the detection of substellar companions (Gladysz et al., 2010; Frazin, 2016; Meeker et al., 2018; Stangalini et al., 2018).

In this paper, we present an improved version of SSD to exploit noise-free photon-counting cameras like MEC, the Microwave Kinetic Inductance Detector (MKID) Exoplanet Camera (Walter et al., 2018; Meeker et al., 2018) on Subaru Telescope’s SCExAO instrument (Lozi et al., 2018). Our approach statistically distinguishes a combination of constant and speckle intensity (both ultimately from the bright star) from incoherent light (from a planet or disk). At small separations, where ADI and SDI are least effective and the scientific questions are most pressing (Mawet et al., 2012), we expect SSD to offer strong improvements in the limiting detection contrast. This paper describes the new photon counting SSD technique and demonstrates its performance on simulated data.

We organize the paper as follows. In Section 2 we detail the simulation of photon lists to emulate the data expected from photon-counting cameras like MEC. Section 3 describes a formal extension of previous SSD analysis techniques with millisecond images using a maximum likelihood algorithm. We simulate the performance for various atmospheric conditions, planet brightnesses, and effective exposure times. Section 4 presents the new photon-counting SSD algorithm that estimates the incoherent light from a companion or disk directly from individual photon arrival timestamps. We demonstrate the algorithm on a simulated telescope image. We discuss the main results in Section 5, and conclude with Section 6.

2. Simulating Photon Arrival Times

2.1. Modeling the Stellar Speckle Intensity

The statistics governing the off-axis intensity distribution of coherent light with a partially developed speckle pattern have been studied at length. Originally derived by Goodman (1975) and verified experimentally by Cagigal & Canales (2001) and Fitzgerald & Graham (2006), the probability of getting an instantaneous intensity $I$ given $I_{c}$ and $I_{s}$ is governed by the Modified Rician (MR) distribution, defined as:

[TABLE]

where $I_{0}[x]$ denotes the zero-order modified Bessel function of the first kind. The parameter $I_{c}$ represents the intensity of the “constant” part of the diffraction pattern, i.e. the PSF of a star without the atmosphere, while $I_{s}$ is the intensity of the seeing halo which manifests as a “speckle” pattern.

Typically the AO system attempts to confine all the speckles into the constant diffraction pattern. A coronagraph can remove or transform this coherent portion out of the image plane. For all of the calculations in this paper we assume that the Strehl ratio, and by extension $I_{c}$ and $I_{s}$, remain constant.

Speckle intensity is correlated temporally. In the limit of Kolmogorov atmospheric turbulence and frozen flow, the speckle decorrelation time can be thought of as the wind crossing time across the telescope pupil (Macintosh et al., 2005). In reality, the turbulence is not Kolmogorovic and the AO loop and telescope vibrations can further complicate the speckle temporal power spectrum density (PSD) (Stangalini et al., 2016). This may result in a faster speckle decorrelation time, or a temporal PSD described by multiple exponential timescales. For the purposes of this paper we characterize the speckle PSD as a simple exponential decay with a characteristic speckle lifetime of $\tau_{s}=0.1$ s to roughly match empirical data in the near infrared (Fitzgerald & Graham, 2006; Meeker et al., 2018; Goebel et al., 2018). The speckle lifetime can change drastically depending on the atmospheric conditions but our results here can be qualitatively understood in those cases by scaling all parameters in time. A core condition for the SSD technique is for many photons to arrive in a single speckle decorrelation time: if this is satisfied, individual speckle fluctuations can be probed.

2.2. Modeling an Incoherent Source

A Poisson source incoherent with the stellar light and parameterized by intensity $I_{p}$ can be injected and the relevant photons uniquely identified in simulation. Any incoherent Poisson sources will be represented by this term, including binary companions, planetary companions, extended sources, dark current, and read noise. In the case that $I_{s}=0$, the MR distribution reduces to a Poisson distribution with intensity $I_{c}$ and will be indistinguishable from $I_{p}$. In the demonstrations that follow, we identify $I_{p}$ as an injected planetary companion.

The intensity fluctuations associated with speckles elevate the noise floor above the typical shot noise, i.e., it is harder to measure the planet’s intensity when it is embedded in a boiling speckle field. If the speckle temporal information is marginalized over, then we show in the Appendix that the total variance of temporally correlated intensities obeying MR statistics combined with photon shot noise is

[TABLE]

The variance of the measured intensity is inversely proportional to the total integration time, $T_{tot}$, as expected. This is the long exposure ($t_{\mathrm{exp}}\gg\tau_{s}$) photon noise limit faced by all PSF subtraction techniques like ADI and SDI. In the limit that $2\tau_{s}I_{s}\ll 1$, this variance reduces to pure shot noise on the number of photons. For high contrast imaging, typical parameters might be $\tau_{s}=0.1$ s, $I_{s}=50$ s*-1*, and $2\tau_{s}I_{s}=10$, in which case the noise from speckle fluctuations will dominate. In Section 4, we show that this noise can be overcome by temporally resolving individual fluctuations in the speckle background.

2.3. Generating Mock Photon Lists

We have developed code111Part of the MKID Pipeline python package available at https://github.com/MazinLab/MKIDPipeline for quickly generating mock photon lists with an optional injected planet, corresponding to the output of a single MKID-like pixel which is single photon counting with low noise. The photon lists obey the following rules:

The underlying intensities are MR distributed, but are correlated in time with $\langle I(t)I(t+\delta t)-\langle I\rangle^{2}\rangle\propto\exp{\left[-\delta t/\tau_{s}\right]}$. 2. 2.

$I_{c}$, $I_{s}$, and $I_{p}$ are independently specified by the user, such that $I_{c}$ and $I_{s}$ govern the MR statistics of the stellar intensity and $I_{p}$ is the mean count rate of a Poisson source. The total intensity should have an expectation value of $\left\langle I\right\rangle=I_{c}+I_{s}+I_{p}$. 3. 3.

Due to the intrinsic dead time in an MKID, photons are removed from the list if they arrive within $\tau_{0}$ of the previous valid photon’s arrival time.

The procedure for generating photon timestamps begins with creating a correlated list of random numbers that follow a Gaussian distribution. The random numbers are transformed to a uniform distribution ranging from 0 to 1, and finally transformed again to a MR distribution. The correlated MR sequence defines the intensity in photon counts for a small ($\ll\tau_{s}$) time bin (we used 200 $\mu$s), and a Poisson draw on that “instantaneous” intensity determines the number of photons that will finally be placed into that bin. The photons are distributed according to a uniform distribution in each bin. For simplicity, we assume neighboring pixels to have uncorrelated photon lists.

Since MKIDs do not have the same dark current or readout noise as conventional semiconductor detectors we do not add any additional noise. However, during high count rates photons can be lost due to a firmware triggering lockout that acts as a non-paralyzable dead time (van Eyken et al., 2015). We use $\tau_{0}=10\leavevmode\nobreak\ \mu$s to match the latest firmware implemented on MEC. This dead time formulation could be used with quasi-photon counting EMCCDs to account for photon pile-up.

3. SSD with Millisecond Images

Past efforts have been successful in extracting $I_{c}$ and $I_{s}$ parameters of astrophysical sources from a series of millisecond images (Fitzgerald & Graham, 2006; Gladysz et al., 2010; Meeker et al., 2018). In these experiments images were acquired with a fixed exposure time and the intensities at spatial locations under test were extracted into a light curve like that shown in the top panel of Figure 1. The distribution of intensities formed a histogram to which a modified Rician function could be fit. The middle panel of Figure 1 shows this process for three exposure times. Since faint companions can masquerade as static speckles, a large $I_{c}/I_{s}$ ratio can be used as a merit function for planet detection in comparison to the rest of the image field (Gladysz et al., 2010; Meeker et al., 2018).

3.1. Maximum Likelihood Model for Discrete Light Curves

In this section we extend previous work to a more formal approach by finding the best-fit values of $I_{c}$, $I_{s}$, and $I_{p}$ with a maximum likelihood algorithm operating on the light curve. This allows for a direct detection of the non-stellar intensity, $I_{p}$, in the form of a point source like a planet or an extended source like a protoplanetary disk. Furthermore, we can naturally introduce photon noise from a low-intensity photon counting regime. By allowing the MR distribution in Equation (1) to suffer a Poisson-Mandel transformation (Cagigal & Canales, 1999; Aime & Soummer, 2004), the discrete stellar intensity distribution becomes

[TABLE]

where $L_{n}$ is the $n^{\mathrm{th}}$ Laguerre polynomial. The units of intensity for this section are number of photons per exposure time, $t_{\mathrm{exp}}$. If a planet (or some other source that is incoherent with the star) exists in the field, the intensity distribution at that location will be the discrete convolution of $p_{\star}$ with the planet’s independent probability distribution $p_{p}$. To simplify, we assume the incoherent source has a Poisson probability distribution, $p_{p}[m|I_{p}]=\exp\left[-I_{p}\right]I_{p}^{m}/m!$ with average intensity $I_{p}$. The likelihood of the $i$th bin of a light curve containing $k$ photons given $I_{c}$, $I_{s}$, and $I_{p}$ is

[TABLE]

which becomes

[TABLE]

The likelihood of the entire light curve is $\mathcal{L}=\prod_{i}\mathcal{L}_{i}$. We use a Newton conjugate-gradient search to find the most likely values of $I_{c}$, $I_{s}$, and $I_{p}$.

The maximum likelihood estimate of the intensity distribution is overplotted onto the light curve histograms for three different exposure times in the middle panel of Figure 1. The bottom panel shows the likelihood functions marginalized over $I_{p}$ for these three exposure times.

3.2. Performance of Millisecond Imaging SSD

To understand the performance of this millisecond imaging SSD algorithm and the extent to which it can improve exoplanet detections, we produce receiver operator characteristic (ROC) curves (Tanner Jr. & Swets, 1954; DeLong et al., 1988; Krzanowski & Hand, 2009; Jensen-Clem et al., 2017). We achieve this by generating an ensemble of mock photon lists with the same nominal values of $I_{c}$, $I_{s}$, and $I_{p}=0$ and calculating maximum likelihood estimates to build up a distribution of $I_{p}$ corresponding to the signal-absent hypothesis. Next we inject a planet with $I_{p}>0$ and calculate maximum likelihood estimates to build a distribution for $I_{p}$ corresponding to the signal-present hypothesis. These distributions are shown in the left column of Figure 2. The cumulative distributions define the false positive and true positive fractions, which in turn define the shape of the ROC curve for a planet with a given intensity.

The right column of Figure 2 demonstrates the performance of the millisecond imaging SSD algorithm for different exposure times. The shape of the discrete intensity distribution can change depending on exposure time, which systematically affects the resulting maximum likelihood estimates (MLE) for $I_{c}$, $I_{s}$, and $I_{p}$ along with their uncertainties. In some cases there is a significant probability for the most likely value of $I_{p}$ to be equal to zero. In the case that $I_{c}\gg I_{s}$, the MLE probability can be multimodal with a significant peak at $\textrm{MLE}\leavevmode\nobreak\ I_{p}\approx\textrm{True}\leavevmode\nobreak\ I_{c}+\textrm{True}\leavevmode\nobreak\ I_{p}$. Both these behaviors appear for the same reason: when there is little modulation by $I_{s}$ of photons associated with $I_{c}$, then $I_{c}$ becomes difficult to distinguish from $I_{p}$. However, the total flux is still accurately recovered.

4. SSD in the Photon Counting Regime

4.1. Maximum Likelihood Model for Photon Arrival Times

We find that there exists an optimal camera frame rate that leads to the most precise probability distribution for $I_{c}$, $I_{s}$, and $I_{p}$. For exposure times that are too long ($t_{\mathrm{exp}}\gg\tau_{s}$) the speckle temporal information is averaged over. For exposure times that are too short ($t_{\mathrm{exp}}\ll\delta t$, the photon inter-arrival time) the correlation between subsequent photon arrivals is lost because frames are interpreted as an unordered set containing only 0 or 1 photons. In order to avoid this pitfall, we develop the posterior probability for $I_{c}$, $I_{s}$, and $I_{p}$ directly from the set of photon inter-arrival times.

We start by considering the (normalized) probability density for the next inter-photon arrival interval, $\delta t$, given a fixed intensity $I$:

[TABLE]

The speckle field intensity is not fixed but varies in time with the modified Rician probability density described in Equation (1). At a fixed point in time, the probability density of the next photon arrival time given $I_{c}$ and $I_{s}$ becomes

[TABLE]

where we have integrated the Poisson probability density from Equation (6) over all possible instantaneous stellar intensities, $I$. Intensities are considered to be in units of photons per second. We assume here that the speckle lifetime is much longer than the time it takes a photon to arrive, $\tau_{s}\gg\delta t$.

With a set of photon inter-arrival times, $\{\delta t_{i}\}$, we want the relative probability that a $\delta t$ is realized. At moments that happen to have higher instantaneous intensities the number of short $\delta t_{i}$’s will be increased. Thus, the relative probability of a photon inter-arrival time being in our data becomes

[TABLE]

Finally, we consider the case that light incoherent with the star, such as from a planet, is in the field. We consider only the simplest case in which this source has constant intensity $I_{p}$ and is governed by Poisson statistics. In this case, the relative probability of $\delta t$ given $I_{c}$, $I_{s}$, $I_{p}$ is

[TABLE]

which can be evaluated analytically. We then find the normalization constant by setting

[TABLE]

Equation (10) accounts for the non-paralyzable detector dead time, $\tau_{0}$, intrinsic to MKIDs by replacing the lower limit of integration with $\tau_{0}$. In Equation (9), we set the likelihood equal to zero for $\delta t<\tau_{0}$. It is convenient to use the change of variables,

[TABLE]

The log-likelihood, $\log(\mathcal{L})[I_{c},I_{s},I_{p}]$, then becomes

[TABLE]

As in Section 3, we assume that AO performance remains stable so that $I_{c}$ and $I_{s}$ remain constant over the course of observations.

Equation (4.1) for the photon counting SSD algorithm replaces Equation (3.1) from Section 3. We use a Newton conjugate-gradient search to find the maximum of the log-likelihood space and recover the best estimates for $I_{c}$, $I_{s}$, and $I_{p}$.

The photon-counting SSD algorithm consistently outperforms the millisecond imaging SSD algorithm from Section 3, which marginalizes over temporal information via the exposure time (see Figure 2). In Figure 3 we show the performance of the photon counting SSD algorithm under conditions given by various combinations of $I_{c}$ and $I_{s}$, and with various planet brightnesses $I_{p}$. The algorithm performs well with a high true positive detection rate even for the case when $I_{c}$ and $I_{s}$ are both large. However, the performance suffers in the case of $I_{c}\gg I_{s}$. As in millisecond imaging SSD, $I_{p}$ and $I_{c}$ become indistinguishable in this limit.

4.2. Maximum A Posteriori Estimation

Prior knowledge of a parameter can improve the estimates of $I_{c}$, $I_{s}$, and $I_{p}$. In our case, we commonly have information on the $I_{c}$ parameter which corresponds to the static or quasi-static speckle point spread function either from a telescope model or measured on a reference star. For such situations, the log likelihood function in Equation (4.1) can be modified with a Gaussian prior as

[TABLE]

where $\tilde{I}_{c}\pm\sigma[\tilde{I}_{c}]$ is the prior on $I_{c}$. The new estimates of $I_{c}$, $I_{s}$, and $I_{p}$ become maximum a posteriori (MAP) estimates.

4.3. Performance on Simulated Telescope Image

We evaluate the performance of the photon-counting SSD algorithm on a simulated 30 s telescope image without a coronagraph. We identify $I_{c}$ as the Airy-ring pattern of a diffraction limited telescope with a circular unobstructed aperature, $I_{s}$ as the seeing halo from atmospheric speckles with a Strehl ratio $\Sigma I_{c}/\Sigma(I_{c}+I_{s})=0.7$, and $I_{p}$ as a series of injected planets at various brightnesses and separations. The faintest planet’s total intensity is 40 photons/s, for a contrast of $5\cdot 10^{-5}$ with the host star, and is separated by 3.5 $\lambda$/D. Assuming a 5% end-to-end throughput on an 8.2 m telescope, the stellar magnitude in the near infrared is approximately $J=10$. We assume a speckle decorrelation time of $\tau_{s}=0.1$ s and detector dead time of $\tau_{0}=10$ $\mu$s.

Each pixel has an independent 30 s photon list generated from the “True” $I_{c}$, $I_{s}$, and $I_{p}$ shown in panels (c), (f), and (i) respectively of Figure 4. Spatial correlations in the photon lists are ignored for simplicity. The average intensity realized for each pixel is shown in Figure 4a. Figure 4b is the average intensity minus the expected light from the star, $(\textrm{True}\leavevmode\nobreak\ I_{c}+\textrm{True}\leavevmode\nobreak\ I_{s})$, which illustrates the best possible long exposure PSF subtraction (compare the background variance to Equation (2)). The MLE $I_{c}$, $I_{s}$, and $I_{p}$ are shown in Figure 4d, g, and j respectively. In Figure 4e, h, and k, we calculate the MAP estimates. We used the True $I_{c}\pm 3\cdot\sqrt{\text{True}\leavevmode\nobreak\ I_{c}}$ as a Gaussian prior on $I_{c}$; in practice one could use a telescope model or a reference PSF. The central black dot with radius 1.22 $\lambda$/D is not a coronagraph but simply obscures the on-axis light for convenience.

Comparing Figure 4d to Figure 4b shows that the MLE $I_{p}$ from the photon counting SSD algorithm recovers the injected planets better than a perfect stellar PSF subtraction (i.e. subtraction of the true $I_{c}$ and $I_{s}$). Figure 4d and e show that the MAP estimate for $I_{p}$ is not significantly better than the MLE $I_{p}$ (although the MAP estimate for $I_{c}$ is more precise). This is surprising because the MAP estimate includes a prior on $I_{c}$ that should help discriminate between $I_{c}$ and $I_{p}$. That this is not the case indicates we can take full advantage of the photon counting SSD algorithm without prior knowledge of the telescope PSF.

For the central pixel of each planet, we used $10^{5}$ independent photon lists to calculate the signal-to-noise ratio $\textrm{S/N}=(\langle I_{p}\rangle-\langle\textrm{Background}\rangle)/(\textrm{std. dev.}\langle I_{p}\rangle)$ where the $\langle\textrm{Background}\rangle$ is estimated by not injecting a planet. These are recorded in Table 1. The long exposure photon noise limit is also recorded in Table 1 where the estimated $I_{p}$ is equal to the total flux minus the light from the star, $(\textrm{True}\leavevmode\nobreak\ I_{c}+\textrm{True}\leavevmode\nobreak\ I_{s})$. While the results are from the simulated ensemble, they match the results from Equation (2). Table 1 is representative of 30 seconds of data, but the S/N will scale with $\sqrt{T_{\rm tot}}$. For a 2 minute exposure, all values in Table 1 should be scaled up by a factor of 2. The S/N was calculated using only the central pixel for convenience but would be larger if the surrounding pixels were considered.

5. Discussion

While it remains impossible to beat the photon shot noise $\sqrt{N}$, Figure 4 shows that the photon counting SSD algorithm can beat the long exposure ($t_{\mathrm{exp}}\gg\tau_{s}$) photon noise limit described by Equation (2). Table 1 quantifies the improvement in the S/N as a factor of 3 in the case of faintest planet ($5\cdot 10^{-5}$) at the nearest separation ($3.5\leavevmode\nobreak\ \lambda/D$). This is possible because the speckle fluctuations are temporally resolved ($t_{\mathrm{exp}}\ll\tau_{s}$) and individual speckles are probed by multiple photons ($\delta t\ll\tau_{s}$). In the case that the fluctuations from stellar speckles dominate the variance of the total intensity ($2I_{s}\tau_{s}\gg 1$), fast, noiseless detectors like MKIDs or EMCCDs are needed to dig beneath the noise. This often occurs in high contrast imaging at small separations ($\lesssim$ seeing radius) and is especially important at $\lesssim 5\leavevmode\nobreak\ \lambda/D$ where ADI and SDI start to lose their effectiveness (depending on spectral coverage, sky rotation, and AO performance).

SSD will benefit ADI and SDI by reducing speckle noise from the data that is fed into those algorithms. ADI processing can be approached the same way as usual, but instead of using raw images one would use the $I_{p}$ images produced with SSD. The modulation of the planet location will be unaffected by SSD. Similarly for SDI, the algorithm would be given $I_{p}$ maps at various wavelengths. This approach would require wavelength information for each detected photon, which is an intrinsic feature of an MKID detector.

Our SSD algorithm does not perform well when $I_{c}\gg I_{s}$ as seen in Figure 3. In this regime there is little modulation of the static speckle intensity $I_{c}$ by the atmospheric speckle field $I_{s}$, and as a result $I_{c}$ starts to become indistinguishable from the Poisson distributed $I_{p}$. This can be greatly mitigated with a coronagraph and possibly active speckle nulling (Martinache et al., 2014) which directly reduce $I_{c}$ in the image.

On a space-borne telescope atmospheric speckles are not a concern, implying that the image-plane intensity will have different temporal behavior from ground-based observatories. While our SSD algorithm relies on the intensity following a MR distribution, in general any distribution can be used so long as it is known. It may also be possible to modulate speckles in a controlled way using onboard AO.

6. Conclusions

In this paper we exploit photon arrival time statistics with a stochastic speckle discrimination (SSD) algorithm to distinguish planets from speckles. We first extend previous work with a formalized maximum likelihood algorithm operating on light curves with fixed, albeit fast, exposure times. We find that the likelihood space can sometimes result in bimodal behavior in the case of $I_{c}\gg I_{s}$. Additionally, the choice in exposure time can systematically skew the MLE of $I_{p}$ and inflate its variance. More generally, with a fixed exposure time, the performance will change as a function of parameters $I_{c}$, and $I_{s}$. This is a problem because $I_{c}$ and $I_{s}$ can change with observing conditions as well as with separation from the host star.

To overcome these difficulties, we have developed a new photon-counting SSD algorithm that calculates the maximum likelihood for $I_{p}$ directly from the individual photon arrival times. With this approach the likelihood space becomes smooth and unimodal and the precision is maximized. The planet detection performance can be better by a factor of 2 than perfect stellar PSF subtraction of a long exposure. This requires fast noiseless detectors like MKIDs.

We have made several simplifying assumptions in our analysis. We take the speckle temporal PSD to be described by the single exponential timescale $\tau_{s}$, we assume $I_{c}$, $I_{s}$, and $I_{p}$ to remain constant, and we assume that the MR distribution accurately describes the off-axis stellar intensity. Finally, we ignore chromaticity. These assumptions represent avenues of exploration for future work. The speckle temporal PSD can be measured and used to more accurately simulate photon lists. Since the instantaneous Strehl can be measured, a future implementation of this algorithm might include that as a priori information in the log likelihood model. This would also inform additional variance on $I_{p}$ apart from the Poisson noise.

SSD algorithms are ultimately constrained on two fronts. First, the photon arrival time $\delta t$ must be much shorter than the speckle decorrelation time $\tau_{s}$: we need many photons to characterize the properties of a materialized speckle. Second, the performance degrades when $I_{c}$ is large but $I_{s}$ is small, because the incoherent planet light masquerades as the static speckles described by $I_{c}$. Fortunately, this can be improved with a coronagraph.

In return, SSD algorithms are most useful when $I_{s}$ is large, which is often inescapable at small inner working angles. Furthermore, the results are not directly dependent on separation (although they are dependent on $I_{c}$ and $I_{s}$, which are larger at small separations). This makes SSD a powerful post-processing technique at small inner working angles where it can complement more established techniques like ADI and SDI. At high speckle intensities, SSD can even outperform the theoretical limits of a perfect implementation of ADI and SDI.

This work was supported by the National Science Foundation Grant 1710385 and NASA ROSES grant NNX15AG23G. We would like to thank Michael Fitzgerald (University of California, Los Angeles) for helpful comments during the early stages of this work.