© The Authors 2023

1 Introduction

Ground-based exoplanet imagers typically consist of an adaptive optics (AO) system and a coronagraph. The AO is normally used to maximize the Strehl ratio in the science image, while the coronagraph minimizes diffracted starlight in the image plane. We refer to the focal-plane region where the coronagraph minimizes diffracted starlight as the “dark hole”. In the case of certain coronagraphs, such as apodizing phase plate (APP) coro-nagraphs, this can be a specific region within the focal plane; whereas, in the case of other coronagraphs such as Lyot coron-agraphs, this can refer to the entire focal plane. Here we define the Strehl ratio as the ratio of the intensity of the central stellar point spread function (PSF) core to the diffraction-limited intensity at the center of the PSF without the coronagraph; our definition, therefore, includes the reduction of the Strehl function due to a pupil-plane coronagraph. We define the raw contrast as the ratio of the average stellar intensity within the dark hole to the intensity of the central stellar PSF core. Within the dark hole, the exoplanet can be imaged with high contrast. The review papers by Ruane et al. (2018), Jovanovic et al. (2018), and Snik et al. (2018) summarize the various tools and techniques used in cutting-edge high-contrast exoplanet imagers. If we directly use the AO to maximize the dark-hole raw contrast rather than to minimize residual phase error, we can combine the AO and a pupil-plane coronagraph into a real-time, adaptive, pupil-plane phase coronagraph, which can achieve much better contrast at the cost of a slightly reduced Strehl ratio (Radhakrishnan et al. 2018). An AO control strategy based on this adaptive coron-agraph concept would require no additional optics and would provide the best contrast given the limited degrees of freedom of the deformable mirror (DM) in the AO system.

The maximum achievable contrast is further limited by differential wavefront aberrations between the science and the wavefront-sensor (WFS) optical paths. These quasi-static aberrations are referred to as non-common path aberrations (NCPA). Information about the NCPA must be obtained from the science image to estimate and compensate for them. The science image is typically integrated over a period of several milliseconds to seconds, during which the speckles caused by residual atmospheric turbulence rapidly evolve. Most NCPA retrieval techniques rely on the assumption that the average phase of the residual turbulence is zero during the science-image integration time. If the AO is used to flatten the wavefront, this is a reasonable assumption to make, as the residual turbulence after AO correction consists of low-amplitude high spatial-frequency components with a zero-mean distribution. This allows algorithms such as those based on phase-diversity phase retrieval (Mugnier et al. 2006) to estimate the NCPA phase from the integrated focal-plane image. These algorithms make use of one or more focal-plane image(s) created with precisely known added phase aberrations in addition to the science image to estimate the NCPA.

In the case of the adaptive coronagraph described above, however, the average phase during the integration time is not zero. This is because the DM adds phase offsets to the wave-front specifically to maximize the contrast within the dark hole. The average phase offset is not zero and results in an inaccurate NCPA phase estimate when using standard phase-diversity approaches, such as those described in Keller et al. (2012) and Por & Keller (2016).

To accurately estimate the NCPA phase using phase-diversity approaches, we have to account for 1) the nonzero average residual phase and 2) the averaging of the PSF over time. In the following sections, we derive an analytical expression for the time-averaged post-coronagraphic PSF in the presence of residual atmospheric turbulence, in analogy to the calculations in Sivaramakrishnan et al. (2002). We then calculate the components of the PSF caused by averaging over many different turbulence phase screens and remove them to debias the corresponding focal-plane images. We then use classical phase diversity on the debiased images. This gives us an accurate estimate of the NCPA phase, which we can feed back into our contrast control to achieve a significantly improved contrast.

Section 2 gives an overview of the various techniques currently in use to estimate and correct NCPA. We then explain the complications arising from the nonzero average phase, for example when using the AO to maximize the contrast in a specific part of the field of view. Section 3 presents the derivation of the analytical expression for the post-coronagraphic long-exposure PSF, and the use of phase diversity to estimate NCPA. Section 4 describes our simulation environment and presents the results of closed-loop contrast control with NCPA correction. Finally, Sect. 5 discusses our results and potential issues that may arise when implementing our new approach in an instrument.

2 Overview of NCPA retrieval algorithms

Mitigating the effects of the NCPA involves measuring the electric field in the focal plane. Once measured, the DM can be driven to minimize the electric field in the dark hole region of the focal plane, using techniques such as energy minimization (Malbet et al. 1995), speckle nulling (Bordé & Traub 2006), electric field conjugation (EFC; Give’on et al. 2007), stroke minimization (Pueyo et al. 2009), nonlinear dark hole control (Paul et al. 2013b), or speckle minimization with a self-coherent camera (Mazoyer et al. 2014; Singh et al. 2019; Galicher et al. 2019). Spatial linear dark field control (Miller et al. 2021; Bos et al. 2021) is another technique that uses the approximately linear response of some bright regions of the focal-plane intensity to the DM action.

