1.

Introduction

With thousands of exoplanets having been discovered primarily via indirect detection methods, including transit photometry and radial velocity measurements, further studies are desired to understand the dynamics of planet formation, interactions with debris disks and aid the discovery of potentially habitable worlds. These studies can be enabled if the challenges of direct imaging are overcome, of which the two primary challenges are the relatively small angular separation between the host star and the exoplanet as well as the small planet-to-star flux ratio. For perspective, a Sun–Earth analog at 10 parsecs will have an angular separation of $\sim 0.1$ arcsec and a flux ratio of $\sim {10}^{-10}$.¹ While coronagraphs and high-order wavefront control algorithms have been demonstrated under laboratory conditions with ground-based instruments,^2,3 atmospheric turbulence, telescope stability, and absorption in certain bandpasses limit ground-based performance. A key milestone for direct imaging will be the launch of the Nancy Grace Roman Space Telescope, which will carry an onboard coronagraph with multiple modes for imaging and spectral characterization.⁴

While acting as a technology demonstration, the Roman Coronagraph is expected to bridge the performance gap between current high-contrast imaging capabilities and what is expected for a future 6-m class Habitable Worlds Observatory (HWO) recommended by the decadal survey to image Earth-like exoplanets.⁵ Predicted to achieve contrasts on the order of ${10}^{-9}$,⁶ the Coronagraph Instrument is particularly focused on demonstrating two distinct coronagraph designs, low-noise electron multiplying detectors, and high-order wavefront sensing and control (HOWFSC) algorithms required to suppress speckles from quasi-static aberrations within the optical train. The primary method intended for the HOWFSC scheme has been the use of model-based electric field conjugation (EFC).⁷ Described in Sec. 2.1, EFC is combined with pairwise probing (PWP) to sense and minimize the focal plane electric field using an instrument model.

Currently, EFC is expected to achieve the desired contrasts for Roman. However, the contrast attained by EFC has been demonstrated to degrade, in some cases by an order of magnitude, due to model errors, particularly those associated with the deformable mirror (DM) translation and rotation.⁸ This is the motivation to study the alternative implicit EFC method⁹ where the Jacobian is constructed entirely from empirical data and inherently includes all perturbations in the optical system. For a space-borne telescope, implicit EFC (iEFC) can be a particularly valuable method as it can mitigate the risk associated with instrument performance in the presence of instrument uncertainties. In addition, iEFC does not require the focal plane electric field to be sensed prior to computing the actuator commands. This simplifies the controller because only a single regularization parameter needs to be optimized, whereas an implementation of EFC with PWP requires two regularization parameters, one for the initial electric field estimation during PWP and one for the conjugation using EFC.⁹

Here, iEFC is extended to use two DMs and proven capable of creating an annular dark hole such that it can be applied to the Roman Coronagraph as well as designs for future missions such as the HWO. The wide-field-of-view shaped pupil coronagraph (SPC-WFOV) mode is chosen because preliminary simulations demonstrated that iEFC is likely not suited to the hybrid Lyot coronagraph (HLC) due to the frequent relinearization of the Jacobian required for the HLC if the DM design patterns are not used¹⁰ or are inaccurate.¹¹ Specifically, the empirical calibration of the iEFC Jacobian would likely make repeated relinearizations unfeasible for the HLC because of the required duration of each Jacobian calibration.

Section 2 provides an overview of the EFC and iEFC algorithms, including how an arbitrary basis of modes can be used for each method. In addition, iEFC is extended to use two DMs and initially demonstrated with a scalar vortex coronagraph (SVC) model. Section 4 then includes a demonstration of EFC for the SPC-WFOV mode that illustrates how instrument perturbations unaccounted for in a model result in contrast degradation and a more complex control scheme. In Sec. 5, the efficacy of iEFC is tested with monochromatic simulations using multiple modal bases to demonstrate which is optimal for the SPC-WFOV instrument. Section 6 then expands the iEFC simulations for a single narrowband wavefront control bandpass and a broader 10% bandpass. Finally, Sec. 7 integrates flux estimates for a reference star and an electron multiplying charge-coupled device (EMCCD) model into the simulations to test the feasibility of iEFC for both bandpasses considered. The results demonstrate that iEFC can be a suitable method for the SPC-WFOV mode and can mitigate the risk of pupil misalignments after the launch of Roman. However, the disadvantage of iEFC will be noise accumulated during calibration and the integration time required to acquire all calibration frames. Given that EFC has been rigorously tested and simulated for the Roman Coronagraph and is expected to achieve the desired contrast goals, iEFC should likely be tested during the latter stages of the Roman mission to further study the efficacy of the method. This would inform future missions, ideally with less obstructed and larger apertures capable of capturing more photons to reduce calibration times and mitigate the effects of calibration noise.

2.

Definitions and Algorithms

For clarity between EFC and iEFC, Table 2 in Appendix A has a list of variables used for the formalism of both methods. Variables in bold font represent column vectors, scalars are all lowercase variables, and matrices and scalar fields are all uppercase variables. The precise shape of all matrices is also included in Table 2.

The images throughout these simulations are presented in units of normalized intensity (NI), defined as

2.1.

EFC

