# Fast Coherent Differential Imaging on Ground-based Telescopes Using the Self-coherent Camera

Benjamin Gérard, Christian Marois, Raphaël Galicher

► To cite this version:

Benjamin Gérard, Christian Marois, Raphaël Galicher. Fast Coherent Differential Imaging on Ground-based Telescopes Using the Self-coherent Camera. The Astronomical Journal, 2018, 156 (3), pp.106. 10.3847/1538-3881/aad23e. hal-02293196

HAL Id: hal-02293196

https://hal.science/hal-02293196v1

Submitted on 22 Aug 2025

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_1.png" /></figure>

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

# Fast Coherent Differential Imaging on Ground-based Telescopes Using the Self-coherent Camera

Benjamin L. Gerard<fnref n="1" />,<fnref n="2" />, Christian Marois<fnref n="1" />,<fnref n="2" />, and Raphaël Galicher<fnref n="3" />

Received 2018 March 23; revised 2018 June 2; accepted 2018 July 6; published 2018 August 20

## Abstract

Direct imaging and spectral characterization of exoplanets using extreme adaptive optics (ExAO) is a key science goal of future extremely large telescopes and space observatories. However, quasi-static wavefront errors will limit the sensitivity of this endeavor. Additional limitations for ground-based telescopes arise from residual AO-corrected atmospheric wavefront errors, generating millisecond-lifetime speckles that average into a halo over a long exposure. A solution to both of these problems is to use the science camera of an ExAO system as a wavefront sensor to perform a fast measurement and correction method to minimize these aberrations as soon as they are detected. We develop the framework for one such method based on the self-coherent camera (SCC) to be applied to ground-based telescopes, called the Fast Atmospheric SCC Technique (FAST). We show that with the use of a specially designed coronagraph and coherent differential imaging algorithm, recording images every few milliseconds allows for a subtraction of atmospheric and static speckles while maintaining an algorithmic exoplanet throughput close to unity. Detailed simulations reach a contrast close to the photon noise limit after 30 s for a 1% bandpass in the H band on both zeroth and fifth magnitude stars. For the latter case, this is about 110 times better in raw contrast than what is currently achieved from ExAO instruments if we extrapolate for an hour of observing time, illustrating that the improvement in sensitivity from this method could play an important role in the future detection and characterization of lower mass exoplanets.

**Key words:** planets and satellites: detection – instrumentation: adaptive optics – techniques: image processing – techniques: interferometric

## 1. Introduction

Direct imaging of exoplanets on ground-based telescopes is currently limited in sensitivity by both dynamic and quasi-static wavefront errors. Dynamic wavefront errors arise from the residual aberrations left uncorrected by an adaptive optics (AO) system, seen by the science camera as short-lived atmospheric speckles to longer-lived dome seeing effects, ultimately averaging into a halo over a long exposure. Quasi-static wavefront errors instead arise from uncalibrated imperfections in the polishing and reflectivity of the optics in an instrument as well as from slower dome seeing effects. These aberrations produce a speckle pattern on the science camera. The pattern is quasi-static in nature due to slower thermal and flexure-related changes in the instrument and/or variable dome conditions, ultimately requiring daily or often hourly calibrations.

In the current paradigm of ground-based exoplanet imaging, optical correction of residual AO atmospheric speckles and quasi-static speckles is respectively addressed as follows:

1. by using the highest performing, “extreme” AO (ExAO) systems that produce a residual wavefront error of typically less than about 100 nm rms in the H band,

2. mainly by using high-quality optics that produce a quasi-static wavefront error of a few nanometers rms (e.g., Marois et al. 2008), although ongoing efforts are in place to calibrate low-order aberrations with the deformable mirror (DM) of an ExAO system (e.g., N’Diaye et al. 2016, Lamb et al. 2018).

Subsequently, two main post-processing strategies are then used to further attenuate the speckle pattern:

1. spectral differential imaging (SDI, Racine et al. 1999; Marois et al. 2000; Sparks & Ford 2002), where images are recorded at different wavelengths either in multiple narrow bands or with an integral field spectrograph (IFS) to distinguish magnification of speckles as a function of wavelength versus an exoplanet’s fixed position, and

2. angular differential imaging (ADI, Marois 2004; Marois et al. 2006), where the telescope’s pupil orientation is fixed while tracking a star such that the relative position angle of an exoplanet rotates while the position angle of the telescope optics remains fixed.

Although these two observing techniques are now often combined with advanced post-processing methods, such as a least-squares, singular value decomposition, and/or approaches involving principal component analysis (e.g., Lafrenière et al. 2007; Marois et al. 2010; Soummer et al. 2012; Devaney & Thiébaut 2017), SDI and ADI are still fundamentally limited by chromaticity (e.g., Marois et al. 2008) and temporal stability (e.g., Marois 2004; Vigan et al. 2010), respectively.

Atmospheric speckles provide a further limitation to temporal stability. Macintosh et al. (2005) studied the statistics of atmospheric speckles in an AO-corrected coronagraphic image. In the absence of quasi-static aberration, if images are photon noise limit so that atmospheric speckle amplitudes are consistently below the photon noise limit in individual exposures, the contrast (see Appendix B for a formal definition) will scale as $t^{-0.5}$, where t is the exposure time. If atmospheric

<footnote n="1">¹ University of Victoria, Department of Physics and Astronomy, 3800 Finnerty Rd, Victoria, V8P 5C2, Canada; bgerard@uvic.ca</footnote>

<footnote n="2">² National Research Council of Canada, Astronomy & Astrophysics Program, 5071 West Saanich Rd, Victoria, V9E 2E7, Canada</footnote>

<footnote n="3">³ Lesia, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, Univ. Paris Diderot, UPMC Univ. Paris 06, Sorbonne Paris Cité, 5 place Jules Janssen, F-92190 Meudon, France</footnote>
speckle amplitudes are instead above the photon noise limit in individual exposures, contrast will remain flat when $t < \tau_{\text{spec}}$, where $\tau_{\text{spec}}$ is the speckle lifetime. If the image is dominated by atmospheric speckles and $t > \tau_{\text{spec}}$, contrast will scale as $(t/\tau_{\text{spec}})^{-0.5}$ while atmospheric speckles average into a halo that is always above the photon noise limit. Macintosh et al. (2005) show that the speckle lifetime scales with the atmospheric clearing time of a telescope aperture, which for 10 m class telescopes is typically of the order of a fraction of a second. This timescale will thus motivate our adoption of a fast (i.e., integrating on timescales less than the speckle lifetime) measurement and subtraction strategy.

A family of solutions that can be used to correct for both quasi-static and residual atmospheric aberration, separate from ADI and/or SDI, is called coherent differential imaging (CDI). CDI uses the science camera of an ExAO system to measure and subtract coherent speckles from the star without removing any incoherent exoplanet light. Over the past decade, a number of CDI techniques have been proposed, developed, and tested in the lab (e.g., Guyon 2004; Baudoz et al. 2006; Bordé & Traub 2006; Give’On et al. 2007; Serabyn et al. 2011; Sauvage et al. 2012). These algorithms are generally designed for stable space telescope systems and require many iterations to converge at a minimum contrast. However, some CDI algorithms have recently been demonstrated on-sky, including:

1. Martinache et al. (2014) present the first on-sky CDI demonstration, using the speckle nulling algorithm originally proposed by Bordé & Traub (2006). This process involved stepping through the unknown phase of each speckle in order to minimize the focal plane intensity in a region of interest, ultimately converging after 20 iterations, where each iteration was a five second exposure. This iterative process limited contrast improvement to the first bright diffraction ring at $2.3 \lambda/D$ (where $\lambda$ is the wavelength of light and $D$ is the telescope diameter), a relatively stable quasi-static feature for these timescales of the order of a minute. Residual atmospheric speckles and quasi-static speckles at higher spatial frequencies were not correctable at these frame rates, ultimately yielding an improvement in contrast by a factor of 3.

2. Bottom et al. (2017) demonstrate a CDI approach using phase-shifting interferometry, where one path in a Mach–Zehnder interferometer is advanced through four piston offsets between 0 and $3\pi/2$ radians. This approach is demonstrated on-sky at Palomar to present the first CDI-based detection of a substellar companion, the brown dwarf HD 49197b. Again, because this algorithm is iterative in nature, measurements and corrections were could only made on timescales longer than a few seconds. Correction of residual atmospheric and quasi-static speckles that generate temporal instability on shorter timescales was not feasible, ultimately yielding an improvement in contrast by a factor of about 1–4, depending on the different post-processing methods that were used.