With regard to measuring the electric field in the focal plane, several promising techniques have been developed. Some of these methods measure the focal-plane electric field by modulating the deformable mirror (Thomas et al. 2010; Kasper et al. 2013) and measuring its response in the focal plane. This results in a disruption of the science-image acquisition. Repeated disruption of the acquisition for focal-plane electric-field sensing reduces the effective integration time that can be spent on a given target before the contrast degrades and the focal-plane electric field has to be measured again. The Fast & Furious algorithm (Keller et al. 2012; Korkiakoski et al. 2014) has the advantage of estimating and compensating for the NCPA in a few iterations without disrupting the science image acquisition. Another method is to add a portion of the starlight to interfere with the residual starlight in the focal plane, such as in the Self Coherent Camera (Baudoz et al. 2005).

Some techniques have been developed to measure the focal-plane electric field without disturbing the science image. For example, Martinache (2013) used kernel phases to estimate the focal-plane electric field, provided that the wavefront errors are small (<1 radian) and that the pupil has some asymmetry. Another technique is millisecond imaging (Rodack et al. 2021), which uses short-exposure images to effectively freeze the turbulent layers in the atmosphere, enabling the use of WFS telemetry to probe the optical system for aberrations. The use of pupil-plane holograms (Por & Keller 2016) allows for several phase probes to be simultaneously applied to holographic copies of the science PSF, which enables phase-diversity-based electric field measurements without modulating the DM. Special coron-agraphic designs can be used to create single focal-plane images with sufficient information to perform phase retrieval, for example as was done by Wilby et al. (2017). Promising work is being done by Bos et al. (2019) with vector-APP coronagraphs in combination with asymmetric apertures to produce science images that can directly be used to estimate wavefront aberrations. This eliminates the need for focal-plane images produced by multiple beams with different diversity phases, thereby maximizing the number of useful photons in the science image. Another approach uses the residual AO speckles themselves as phase probes to estimate the NCPA, as in phase-sorting interferome-try (Codona et al. 2008; Codona & Kenworthy 2013), instead of using a separate phase diversity image. However, this requires science images with integration times that are not much longer than the WFS integration time.

One of the most straightforward methods to estimate the NCPA in a coronagraphic high-contrast imaging system is the COronagraphic Focal-plane wave-Front Estimation for Exoplanet detection (COFFEE; Sauvage et al. 2012; Paul et al. 2013a,b). As with other phase-diversity-based methods, COFFEE makes use of two (or more) focal-plane images generated downstream of a coronagraph to estimate the NCPA. COFFEE has been shown to accurately estimate and compensate for higher-order quasi-static aberrations in the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument at the Very Large Telescope (VLT; Paul et al. 2014b). COFFEE can be used to simultaneously estimate both the phase and amplitude aberrations (Paul et al. 2014a; Herscovici-Schiller et al. 2018).With the exception of COFFEE, the other NCPA retrieval methods described above do not take the residual turbulence into account and do not account for a non-zero residual phase arising from a different AO control strategy. Herscovici-Schiller et al. (2019) have shown that COFFEE can also be used when residual wavefront aberrations are present due to partial AO correction of atmospheric turbulence, and the focal-plane images are integrated over evolving residual turbulence. However, Herscovici-Schiller et al. (2019) make the assumption that the turbulence is stationary and ergodic, and do not take nonzero average residual turbulence into account. In this article, we describe a method to accurately measure and compensate for the NCPA in the case of a pupil-plane coronagraphic imaging system. Our method makes no assumptions about the residual turbulence. We use the pupil-plane measurements provided by the WFS to calculate the effects of the residual wavefront aberrations and their averaging over time, and remove these effects from the integrated science image.

3 NCPA phase retrieval in nonzero average phase scenarios

The techniques described above that estimate NCPA from longexposure coronagraphic images make the assumption that the average phase of the residual turbulence is zero over the integration time of the science image. This is a reasonable assumption to make when driving the AO to minimize phase aberrations in a long-exposure science image. However, when driving the AO to directly optimize the contrast, the DM adds phase offsets to the wavefront to keep the contrast high, and the average effect of these offsets is not zero. The same issue occurs when the focal plane images are integrated over a relatively short time, that is not long enough for the phase residuals to average out sufficiently.

In the following sections, we present a technique to estimate the NCPA from long-exposure coronagraphic images without making any assumptions about the residual turbulence. Our approach is tailored to imaging systems with pupil-plane coro-nagraphs, such as the APP coronagraph (Codona et al. 2006).

3.1 Derivation of analytical expression for the PSF