Following the formalism of EFC from Give’on et al.,⁷ the fundamental concept is that the irradiance in the image plane can be minimized by computing DM commands that destructively interfere with the measured electric field.⁷ This is done using a model of the instrument to estimate the electric field response of DM modes in the image plane. For a coronagraph, this electric field may be written as some linear operation $C[]$ acting on the pupil plane wavefront $A(x,y)e^{i\mathrm{\Phi }(x,y)}$ where $A(x,y)$ and $\mathrm{\Phi }(x,y)$ are the scalar fields corresponding to the pupil plane amplitude and phase aberrations, respectively.

Eq. (2)

$$E_{\mathrm{im}}(x,y)=C[A(x,y)e^{i\mathrm{\Phi }(x,y)}].$$

With coronagraphs comprising multiple relays between pupil planes and focal planes, $C[]$ is often a series of Fourier transforms propagating the wavefront between the transverse planes. For a system with a single DM, the electric field can be written as

Eq. (3)

$$E_{\mathrm{im}}(x,y)=C[A(x,y)e^{i\mathrm{\Phi }(x,y)}{e}^{i{\mathrm{\Phi }}_{DM}(x,y)}].$$

Here, ${\mathrm{\Phi }}_{\mathrm{DM}}(x,y)$ is the phase induced by the surface of the DM expressed as a sum of weighted influence functions $F_n(x,y)$ with the expression

Eq. (4)

$${\mathrm{\Phi }}_{\mathrm{DM}}(x,y)=\frac{4\pi }{\lambda }\sum _{n=1}^{N_{\mathrm{acts}}}{a}_n{F}_n(x,y).$$

Assuming that the DM surface imparts a small phase, the DM phasor can be approximated as $e^{i{\mathrm{\Phi }}_{\mathrm{DM}}(x,y)}\approx 1+i{\mathrm{\Phi }}_{\mathrm{DM}}(x,y)$, allowing for the following simplification to the image plane electric field

Eq. (5)

$$E_{\mathrm{im}}(x,y)=C[A(x,y)e^{i\mathrm{\Phi }(x,y)}]+iC[{\mathrm{\Phi }}_{\mathrm{DM}}(x,y)].$$

Once sampled by the science camera detector, $E_{\mathrm{im}}(x,y)$ is written as the concatenated vector ${\mathbf{E}}_{\mathbf{im}}$ containing the real and imaginary components of the electric field at each pixel in the dark hole. The individual terms for the electric field contributed by the system aberrations and the DM influence functions are now written as

Eq. (6)

$$C[A(x,y)e^{i\mathrm{\Phi }(x,y)}]={\mathbf{E}}_{\mathbf{ab}},$$

and

Eq. (7)

$$iC[{\mathrm{\Phi }}_{\mathrm{DM}}(x,y)]=G_1{\mathbf{A}}_{\mathbf{1}},$$

respectively. Here, $G_1$ is the Jacobian that transforms the actuator heights of DM1 (${\mathbf{A}}_{\mathbf{1}}$) into the electric field contributed by the DM at the focal plane. With this, the image plane electric field is written as the following system of linear equations.

Eq. (8)

$${\mathbf{E}}_{\mathbf{im}}={\mathbf{E}}_{\mathbf{ab}}+G_1{\mathbf{A}}_{\mathbf{1}}.$$

As presented by Give’on et al.,⁷ when a second DM is used outside a pupil plane, this system of linear equations is expanded to include the term $G_2{\mathbf{A}}_{\mathbf{2}}$. While similar to the Jacobian for the pupil plane DM, the Jacobian $G_2$ includes the angular spectrum propagation effects from the plane of the second DM back to the pupil plane. As a system of linear equations, the two Jacobians can be concatenated to capture the influence of all DM actuators such that the final equation for the image plane electric field is simplified to

Eq. (9)

$$\begin{gathered} {\mathbf{E}}_{\mathbf{im}}={\mathbf{E}}_{\mathbf{ab}}+G_1{\mathbf{A}}_{\mathbf{1}}+G_2{\mathbf{A}}_{\mathbf{2}}={\mathbf{E}}_{\mathbf{ab}}+G_{\mathrm{EFC}}\mathbf{A}, \\ \mathrm{where}{G}_{\mathrm{EFC}}=[G_1,G_2]\hspace{0.17em}\mathrm{and}\hspace{0.17em}\mathbf{A}=\left[\begin{array}{c}{\mathbf{A}}_{\mathbf{1}}\\ {\mathbf{A}}_{\mathbf{2}}\end{array}\right]. \end{gathered}$$

In the case where an arbitrary basis of DM modes is used instead of individual actuators, this can be written more generally in terms of the modal coefficients ${\mathbf{m}}_{\mathbf{c}}$. The basis of DM modes chosen are referred to as the calibration modes that are contained in the matrix $M_{\mathrm{modes}}$. Note that the shape of the Jacobian will depend on the amount of calibration modes used, and the equation for the image plane electric field is then

Eq. (10)

$${\mathbf{E}}_{\mathbf{im}}={\mathbf{E}}_{\mathbf{ab}}+G_{\mathrm{EFC}}{\mathbf{m}}_{\mathbf{c}}.$$