In this paper, we develop the framework for another CDI approach, called the self-coherent camera (SCC, Baudoz et al. 2006; Galicher et al. 2010), to be applied to ground-based telescopes in order to correct both atmospheric and quasi-static aberration, and hereafter referred to as “FAST” (for “Fast Atmospheric SCC Technique”). In principle, our FAST

approach is no longer limited by temporal stability as long as images are acquired faster than the dominant atmospheric and/or quasi-static speckle lifetime. In Section 2 we present the foundation for our new FAST solution, including an optimized focal plane mask (FPM), called the SCC FPM (Section 2.1), long exposure simulations using the SCC FPM that motivate the potential gain of a FAST strategy (Section 2.2), and a post-processing strategy to measure and subtract coherent stellar speckles (Section 2.3). In Section 3 we present detailed numerical simulations of FAST performance, and then we provide concluding remarks in Section 4. The appendices include a review of the main principles of the SCC (Appendix A, including post-processing algorithms in Appendix A.2 and photon noise propagation in Appendix A.3) and the parameters and assumptions used in our simulations (Appendix B). At this time we do not yet consider chromaticity/broadband solutions, instrument-specific simulations, DM control, or a laboratory demonstration.

## 2. The FAST Solution

In this section we will both present and justify the need for a FAST solution, which includes our new proposed SCC FPM and a post-processing subtraction strategy designed for millisecond-timescale exposures. Simulations using the SCC FPM for a 30 s exposure will also illustrate the limitations of typical ground-based CDI observations due to residual atmospheric turbulence.

### 2.1. The SCC FPM

The conventional SCC design (see Appendix A) will only work on long exposures to measure quasi-static aberration. This restriction is a result of the amount of light transmitted through the pinhole in the Lyot plane. All previous SCC papers have focused on correction of quasi-static aberration of the order of a few nanometers rms because of this limitation; here, we now address wavefront control with the SCC more generically, by measuring and correcting for speckles on all timescales, including atmospheric speckles.

Assuming that we want to detect and subtract speckles as soon as they are visible and before they disappear, when a few photons are recorded for a speckle we ideally want the modulation amplitude to be at least as large as the speckle amplitude that we want to remove, requiring that the envelope of the pinhole point-spread function (PSF) match the amplitude and power law of the stellar speckles. If this requirement is not met, fringes on a speckle of interest will be detected at a relatively lower signal-to-noise ratio (S/N) than the speckle, which will generate a noisier/biased calibration. Integrating for longer will increase the fringe S/N to mitigate this problem, but this integration time must always in principle be shorter than the speckle lifetime; otherwise a sufficient S/N measurement of a speckle's fringes can never be made before the speckle disappears.

We ultimately found that the most successful way to get enough light through the pinhole was to design a new coronagraph, called the SCC FPM, which is outlined below. We replace the amplitude mask of the FPM in a simple Lyot coronagraph and instead apply a tilt to separations between zero and $3 \lambda/D$, creating an off-axis pupil copy in the Lyot plane. This FPM acts as a spatial filter for the downstream off-axis pupil copy, filtering only low-order modes ($\le 3$ cycles/
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: instead of using an amplitude FPM, a phase shift is applied in the focal plane, called the SCC FPM. Middle: the Lyot plane intensity after using the SCC FPM. The main pupil is on the left, and the off-axis pupil created from the SCC FPM is on the right. Right: an enlarged region of the blue rectangle in the middle panel. The dashed purple circle shows the theoretical minimum size of this off-axis pupil copy coming from the SCC FPM of diameter 6 λ/D compared to our measured FWHM (the dotted blue circle), illustrating that we are within a factor of about 2 of this limit. The solid black circle shows the size of the SCC pinhole for comparison; it is still well below the diffraction limit (the FPM inner working angle for a diffraction limit corresponding to the maximum pinhole size is 9.3 λ/D). The integrated intensity in the Lyot pinhole with this FPM is increased by a factor of about 3.2 × 10⁵ compared to the normal (unapodized) Lyot coronagraph.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) The SCC image, with photon noise, for a 1 ms exposure of an $m_H = 0$ star using the SCC FPM design shown in Figure 1 and the wavefront parameters described in Appendix B. (b) The MTF of (a), illustrating that fringes (corresponding to the higher spatial frequency sidelobes in this image) are detected above the background photon noise floor.</figcaption></figure>

pupil). After applying a grid search of amplitude and width parameters for various additional low-order aberrations, we found the best modulation amplitude results when using a symmetric, two-dimensional Gaussian of amplitude 2.9 µm and FWHM 4.7 λ/D, shown in combination with the FPM tilt in Figure 1 along with the corresponding Lyot plane intensity. The right panel of Figure 1 shows the diffraction limit of an FPM of diameter 6 λ/D (i.e., for a telescope with diameter D/6), and illustrates that our best solution is within a factor of 2 of this limit. Although our SCC FPM solution redistributes more light toward the center of the off-axis pupil, we are not actually “shrinking” this pupil; this would require either a smaller telescope and/or a mask that completely removes the smallest spatial scales of the complex electric field in the focal plane.

A significant amount of light is still blocked from going through the pinhole. The pinhole cannot be made larger due to the limitations of Equation (9) (which requires that the 1.22 λ/d minimum of the pinhole PSF lie outside the AO control region, where d is the pinhole diameter remapped to the entrance pupil). This maximum pinhole size is also impossible to reach with our SCC FPM design of diameter 6 λ/D, as illustrated in the right panel of Figure 1; the diffraction limit is D/6 whereas the largest pinhole size is (from Equation (9))

$D/18.5$ (i.e., reaching this diffraction limit would require an FPM of diameter 18.5 λ/D, or an inner working angle of 9.3 λ/D). However, despite these limitations, by using the SCC FPM design in Figure 1 a sufficient amount of light is transmitted through the pinhole to detect fringes for exposures of the order of a few milliseconds; the image and modulation transfer function (MTF) for a 1 ms exposure of an $m_H = 0$ star using the new SCC FPM design are shown in Figure 2, illustrating that the MTF sidelobes of higher spatial frequency are detected above the photon noise floor.

## 2.2. Long Exposure Simulations and Analysis

In this section we argue that, even with the advent of our new SCC FPM design, the typical CDI strategy of taking long exposures to measure and correct for quasi-static aberration is still limited by the atmosphere. Using the simulation parameters and SCC FPM design presented in Appendix B and Section 2.1, respectively, Figure 3 illustrates the sensitivity limitations for a 30 s exposure of an $m_H = 0$ star; it shows (a) the recorded image, (b) a perfect subtraction of everything that is fringed in image (a) (i.e., using Equation (A8), assuming a noiseless, simultaneous measurement of the pinhole PSF during that 30 s exposure), (c) a high-pass filter of this
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) A recorded image for a 30 s exposure of an $m_H = 0$ star. (b) A subtraction of the left image using Equation (A8), assuming a perfect simultaneous measurement of the pinhole PSF during that 30 s exposure. (c) A high-pass filter of (b); a simulated exoplanet at $(x, y) = (10, 0)$ $\lambda/D$ at a flux-normalized contrast of $10^{-5}$ is identified by the blue circle. (d) A contrast curve of images (a)–(c) compared to the photon noise limit (labeled “phlim”). (e) Contrast at $10 \lambda/D$ vs. time after applying the same subtraction algorithm as used on image (c) to stacked images at all times in our simulation for both zeroth and fifth magnitude stars; the photon noise limit is also shown for comparison to each star.</figcaption></figure>

subtracted image using a $2 \times 2 \lambda/D$ median boxcar filter, (d) a contrast curve of these three images compared to the input image and photon noise limit, and (e) the contrast at $10 \lambda/D$ versus time of subtracted, high-pass-filtered images compared to the photon noise limit for both $m_H = 0$ and $m_H = 5$ stars.

We find that the contrast in Figure 3(b) is about 50% of the contrast in the input image (a), a non-negligible component of the image that is “invisible” to the SCC (or any CDI algorithm) because fringed atmospheric speckles are blurred out over the long exposure. The spatial scales in the residual atmospheric halo of Figure 3, labeled “output,” suggest that a high-pass filter of this image will improve the contrast, and it is thus shown to the right as “$F_{HP}{output}$.” Indeed, Figure 3(d) shows that high-pass filtering does improve the contrast by up to a factor of about 5 throughout the AO control region, although this is still a factor of about 3–20 above the photon noise limit. This contrast curve also illustrates that applying the same high-pass filter to the normal photon noise limit, derived from Appendix A.3, has a negligible effect.

