Precision speckle pattern reconstruction for high contrast imaging

In High Contrast Imaging, a large instrumental, technological and algorithmic effort is made to reduce residual speckle noise and improve the detection capabilities. In this work, we explore the potential of using a precise physical description of speckle images, in conjunction with the optimal detection statistic to perform High Contrast Imaging. Our method uses short-exposure speckle images, reconstructing the Point Spread Function (PSF) of each image with phase retrieval algorithms. Using the reconstructed PSF's we calculate the optimal detection statistic for all images. We analyze the arising bias due to the use of a reconstructed PSF and correct for it completely up to its accumulation over $10^4$ images. We measure in simulations the method's sensitivity loss due to overfitting in the reconstruction process and get to an estimated 5$\sigma$ detection limit of $5\times 10^{-7}$ flux ratio at angular separations of $0.1 -0.5^{\prime\prime}$ for a $1h$ observation of Sirius A with a 2m-telescope.


INTRODUCTION
High Contrast Imaging (HCI), the method to detect exoplanets in which the exoplanet is seen directly as an additional source near its host star, is an important method for current and future characterization of exoplanets atmospheres (Currie et al. (2011), Konopacky et al. (2013)) as it measures light from the exoplanet itself.
HCI requires separating the faint exoplanet from its bright nearby host star, which is highly dependent on the telescope's angular resolution.For ground-based observatories, the atmosphere presents an additional constraint on the telescope's angular resolution.The turbulent flow of air in the atmosphere creates index of refraction variations, that create rapidly changing phase aberrations.The phase aberrations degrade the telescope's angular resolution, creating speckle images for short-exposure imaging and a broad seeing disk image for long-exposure.
A way to overcome atmospheric seeing is the use of Adaptive Optics (AO) systems that sense and correct, in real-time, the atmospheric phase aberrations.Such systems bring the angular resolution from the seeing limit close to the telescope's diffraction scale and are used by the leading HCI instruments on state-of-the-art telescopes (VLT-SPHERE Beuzit et al. (2019), Subaru-CHARIS Groff et al. (2015) and MEC Walter et al. (2020), Gemini-GPI Nielsen et al. (2019)).
Analyzing the idealized case of HCI with perfectly known, yet uncorrected, Point Spread Function (PSF), we obtain that even a modest-sized telescope in the seeing limited case can reach a fantastic contrast, comparable to the best performance of state-of-the-art facilities, as can be seen in the dashed lines in Figure 1.
Motivated by this computation, in this work, we propose a method to perform HCI using short exposure speckle images, our scheme is illustrated in Figure 2, the measurements are both imaging and wavefront sensing, a sequence of algorithms is used to reconstruct the atmospheric phase aberrations from the measurements 4.1, the optimal detection statistic is then calculated 2.3 and Parametric Bootstrap method is used to correct the bias arising in the procedure 4.2.Analysis of the expected performance is presented in 2.3 and 4.3 and the results presented in Figure 1. Figure 1.5σ Detection limit at different separations.The detection limit is calculated for a 2m-telescope, at seeing condition r0 = 15cm, for total of 10 14 photons (corresponds to observing Sirius A for ∼ 1h at 500-700nm with efficiency of ∼ 70%).The colored-dotted lines are the optimal detection limit described in Section 2, the colored-solid include overfitting losses described in Section 4, the black line shows the performance of GPI from Ruffio et al. (2017).The plot was created using the package by Bailey (2021).Schematically in our method, we take short exposure images, reconstruct the phase aberrations for each image, calculate the statistic using the image and reconstructed phase and sample to correct for its bias.

HIGH CONTRAST IMAGING THROUGH THE ATMOSPHERE
When imaging from the ground, the atmosphere presents a difficulty.As air flows turbulently, differ-ent regions contain air of different velocities and densities that cause spatial variations in the index of refraction.As light passes through this medium with an in-homogeneous index of refraction, optical path differences accumulate, and phase aberrations are presented on the pupil plane of our telescope.Phase aberrations degrade the PSF of our telescope, and the atmospheric phase aberrations change rapidly, the instantaneous degraded PSF is typically comprised of diffraction-limited speckles, called speckle image, and the time-averaged called the seeing disk.