Figure 1 illustrates the typical practices for how EFC is calibrated and performed. Prior to the EFC control loop, the Jacobian is computed using the central difference approximation in Eq. (11) to estimate the response of the electric field in the image plane for each calibration mode. To do so, the positive and negative calibration modes are applied in a Fourier-based optical model with $a_{\mathrm{mode}}$ being the chosen calibration amplitude to calculate the fields referenced as ${\mathbf{E}}_{\mathbf{P}}$ and ${\mathbf{E}}_{\mathbf{N}}$, respectively. While it is possible to estimate the response using the nominal electric field and only applying a single positive or negative calibration mode, the central difference approximation decreases the error from nonlinear terms.¹³ Each response was stored in the Jacobian as a single vector ${\mathbf{R}}_{\mathbf{i},\mathbf{EF}}$ where $i$ denotes the index of the calibration mode.

Eq. (11)

$${\mathbf{R}}_{\mathbf{i},\mathbf{EF}}=\frac{{\mathbf{E}}_{\mathbf{P}}-{\mathbf{E}}_{\mathbf{N}}}{2a_{\mathrm{mode}}},$$

Eq. (12)

$$G_{\mathrm{EFC}}=[\begin{array}{cccccc}{\mathbf{R}}_{\mathbf{1},\mathbf{EF}}& {\mathbf{R}}_{\mathbf{2},\mathbf{EF}}& \dots & {\mathbf{R}}_{\mathbf{i},\mathbf{EF}}& \dots & {\mathbf{R}}_{{\mathbf{N}}_{\mathbf{modes}},\mathbf{EF}}\end{array}].$$

Fig. 1

Functional steps for model-based EFC. The Jacobian $G_{\text{EFC}}$ is constructed by applying positive and negative calibration modes within the model and calculating the electric field at the focal plane. The response/derivative is computed from the difference normalized by the amplitude of the mode applied ($a_{\mathrm{mode}}$). The EFC controller uses the pseudo-inverted Jacobian to compute the update to the DM actuators using the estimated electric field ${\mathbf{E}}_{\mathbf{ab}}$. All boxes in yellow indicate a step that is performed computationally.

During each iteration of EFC, the solution for modal coefficients to minimize the electric field is commonly found by minimizing a cost function such as

Eq. (13)

$$J=|{\mathbf{E}}_{\mathbf{ab}}+G_{\mathrm{EFC}}{\mathbf{m}}_{\mathbf{c}}{|}^2+\lambda |{\mathbf{m}}_{\mathbf{c}}{|}^2,$$

where $\lambda$ is a regularization parameter penalizing the use of large actuator strokes. Throughout the simulations presented, the “$\beta$ regularization” technique described by Sidick et al.¹⁴ is used for both EFC and iEFC simulations. With this technique, the regularization parameter chosen by the user is $\beta$. Once specified, the control matrix ($M_{\mathrm{control}}$) is computed numerically using the equation

Eq. (14)

$$M_{\mathrm{control}}=(G_{\mathrm{EFC}}^{⊺}{G}_{\mathrm{EFC}}+{\alpha }^2{10}^{\beta }I{)}^{-1}{G}_{\mathrm{EFC}}^{⊺}.$$

Here, $\alpha$ is the maximum singular value of the square matrix $G_{\mathrm{EFC}}^{⊺}{G}_{\mathrm{EFC}}$. The modal coefficients that minimize the cost function are calculated with the matrix-vector product $-M_{\mathrm{control}}{\mathbf{E}}_{\mathbf{ab}}$. These modal coefficients can be transformed into actuator heights using $M_{\mathrm{modes}}$ as demonstrated in the following equation:

Eq. (15)

$$\mathbf{A}=\left[\begin{array}{c}{\mathbf{A}}_{\mathbf{1}}\\ {\mathbf{A}}_{\mathbf{2}}\end{array}\right]=M_{\mathrm{modes}}{\mathbf{m}}_{\mathbf{c}}=-M_{\mathrm{modes}}{M}_{\mathrm{control}}{\mathbf{E}}_{\mathbf{ab}}.$$

A critical component to EFC is how the electric field in the focal plane is estimated prior to the actuator command being calculated. While the self-coherent camera is a method of estimating the focal plane electric field using interference from a reference channel,¹⁵ the common technique proposed for Roman is PWP.¹⁶ This method uses at least two linearly independent DM probe commands to generate phase diversity within the focal plane. For each probe command, two images are recorded, one with a positive probe and one with a negative probe. The difference images of each probe command are then used to estimate the electric field in each pixel of the desired dark hole by solving Eq. (16) for ${\mathbf{E}}_{\mathbf{ab}}$ where $\mathit{\delta }$ is the concatenated vector of difference images and $M_{\mathrm{probes}}$ is a matrix with the estimated electric field contributed by the DM probes. Solving this system of linear equations requires another pseudo-inverse in which the regularization parameter ${\lambda }_{\mathrm{PWP}}$ is used to make the algorithm robust against noise.

Eq. (16)

$$\mathit{\delta }=M_{\mathrm{probes}}{\mathbf{E}}_{\mathbf{ab}},$$

Eq. (17)

$${\mathbf{E}}_{\mathbf{ab}}=(M_{\mathrm{probes}}^{⊺}{M}_{\mathrm{probes}}+{\lambda }_{\mathrm{PWP}}I{)}^{-1}{M}_{\mathrm{probes}}^{⊺}\mathit{\delta }$$