As discussed in Section 1, Figure 3(e) illustrates the conclusions from Macintosh et al. (2005). For a fifth magnitude star, $F_{HP}{output}$ is dominated by photon noise and so contrast continually improves as $t^{-0.5}$. For the zeroth magnitude star, $F_{HP}{output}$ is flat and dominated by atmospheric speckles at $t < \tau_{spec} \approx 200$ ms and averaging as $(t/\tau_{spec})^{-0.5}$ at $t > \tau_{spec}$. More generally, we note that although the spatial and temporal statistics of residual atmospheric turbulence will change with varying observing conditions (causing a variable gap in

contrast), the main point of Figure 3 is to illustrate the following averaging properties for atmospheric speckles during long exposures:

1. Atmospheric speckles leave a residual halo that cannot be measured optically by the SCC or any CDI algorithm.

2. On bright stars, this ultimately limits the achievable contrast, even after post-processing, to a non-negligible factor above the photon noise limit.

3. On dim stars, a negligible number of photons are recorded per speckle lifetime, and so this effect is either at or buried within the photon noise limit.

Points 1 and 2 above illustrate a fundamental limitation on achievable raw contrast with any CDI algorithm during a long exposure (e.g., Sauvage et al. 2012). Thus, our next FAST approach, discussed below in Section 3, is designed to minimize this effect by running fast exposures every few milliseconds, subtracting these recorded images, and then stacking the subtracted residuals to obtain contrast improvement proportional to $t^{-0.5}$.

## 2.3. Estimation of the Pinhole PSF

The main principle of the SCC (Baudoz et al. 2006; Galicher et al. 2010, and subsequent papers) is summarized in Appendix A. From this section and previous work, we note that two different Fourier filtering algorithms on the recorded
SCC image, I, provide

$$|I| \equiv |F^{-1}\{F\{I\}m_1\}| = |A_S| |A_R|, \quad \text{and} \quad (1)$$

$$F^{-1}\{F\{I\}m_2\} = |A_S|^2 + |A_R|^2 + |A_P|^2, \quad (2)$$

where $m_1$ and $m_2$ are different binary masks, and $A_S$, $A_R$, and $A_P$ are the electric fields in the detector plane from the stellar speckles, reference pinhole, and exoplanet, respectively.

From these two equations with three unknowns, we cannot retrieve the exoplanet image $|A_P|^2$. Baudoz et al. (2012) approach this problem by recording a calibration image of $I_R \equiv |A_R|^2$ (i.e., the pinhole PSF), using an internal source and occulting the main beam of the Lyot stop. This is a viable approach to reconstructing the speckle pattern $|A_S|^2$ induced by small and/or static aberrations. However, it may not be so for measuring atmospheric speckles in exposures of the order of a millisecond unless the electric field going through the pinhole is relatively stable across multiple realizations of a translating atmospheric phase screen. We discuss this calibrated reference imaging approach later in Section 2.3.2. In Section 2.3.1, we discuss the approach of direct, on-sky measurement of the pinhole PSF. Finally, in Section 2.3.3 we present a filtering algorithm, using only the recorded on-sky SCC image to reconstruct the pinhole PSF.

### 2.3.1. Direct Measurement

Simultaneously obtaining a high-S/N image of the pinhole PSF along with a noisy target image should result in a contrast improvement proportional to $t^{-0.5}$ (see Appendices A.2 and A.3 for a further discussion). However, in reality this measurement of the pinhole PSF will also be affected by photon noise due to the finite amount of light transmitted through the SCC pinhole. Using the SCC FPM, we simulate a separate pinhole PSF measurement assuming use of a 50/50 beam splitter placed on the pinhole: 50% of the light going through the pinhole is transmitted through the beam splitter and focused onto the science camera, while the other 50% is redirected by a small angle (i.e., just downstream of the pinhole) and focused onto the same camera. This would also require a field stop in an upstream focal plane (either at or downstream from the FPM plane) to block any stellar background that would otherwise prevent any pinhole PSF light from being recorded with a reasonable S/N. The beam splitter will induce quasi-static and chromatic aberrations on the reference beam, but we do not simulate them in this paper because their effects should be of second order because the beam splitter is downstream from the FPM and can be super-polished to extremely high precision.

With the setup discussed above, the pinhole PSF for a one millisecond exposure of an $m_H = 0$ star is shown in Figure 4(a). Although this image is indeed recording as much as four photons per pixel in a one millisecond exposure, photon noise still clearly dominates compared to the image without photon noise in Figure 4(e). However, looking at the MTF of Figure 4(a) in Figure 4(b) shows that only the central region of Fourier space has physical representation of spatial scales larger than $\lambda/d$ (i.e., the central region of the MTF), while all other higher spatial frequencies in the MTF are only from recorded photon noise. Thus, using the Fourier filter $m_3$ in Figure 4(c) allows the noisy image in Figure 4(a) to be converted to the Fourier-filtered image in Figure 4(d). The ratio

of Figure 4(d) Figure 4(e), illustrated in Figure 4(f), shows that this reconstruction, even for an extremely low number of recorded photons, is still accurate to about ±10%. As long as this offset is random on a frame-by-frame basis (i.e., generated from photon noise), stacking the output of a subtracted image using Equations (1) and (2) (i.e., where a filtering algorithm as in Figure 4(d) is used to estimate the pinhole PSF) will continue to average randomly (see Appendix A.3 for a further discussion). Only a static offset generating speckles that are consistently over- or under-subtracted will create static residuals after stacking. This same principle applies to the other two-pinhole PSF reconstruction methods in Sections 2.3.2 and 2.3.3. Intuitively, there is no reason to believe here that any static offset would be present in Figure 4(e), since these residuals are generated from photon noise present at the same low spatial frequencies as the pinhole PSF, and low-order diffraction and/or atmospheric effects from the FPM are unlikely to cause such static offsets (but see Section 2.4).

### 2.3.2. Calibration Algorithm

Similar to the approach in Baudoz et al. (2012), a daytime pinhole PSF can be recorded with very high S/N by using an internal source and masking the main beam in the Lyot plane. However, this recorded calibration pinhole PSF will be without

1. an intrinsic flux normalization, unlike the direct measurement in Section 2.3.1, and

2. effects from residual AO atmospheric turbulence.

The implications of point 2 above are that additional low-order aberrations from the atmospheric phase screen (i.e., at spatial frequencies of less than 3 cycles/pupil for an SCC FPM of diameter $6 \lambda/D$) that are not present in the calibration pinhole PSF will ultimately limit the contrast achievable with this method. This limitation arises from non-centro-symmetric features of the on-sky pinhole PSF relative to the static pinhole PSF (discussed further in Section 2.4). Furthermore, we simulate the differences in static phase screens between a daytime calibration and on-sky measurement by using a different specific wavefront map of static phase and amplitude aberration for the calibration pinhole PSF versus the on-sky static component, where both are still normalized to the same rms values as described in Appendix B but are fully decorrelated from one another. Then, assuming we have properly calibrated the pinhole PSF with this approach, Equations (1) and (2) allow a subtraction of $|A_S|^2$ and $|A_R|^2$ without subtracting $|A_P|^2$.

### 2.3.3. Reconstruction Algorithm

In this section we will revisit the underdetermined system of Equations (1) and (2) (i.e., two equations and three unknowns). If $|A_P|^2$ could be removed from Equation (2), we would have a determined system. Noting that the spatial scale of $|A_R|^2$ is limited to $\lambda/d$, running a low-pass filter on Equations (1) and (2) to retain spatial scales larger than or equal to $\lambda/d$ should have the effect of “filtering the exoplanet light out of the image,” thus providing two equations with two unknowns and a solution for $|A_R|$. These steps are outlined below, where $F_{LP}\{$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) A recorded pinhole PSF, with photon noise, for a 1 ms integration on an $m_H = 0$ star using the SCC FPM. (b) The MTF of (a). (c) A binary mask used to Fourier filter (a) so that photon noise is suppressed on spatial scales below $\lambda/d$. (d) The result of the result of Fourier filtering (a): filter $I_R = F^{-1}\{F\{\text{noisy } I_R\}m_3\}$. (e) A recorded pinhole PSF where no photon noise is simulated. (f) The ratio of (d) to (e).</figcaption></figure>

is a linear low-pass filter operator:

$$ I_{m,lf} \equiv F_{LP}\{|I|\} = F_{LP}\{|A_S| |A_R|\} \\ \approx F_{LP}\{|A_S|\} |A_R|, \quad \text{and} \quad (3) $$