Atmospheric seeing
The PSF is a characteristic of the optical system that can be calculated for monochromatic illumination as the Fourier transform of the aperture function (B) and the phase aberration (φ) (1) In the case of perfect imaging from a circular aperture, the PSF will take the shape of an Airy disk with an angular scale of the diffraction limit λ/D.The phase aberrations caused by the atmosphere degrade the PSF significantly and are a topic of many theoretical and observational studies.As a model for the atmospheric phase aberrations, we will use the classical Kolmogorov phase aberration structure function with a typical scale r 0 called the Fried parameter The Fried parameter can also be understood as the size of the telescope that transitions from being diffraction limited to seeing limited.
To model imaging in a finite wavelength band, we integrate the single wavelength PSFs incoherently and include the two leading effects that change the PSF as a function of the wavelength: where E λ is a term representing the spectrum of the source, and S λ/λ0 an operator that stretches the coordinates as the diffraction scale increases.

Statistical model
To push the limits of direct imaging we have to stay as close as possible to optimality, therefore we will start by mathematically formulating the high contrast imaging detection problem rigorously.
Directly imaging exoplanets can be formulated as a hypothesis testing question in which we want to distinguish between the null hypothesis of the images made of a single star with some flux f and the alternative hypothesis of an additional source with flux fraction at relative position q The mapping between the true sky T and the measured images is by the PSF, with additive noise we will assume has a shot-noise component related to the source, shot-noise related to the sky background and detector read noise, taking the Gaussian approximation of the Poisson distribution (justified as we have a high number of photons per pixel) where N is the normal distribution, b is the noise from sky background and read-noise, ⊗ is the convolution operator.

Detection statistic
To distinguish between the two hypotheses we will use the following statistic (with further discussion in Appendix A) where ← − P is the reverse of P , and φ is the Maximum-A-Posteriori estimator of the atmospheric phase aberrations phi.We can read this statistic as a filter matching the image's deviation from a point source, and subtracting the location of the primary source as the unknown flux would create a flux deficit at its location.And the variance associated with the image's variance is We can estimate the 5σ detection threshold of this statistic for the ideal case in which we know φ exactly min : 5 by simulating atmospheric phase screens from the power spectrum calculated by Noll (1976), the resulting contrast thresholds as a function of separation for different bandwidths are shown in the colored dashed lines in Figure 1.

SIMULATIONS
To test our method we used end-to-end numerical simulations of the proposed measurement instruments.
We generate atmospheric aberrations by sampling the Fourier modes of the Kolmogorov spectrum as derived by Noll (1976).Images are calculated with the assumed true sky image and PSF as calculated in Equation 3, and then sampled according to Poisson with additional Gaussian noise.The simulation of the wavefront sensor neglects chromatic behavior, as only the centroid position will be used.
Reference values for simulation parameters are listed in Table 1.

METHODS
In order to use this statistic, we first need to show a method to calculate φ, and even though we know how the statistic distributes given the correct φ, we need to examine carefully how it distributes when using φ.

Recovering the atmospheric phase aberrations
In order to calculate the statistic from Equation 7we have to find argmax{P (I|H 0 , P φ ) P (φ)} , the Maximum-a-Posteriori (MAP) estimator for the atmospheric phase aberrations.
The general problem of phase retrieval, recovering a phase from its Fourier modulus in the presence of noise is known to be hard and a topic of many studies (for example the reviews Fienup (2013), Shechtman et al. (2015)).Therefore we employ a simple yet powerful sequence of algorithms.Starting from direct measurement of the phase using a Shack-Hartmann WaveFront Sensor (SHWFS, Platt & Shack (2001)), then improving that estimator using the Gerchberg-Saxton (GS, Gerchberg & Saxton (1972)) algorithm and finally using Gradient Descent to ensure we get the MAP.
This procedure converges to the correct MAP for most instances in the case of enough photons (∼ 10 6 per image), reasonable atmosphere (r 0 = 15cm@D = 2m), but only for imaging in a single wavelength, further discussion, and elaborate results are presented in Appendix B.