For our derivation, we assume that the science focal-plane image is integrated over a certain time δt and that we perfectly know the amplitude transmission in the pupil plane. We use two focal-plane images: one in focus and one with a known defocus that is added before the coronagraph. We chose this method as a proof of concept; however, it is not an optimum method of NCPA estimation because half the photons are not used in the science focal plane. We chose an in-focus and an out-of-focus beam to provide the phase diversity for the sake of simplicity; interested readers can refer to Lee et al. (1999), who explored optimal phase diversity, although in a noncoronagraphic context. In Sect. 5 we discuss alternative methods.

When estimating the phase from focal-plane images, we have to account for the fact that the residual phase averaged over the integration time is not zero. This nonzero average phase combined with the phase residuals must be accounted for in the in-focus and out-of-focus images before phase diversity can be applied. We show that the effects of residual turbulence do not allow us to simply assume that standard phase diversity estimates the sum of the actual NCPA and the average residual phase.

In the following, we derive an estimate for the turbulence-induced contributions to the coronagraphic PSF by extending the work done by Sivaramakrishnan et al. (2002); the mathematical details are virtually identical, but the interpretation is different. The coronagraph is assumed to be a shaped pupil or an apodiz-ing phase plate. We continuously record the phases ϕ_t measured by the WFS N times during the integration time of the science image. In contrast to Sivaramakrishnan et al. (2002) we cannot assume that the time-averaged phase is zero. We therefore define a quasi-static, complex aperture that includes the average phase as well as the amplitude and phase of the pupil coronagraph. The average of the recorded WFS phases is . The quasi-static complex aperture is then given by (1)

We note that A_i is the aperture illumination function, which is typically 1 within the aperture and 0 outside; A_c and ϕ_c are the amplitude and phase additions by the coronagraph, respectively; and ϕ_d is the known defocus added to the phase in the diversity image. For an apodizing phase-plate coronagraph A_c = 1, ϕ_c ≠ 0, and for a shaped pupil A_c is 0 or 1 depending on the position in the pupil and ϕ_c = 0.

We can now closely follow Sivaramakrishnan et al. (2002) by replacing their static, real aperture with our quasi-static, complex aperture A. We also need to subtract the average phase ϕ_avg (over the science image integration time) from each phase ϕ_t since the average phase is already included in the quasi-static complex aperture. We refer to these mean-subtracted phases as ψ_t where (2)

At a given time t, the instantaneous pupil-plane electric field after the coronagraph is given by (3)

The time-averaged focal-plane image can then be approximated as the sum of the quasi-static coronagraphic PSF for all ϕ_t = 0, a background “halo”erm P_halo, and a term corresponding to the speckles pinned to the bright rings of the PSF P_spec: (4)

Here, P_sci is the integrated science image. Using the Fraunhofer propagation operator (Goodman 2005), which is proportional to the Fourier transform operator ℱ, we can write the three terms as (5) (6) (7)

We note that ℜ indicates the real part of the complex value and the asterisk refers to the complex conjugate.

For the phase diversity image, we perform exactly the same calculations but using the quasi-static complex aperture for the out-of-focus image given by Eq. (1). This leads to the corresponding out-of-focus phase diversity post-coronagraphic PSF P_div for ϕ_t = 0 as well as the out-of-focus speckle and halo terms. This analytical expression is similar to the one derived in Herscovici-Schiller et al. (2017), but we found the expression in Sivaramakrishnan et al. (2002) easier to implement and interpret.

3.2 Phase-diversity-based phase retrieval

Throughout the integration time of the science camera, we record the phase measurements ϕ_t provided by the WFS. This allows us to calculate ψ_t, and hence P_halo and P_spec for both the in-focus and out-of-focus images. We then remove the halo and speckle terms from the in-focus image P_sci and the out-of-focus image P_div using Eq. (4): (8) (9)

The resulting debiased images are P₀ and ^divP₀, respectively. We then estimate the NCPA ϕ_ncp in terms of coefficients of a truncated mode basis similar to COFFEE (Sauvage et al. 2012), (10)

where z_j are the first M Zernike polynomials and a_j are their corresponding coefficients. In our simulations, we chose M to be 10, since we make the assumption that the NCPA are comprised of only low-order Zernike modes (up to the third radial order in this case).We exclude the piston term z₀ in our reconstruction because it has no effect on the focal-plane intensities. We model the images formed with the estimated aberration as follows: (11) (12)

The phase retrieval routine minimizes the difference between the debiased images and the modeled in-focus and out-of-focus images. We specifically minimize (13)