$$ \begin{aligned} \operatorname{im}_{\text{lf}} &\equiv F_{\text{LP}}\{F^{-1}\{F\{I\}m_2\}\} = F_{\text{LP}}\{|A_S|^2 + |A_R|^2 + |A_P|^2\} \\ &\approx [F_{\text{LP}}\{|A_S|\}|^2 + [F_{\text{LP}}\{|A_P|\}|^2 + |A_R|^2 \\ &\approx [F_{\text{LP}}\{|A_S|\}|^2 + |A_R|^2. \end{aligned} \quad (4) $$

We found that using a 3.9 × 3.9 $\lambda/D$ median boxcar filter for $F_{LP}\{\}$ provided both a good contrast and algorithmic exoplanet throughput (see Appendix B for our formal definition), although these metrics were not optimized with any tolerance requirements.

We are able to remove the $F_{LP}\{\}$ from $|A_R|$ because a 3.9 $\lambda/D$ median boxcar kernel has a negligible impact on the smallest $\lambda/d$ scales ($\lambda/d = 37 \lambda/D$ in our simulation) of $|A_R|$. Additionally, in Equation (4) we assume that $[F_{LP}\{|A_P|\}|^2 \approx 0$ because the central $\lambda/D$-sized core of the exoplanet PSF is removed by the low-pass filter (only about 2% of the peak flux in an exoplanet core remains after applying the low-pass filter), and we ignore the higher-order halo terms of the exoplanet PSF in the context of point source detection (Perrin et al. 2003). However, we do not make the same argument for the low-pass-filtered starlight term, $F_{LP}\{|A_S|^2\}$; because this term is of similar magnitude to the $|A_R|^2$ term, ignoring the low spatial frequency information of the stellar PSF would otherwise bias our estimate of $|A_R|^2$. Ultimately, the exoplanet PSF halo term

is absorbed into $F_{LP}\{|A_S|^2\}$ and/or $|A_R|^2$, which will eventually cause static biases in the regime of high exoplanet S/N. Thus, combining Equations (3) and (4) (now two equations and two unknowns), we find that

$$ |A_R|^2 = \frac{\operatorname{im}_{\text{lf}} - \sqrt{\operatorname{im}_{\text{lf}}^2 - 4 I_{m,lf}^2}}{2}. \quad (5) $$

Finally, combining Equations (1), (2), and (5) then allows us to subtract all of the terms in Equation (2) except for $|A_P|^2$.

## 2.4. Flux Normalization and Static Limitations

In this section we present a flux normalization procedure and discuss static limitations of the subtraction algorithms presented in Section 2.3. For all algorithms, we found additional static systematic offsets that ultimately prevented contrast improvement below a certain level. For this reason, we implemented two additional steps to flux-normalize the reconstructed image without subtracting any additional exoplanet flux:

1. Annular flux normalization by summation:

$$ \begin{aligned} \operatorname{im}_{\text{subt}}(r) &= F^{-1}\{F\{I_{\text{noisy}}\}m_2\} - c_{\text{norm}}(r) \\ &\quad \times [|F^{-1}\{F\{I_{\text{noisy}}\}m_1\}|^2 / \tilde{I}_R + \tilde{I}_R], \quad \text{where} \\ c_{\text{norm}}(r) &= \sum_r [F^{-1}\{F\{I_{\text{noisy}}\}m_2\}] / \sum_r \\ &\quad \times [|F^{-1}\{F\{I_{\text{noisy}}\}m_1\}|^2 / \tilde{I}_R + \tilde{I}_R], \end{aligned} \quad (6) $$
where $I_{noisy}$ represents the recorded noisy SCC image, $\tilde{I}_R$ represents the estimated/calibrated pinhole PSF via any of the algorithms presented so far, and $\Sigma_r$ denotes a sum of all the pixels within $0.3 \lambda/D$ of a given separation $r$. Equation (6) is equivalent to forcing the mean of the output subtracted image, $im_{subt}$, to 0. Note that, even though the target image, $F^{-1}\{F\{I_{noisy}\}m_2\}$, includes exoplanet light and therefore biases the normalization coefficient at the separation of the exoplanet, we found about 0%–15% losses in algorithmic exoplanet throughput from this effect for our simulated exoplanets at flux-normalized contrasts levels of $10^{-5}$ and $5 \times 10^{-5}$ (see Section 3), and otherwise found significant improvement of normalization by summation versus, for example, robust standard deviation.

2. High-pass filtering the already-flux-normalized image from Equation (6), but only after stacking images to the desired integration time. The annular normalization approach in step 1 will not remove symmetric, low-order residual aberrations that can build up over time because of small calibration errors, e.g., as discussed in Section 2.3.2. Thus, we use an aggressive $2.1 \times 2.1 \lambda/D$ median boxcar filter to generate a low-pass-filtered version of the flux-normalized subtracted image, $F_{LP}\{im_{subt}\}$, and then high-pass-filter the image via

$$F_{HP}\{im_{subt}\} = im_{subt} - F_{LP}\{im_{subt}\}. \quad (7)$$

Equation (7) is the final form of stacked images that are displayed and shown in contrast curves of this paper for the calibration and reconstruction algorithms. We found that annular normalization degraded contrast in the case of the direct pinhole PSF measurement algorithm, and so for this algorithm we use only a high-pass filter. However, interestingly, we did find static limitations for the direct pinhole PSF measurement algorithm if we did not use a high-pass filter, revealing the presence of static low-order aberrations as a result of our Fourier filtering algorithm illustrated in Figure 4. For an analogous comparison we also subsequently high-pass filter the photon noise limit images (see Appendix A.3 but also Figure 3, which shows this has a negligible effect).

Even in the absence of photon noise, there is still a fundamental limit on the performance of the calibration and reconstruction algorithms, shown in Figure 5. The labels “algo,” “cal,” and “filt” represent the reconstruction, calibration, and direct pinhole PSF algorithms presented in Sections 2.3.3, 2.3.2, and 2.3.1, respectively. Without photon noise, the direct pinhole PSF is limited by numerical noise at the $10^{-18}$ level (see Appendix A.3), and so this subtraction algorithm does not suffer from the same level of performance degradation in terms of achievable contrast improvement as the other two do.

Figure 5(a) shows the factor by which contrast improves compared to an input noiseless target image (but still with atmospheric and static wavefront components as described in Appendix B) and thus illustrates that we never expect to improve contrasts in a single subtraction by better factors of about 20 and 6 at $10 \lambda/D$ with the algo and cal algorithms, respectively. Because no photon noise is simulated here, these limits are related to imperfections in the estimate of the pinhole PSF, which ultimately cause an asymmetric speckle pattern in the image that biases the normalization coefficient from

Equation (6). Specifically, this bias arises from the term $[|F^{-1}\{F\{im_{noisy}\}m_1\}|^2/\tilde{I}_R + \tilde{I}_R]$ generated for the pinhole PSF calibration/reconstruction, since a combination of phase and amplitude aberration in the pupil plane produces a non-centrosymmetric speckle pattern (see Appendix B for a further discussion). This asymmetry will limit the contrast achievable in Equation (6) by using a single normalization coefficient at a fixed separation. However, over a long exposure the limits from Figure 5(a) may not be the same; a large number of uncorrelated residual atmospheric speckle patterns may average to allow greater improvement to be reached than the limits for a realization with a single atmospheric phase screen.

The predicted contrast limits for a long exposure are illustrated in Figure 5(b), showing the subtraction of a target image with only static aberration (i.e., no atmospheric or photon noise component is simulated). The atmospheric averaging effect described above should in principle not reach below the limits in Figure 5(b). Thus, once photon noise is added, which should in principle not create additional static effects, the contrast improvement proportional to $t^{-0.5}$ from stacking images is expected to flatten out around these levels.

Figure 5 also shows a variation in contrast improvement for the calibration algorithm depending on the specific phase screen realization used to generate the pinhole PSF. This illustrates the sensitivity of this technique to pupil modes of any spatial frequency that vary in each phase screen realization (although tip and tilt are removed by a least-squares subtraction from each static phase screen; see Appendix B for further details) and/or to low-order modes diffracted from the FPM into the off-axis Lyot stop pinhole; one instance of this effect is clearly seen in Figure 5(d), which shows residual low-order features that span the full AO control region (i.e., every individual subtracted cal image shows a different residual low-order pattern). This sensitivity can be either rewarding or degrading, varying contrast factors of $\pm 10$ in different parts of the image, and thus warrants a further study of how to favorably optimize this variation (e.g., using a deconvolution procedure to remove the dominant modes that are causing the contrast degradation), but that is beyond the scope of this paper.

### 3. FAST Simulations and Analysis

With the framework developed in Section 2, we consider the FAST approach to measuring and subtracting both quasi-static and atmospheric speckles over short exposures, in contrast to the long-exposure strategy in Section 2.2. Figure 6 shows the results of recording and subtracting an 8 ms exposure of a zeroth magnitude star and a 30 ms exposure of a fifth magnitude star according to each method described in Section 2.3. Figure 6 illustrates that we are able to reach close to the photon noise limit in a relatively short exposure for both a zeroth and a fifth magnitude star, although this limit is only a factor of about 10 times below the input target image. Assuming the use of a photon counting camera, photon noise will dominate in these short exposures. As already discussed in Section 1 and illustrated in Figure 3(e), the speckle lifetime in our simulations is about 200 ms. Thus, the choice of a 30 ms integration time for the $m_H = 5$ simulations does not risk the effect of AO halo averaging over multiple speckle lifetimes. With this concept in mind of running fast enough to “freeze” the atmospheric speckles, we also choose to use 8 ms exposures on the zeroth magnitude simulation. Although a frame rate of 1 kHz also reaches the same photon noise limit
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Performance limitations of the calibration (shown as “cal”) and pinhole PSF reconstruction (“algo”) algorithms, in the absence of photon noise. (a) Contrast improvement factor vs. radial separation for the algo and cal algorithms (i.e., a ratio of contrast curves between the input target image and output subtracted image). The target image includes both atmospheric and static wavefront components (see Appendix B for additional details). The cal algorithm is shown as the median and interquartile range (IQR) from 100 different iterations that each use a different static phase screen to generate the calibration pinhole PSF (but still with the same rms amplitude). (b) Contrast curves for the algo and cal algorithms applied to an SCC image with only static aberration, predicting the ultimate limits for either algorithm during a long exposure. (c) and (d) The algo and cal images, respectively, used to generate the contrast improvement curves in (a); a single random realization is selected from the 100 different cal iterations.</figcaption></figure>

after 30 s, the former frame rate of 125 Hz provides a larger contrast improvement in individual subtractions, as illustrated by comparing the two panels of Figure 6.

With the use of this subtraction procedure on individual exposures, Figures 7 and 8 show the results of averaging out to a cumulative exposure time of 30 s. An exoplanet is simulated at $(x, y) = (+9, 0)$ $\lambda/D$ with flux-normalized contrasts of $10^{-5}$ and $5 \times 10^{-5}$ for the zeroth and fifth magnitude stars, respectively. The respective algorithmic exoplanet throughputs for the algo, cal, and filt algorithms are 1.18, 1.0, and 0.8 for the zeroth magnitude star and 0.88, 0.84, and 0.77 for the fifth magnitude star.

As predicted in Section 2.4 (Figure 5), the contrast curves for the pinhole PSF calibration and reconstruction algorithms flatten out, likely due to static limitations from a combination of how well the pinhole PSF is estimated and how the

unrealistic assumptions of azimuthal symmetry in the images impact our flux normalization procedure. The calibration algorithm for the $m_H = 5$ simulation does not yet flatten out because the photon noise limit improves in contrast by only a factor of about 10 over 30 s, still above the limits illustrated for the $m_H = 0$ simulation. The predicted limits from Section 2.4, Figure 5(b)—the best achievable contrasts for algo and cal of $4 \times 10^{-6}$ and $2 \times 10^{-5}$, respectively—are mostly consistent with the results from Figure 7—measured 30 s contrasts for algo and cal of $1 \times 10^{-5}$ and $6 \times 10^{-6}$, respectively. Interestingly, the addition of atmospheric speckles and/or photon noise decreases the limit of the reconstruction algorithm by a factor of 2.5 but increases that of the calibration algorithm by a factor of 3.3.

Figure 7 shows that contrast in the direct pinhole PSF measurement algorithm scales as $t^{-0.5}$ out to 30 s in both
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Contrast curves showing the recorded input image, output subtracted images, and photon noise limit for (a) an 8 ms exposure on a zeroth magnitude star, and (b) a 30 ms exposure on a fifth magnitude star. The legend in (b) applies to both panels. The labels “algo,” “cal,” and “filt” represent the reconstruction, calibration, and direct pinhole PSF algorithms presented in Sections 2.3.3, 2.3.2, and 2.3.1, respectively.</figcaption></figure>

simulations. This is an important result, and shows that in principle a direct, simultaneous measurement of the full SCC image and pinhole PSF, even if each exposure can record only a few photons, is not limited by static aberration. For example, extrapolating the $m_H = 0$ results, 84 hr of telescope time would reach a $5\sigma$ contrast of $10^{-8}$ at $10\lambda/D$ for a 1% bandpass.

Running simulations longer than 30 s would not provide additional insight into the physical limitations considered in this paper, since the limiting and continuous effects observed in the upper left panel of Figure 7 provide sufficient information to predict the behavior out to longer exposure times for either the zeroth or fifth magnitude case. These effects are also predicted and discussed in a noiseless analysis of the limitations from static aberration in Section 2.4. For the fifth magnitude case, the algo and cal algorithms will flatten out at around the same contrast improvement factors as the zeroth magnitude case (approximately 12 and 20, respectively; we can already begin to see this flattening for the algo algorithm after 30 s in the lower left panel of Figure 7), while the filt algorithm will continuously improve within a factor of a few from the photon noise limit all the way out to the end of the simulation. The same behavior applies to the zeroth magnitude case, for which all three subtraction algorithms clearly predict continuously flat performance (for algo and cal) or a performance improving proportional to $t^{-0.5}$ (for filt) out to any exposure time.

## 4. Summary and Conclusions

In this paper we have presented a new novel approach to CDI on ground-based telescopes using the SCC, called the Fast Atmospheric SCC Technique (FAST). Our main findings are as follows.

1. We design a new FPM specialized for our FAST approach that allows fringes to be detected in exposures that are of the order of a few milliseconds (Section 2.1).

2. We illustrate that the standard approach to ground-based CDI (i.e., taking a long exposure to average out atmospheric speckles and measure only the remaining static speckles) is limited in achievable contrast by

residual atmospheric speckles that are above the photon noise limit but not measurable/correctable during long exposures (Section 2.2).

3. Developing on the framework of Baudoz et al. (2012), we show that a post-processing algorithm can be used to reconstruct and subtract both atmospheric and static components of an on-sky SCC image without subtracting an incoherent, off-axis exoplanet, but that this process benefits from additional knowledge of the SCC pinhole PSF (Section 2.3).

4. We propose three different algorithms to estimate the missing pinhole PSF information:

(a) a direct, “live” simultaneous measurement of the pinhole PSF along with the SCC image (Section 2.3.1),

(b) use of a daytime pinhole PSF calibration image (Section 2.3.2), and

(c) a pinhole PSF reconstruction algorithm that uses only the recorded SCC image and a series of filtering algorithms (Section 2.3.3).

5. In Section 3, with the new SCC FPM design from Section 2.1 and post-processing algorithms from Section 2.3, we illustrate that obtaining fast exposures of the order of a few milliseconds to freeze the atmosphere is a solution to the long-exposure limitations from Section 2.2. With this approach, we find that the direct pinhole PSF measurement algorithm continually improves contrast with time close to the photon noise limit all the way out to the end of our 30 s simulations for both a zeroth and a fifth magnitude star. This result suggests that high-contrast imaging instruments using this technique would no longer be limited by static or atmospheric aberration, and that “raw” contrast would improve simply by integrating longer.

This paper sets the foundation for fast focal-plane wavefront sensing with the SCC. Using the FAST framework presented in this paper, future papers will investigate FAST chromaticity, DM control, instrument-specific simulations, and a laboratory demonstration. Addressing these topics is essential before a FAST strategy can be implemented on-sky; we are currently considering FAST applications to a future upgrade of the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Simulation results for a 30 s stacked exposure on zeroth and fifth magnitude stars (top and bottom rows, respectively). The left and right columns show the contrast at 10 $\lambda/D$ vs. time and full contrast curves after 30 s, respectively. The legend in the lower left panel applies to all panels.</figcaption></figure>

Gemini Planet Imager (GPI) and/or Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) instrument(s).

We gratefully acknowledge research support of the Natural Sciences and Engineering Council (NSERC) of Canada through the Postgraduate Scholarships-Doctoral (PGS-D) program and the Technologies for Exo-Planetary Science Collaborative Research and Training Experience (TEPS CREATE) program. The authors wish to thank B. Macintosh, O. Guyon, J. Lozi, P. Pathak, and A. Sahoo for helpful discussions that improved this manuscript. The authors thank the anonymous referee for his or her comments and suggestions that have significantly improved this manuscript.

## Appendix A
The SCC

### A.1. Principle

The SCC was originally developed by Baudoz et al. (2006) and then modified by Galicher et al. (2010) specifically to be integrated with high-contrast imaging, utilizing the diffraction effects of a coronagraph. The optical effect of adding a hard-edged FPM, either in amplitude or in phase, to suppress the core of the stellar point-spread function (PSF) is to diffract light to higher spatial frequencies in the downstream pupil plane

than are present in the original entrance pupil. A Lyot stop is then used in this downstream pupil plane in order to suppress such diffraction effects, which would ultimately be seen in the coronagraphic image. However, the main design difference of the SCC is the additional use of an off-axis pinhole in the Lyot stop. The light transmitted through this pinhole will then interfere with the main Lyot beam in the downstream focal plane, and a camera placed in this focal plane with the necessary sampling will resolve fringes in the image. These recorded fringes retain both phase and amplitude information from the complex electric field of a star. However, no information from the exoplanet is encoded as long as its light does not lie on a discontinuity of the FPM (e.g., on the edge of an amplitude/phase mask). This concept is illustrated in Figure 9. With this setup, a recorded monochromatic image at wavelength $\lambda$ is given by (Baudoz et al. 2006)

$$I(\alpha) = |A_S(\alpha)|^2 + |A_P(\alpha)|^2 + |A_R(\alpha)|^2 + 2 \operatorname{Re} \left\{ A_S(\alpha) A_R(\alpha)^* e^{\left(\frac{2i\pi\alpha\xi_0}{\lambda}\right)} \right\}, \quad (8)$$

where $A_S(\alpha)$, $A_P(\alpha)$, and $A_R(\alpha)$ are the focal plane complex electric field components of the star, planet, and light from the Lyot stop pinhole, respectively, $\xi_0$ is the separation between the optical axis and the center of the pinhole in the plane of the Lyot
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The stacked images for zeroth and fifth magnitude stars (top and bottom rows, respectively) after 30 s from each PSF subtraction method presented in this paper: the algo, cal, and filt algorithms (left, middle, and right column, respectively). An exoplanet is simulated in each image at flux-normalized contrasts of $10^{-5}$ and $5 \times 10^{-5}$ for the zeroth and fifth magnitude stars, respectively.</figcaption></figure>

stop, and $\alpha$ defines the $(x, y)$ position in the image. The modulation term in Equation (8), (2)Re{$A_S(\alpha)A_R(\alpha)^* e^{(2i\pi\alpha\xi_0)/\lambda}$}, is a cosine term with a sub-$\lambda/D$ ($D$ is the diameter of the telescope) spatial scale defined by the separation in the Lyot plane between the main beam and the pinhole. The use of $\alpha$ will be removed from subsequent notation in this paper for simplicity, and similarly the modulation term will be replaced by $M(|A_R| |A_S|)$, since in this paper we will only be considering speckle subtraction by image processing instead of DM control.

In order to prevent $|A_R|$ from reaching a minimum inside the region of the image correctable by a DM, which would attenuate most of the fringes at that separation and prevent a measurement/correction at sufficient S/N, the maximum pinhole diameter, $d$, is (Mazoyer et al. 2014)

$$d \le 1.22\sqrt{2} D/N_{\text{act}}, \qquad (9)$$

where $N_{\text{act}}$ is the number of DM actuators across the pupil. This region of the image correctable by a DM will hereafter be referred to as the “AO control region.” In general, Equation (9) is only a limitation for the area of interest that we want to correct with a DM; a subtraction via image processing may require either a smaller or a larger pinhole (to correct for a larger or smaller area in the focal plane, respectively, than the AO control region).

## A.2. Post-processing Algorithms

The modulation term from Equation (8) can be algorithmically isolated in spatial frequency to prevent any loss of exoplanet light. The full processing algorithm to obtain this information from the recorded image is described in

Galicher et al. (2010, and references therein) and illustrated in Figure 10(a) and briefly summarized below:

1. taking the Fourier transform of the recorded image to obtain the complex-valued optical transfer function (OTF),

2. multiplying the OTF by a binary mask, $m_1$, to isolate an OTF sidelobe, generated from spatial scales smaller than $\lambda/D$ in the image, that contains only the fringed starlight. In order for this sidelobe to not overlap in the OTF plane with the main beam, the center of the pinhole in the Lyot plane must be placed at minimum separation, $\xi_0$, from the center of the pupil of (Galicher et al. 2010):

$$\xi_0 \ge 1.5 D + d/2 \qquad (10)$$

3. shifting the unmasked complex array to the center of the image, and then

4. applying an inverse Fourier transform back to the image plane.

This post-processing algorithm, or Fourier sidelobe isolation algorithm, is also the standard wavefront sensing algorithm used for DM control and is crucial to preserving optical exoplanet throughput. As long as no exoplanet light is sent through the reference pinhole, all of the exoplanet information will reside in the main beam of the OTF, whose amplitude is shown in Figure 10; by isolating the OTF sidelobe, we are directly measuring the electric field of the star without confusion/contamination from the exoplanet.

The product produced from step 4 above, referred to in the literature as $I_-$, is an algorithmic solution to isolate the complex amplitude of the modulation term in Equation (8). By masking only one sidelobe in the OTF, this removes both the factor of two in the modulation term and the modulated component of
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Upper left: an FPM is used to occult the PSF core. Lower left: this effect creates diffracted light in the downstream pupil plane at higher spatial frequencies than are present in the original entrance pupil. This effect occurs from any FPM, either in phase and/or in amplitude. The contrast scale is different in this image from the other three; the central pupil is saturated to illustrate the relative intensity of this diffraction effect at higher spatial frequencies, of the order of 10⁶ times dimmer compared to the central pupil. Lower right: a Lyot stop is used to suppress these diffraction effects, but a small off-axis pinhole transmits some of this diffracted light. Upper right: intensity in the downstream image plane after this Lyot stop is applied, showing sub-λ/D scale fringes from the interference between the pinhole and main Lyot beam. These fringes provide a direct measurement of the complex electric field in the image plane. The relative intensity in the Lyot pinhole is enhanced by a factor of 2.5 × 10⁵ in the upper and lower right images for illustration purposes. All images show the intensity of the complex electric field and are normalized to the peak intensity of the image core if no FPM was present.</figcaption></figure>

that term, M, since all other baselines are also removed (shifting the OTF sidelobe to the center of the Fourier plane after applying a binary mask will only change the phase of I_ in the focal plane, not the amplitude). Thus, taking the absolute value of I_ gives

$$|I_{-}| = |F^{-1}\{F\{I\}m_{1}\}| = |A_{S}| |A_{R}|. \qquad (11)$$

The phase term from this filtering algorithm is used for DM control in a least-squares approach by Mazoyer et al. (2014) and subsequent papers, but it will not be used in this paper for our PSF subtraction approach.

Another image term can be separately reconstructed by placing a different binary mask in the OTF plane, illustrated in Figure 10(b): removing the modulation term by isolating the main beam with an amplitude mask, m_2, then doing an inverse Fourier transform and taking the absolute value provides the unmodulated terms of the image:

$$F^{-1}\{F\{I\}m_{2}\} = |A_{S}|^{2} + |A_{R}|^{2} + |A_{P}|^{2}. \qquad (12)$$

No absolute value or real component operator is needed on the output of Equation (12); $F^{-1}\{F\{I\}m_2\}$ is entirely a real image, since the Fourier filtering mask $m_2$ is centro-symmetric. As discussed in Section 2.3, if there was a way to estimate the pinhole PSF, this would then allow a reconstruction of $|A_S|$ using Equation (11) and thus a full subtraction of all terms in Equation (12) other than $|A_P|^2$. We discuss methods of pinhole PSF estimation in Section 2.3. Here, assuming we have obtained a noiseless image of the pinhole PSF, from the same wavefront realization as Equation (12), the subtracted image,

im\_subt, is then

$$ \begin{aligned} \mathrm{im}_{\mathrm{subt}} &= F^{-1}\{F\{I\}m_2\} - \left(\frac{|I_-|^2}{I_R} + I_R\right) \\ &= (|A_S|^2 + |A_R|^2 + |A_P|^2) - \left(\frac{|I_-|^2}{I_R} + I_R\right) \\ &= |A_P|^2 + \left(|A_S|^2 - \frac{|I_-|^2}{I_R}\right). \end{aligned} \qquad (13) $$

Using the simulation parameters described in Appendix B and a single atmospheric and static phase screen realization, we obtain a starting $1\sigma$ contrast across the full $32 \times 32 \lambda/D$ AO control region of $6.3 \times 10^{-4}$. Using this noiseless target image to calculate $|I_-|$, and then additionally using the noiseless pinhole PSF, the resulting subtracted image from Equation (13) reaches a $1\sigma$ contrast of $2.0 \times 10^{-18}$ across the full AO control region (when no companion is simulated); this limit is from unphysical numerical noise in our simulations, which illustrates that we will not be anywhere close to seeing these effects once we consider the effects of photon noise, and that instead when analyzing noisy images we are isolating the physical effects of photon noise propagation through our coronagraphic image reconstruction algorithm. Thus, with no photon noise,

$$ \begin{aligned} & |A_S|^2 - \frac{|I_-|^2}{I_R} = 0, \text{ or} \\ & \frac{|F^{-1}\{F\{I\}m_1\}|^2}{I_R} = |A_S|^2. \end{aligned} \qquad (14) $$
However, in reality photon noise sets a fundamental limit on the contrast we can achieve in Equation (13). These effects will be discussed next in Section A.3.

### A.3. Photon Noise

To understand how photon noise is propagated through the subtraction algorithm in Equation (13), independent of the pinhole PSF estimation, we use the same noiseless $I_R$ but in addition a noisy target image, $im_{noisy}$, to estimate the modulation amplitude, $|L|_{\text{lest}}$. By decomposing the noisy image into a noiseless component and a noise-only component (i.e., $im_{noisy} = im_{noiseless} + [im_{noisy} - im_{noiseless}]$) the subtracted noisy image is then

$$ 
\begin{aligned}
im_{\text{subt, noisy}} &= F^{-1}\{F\{im_{\text{noisy}}\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2}{I_R} + I_R \right) \\
&= F^{-1}\{F\{im_{\text{noiseless}} + [im_{\text{noisy}} - im]\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2}{I_R} + I_R \right) \\
&= F^{-1}\{F\{im_{\text{noiseless}}\} m_2\} \\
&\quad + F^{-1}\{F\{im_{\text{noisy}} - im_{\text{noiseless}}\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2}{I_R} + I_R \right) \\
&= (|A_S|^2 + |A_R|^2 + |A_P|^2) \\
&\quad + F^{-1}\{F\{im_{\text{noisy}} - im_{\text{noiseless}}\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2}{I_R} + |A_R|^2 \right) \\
&= |A_P|^2 + F^{-1}\{F\{im_{\text{noisy}} - im_{\text{noiseless}}\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2 - |F^{-1}\{F\{im_{\text{noiseless}}\} m_1\}|^2}{I_R} \right) \\
&= |A_P|^2 + F^{-1}\{F\{im_{\text{noisy}} - im_{\text{noiseless}}\} m_2\} \\
&\quad - \left( \frac{|F^{-1}\{F\{im_{\text{noisy}}\} m_1\}|^2 - |F^{-1}\{F\{im_{\text{noiseless}}\} m_1\}|^2}{I_R} \right) \\
&\equiv |A_P|^2 + T_1 - T_2,
\end{aligned} 
\quad (15) $$

where the second to last line of Equation (A8) uses the numerical result from Equation (A7). Equation (A8) illustrates that the subtraction algorithm presented above (when using a simultaneous noiseless pinhole PSF and if no exoplanet is simulated) is the photon noise limit and that two different terms contribute to this limit:

1. $T_1$ is the photon noise limit from a noisy SCC image but Fourier-filtered to remove fringes according to Figure 10(b) so that no spatial scales exist in the image that are smaller than $\lambda/D$, and

2. $T_2$ is the photon noise limit of the reconstructed $|A_S|^2$ (i.e., starlight) term from a noisy SCC image but Fourier-filtered to isolate the modulation amplitude according to Figure 10(a).

Images and contrast curves for $im_{noisy}$, $im_{subt, noisy}$, $T_1$, and $T_2$ are shown in Figures 11(a) and (b), respectively, for a 1 ms

exposure on an $m_H = 0$ star along with the default parameters from Appendix B. A normal amplitude mask FPM and Lyot stop do not diffract enough light into the off-axis Lyot pinhole to see any fringes above the photon noise for a millisecond exposure, and so, as in Figure 9, we unphysically enhance the pinhole intensity in Figure 11 by a factor of $2.5 \times 10^5$. A physical solution to this problem will be presented in Section 2.1.

Figure 11 illustrates that photon noise in exposures of the order of a few milliseconds, even on bright stars, sets a fundamental limit on the contrast achievable with this coronagraphic image reconstruction technique. Comparing the two terms that contribute to this photon noise limit, the noise from $T_2$ dominates over the noise from $T_1$, illustrating that higher spatial frequencies of the image at the positions of the MTF sidelobe are more affected by photon noise (i.e., in terms of contrast) than the central beam of the MTF. This greater impact of photon noise at higher spatial frequencies can be explained by comparing the “S/N” of the MTF sidelobe versus the central beam in Figure 11(a): although the MTF sidelobes are clearly detected above the background photon noise, the central beam is still significantly brighter than the sidelobes and thus detected at a relatively higher “S/N.” Accordingly, at this higher spatial frequency, pixel-to-pixel noise in the image ultimately propagates through the $|L|$ reconstruction algorithm to $\lambda/D$ spatial scales in $|L|_{\text{lest}}$, causing $\lambda/D$ scale residuals in $T_2$ that dominate over $T_1$ instead of just pixel-to-pixel dominated photon noise. Ultimately, the relative impact of $T_1$ versus $T_2$ will depend directly on the modulation amplitude and thus the intensity through the off-axis pinhole of the Lyot stop (i.e., if the MTF sidelobes in Figure 11(a) were brighter, or detected at a higher “S/N,” $T_2$ would make a relatively lower contribution to the overall $T_1 - T_2$ photon noise limit). At the moment this pinhole light is enhanced by a factor of $2.5 \times 10^5$ (in intensity) for illustration, but this is indeed set to a similar value to that provided by the physical solution presented in Section 2.1. Thus, the illustration here showing that the noise from $T_2$ dominates over $T_1$ will also hold throughout this paper.

Figure 11(b) and Equation (A8) show that a direct measurement of the pinhole PSF at infinite S/N would allow subtraction to the photon noise limit (by definition), since $im_{subt, noisy} = T_1 - T_2$ (if no exoplanet is simulated). With a noiseless pinhole PSF, the photon noise limit is generated entirely from Poisson noise and linear operators, suggesting that the $im_{subt, noisy}$ outputs of each new recorded image should be completely uncorrelated from one another and that a continuous stacking of these subtracted millisecond images would improve contrast in proportion to $t^{-0.5}$. Thus, $T_1 - T_2$ or its filtered version (see Section 2.4) is shown in plots as “the photon noise limit” throughout this paper. Regarding additional noise propagated through our subtraction algorithms (Section 2.3) as the result of an imperfect estimate of the pinhole PSF, we can see by looking at the second to last line of Equation (A8) that as long as our estimate of the pinhole PSF is

1. symmetric (in pixel value distribution) around the true noiseless $I_R$ value, and

2. uncorrelated from frame to frame each time a new noisy image is generated,
<figure><img src="image_12.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

**Figure 10.** Two different Fourier filtering algorithms designed to isolate different SCC image terms in spatial frequency. (a) From Baudoz et al. (2006), a Fourier transform of the image (the optical transfer function, or OTF) yields (in amplitude) a main beam and two sidelobes of higher spatial frequency, the latter of which represents the electric field of the SCC image modulation amplitude and can be isolated as illustrated here. (b) A separate but similar Fourier filtering algorithm to isolate the unfringed image components.

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

**Figure 11.** (a) Assuming perfect knowledge of the pinhole PSF, various terms from Equation (A8) are shown for a 1 ms exposure, with photon noise, on an $m_H = 0$ star, with the pinhole intensity enhanced by a factor of $2.5 \times 10^5$, and all other parameters are the default values from Appendix B: the recorded noisy image (upper left), the subtracted image/photon noise limit (upper middle), the photon noise limit from a noisy SCC image but Fourier-filtered to remove fringes so that no spatial scales exist that are smaller than $\lambda/D$ ($T_1$, upper right), and the photon noise limit of the reconstructed $|A_s|^2$ term from a noisy SCC image but Fourier-filtered to isolate the modulation amplitude ($-T_2$, lower right). For all these images the FOV is $32 \times 32 \lambda/D$ and a color bar is shown on the right to illustrate the number of photons recorded and remaining in the subtracted images. The MTF of the upper left image is also shown in the lower left, illustrating the relatively higher impact of photon noise on the sidelobes (i.e., modulation amplitude) compared to the central beam. (b) The contrast curves for the images in panel (a).

the contrast should still improve in proportion to $t^{-0.5}$ simply by stacking subtracted images, although at a slightly worse contrast from the photon noise limit. However, as in Sections 2.4 and 3, in some cases one or both of these conditions may not be met, ultimately limiting the contrast achievable even after an infinite exposure time.

## Appendix B
### Simulation Parameters and Assumptions

Throughout this paper, unless explicitly stated, simulations are run at 1.6 µm, using an 8.2 m telescope, an entrance pupil with secondary obscuration (but no spiders) analogous to the Gemini South telescope, and a 1% bandpass filter. Mazoyer et al. (2014) show that chromaticity effects on the SCC performance should be negligible over this narrow bandwidth and have similar performance to the monochromatic case, and so here we only use the 1% bandpass filter for photon counting purposes in an otherwise achromatic simulation. For a Lyot coronagraph we use an amplitude FPM of diameter $5 \lambda/D$, and a circular Lyot

stop, 4% undersized relative to the entrance pupil; we did not use an apodizer in our simulations. The term “algorithmic exoplanet throughput” refers to the ratio of the peak exoplanet flux between the output and input of a PSF subtraction algorithm, with a recorded image as the input and a PSF-subtracted image (i.e., using image processing) as the output; this expression does not refer to optical/instrumental throughput of exoplanet light. We use the Fraunhofer approximation to propagate the electromagnetic field between pupil and focal planes. We assume the use of a $32 \times 32$ actuator square DM across the pupil to define the AO control region as a $32 \times 32 \lambda/D$ box around the optical axis of the image plane. We use a pinhole along the $+x$ axis of the Lyot plane with a relative size (Equation (9)) of $d/D = 1.22\sqrt{2}/32 \approx 5.4\%$ and relative separation (Equation (10)) of $\xi_0 = 1.6D + d/2$ (we use a factor of 1.6 instead of 1.5 so that the MTF sidelobes are “comfortably” not overlapping). We use a beam ratio (image size divided by pupil diameter, or number of pixels per $\lambda/D$) of 5.36. When simulating exposures to estimate photon noise, we calculate flux
in number of photons per second for the above parameters from Bessell et al. (1998) and then assume a transmission through the atmosphere of 90%, transmission through the telescope and instrument of 20%, and detector quantum efficiency of 80%. In this paper, $5\sigma$ contrast curves are calculated from an image as follows.

1. Normalize the image: divide the propagated coronagraphic image by the peak pixel value of the non-coronagraphic image, obtained from the propagated focal-plane intensity before applying an FPM but with a Lyot stop placed in the entrance pupil.

2. Calculate the standard deviation of the image output from step 1 in an annulus $0.6 \lambda/D$ wide centered at radial separation $\rho$, and then multiply the result by 5 to obtain a "$5\sigma$ contrast."

3. Repeat step 2 at each radial separation, increasing by increments of $0.2 \lambda/D$ between 3 and $16 \lambda/D$, and then plot each contrast value versus each $\rho$ value to obtain a "$5\sigma$ contrast curve."

To simulate quasi-static aberration, we tune the spatial power spectral density (PSD) amplitude and power law of both phase and amplitude aberration to match a raw GPI contrast curve of a single IFS slice from a stable, 30 s exposure image; in this sequence, consecutive exposures are about 90% correlated, suggesting that we are seeing mostly static wavefront error and a small amount of residual atmosphere. Additionally, GPI images show no/very little symmetry; as in Marois (2004), pure phase or pure amplitude aberration should be approximately symmetric in a coronagraphic image for a small amount of aberration (i.e., so that a first-order Taylor expansion of the PSF from Perrin et al. 2003 is a valid approximation), and so the absence of symmetry implies an approximately equal contribution of phase and amplitude aberration to generating quasi-static speckles. With this in mind, we use a 25 nm rms, -1.5 power-law phase aberration and a 1% rms, -2 power-law (in intensity) amplitude aberration to generate a static phase and amplitude aberration, from both of which tip and tilt are subsequently removed by a least-squares subtraction.

To approximate an ExAO system, we use a 100 nm rms phase screen conjugated to the entrance pupil, or by the Maréchal approximation a Strehl ratio of 0.87, with a -2 power law, which is typical for AO-corrected residual turbulence (J. P. Véran 2018, private communication), subsequently removing tip and tilt to simulate the effect of a low-order wavefront sensor and remove lower order pointing effects. We do not consider the effects of residual uncorrected aberrations (i.e., the phase screen PSD amplitude and shape should be significantly larger, of the order of tens of microns rms, and closer to -11/3, respectively, beyond the spatial frequencies corresponding to the DM control radius) and therefore will not analyze performance beyond the control radius of $16 \lambda/D$. We use the framework from Srinath et al. (2015), where the $\alpha$ parameter, a dimensionless number

between 0 and 1, is used to simulate atmospheric "boiling," or deviating random effects from pure frozen flow; e.g., $\alpha = 0.95$ means that 95% of the current phase screen (or $\sqrt{100^2(0.95)} \approx 97$ nm rms) will be translated to the next phase screen realization via a Fourier shift based on the time interval, telescope diameter, and wind speed and direction, whereas the other $\sqrt{100^2(1-0.95)} \approx 22$ nm rms component of the next phase screen realization will be randomly generated, but with the same -2 power law. In this paper we use a wind speed of $10 \text{ m s}^{-1}$ in the $+x$ pupil direction, 1 ms time intervals, and $\alpha = 0.95$ (i.e., 5% random turbulence is added every millisecond).

## ORCID iDs

Benjamin L. Gerard https://orcid.org/0000-0003-3978-9195

Christian Marois https://orcid.org/0000-0002-4164-4182

## References

Baudoz, P., Boccaletti, A., Baudrand, J., & Rouan, D. 2006, in IAU Coll. 200: Direct Imaging of Exoplanets: Science & Techniques, ed. C. Aime & F. Vakili (Cambridge: Cambridge Univ. Press), 553

Baudoz, P., Mazoyer, J., Mas, M., Galicher, R., & Rousset, G. 2012, Proc. SPIE, 8446, 84468C

Bessell, M. S., Castelli, F., & Plez, B. 1998, A&A, 333, 231

Bordé, P. J., & Traub, W. A. 2006, ApJ, 638, 488

Bottom, M., Wallace, J. K., Bartos, R. D., Shelton, J. C., & Serabyn, E. 2017, MNRAS, 464, 2937

Devaney, N., & Thiebaut, É 2017, MNRAS, 472, 3734

Galicher, R., Baudoz, P., Rousset, G., Totems, J., & Mas, M. 2010, A&A, 509, A31

Give'On, A., Belikov, R., Shaklan, S., & Kasdin, J. 2007, OExpr, 15, 12338

Guyon, O. 2004, ApJ, 615, 562

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

Lamb, M., Norton, A., Macintosh, B., et al. 2018, Proc. SPIE, 10703, 107035M

Macintosh, B., Poyneer, L., Sivaramakrishnan, A., & Marois, C. 2005, Proc. SPIE, 5903, 170

Marois, C. 2004, PhD dissertation, Université de Montréal

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., Macintosh, B., Soummer, R., Poyneer, L., & Bauman, B. 2008, Proc. SPIE, 7015, 70151T

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

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

Mazoyer, J., Baudoz, P., Galicher, R., & Rousset, G. 2014, A&A, 564, L1

N'Diaye, M., Vigan, A., Dohlen, K., et al. 2016, A&A, 592, A79

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

Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, PASP, 111, 587

Sauvage, J.-F., Mugnier, L., Paul, B., & Villecroze, R. 2012, OptL, 37, 4808

Serabyn, E., Wallace, J. K., & Mawet, D. 2011, ApOpt, 50, 5453

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

Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543

Srinath, S., Poyneer, L. A., Rudy, A. R., & Ammons, S. M. 2015, OExpr, 23, 33335

Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, 407, 71