Detection in the presence of learned PSF -H 0
When applying the statistic from Equation 7on learned PSF we have to deal with the effect of its inevitable errors.From simulations we learn, as expected, that our phase estimator can be modeled as distributing normally with some bias µ and covariance Σ around the correct atmospheric phase aberration φ ∼ N (φ + µ, Σ) . (10) From Equations 3 and 7 we can see that our statistic is a non-linear function of the phase estimator.Therefore, our estimator for the statistic, calculated using our estimator of the phase, is slightly biased but as it is a bias shared by all images, when we accumulate the statistic over many images it accumulates too and put an upper limit on the number of images we can use once it becomes significant.As the notation suggests, this bias depends on φ, and for simplicity, we will assume it varies slowly relative to our phase estimation error and therefore can be expanded in a series and its expectation value We extend the method of Parametric Bootstrap Dekking (2005) to estimate this bias.By sampling images and solving their corresponding phases based on the recovered phase estimator we can sample the distribution in Equation 10φ1 and calculate the following linear combination that has the expectation value as the bias we want to correct (detailed calculation in Appendix C) In practice, we calculate this term a few times to get better convergence to the mean and have a handle on its variance, which is empirically smaller than the inherent variance of the statistic, that is calculated in Equation 8, and can be reduced as 1/ √ N with more simulations.

Signal loss due to over fitting
After we made sure we are recovering null detections for images sampled from H 0 , we need to test the expected performance of the method for signals from H 1 .
An inevitable tension is present in our phase aberration estimation process.The posterior probability of the phase aberration, φ, is dominated by the deviation of its corresponding PSF, P φ , from the image.This will lead, in the presence of an additional source, to overfitting of the secondary as part of the primary.Detecting using an overfitted φ subtracts part of the secondary and therefore loses some of the signal available with perfect knowledge of φ.To quantify this effect, we perform injection-recovery simulations for different separations and different brightness.
In general, we expect the recovered SNR to be some function of the injected SNR and the separation from the star at which the signal is injected.In the weak signal limit at which we work, the injected SNR per image is small SNR inj 1, and we expect the recovered SNR to be proportional to the injected SNR, as can be seen in Figure 4a, and the dependence on the location to be some smooth function that takes into account the assumed phase aberrations spectrum and the degeneracy between the phase aberrations and the additional source, as can be seen in Figure 4b.
We use the values of α(|q|) from Figure 4 to convert the optimal dotted lines in Figure 1 to the more realistic solid lines.In the future, we will investigate using several wavelength bands in tandem to mitigate this loss, using the a-chromatic position of the secondary and the chromatic effect of the phase aberrations.