where a has the individual mode amplitudes a_c as its components; also, P₀ and P_0div are the cleaned in-focus and out-of-focus images, respectively, obtained from reading out the focal-plane camera and subtracting the effects of the average phase and the residual turbulence as described above. Furthermore, and are the estimated in-focus and out-of-focus images. Similar to the approach in COFFEE, we use a limited memory variable metric Broyden-Fletcher-Goldfarb-Shanno (BFGS) method to perform the minimization of J(a) (Press et al. 2007). We start with an initial estimate of zero for all a_j. The resulting optimum coefficients describe the NCPA in terms of Zernike polynomials.

One might expect that one could directly estimate ϕ_avg + ϕ_NCPA instead of incorporating ϕ_avg into the PSF models and taking it into account when removing the effects of turbulence. However, the structure of these two terms is vastly different. On the one hand, ϕ_NCPA is dominated by low-order optical aberrations that can be described by a few free parameters. On the other hand, ϕ_avg is dominated by high-order wavefront modes due to residual turbulence and contrast control and therefore has a large number of degrees of freedom. Since we do know ϕ_avg, we include it in the model instead of trying to retrieve it along with ϕ_NCPA.

Fig. 1 APP phase addition and corresponding PSF. Left: phase added to the wavefront by the APP. Axes are in units of meters on the DM, where the DM diameter is 0.0105 m. The color scale ranges from –π to π radians. Right: instantaneous PSF produced with a flat incident wavefront and this APP. The wavelength is 532 nm. The axes are in units of λ/D. The image is normalized by dividing by the maximum pixel intensity, and the color scale is in log10 units.

Fig. 2 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to minimize phase errors. There are no NCPAs in this simulation. The dashed gray lines are purely for comparison, as they correspond to time stamps at which NCPAs are estimated in later simulations. The dashed blue line indicates the average Strehl ratio and the dashed red line indicates the average raw contrast for the case when there are no NCPAs and the AO is driven to flatten the wavefront.

Fig. 3 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. There are no NCPAs in this simulation. The setup is the same as in Fig. 2, but in this case the DM is driven to maximize raw contrast at the dark hole instead of flatten the wavefront. The dashed gray lines are purely for comparison, as they correspond to time stamps at which NCPAs are estimated in later simulations. The dashed blue line indicates the average Strehl ratio and the dashed red line indicates the average raw contrast for the case when there are no NCPAs and the AO is driven to maximize the raw contrast.

Fig. 4 Optical and control setup used for simulations, described in Sect. 4. The wavefront is propagated through atmospheric turbulence, followed by a DM. The beam is then split, with half the light going to a WFS and the other half propagating through the APP coronagraph to the science camera. A phase aberration is added to the latter beam to simulate NCPA. The DM is controlled using a nonlinear control system as described in Radhakrishnan et al. (2018) to maximize focal-plane contrast rather than flatten the wavefront. The controller is driven based on inputs from the WFS as well as information from the focal plane.

Fig. 5 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to minimize phase errors. The imaging camera is read out every 200 ms as indicated by the vertical, dashed, gray lines. The science image and diversity image are used to directly estimate the NCPA, without first removing the halo and speckle terms from the images. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

4 Simulations

4.1 Simulation environment

The simulation environment is developed in hcipy, a Python framework for performing optical propagation simulations for high-contrast imaging applications (Por et al. 2018). The setup consists of a perfect telescope with an unobstructed circular aperture 4 m in diameter, a 10.5 mm diameter DM with 25 actuators across the diameter, a coronagraph, a perfect wavefront sensor with 35 subapertures across the diameter of the pupil plane, and a perfect imaging camera with 8 pixels per diffraction width. The DM influence functions have a Gaussian shape with the full width at half maximum of each actuator being equal to 1/25th of the diameter of the DM. Downstream of the DM is an APP coronagraph (Kenworthy et al. 2007). We chose to use an APP because as a pupil-plane optic, it is largely insensitive to tip and tilt errors as opposed to most focal-plane coronagraphs, and because it only modifies the phase of the wavefront. After the coronagraph, we split the beam and add a defocus of 0.4rad (RMS) to one of the beams. The beams are simultaneously imaged onto the science camera. We refer to the in-focus image as the science image and the out-of-focus image as the diversity image.

The APP coronagraph we use in our simulations creates a rectangular dark hole on one side of the PSF core. This dark hole extends from 2 to 8λ/D in the radial direction and from –3 to 3λ/D in the azimuthal direction. With a perfectly flat incident wavefront, the Strehl ratio of the instantaneous PSF is approximately 75% (due to the phase introduced by the APP coronagraph) and the contrast is approximately 2.1 × 10^–7. The phase addition to the wavefront by the APP and the corresponding PSF are shown in Fig. 1. We consider two sources of aberrations in our simulations - phase fluctuations due to residual, uncorrected atmospheric turbulence and static, low-order, phase-only NCPA introduced by the instrument. We simulate atmospheric turbulence with the von Karman statistical model (von Karman 1948), assuming a Fried parameter of 20 cm at 500 nm over the 4-m circular aperture, and an outer scale of 40 m. This is used to generate an evolving phase screen at a wavelength of 532 nm, which is propagated through the optical system. We assume frozen-flow atmospheric turbulence and use a wind speed of 8 m s^–1. The WFS is read out at a simulated frame rate of 1 kHz. The slower, imaging science camera is read out every 200 ms of simulated time. During the science exposure, the atmosphere evolves and the DM directly optimizes the focal-plane contrast based on the WFS measurements as described in Radhakrishnan et al. (2018).