Various other implementations of PWP have been developed, including extensions to broadband estimation¹⁷ and implementations of a Kalman filter for improved estimation,¹⁸ but all rely on the difference images from pairwise probes.

2.2.

iEFC

As described by Haffert et al.,⁹ iEFC includes no explicit calculation of the electric field within each iteration because the electric field is not the desired control variable, rather, difference images of probes are the control variable. The formalism of PWP demonstrates that difference images of two or more linearly independent probes are related to the focal plane electric field by a linear transformation making them a linear proxy for the electric field. Therefore, minimizing the difference images also minimizes the electric field.

Because the control variable for iEFC is measurable, a Fourier optics model of the instrument is not required to compute the Jacobian. Instead, iEFC is calibrated by measuring the response ${\mathbf{R}}_{\mathbf{i},\delta }$ for each calibration mode. Here, the response measured for a single mode is another central difference of the measurements for both the positive and negative calibration modes. As explained in Eq. (19), the individual measurements themselves are the difference images ${\mathit{\delta }}_{+j}$ and ${\mathit{\delta }}_{-j}$. The value of $j$ references which probe the images are acquired for, $+$ and $-$ reference the sign of the calibration mode, while $P$ and $N$ reference the sign of the probe command. Each measured response is now a double difference that is stored in the iEFC Jacobian $G_{\mathrm{IEFC}}$.

Eq. (18)

$${\mathit{\delta }}_{+j}=\frac{{\mathbf{I}}_{+jP}-{\mathbf{I}}_{+jN}}{2a_{\mathrm{probe}}}\mathrm{and}{\mathit{\delta }}_{-j}=\frac{{\mathbf{I}}_{-jP}-{\mathbf{I}}_{-jN}}{2a_{\mathrm{probe}}},$$

Eq. (19)

$${\mathit{\delta }}_{+}=\left[\begin{array}{c}{\mathit{\delta }}_{+1}\\ {\mathit{\delta }}_{+2}\end{array}\right]\mathrm{and}{\mathit{\delta }}_{-}=\left[\begin{array}{c}{\mathit{\delta }}_{-1}\\ {\mathit{\delta }}_{-2}\end{array}\right],$$

Eq. (20)

$${\mathbf{R}}_{\mathbf{i},\delta }=\frac{{\mathit{\delta }}_{+}-{\mathit{\delta }}_{-}}{2a_{\mathrm{mode}},}$$

Eq. (21)

$$G_{\mathrm{IEFC}}=[\begin{array}{cccccc}{\mathbf{R}}_{\mathbf{1},\delta }& {\mathbf{R}}_{\mathbf{2},\delta }& \dots & {\mathbf{R}}_{\mathbf{i},\delta }& \dots & {\mathbf{R}}_{{\mathbf{N}}_{\mathbf{modes}},\delta }\end{array}].$$

Figure 2 illustrates the procedures for calibrating iEFC and closing the loop. Similar to EFC, the cost function can be defined such that the difference images will be minimized. This cost function can be used to then compute a control matrix for $M_{\mathrm{control}}$. Each iteration of iEFC then consists of measuring the difference images $\mathit{\delta }$ and solving for the modal coefficients by multiplying by $M_{\mathrm{control}}$.

Eq. (22)

$$J=|\mathit{\delta }+G_{\mathrm{IEFC}}{\mathbf{m}}_{\mathbf{c}}{|}^2+\lambda |{\mathbf{m}}_{\mathbf{c}}{|}^2.$$

Fig. 2

Critical steps for generating a Jacobian for iEFC by calibrating a chosen set of DM modes using difference images of probes concatenated into the vectors ${\mathit{\delta }}_{+}$ and ${\mathit{\delta }}_{-}$. By inverting the Jacobian and measuring $\mathit{\delta }$ at each iteration, the difference images are minimized to reduce the electric field amplitude. Note that each measurement of $\mathit{\delta }$ is also normalized by the amplitude of the probe ($a_{\mathrm{probe}}$).

The two-DM iEFC controller is initially tested using an SVC model created with the High Contrast Imaging in Python¹⁹ software to verify the capability of creating an annular dark hole. The model uses a charge 6 vortex with a 90% diameter Lyot stop. The wavelength for the simulation is 650 nm with two $34\times 34$ actuator DMs. The DM pupil is 10 mm in diameter with a 300-mm separation between the DMs (roughly a Fresnel number of 512). Here, the control region spans from $2\lambda /D$ to $10\lambda /D$. The calibration modes chosen are a basis of Fourier modes that span frequencies from $1.5\lambda /D$ to $12\lambda /D$. The probes used are two single actuator pokes of adjacent actuators that are indicated by white circles in Fig. 3 on DM1 (bottom left). The choice of Fourier modes, sampling, and probes will be expanded upon in Sec. 5. While this result is for an ideal monochromatic vortex coronagraph, the application of iEFC with two DMs is found to be suitable for annular dark holes.

Fig. 3

Annular dark hole created using iEFC with an SVC model where the final mean contrast is $4.94\times {10}^{-10}$ after 30 iterations. The two white circles on the first DM indicate the actuators that are used as probes.

3.

Roman SPC Physical Optics Model