CONCLUSION AND OUTLOOK
We presented a method to reconstruct the atmospheric phase aberrations using phase retrieval algorithms, and use its corresponding PSFs to calculate the optimal statistic to detect the presence of a secondary faint source in a sequence of short-exposure speckle images.We presented a method to correct the arising  bias in our statistic, showing its complete correction for stacks of up to 10 4 images, and a pathway to increase it by expanding Equation 12 beyond second order, which is crucial for the use of the detection method.
The leading idea for our method is to optimize the PSF of the image through parameters of a physically constrained model, use the optimal statistic, and correct for its arising bias or other artifacts.The analysis presented in this work treats the case of imaging with no AO, but the method isn't limited to such observations.A probable avenue for improvement is to use AO to some extent which will reduce the photon noise of the host star, and concentrate the light of the planet to a diffraction-limited spot, thus enabling even better detection of the secondary source.
To advance this method to achieve high contrast on sky, further work is required, including: • Better phase retrieval scheme needs to be devised, to enable the retrieval from images taken in a narrow wavelength band.
• Demonstration of the method in-lab is essential and will ensure we can precisely describe speckle images.
• Generalizing the method for simultaneous imaging in a few bands, which is important for collecting enough light, and might make the phase retrieval more regular and therefore easier.
constant slope the image is an off-centered point, with the tip-tilt of the point related linearly to the slope of the wavefront.
In our simulation, we used an array of 26 × 26 sub-apertures which for our reference seeing conditions and telescope size resulted in each sub-aperture being of the size ∼ r 0 /2 × r 0 /2, to recover the phase aberrations from the WFS we used a maximum-a-posteriori estimator in the Gaussian regime and using a measured design matrix, as done by Clare (2004).An obvious, order-unity difference is visible between the two PSFs, therefore the measured phase is not accurate enough to be used for detection.
The phase aberrations recovered by the WFS have significant errors, with the PSF calculated using them having order-unity deviations from the input PSF as can be seen in Figure 5, therefore cannot be used to calculate the statistic.To improve the phase estimator we use further optimization techniques.

B.2. Gerchberg-Saxton algorithm
The Gerchberg-Saxton algorithm is an algorithm that was designed and proved to solve the phase retrieval problem under the l 2 norm by Gerchberg & Saxton (1972).We started the algorithm with the initial guess of the pupil function and the phase measured by the WFS, from that the algorithm transforms the guessed field to the image plane, constraining the amplitude by the measured image, then returning to the pupil plane constraining the pupil function again and repeating until convergence or some maximal number of iterations.
To quantify the convergence of the algorithm we compare an observable, likelihood-based criteria and an unobservable, phase error criteria, we simulated images in the single wavelength case for different seeing and flux conditions we examine the performance of the algorithm in the observable criteria and its agreement with the unobservable one.Table 2. Success rate percent for different seeing conditions and photons per speckle image (in the used sub-band).The success rate was determined according to the likelihood threshold for 10 3 attempts.
In Table 2 we show the fraction of phase aberrations that are successful in converging to a solution that passes the likelihood threshold, we see that the algorithm converges more for smaller phase aberrations or flux.However, in Table 3 we see the phase error of runs that converged successfully in the likelihood sense is high for the small flux case, which we can understand as not constraining enough.To ensure the gradient-based method will work we want the phase error to be smaller than a radian so the problem will be close to quadratic or some low-order polynomial.In future works we will expand this optimization step to allow convergence for finite bandwidth imaging.

B.3. Gradient descent
To go from the l 2 solution to the MAP estimator we optimize a function proportional to the log-posterior L(φ) := x,y (I − f P φ ) using gradient descent with the method to rapidly calculate the gradient as done by Fienup (1999), this method takes advantage of the Fourier transform in Equation 3 to calculate it with only 2 FFT operations

Figure 2 .
Figure2.Schematic diagram of the HCI method.Schematically in our method, we take short exposure images, reconstruct the phase aberrations for each image, calculate the statistic using the image and reconstructed phase and sample to correct for its bias.

Figure 3 .
Figure 3. Detection statistic (a) without and (b) with bias correction.The images were simulated for a 2m-telescope, at seeing condition r0 = 15cm, for a bandwidth of 10%, and the statistic was accumulated over 10 4 images.

Figure 4 .
Figure 4. Injection-Recovery simulation.(a) The recovered ratio for different injected SNR per image, injection at |q| = 0.2 , we can see a plateau as we go for weaker signals, (b) recovery ratio for different separations at SNRinj = 1 per image.

Figure 5 .
Figure 5.Comparison between a (a) sampled PSF and the (b) the PSF corresponding to the phase measured by the SHWFS.An obvious, order-unity difference is visible between the two PSFs, therefore the measured phase is not accurate enough to be used for detection.

Table 1 .
Reference parameters for simulations.

Table 3 .
Phase root-mean-square for different seeing conditions and photons per image.The phase root-mean-square was averaged for runs that passed the likelihood threshold.