In all our simulations, we simulate the NCPA by assuming that the aberrations are phase-only, are static within the time frames we are simulating, comprise of only astigmatism and coma modes, and the total phase aberration is capped at 0.5 rad RMS. We ran numerical simulations with 100 different realizations of the NCPA. In the following sections, we present the results for one of these realizations. We obtained similar results with all 100 realizations of the NCPA, with respect to the improvement in contrast and in Strehl ratio.

Fig. 6 Actual and estimated NCPA with wavefront flattening. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using phase-diversity-based phase retrieval on the science and diversity images without clean up. In this case the DM is driven to flatten the wavefront, i.e., minimize phase aberrations. Bottom: residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

Fig. 7 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup is the same as in the previous figure. The imaging camera is read out every 200 ms as indicated by the vertical, dashed, gray lines. The science image and diversity image are first cleaned by removing the halo and speckle terms and are then used to estimate the NCPA. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

Fig. 8 Actual and estimated NCPA with wavefront flattening and correcting for halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using phase-diversity-based phase retrieval on the science and diversity images after first removing the halo and speckle terms from them. As in the previous figure, the DM is driven to flatten the wavefront, i.e., minimize phase aberrations. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

Fig. 9 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to maximize the focal-plane contrast. NCPAs are estimated from the focal plane images every 200 ms (vertical, dashed, gray lines) without removing the halo and speckle terms. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

4.2 Simulation without NCPAs

For illustrative purposes, we first present the results of numerical simulations with this optical setup without any NCPAs. Figure 2 shows the evolution of the instantaneous Strehl ratio and the contrast over 700 ms of simulation time, when the DM is driven to flatten the wavefront. The science image in this case has a Strehl ratio of ≈64% and a raw contrast of 1.7 × 10^–4 at the dark hole. Figure 3 similarly shows the evolution of the instantaneous Strehl ratio and contrast when the DM is driven to maximize focal-plane dark-hole raw contrast rather than flatten the wavefront, as described in Radhakrishnan et al. (2018). The initially rapid improvement in contrast in this case is due to the convergence to a higher contrast while the contrast optimization loop closes. In this case, the resulting integrated science image has a Strehl ratio of ^63% and the dark hole has a raw contrast of 2.2 × 10^–5.

4.3 Including NCPA estimation in the adaptive coronagraph

The AO control system uses real-time measurements from the WFS to drive the DM. This contrast-control system drives the DM to directly optimize the estimated focal-plane contrast in the dark hole created by the coronagraph rather than minimize the wavefront error. We refer to this combination of AO and coronagraph working together as an adaptive coronagraph. The adaptive-coronagraph control system used here also includes the estimates from the NCPA phase retrieval routine, as shown in Fig. 4. The modified control system comprises two loops: a fast loop where the residual turbulence is measured by the WFS and the DM is driven to maximize focal-plane contrast for that particular turbulence, and a slow loop where the NCPA phase is measured and periodically fed to the AO controller to be included in the contrast optimizer. After every science and diversity camera exposure, we estimate the NCPA phase using the phase-diversity approach described in Sect. 3 and include this NCPA estimate in the contrast controller. The combination of the two loops ensures that both the dynamic aberrations due to turbulence and the quasi-static NCPA are accounted for by the controller when maximizing the contrast. This ensures the best possible contrast within the focal-plane dark hole for a given AO system.

Fig. 10 Actual and estimated NCPA with contrast control and not including halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when directly using phase-diversity-based phase retrieval on the science and diversity images without clean up. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

Fig. 11 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP coronagraph, and an imaging camera. The AO operates at 1 kHz and is driven to maximize the dark-hole contrast. NCPAs are estimated every 200 ms (indicated by the vertical, dashed, gray lines) using the technique described in Sect. 3. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

4.4 Inaccurate NCPA estimation