The HOWFSC simulations performed here make use of an end-to-end physical optics model for the Roman Coronagraph. The model was created using POPPY²⁰ as the backend physical optics propagation software, but the prescription and data are the same as those in the roman_phasec_proper (v1.4) package.²¹ The motivation for choosing POPPY was to leverage end-to-end simulations entirely on a graphics processing unit (GPU) as CuPy²² was recently added as a computation feature in POPPY. The GPU computations greatly reduce simulation times and enable a more extensive investigation given the large number of computations required to simulate an iEFC Jacobian. When compared with the results of the PROPER model, the mean NI of the POPPY model was found to be within a percent for an unocculted point-spread function and within a few percent for a coronagraphic image. More comparisons of the speckles and image morphology have been done in previous works.²³

Within the physical optics model, each optic includes a surface error map that is applied to the wavefront during propagation. In all EFC and iEFC simulations, the initial state of the Fresnel model includes a DM1 flatmap correcting the pupil aberrations that improves the initial contrast for HOWFSC by more than an order of magnitude. Figure 4 demonstrates the initial images from the model prior to any HOWFSC.

Fig. 4

Initial state of the physical optics model prior to any HOWFSC. The image of the unocculted on-axis source (a) is the reference image to which coronagraph frames are normalized. Panel (b) is the initial coronagraphic image without correction by a DM1 flatmap followed by the image on panel (c) using the DM1 flatmap. This flatmap is used at the start of all simulations such that HOWFSC requires less stroke and attains better contrasts.

4.

Motivation for iEFC

A risk for the Roman Coronagraph will be how well can the model be calibrated to the instrument once the instrument is in orbit. As such, many efforts have been made to investigate the impact of model errors. Previous simulations and experiments by Sidick et al.¹⁴ and Marx et al.²⁴ have demonstrated that using a $\beta$ scheduling approach can mitigate the impact of model errors because periodically relaxing the regularization parameter $\beta$ can improve the correction of poorly controlled modes even though overall contrast can degrade. By making the $\beta$ value more strict again, a better overall contrast can be attained. Effectively, this means that model errors can be mitigated given a larger number of iterations to perform HOWFSC.

Figure 5 presents the results of a single case study investigating the performance of the SPC-WFOV mode in the presence of an unknown lateral shear of the shaped pupil mask (SPM). First, EFC is performed on the SPC-WFOV mode with no model errors. In this scenario, the end-to-end model is used to compute the Jacobian, and the same model is then used within the EFC controller. Note that PWP is ignored by using the model to directly compute the electric field from which the actuator commands are calculated, so perfect electric field estimation is simulated. With no model errors, EFC rapidly converges to $1.43\times {10}^{-10}$ after just 15 iterations.

Fig. 5

When no model errors are simulated, EFC converged to the result in panel (a) after just 15 iterations. When introducing the SPM shear of $\sim 85\mu \mathrm{m}$, the result in panel (b) was attained after 60 iterations. Panel (c) illustrates the mean contrast per iteration for both simulated scenarios. Panel (d) presents the $\beta$ regularization value used during each iteration, and the green lines represent each iteration at which the Jacobian was recomputed using the nonlinear coronagraph model.

For the scenario with an unknown lateral shear of the SPM, the model with no errors is again used to compute the initial Jacobian. However, a second model where the SPM has been shifted by 0.5% of the pupil diameter in the $y$-direction is used within the EFC controller to compute the electric field and apply the DM commands. In physical units, this shift is $\sim 85\mu \mathrm{m}$, which is $\sim 0.25$ actuators when projected onto the DM pupil. Because of this unaccounted-for error, $\beta$ scheduling and recomputation of the Jacobian become significant within the EFC controller to obtain adequate convergence. How often to switch regularization along with the values to switch between will depend on the coronagraph and the amount of model error. Here, using a strict $\beta$ value of $-2.5$ and a relaxed $\beta$ value of $-5$ while alternating every three iterations was found to be an effective schedule after trial and error. In addition, a new Jacobian is recomputed with the nonlinear coronagraph model every 15 iterations because as DM commands accumulate, the previous Jacobian becomes stale, and EFC would diverge regardless of the $\beta$ schedule. After every relinearization, the $\beta$ schedule is restarted at the strict value of $-2.5$. Relinearizing more frequently at every 12 iterations was also simulated but did not yield significant improvements in contrast or convergence.

While this use of $\beta$ scheduling and recomputation of the Jacobian can mitigate the effect of model errors for EFC, the controller requires more iterations to converge and becomes more complicated with $\beta$ values needing to be tuned according to the coronagraph and model. The multiple relinearizations also increase the computational requirements of the controller, which can become an issue for a HWO since computations will be done on-board instead of using a ground-in-the-loop scheme.²⁵ Additional methods such as the expectation-maximization technique detailed by Sun et al.²⁶ have been proposed to recover the state of a system and improve the model Jacobian, but this adds additional complexity to the controller that is not implemented for comparison here.

5.

Monochromatic iEFC Results

Using the SPC-WFOV mode at 825 nm, the iEFC Jacobian is generated for three types of modes: Fourier modes, individual actuators, and Hadamard modes. Figure 6 depicts examples of both Fourier and Hadamard modes. The Fourier modes are determined by discretizing the focal plane with a chosen spatial frequency sampling and generating cosine and sine commands corresponding to each of the desired spatial frequencies. The Hadamard modes on the other hand are generated based on the total number of actuators on a DM.