In this section, we illustrate the difference between NCPA estimation in the case where the average phase is zero versus the case where the average phase is not zero. Figure 5 shows the evolution of instantaneous Strehl and contrast when the DM is driven to flatten the wavefront without accounting for the effects of the residual turbulence. Every 200 ms the science camera is read out and the NCPA phase is estimated from the science image and diversity image, without taking into account the residual turbulence effects including correcting for the pinned speckles and the halo. The DM is updated to minimize the estimated NCPA phase as far as possible. We use a straightforward phase retrieval algorithm that uses the in-focus and out-of-focus images to estimate the phase, following Eq. (13). Each time we run the phase retrieval routine, we estimate the residual NCPA after DM correction and, within a few science camera read-out cycles, we are able to minimize the phase aberration. Since the DM is driven to flatten the wavefront, the average phase during the integration time of the focal-plane images is close to negligible, and the NCPA phase is estimated with high accuracy, as shown in Fig. 6. Removing the halo and speckle terms from the integrated science and diversity images before using them to estimate the NCPA results in a slightly more accurate estimate. The evolution of the instantaneous Strehl ratio and contrast in this case is shown in Fig. 7. Figure 8 shows the actual and estimated NCPA in this case. To illustrate the problem with NCPA estimation when the average phase is not zero, we first retrieve the phase from the focal-plane images without removing the speckle and halo terms. We run the simulation for 700 ms, running the slow loop every 200 ms. As is apparent in Fig. 9, the contrast initially improves by a little more than an order of magnitude and then at every subsequent estimation of the NCPA it degrades. The inclusion of the NCPA in the control loop results in the DM compensating for the NCPA while optimizing the focal plane contrast. In this case the phase additions by the DM do not minimize the NCPA phase, but rather minimize diffracted starlight within the dark hole in the presence of this NCPA. This results in a greater deviation from zero average phase and hence a stronger effect on the focal-plane images. Using the focal-plane images directly without including the effects of the nonzero average phase therefore results in an inaccurate NCPA phase estimate as shown in Fig. 10. The inaccurate phase estimate degrades contrast when included in the adaptive coronagraph controller, thus closed-loop control even becomes unstable.

Fig. 12 Actual and estimated NCPA with contrast control and including the halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using the cleaned focal-plane images. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

Fig. 13 Evolution of the instantaneous Strehl ratio and the average dark-hole contrast with a 5% broad band imaging. The AO is driven to maximize focal-plane contrast. NCPAs are estimated every 200 ms (indicated by the vertical, dashed, gray lines) using the technique described in Sect. 3. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

Fig. 14 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time when a planet that is five orders of magnitude fainter than the host star is present in the dark hole.

4.5 Accurate NCPA estimation

The results shown in the rest of this paper come from a single simulation run that covers 700 ms of simulated time. The evolving atmospheric turbulence results in dynamic phase aberrations that are measured by the WFS and recorded every 1ms. These measurements include the effects of the DM, which is driven to optimize focal-plane contrast instead of flattening the wave-front. Every 200 ms, the science camera is read out, and the recorded WFS measurements are used to remove the effects of the nonzero average residual phase on the science image and diversity image. The resulting images are used to estimate the NCPA as described in Sect. 3. Figure 11 shows that the contrast improves by almost an order of magnitude once the NCPA is estimated and accounted for by the controller. Similarly, the Strehl ratio increases by almost 50%. We note that the maximum achievable Strehl ratio is around 75% because of the APP. The dashed vertical lines in the figure represent the instances at which the science camera is read out and the NCPA is estimated and compensated for. Figure 12 shows the actual NCPA, the estimated NCPA after 700 ms of simulation time, and the residual NCPA.

Fig. 15 Science camera readouts with planet included. Top-left to top-right to bottom in sequence: in-focus image at each readout of the science camera. The planet is five orders of magnitude fainter than the host star, and is visible within the dark hole. It is indicated here with an arrow. The wavelength is 532 nm. The axes are in units of A/D. The image is normalized by dividing by the maximum pixel intensity, and the color scale is in log10 units.

4.6 NCPA estimation in a broadband setting

So far we have demonstrated NCPA estimation with monochromatic wavelengths, but the same algorithm can also be applied in the case of a finite bandwidth, with few modifications. Figure 13 shows the results of a numerical simulation in which we increased the bandwidth to 5% centered at 532 nm. We assume here that the turbulence phase and the NCPA phase do not change in structure with a change in wavelength. In our model we use the central wavelength of 532 nm to calculate the halo and speckle terms. In the case of the APP, the PSF scales with wavelength. With a broader wavelength band, this effect must be taken into account in the model, for example by taking the weighted average of scaled versions of the PSF. However, in the case of a 5% bandwidth, using just the central wavelength in the model produces satisfactory results. We see the results begin to deteriorate once the bandwidth is increased to 20%.

4.7 NCPA estimation when an exoplanet is present

Another case of interest is when an exoplanet is present within the dark hole. Since the speckles of diffracted starlight are typically two to three orders of magnitude brighter than an exoplanet (Traub & Oppenheimer 2010), the same procedure can be used to estimate the NCPA without the exoplanet light causing a difference between the actual and modeled focal-plane images. We ran a numerical simulation with an exoplanet at 10⁻⁵ times the intensity of the host star. In this case, the simulated exoplanet was as bright as the speckles themselves. We were still able to accurately estimate and account for the NCPA even with the planet present. Figure 14 shows the evolution of the instantaneous Strehl and the contrast during the simulation, and Fig. 15 shows the in-focus image at each readout of the science camera.

5 Conclusions and discussion

Residual atmospheric turbulence in a high-contrast imaging system can have a significant impact on NCPA estimation from focal-plane images if the average residual turbulence phase is not zero during the image integration time. This can occur if the AO system is not driven to flatten the wavefront, but controlled in some other way, for example to maximize the focal-plane contrast (Radhakrishnan et al. 2018). Another situation where this may occur is when estimating NCPAs from focal-plane images that have not been integrated for a long enough duration for the average residual phase to be sufficiently close to zero. Most techniques that involve phase-diversity methods of NCPA estimation result in inaccurate estimates if directly applied to the focal-plane images. Here we have shown that it is possible to model the effects of the nonzero average phase on the focal-plane images and correct for them. This paves the way for accurate NCPA estimation in nonzero-average phase situations.

Estimating the NCPA in this way comes with certain computational demands. The data measured by the WFS must be recorded and stored during the integration time of the science camera. The calculation of the speckle term involves calculating the sum of squares of the measured wavefront phases. This can be calculated in parallel with the WFS readouts. A running average of the post-WFS phase can also be calculated in parallel with the WFS. However, the calculation of the halo term and the complex aperture can only be performed when all the WFS measurements made during the integration time of the focal-plane images are available, and therefore can only be performed after reading out the science camera. Between camera readouts, we calculate the speckle and halo terms, use them to clean the focal-plane images, and then use the L-BFGS-B routine to estimate the NCPA phase from the cleaned images. These calculations must therefore be completed within the integration time of the science camera, which typically ranges from hundreds of milliseconds to seconds. When running the simulations on an Intel Core i5 6440HQ CPU, the calculations required to clean the speckle and halo terms took 13 s of CPU time on average, and the L-BFGS-B optimization to fit the estimated NCPA took 6 s of CPU time on average. Running the simulations in the Google Colaboratory¹ environment reduced these times to 10 s and 3 s, respectively.

Cleaning the integrated camera images is the most time-consuming step since it involves the Fraunhofer propagation of a large amount of data measured by the WFS during the integration time. Instead of using the full Fraunhofer propagator, we could potentially speed up this step by calculating a linear transformation that maps the pupil plane to the central airy core and the dark hole of the focal plane and use this instead of the Fraunhofer propagator. We would then clean only a subset of the focal plane images, which might result in a trade-off between accuracy of the estimated NCPA and computation time. Careful hyperparameter tuning of the optimizer, or maybe even using a different optimizer might also improve computation time. If the NCPA does not change significantly over the camera readout time, we can run the optimizer for more iterations, or alternatively run the optimizer for just one or two iterations per camera read-out. These potential improvements to speed are outside the scope of this paper.

Modern large and extremely large ground-based telescopes will be equipped with multiple graphics processing units (GPUs) for real-time AO control. Software packages such as CACAO (Guyon et al. 2018) and COMPASS (Ferreira et al. 2018) enable such multi-GPU systems to process WFS data and drive the AO at kilohertz frequencies. Estimating the NCPA and driving the AO in the slow loop as shown in Fig. 4 is therefore well within the capability of modern AO control systems.

The type of phase diversity methods discussed here have the disadvantage that half of the light is going to the out-of-focus image that does not contribute to the science image. Therefore, more research on techniques for sensing the focal-plane electric field is needed. Two promising avenues are 1) specially designed vector APP coronagraphs used in combination with asymmetric apertures (Bos et al. 2019) to enable NCPA estimation from a single focal-plane image, and 2) Fast & Furious based techniques which can estimate the focal-plane electric field in just two to three iterations in closed loop (Keller et al. 2012; Korkiakoski et al. 2014; Bos et al. 2020).

Acknowledgements

The work presented in this paper was supported by the Netherlands Organisation for Applied Scientific Research (TNO) and Leiden Observatory. The authors thank Sebastiaan Haffert for his help with the “hcipy” (Por et al. 2018) Python package for high contrast imaging simulations, and David Doelman for the APP coronagraph design used in the simulations (Doelman et al. 2017).

References

All Figures