Fig. 6

Examples of Fourier and Hadamard modes. The Fourier mode is a particular cosine corresponding to a specific spatial frequency that will be used within the control loop. The Hadamard modes are binary commands that span the entire vector space of DM commands.

For each modal basis, the same probes are used to generate the difference images of each calibration mode. These are analogous to the same probes that could be used for PWP. Initially, single actuator pokes are used for the probes as these were also used by Haffert et al.⁹ However, switching to the Fourier probes depicted in Fig. 7 improved the simulations with noise in Sec. 7 because the flux in the control region is increased, allowing for higher signal-to-noise ratio (SNR) and reduced calibration times. For consistency, these monochromatic simulations also implement the same Fourier probes where each DM command is a weighted sum of cosine and sine modes that span the desired control region. This helps to generate a relatively uniform response in the dark hole. The total sum for each probe is then shifted to a region of the DM where the actuators are not obstructed by the SPM.

Fig. 7

Each probe command is a weighted sum of cosine and sine Fourier modes that have been shifted to a region of actuators not attenuated by the SPM. The first probe is solely a superposition of cosine Fourier modes (a), the second is an equal superposition of cosine and sine modes (b), while the final probe is only a superposition of sine modes (c). The difference images of each probe are measured to generate a response or compute modal coefficients.

Overall, four sets of calibration modes were tested; two sets of Fourier modes, one set of single actuator pokes, and one set of Hadamard modes. For each set, the calibration amplitude ($a_{\mathrm{mode}}$) was set to 5 nm, and the probe amplitude ($a_{\mathrm{probe}}$) was set to 20 nm. These values are chosen based on the typical range of values used for testbed experiments with iEFC. After each calibration, iEFC is performed for 30 iterations because that is the desired value for the Roman Coronagraph HOWFSC plan.¹¹

When it came to the Fourier modes, two sets were used because the initial set only considered spatial frequencies within the desired control region. This meant that the only controllable frequencies spanned from $5.4\lambda /D$ to $20.6\lambda /D$ with a sampling of $1\lambda /D$. The result after iEFC yielded a contrast of $4.47\times {10}^{-10}$. However, using an expanded set of Fourier modes that spanned from $1\lambda /D$ to $23\lambda /D$, the final contrast reached is $2.08\times {10}^{-10}$. Figure 8 depicts the two results using both the constrained and extended Fourier modes. In total, the constrained set contained 2552 modes, and the extended set contained 3480 modes.

Fig. 8

The iEFC solutions using Fourier modes constrained to the desired control region reach a contrast of $4.47\times {10}^{-10}$ while the solutions with the extended set of Fourier modes reach $2.08\times {10}^{-10}$. This demonstrates that control of higher and lower spatial frequencies is essential for the HOWFSC method to perform optimally.

For the single actuator modes, the total number of actuators is culled to only the ones within the extent of the geometric pupil, which is $\sim 46.3$ actuators across the diameter and totals 1804 actuators on each DM, or 3608 total modes. Using this basis results in a slight improvement over the extended Fourier modes as the final contrast reached is $1.85\times {10}^{-10}$. Finally, the Hadamard modes are selected by generating a set for each DM. As Hadamard modes are generated in powers of two and a single DM is assumed to use 1804 actuators, a single DM uses 2048 Hadamard modes. Calibrating both DMs meant using a total of 4096 modes, and the final mean contrast was $1.1\times {10}^{-10}$. Depicted in Fig. 9 are the final images for both the single actuator modes and the Hadamard modes.

Fig. 9

Using individual actuator modes, iEFC converged to a contrast of $1.85\times {10}^{-10}$. With Hadamard modes, the final contrast was slightly improved to a value of $1.11\times {10}^{-10}$.

For all simulations of iEFC, the control loop used the same $\beta$ schedule for one-to-one comparisons after tuning the values of the scheduler based on what yielded adequate convergence within 30 iterations. Figure 10 directly compares the performance of iEFC using all calibration modes considered along with the $\beta$ schedule used. Similar to EFC, a strict value of -2.5 was found to be sufficient to converge without degrading contrast while the relaxed value oscillates from −6 to −5. The final $\beta$ used is −4 because that was the most relaxed value that could be used on the final iterations while converging and not requiring any more iterations with a value of −2.5.

Fig. 10

Mean contrasts and the respective $\beta$ values per iteration for each modal basis tested. While the margins between the various modal bases are narrow, it remains notable that the best contrast achieved was $1.11\times {10}^{-10}$ with Hadamard modes.

Although all results with iEFC attained similar convergence and final contrast values, the best result was with Hadamard modes. Further tuning of the $\beta$ values for each modal basis could improve the results individually, but for the one-to-one comparison done here, Hadamard modes proved the most effective. As such, this basis is chosen for the next simulation which includes the same lateral SPM shift that was introduced in the perturbed EFC simulation. The significant difference with iEFC is that the error is present during the calibration given that all empirical data would inherently capture the perturbation. For this simulation, the same $\beta$ regularization schedule was used, and the final solutions depicted in Fig. 11 achieved a slightly worse final contrast with a value of $1.33\times {10}^{-10}$. While the ideal scenario for EFC achieved the same contrast as the best case of iEFC within 15 iterations, this simulation demonstrates that iEFC is a preferable alternative in the scenario where the instrument is not well calibrated because EFC required 60 iterations with the increased computational complexity of relinearizing the Jacobian regularly.

Fig. 11

After simulating iEFC with the same SPM shear introduced in the EFC simulations, a much smaller contrast degradation is noticed with the final contrast now being $1.33\times {10}^{-10}$.

6.

Broadband iEFC Examples

Given that the Roman Coronagraph will have both narrowband filters chosen for wavefront control purposes and broader bandpasses for science observations, additional simulations have been performed for both a wavefront control filter and a science bandpass. Here, filter 4b is chosen as the wavefront control filter with a central wavelength of 825 nm and a full-width at half-maximum of 3.6%.¹⁶ For the science filter, a bandpass of 10% is assumed.

To simulate a single broadband image, multiple wavefronts are propagated in the physical optics model with wavelengths spanning the chosen bandpass. The final image is an incoherent sum of the wavefronts at the image plane. For the narrowband simulations, three equally spaced wavelengths are selected to span the bandpass while five wavelengths are selected for the 10% bandpass. Each wavefront is weighted equally in the final incoherent sum, so no variations in the filter transmittance are being considered.

After calibrating iEFC for each bandpass, the control loop is repeated with results depicted in Fig. 12. Given that Hadamard modes had the best performance in the monochromatic simulations, that was the chosen modal basis here. After iEFC is performed, the 3.6% bandpass reaches a contrast of $1.51\times {10}^{-9}$, while the 10% bandpass reaches a contrast of $6.25\times {10}^{-9}$. Importantly, these are the results after just 15 iterations because convergence was significantly slower after reaching these contrast values.

Fig. 12

Final image using the single narrowband filter (a) illustrates that iEFC reached $1.51\times {10}^{-9}$ contrast, while the image using the 10% bandpass (b) has a contrast of $6.25\times {10}^{-9}$. As depicted in the plot of contrast per iteration (c), both converged in 15 iterations. More iterations were attempted while varying the regularization, but better results likely require additional calibrations of iEFC.

This demonstrates that even in the full 10% bandpass, iEFC can be capable of reaching below ${10}^{-8}$ contrast with the SPC-WFOV mode. Additional calibrations of iEFC after converging for a given bandpass could improve broadband iEFC results, but as explained in Sec. 7, iEFC will likely require lengthy calibration times making recalibration impractical.

7.

iEFC Simulations with Noise

Due to the empirical calibration of iEFC, noise will inherently be included in the Jacobian. To consider the impact of noise, these final simulations consider the same narrowband and broadband cases from Sec. 6, but flux estimates for a given reference star along with a model of Roman’s EMCCD are used to generate more realistic images with photon and detector noise.