Fig. 1 APP phase addition and corresponding PSF. Left: phase added to the wavefront by the APP. Axes are in units of meters on the DM, where the DM diameter is 0.0105 m. The color scale ranges from –π to π radians. Right: instantaneous PSF produced with a flat incident wavefront and this APP. The wavelength is 532 nm. The axes are in units of λ/D. The image is normalized by dividing by the maximum pixel intensity, and the color scale is in log10 units.

In the text

Fig. 2 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to minimize phase errors. There are no NCPAs in this simulation. The dashed gray lines are purely for comparison, as they correspond to time stamps at which NCPAs are estimated in later simulations. The dashed blue line indicates the average Strehl ratio and the dashed red line indicates the average raw contrast for the case when there are no NCPAs and the AO is driven to flatten the wavefront.

In the text

Fig. 3 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. There are no NCPAs in this simulation. The setup is the same as in Fig. 2, but in this case the DM is driven to maximize raw contrast at the dark hole instead of flatten the wavefront. The dashed gray lines are purely for comparison, as they correspond to time stamps at which NCPAs are estimated in later simulations. The dashed blue line indicates the average Strehl ratio and the dashed red line indicates the average raw contrast for the case when there are no NCPAs and the AO is driven to maximize the raw contrast.

In the text

Fig. 4 Optical and control setup used for simulations, described in Sect. 4. The wavefront is propagated through atmospheric turbulence, followed by a DM. The beam is then split, with half the light going to a WFS and the other half propagating through the APP coronagraph to the science camera. A phase aberration is added to the latter beam to simulate NCPA. The DM is controlled using a nonlinear control system as described in Radhakrishnan et al. (2018) to maximize focal-plane contrast rather than flatten the wavefront. The controller is driven based on inputs from the WFS as well as information from the focal plane.

In the text

Fig. 5 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to minimize phase errors. The imaging camera is read out every 200 ms as indicated by the vertical, dashed, gray lines. The science image and diversity image are used to directly estimate the NCPA, without first removing the halo and speckle terms from the images. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

In the text

Fig. 6 Actual and estimated NCPA with wavefront flattening. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using phase-diversity-based phase retrieval on the science and diversity images without clean up. In this case the DM is driven to flatten the wavefront, i.e., minimize phase aberrations. Bottom: residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

In the text

Fig. 7 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup is the same as in the previous figure. The imaging camera is read out every 200 ms as indicated by the vertical, dashed, gray lines. The science image and diversity image are first cleaned by removing the halo and speckle terms and are then used to estimate the NCPA. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

In the text

Fig. 8 Actual and estimated NCPA with wavefront flattening and correcting for halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using phase-diversity-based phase retrieval on the science and diversity images after first removing the halo and speckle terms from them. As in the previous figure, the DM is driven to flatten the wavefront, i.e., minimize phase aberrations. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

In the text

Fig. 9 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP corona-graph, and an imaging camera. The AO operates at 1 kHz and is driven to maximize the focal-plane contrast. NCPAs are estimated from the focal plane images every 200 ms (vertical, dashed, gray lines) without removing the halo and speckle terms. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

In the text

Fig. 10 Actual and estimated NCPA with contrast control and not including halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when directly using phase-diversity-based phase retrieval on the science and diversity images without clean up. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.

In the text

Fig. 11 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time. The setup consists of an AO system, an APP coronagraph, and an imaging camera. The AO operates at 1 kHz and is driven to maximize the dark-hole contrast. NCPAs are estimated every 200 ms (indicated by the vertical, dashed, gray lines) using the technique described in Sect. 3. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

In the text

	Fig. 12 Actual and estimated NCPA with contrast control and including the halo and speckle terms. Top-left: actual NCPA phase aberration. Top-right: estimated NCPA phase after 600 ms when using the cleaned focal-plane images. Bottom: Residual phase aberration after 600 ms. The residual aberration has been magnified by ten to be more visible in the figure.
In the text

Fig. 13 Evolution of the instantaneous Strehl ratio and the average dark-hole contrast with a 5% broad band imaging. The AO is driven to maximize focal-plane contrast. NCPAs are estimated every 200 ms (indicated by the vertical, dashed, gray lines) using the technique described in Sect. 3. The dashed blue and red lines indicate the average Strehl ratio and contrast, respectively, for the case when there are no NCPAs.

In the text

	Fig. 14 Instantaneous Strehl ratio and the average dark-hole contrast as a function of time when a planet that is five orders of magnitude fainter than the host star is present in the dark hole.
In the text

Fig. 15 Science camera readouts with planet included. Top-left to top-right to bottom in sequence: in-focus image at each readout of the science camera. The planet is five orders of magnitude fainter than the host star, and is visible within the dark hole. It is indicated here with an arrow. The wavelength is 532 nm. The axes are in units of A/D. The image is normalized by dividing by the maximum pixel intensity, and the color scale is in log10 units.

In the text