Here, $\zeta$ Puppis is chosen as the reference star being utilized for wavefront control as it was in the Roman observing scenario simulations.²⁷ As described in Appendix B, the flux for $\zeta$ Puppis is computed at each wavelength that will be propagated for the chosen bandpass. The amplitude of each wavefront is then weighted by the square root of the flux estimate for the respective wavelength and propagated to compute an image in units of photons/second/pixel. An additional factor of 0.5 is applied to the flux at the image plane to approximate the throughput losses from reflective surfaces throughout the optical train (note that this is based on a reflectivity of $\sim 97.5\%$ for 25 reflective optics). Next, the image plane flux is input into the emccd_detect ( https://github.com/roman-corgi/emccd_detect) software to simulate an image with electron multiplying (EM) gain and detector noise. The additional parameters chosen to model the EMCCD are provided in Table 1.

Table 1

Parameters used to model the EMCCD. The bias is chosen based on the current best estimate of read noise for the CGI EMCCD.

Again, only Hadamard modes are considered to provide a direct comparison with the noiseless simulations. However, additional methods were implemented during the calibration to prevent both low SNR responses and responses with saturated pixels. The first of these methods is to implement varying calibration amplitudes for modes with a response that is distributed across more pixels in the image. Specifically, the response of each Hadamard mode is first estimated with the magnitude of the mode’s Fourier transform. The maximum value of each estimated response is recorded such that the calibration amplitude can be scaled according to the ratio of the maximum for a given mode and the largest maximum recorded. Therefore, the mode with the largest maximum value in its response uses a calibration amplitude of 5 nm, but all others are multiplied by a scaling factor greater than 1.

The second method used to calibrate both bandpasses is to average and stack frames with varying exposure times. This effectively increases the dynamic range of the images because pixels with large flux values do not get saturated at the lowest exposure time, but pixels with low flux values acquire a higher SNR when using the longer exposure times. For both the 3.6% and 10% bandpasses simulated, exposure times of 0.01, 0.05, and 0.1 s are stacked to generate the final image. The number of frames for each respective exposure time are 10, 4, and 2. The total integration time for each calibration frame is then 0.5 s, so to calibrate all 4096 Hadamard modes, approximately 6.8 h of data acquisition is required. Note that to perform full broadband control using all three narrowband wavefront control filters will require $3\times$ the integration time, a total of 20.4 h, but this is not simulated here. The EM gain used for all frames in the 3.6% bandpass is 250, and the EM gain for frames in the 10% bandpass is 100.

Figure 13 depicts a useful visualization for the quality of the calibration. The images in Fig. 13 illustrates the normalized root mean square (RMS) response of each actuator on DM1. In the case without any noise, the obstructed actuators have a considerably attenuated response, while the two scenarios with noise have a relatively constant bias that is a direct result from noise during calibration. When using the 10% bandpass, the bias of the obstructed actuators is slightly less due to the higher SNR achieved by capturing more photons.

Fig. 13

After calibrating iEFC while including noise, the RMS response of DM1 actuators provides a useful visualization of the quality of the calibration. Panel (a) is the RMS response with a noiseless calibration, while panels (b) and (c) show the RMS responses after calibrating the 3.6% and 10% bandpasses with noise. The constant bias of the obstructed actuators is a direct result of the noise from calibration. Here, the simulation of the 10% bandpass has a slightly lower bias than the 3.6% bandpass because the additional photons improved the SNR.

Because of the noise accumulated during calibration, we see additional degradation in the performance of iEFC, as presented in Fig. 14. With the narrowband filter, the contrast converged to $9.68\times {10}^{-9}$, while the 10% bandpass reached $1.76\times {10}^{-8}$. Each simulation only used 30 iterations because the convergence of iEFC became stagnant even when the regularization value was varied.

Fig. 14

With noise implemented, iEFC reached $9.68\times {10}^{-9}$ using the 3.6% bandpass (a), while $1.76\times {10}^{-8}$ was reached with the 10% bandpass (b). As presented in (c), 30 iterations were used for each simulation because the noise from the calibration degraded the convergence for each bandpass.

These simulations demonstrate the disadvantage of iEFC which is the noise accumulated during calibration. Possible mitigation strategies could be a smarter choice of EMCCD parameters, including the use of higher gain values and photon counting to reduce detector noise, using a brighter reference star (but this will depend on the desired target star and stability of the optical assembly), using longer integration times for calibration, or choosing a more optimal set of modes for the controller. In any case, additional techniques should be investigated to mitigate the impact of noise in the calibration of iEFC if it is to be implemented by Roman.

8.

Conclusion

The iEFC algorithm is an exciting method for the high-contrast imaging community. Like EFC, iEFC is naturally extended to coronagraphs with two DMs and is capable of digging annular dark holes. This makes iEFC a potential method for the Roman Coronagraph and future coronagraphs for the HWO. Here, the efficacy and potential benefits of iEFC have been investigated specifically for the SPC-WFOV mode due to the large linear regime of the DMs. The simulations demonstrate that contrasts comparable to those of EFC can be achieved in the ideal monochromatic scenario and that Hadamard modes are the most effective for the control loop. More significantly, iEFC would reduce the risk of unknown system perturbations that would otherwise require a more complicated $\beta$ scheduling method, more iterations, and greater computational complexity due to relinearizations of the Jacobian. By inherently including the state of the instrument during calibration, iEFC does not need to be relinearized for coronagraphs with large linear regimes such as the SPC-WFOV mode or a vortex coronagraph as simulations here demonstrate that adequate convergence can be obtained within 30 iterations.

Additional simulations demonstrate that iEFC can be calibrated and performed with both narrowband and broadband filters, but the performance of iEFC degrades due to the incoherent superposition of wavelengths. Further simulations including the EMCCD model predict that contrasts of ${10}^{-8}$ can be achieved within 30 iterations after calibrating for an estimated time of 6.8 h when using $\zeta$ Puppis as the reference star. Further investigations including more optimal EMCCD parameters, smarter choice of calibration modes, and improved calibration techniques could reduce the calibration time and improve the final contrast, but other factors such as system drift during calibration should also be considered in the future. Given the simulation results for the time being, the SPC-WFOV mode is a suitable instrument to perform iEFC tests during the latter stages of the Roman mission to inform the WFSC strategies for future coronagraphs such as those considered for the HWO.

9.

Appendix A: Definitions

The variables used for the formalism of EFC and iEFC are summarized here. The shapes of the matrices and vectors are also included to clarify how the algorithms are implemented (Table 2).

Table 2

Significant variables used for the descriptions of EFC and iEFC.

10.

Appendix B: Blackbody Flux Calculations

To include estimates of flux, the blackbody equation

Here, $R$ is the radius of the star, and $d$ is the distance to the star, both of which are included in Table 3. The spectral photon flux is calculated by converting the spectral irradiance using $E=hc/\lambda$. Finally, the spectral flux is integrated over several sub-bandpasses to estimate the flux for each wavelength that is propagated. Figure 15 presents an example of the spectral photon flux computed for $\zeta$ Puppis and seven sub-bandpasses that would be integrated to estimate the photon flux for those seven wavelengths.

Table 3

Parameters used to calculate the spectral radiance of ζ Puppis along with the final solid angle used to calculate flux.

Fig. 15

Example of the spectral flux computed for $\zeta$ Puppis using the blackbody equation. The highlighted regions illustrate the sub-bandpasses integrated to estimate flux values for individual wavelengths centered at each sub-bandpass. Note that the peak of the blackbody is located at $\sim 200\mathrm{nm}$ due to the temperature of the star.