High-order coronagraphic wavefront control with algorithmic differentiation: first experimental demonstration

Abstract. Future space-based coronagraphs will rely critically on focal-plane wavefront sensing and control with deformable mirrors (DMs) to reach deep contrast by mitigating optical aberrations in the primary beam path. Until now, most focal-plane wavefront control algorithms have been formulated in terms of Jacobian matrices, which encode the predicted effect of each DM actuator on the focal-plane electric field. A disadvantage of these methods is that Jacobian matrices can be cumbersome to compute and manipulate, particularly when the number of DM actuators is large. Recently, we proposed a new class of focal-plane wavefront control algorithms that utilize gradient-based optimization with algorithmic differentiation to compute wavefront control solutions while avoiding the explicit computation and manipulation of Jacobian matrices entirely. In simulations using a coronagraph design for the proposed Large UV/Optical/Infrared Surveyor, we showed that our approach reduces overall CPU time and memory consumption compared to a Jacobian-based algorithm. Here, we expand on these results by implementing the proposed algorithm on the High-contrast Imager for Complex Aperture Telescopes tested at the Space Telescope Science Institute and present initial experimental results, demonstrating contrast suppression capabilities equivalent to Jacobian-based methods.


Introduction
Future space coronagraphs attempting to image and characterize Earth-like planets around nearby solar-type stars will rely critically on closed-loop wavefront sensing and control (WFS&C) using deformable mirrors (DMs) to mitigate optical aberrations in the primary beam path.These aberrations, primarily mid-spatial frequency wavefront errors and optical misalignments in the telescope and coronagraph optics, give rise to a speckle floor that is coherent with the star and evolves slowly over time in response to miniscule drifts in the thermal and mechanical state of the observatory.If uncorrected, the speckle floor overwhelms the faint image of the orbiting planet, which is expected to be 10 10 times fainter than the host star at 0.1 arcseconds of separation or less at visible wavelengths. 1 The current state-of-the-art wavefront control algorithms, stroke minimization (SM) 2 and electric field conjugation (EFC), 3 compute DM command solutions using a first-order approximation of the focal-plane electric field in an optimal control framework.In both cases, the optimal DM update is written down in closed form as the solution of a linear system of equations constructed from a Jacobian matrix that describes the influence of each DM actuator on the focalplane electric field.
Until recently, little attention has been paid to the computational demands of SM and EFC.The complexity of SM and EFC is dominated by the cost of computing and manipulating the Jacobian matrix, whose size is proportional to the product of the number of DM actuators and number of dark-zone pixels.In broadband imaging, this is compounded by the requirement to compute a separate Jacobian for each controlled wavelength.The Jacobian is most often modelbased, in which case an optical diffraction model of the coronagraph is evaluated repeatedly to predict the focal-plane influence of each actuator.Moreover, the Jacobian is a linearization of the true, nonlinear behavior of the DMs and must be recalculated periodically as the state of the DMs evolves over time.
As direct imaging missions demand DMs with ever-higher actuator density to enable wider and wider search areas, computational aspects will inevitably become a point of concern from a systems engineering standpoint.In on-orbit WFS&C, all sensing and control computations are processed by the flight computer; conversely, in ground-in-the-loop scenarios, raw data is communicated to a ground-based computing node that calculates the DM correction and relays the correction back to the observatory.Though each approach has tradeoffs, a major advantage of on-orbit WFS&C is the ability to update DM commands more frequently without relying on the continuous availability of a communication link with the ground station.On-orbit WFS&C can help to relax observatory-level wavefront stability requirements by enabling the high-contrast dark zone to be maintained over shorter time intervals.Successful deployment of an on-orbit architecture is predicated on the availability of sufficient computational resources; however, radiation-hard, space-qualified computing hardware lags behind conventional hardware by decades and is extremely resource-limited.This poses a substantial computation capability gap if the current algorithms are expected to be deployed on-orbit on a future direct imaging mission, such as NASA's planned Habitable Worlds Observatory (HWO).
As a case in point, the current state-of-the-art testbed for coronagraph laboratory demonstrations, the decadal survey testbed (DST) at NASA's Jet Propulsion Laboratory, has successfully demonstrated 3.82 × 10 −10 instrumental contrast over a 10% fractional bandpass within an annular dark zone extending from 3 λ 0 ∕D to 8 λ 0 ∕D, where λ 0 is the central imaging wavelength and D is the diameter of the telescope aperture, using two DMs each with 48 × 48 actuators. 4he baseline requirements for an HWO-like mission will be considerably steeper, using the Large UV/Optical/Infrared Surveyor (LUVOIR) and Habitable Exoplanet Observatory (HabEx) flagship mission concepts formulated for the Astro2020 Decadal Survey 5,6 as representative reference designs.The HabEx design features two 64 × 64 DMs, and a search area with a maximal outer radius of 32 λ 0 ∕D over a 20% bandpass-nearly double the total actuator count, a factor of 4 increase in search radius, and a factor of two increase in control bandwidth.Meanwhile, the LUVOIR Architecture "A" reference design includes a pair of 128 × 128 DMs and a dark zone with a 64 λ 0 ∕D maximal outer radius over a 10% fractional bandpass or wider.Because the number of detector pixels in the dark zone scales with the area of the dark zone rather than its radius, these parameters correspond to an increase by a factor of 32 (Habex) and 512 (LUVOIR-A) of the worst-case Jacobian dimensionality, for each control wavelength, compared to DST.
In recent work, we formulated an alternative wavefront control framework that iteratively compute DM updates using gradient-based optimization techniques, which eliminates the need to calculate and manipulate the Jacobian. 7Our approach is based in part on a numerical technique known as algorithmic differentiation (AD), 8 to calculate the gradients of the cost function for optimal control accurately and efficiently.We described an AD-based counterpart to SM,

Notation
In this paper, our principal concern is with algorithms that operate on discrete vector-valued quantities, which are represented in boldface.Many of these quantities vary as a function of control iteration, and are denoted with the subscript k.These may be truly discrete, such as the vector a k of DM actuator command updates, or may represent discretizations of functions of continuous spatial variables, such as electric fields.We denote x as a column vector, its transpose x T as a row vector, and kxk as its Euclidean length.For complex-valued quantities, † denotes the Hermitian transpose.Scalar quantities are denoted in ordinary (i.e., non-boldface) typographic weight.
By convention, the derivative of a scalar with respect to a column vector is a row vector, i.e.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; 3 7 6 where x½n is the n'th element of x.Consequently, the derivative of a column vector with respect to another column vector is a matrix of row vectors: ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 7 ; 3 2 Wavefront Control Using Algorithmic Differentiation The goal of the WFS&C loop in coronagraphy is to drive starlight within the dark zone toward zero over time so that a faint orbiting companion becomes detectable against the reduced starlight background.Each iteration of the WFS&C loop, indexed by the integer k, consists of two steps: (1) an estimation step, in which an estimate Êab;k of the true aberrated electric field E ab;k is formed within the dark zone from focal-plane intensity measurements, and (2) a control step, in which the DM correction is updated to reduce the energy in E ab;k , as shown in Fig. 1.In this paper, we focus principally on the control step.Modern model-based wavefront control algorithms find the DM correction update a k by minimizing some cost function J k ða k ; Êab;k Þ with respect to a k .Usually, J k is constructed to trade off between minimizing starlight and minimizing the size of the correction, which helps to regularize the problem and stabilize the solution.In general, the true relationship between the starlight in the dark zone and the DM correction is highly nonlinear and nonconvex, owing to the fact that the DMs impart phase-only corrections of the form expfiϕ DM g at or near the coronagraph entrance pupil.However, when the optical aberrations are small, we can approximate the true electric field in the coronagraph entrance pupil with a first-order Taylor series expansion about a small update to the DM commands.In this case, the corrected electric field in the dark zone has the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 5 0 1 where Êab;k is the estimate of the aberrated dark-zone electric field produced by the estimation step, and E DM;k is the change in electric field resulting from the update to the DM correction.We can also write E DZ;k in the form ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 4 ; 4 3 5 where G k ≜ ∂E DZ;k ∕∂a k is the Jacobian matrix with dimensions N pix × N act , N pix is the number of dark-zone pixels, and N act is the total number of controllable DM actuators.The intensity from the corrected electric field, integrated over the dark zone, can be written in terms of the Jacobian as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 4 ; 3 6 0 This is a quadratic function of a k , meaning that under this approximation, there exists a unique, optimal DM correction that minimizes the dark-zone starlight.
SM and EFC utilize the relationship in Eq. ( 5) to derive closed-form expressions for this optimal correction in terms of G k that can be evaluated by solving a linear system of dimension N act × N act , as illustrated in Fig. 2. As we discuss in Appendix B, this is equivalent to minimizing the cost function using Newton's method with an exact Hessian matrix.Alternatively, it is possible to find the DM correction by minimizing the cost function J k with respect to a k iteratively, rather than analytically, using gradient-based nonlinear optimization as shown in Fig. 3.
Fig. 2 In EFC and SM, the Jacobian matrix is precomputed using a computer model of the coronagraph to predict the effect of an update to each of the N act DM actuators individually on the darkzone electric field, here represented by the Kronecker δ functions δ n , where δ n ½i ¼ 1 if i ¼ n and zero otherwise.The Jacobian G k and the estimate of the aberrated dark-zone electric field Êab;k together are used to construct a linear system of equations whose solution is the desired DM update a Ã k .We show in Appendix B that this is equivalent to minimizing the wavefront control cost function using Newton's method.
Fig. 1 Closed-loop coronagraphic WFS&C uses images from the science camera, rather than an external wavefront sensing instrument, to estimate the electric field from the host star within the dark zone and drive it toward zero.Model-based WFS&C algorithms use a numerical model of the coronagraph to solve an inverse problem for the unknown electric field and corresponding DM correction, respectively.
To do so eliminates the need to calculate the Jacobian matrix, but requires a way of calculating the gradient vector ∂J k ∕∂a k .Reverse-mode algorithmic differentiation (RMAD) provides a way of doing so that is both computationally efficient and accurate, in the sense that the derivatives computed by RMAD are accurate to machine precision and do not utilize finite-difference approximations. 8The basic principle of RMAD is that any function that can be written down as a sequence of differentiable operations, called the forward model, can be transformed mechanistically to construct a related function, called the adjoint model, that evaluates the derivative of the forward model with respect to any of the intermediate variables encountered during its evaluation, as well as its inputs.Figure 4 illustrates this procedure for the wavefront control cost function J k .We refer the reader to our earlier work for a more comprehensive discussion of the principles of RMAD and its application to wavefront control. 71 Stroke Minimization: From Lagrange Multipliers to Penalty Method SM finds the smallest DM correction that achieves a desired level of stellar intensity integrated over the dark zone, denoted by I target;k . 2 It is solved by finding the stationary point of the Lagrangian function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 8 5 Fig. 3 Our AD-based wavefront controllers AD-PSM and AD-EFC use RMAD to differentiate the wavefront control cost function J k with respect to the DM correction update a k , yielding the gradient vector ∂J k ∕∂a k evaluated at the current iterate a n k .A nonlinear optimization algorithm calculates a new iterate a nþ1 k that reduces the value of the cost function, i.e., J k ða nþ1 k Þ ≤ J k ða n k Þ.This procedure is repeated until the gradient becomes sufficiently small, indicating that the solution a Ã k is at or near a local minimum of the cost function.A starting guess for the solution a 0 k as well as the aberrated electric field Êab;k are the input parameters.
Fig. 4 The forward model for the wavefront control problem maps DM command updates a k to a scalar cost function J k .The DM command vector is split into independent command vectors a 1;k and a 2;k for the pupil-plane and out-of-pupil DM, respectively.These are mapped onto DM surfaces s 1;k and s 2;k using the model in Appendix C, and propagated through an end-to-end coronagraph model to predict the resulting change in dark-zone electric field E DM;k .Reverse-mode AD transforms the forward model into an adjoint model that backpropagates the partial derivatives of J k with respect to each intermediate variable a 1;k , a 2;k , s 1;k , s 2;k , and E DM;k in reverse order, starting from the output on the right.The derivatives with respect to the individual DM command vectors a 1;k and a 2;k are concatenated to form the full gradient vector for optimization.
i.e., a point such that ∂L SM;k ∕∂a 0 k ¼ 0, where a k is the DM actuator command update, E DZ;k is the corrected electric field in the dark zone, a 0 k ≜ ½a T k μ k T , and μ k is the Lagrange multiplier.Because this stationary point is a saddle point, it cannot be reached by minimizing L SM;k with respect to a 0 k .Instead, one chooses a fixed starting value for the Lagrange multiplier, μ 0 k , and minimizes L SM;k with respect to a k to find a corresponding DM solution a 0 k .If the constraint kE DZ;k ða 0 k Þk 2 ≤ I target;k is not satisfied, a larger value μ 1 k > μ 0 k is selected and this procedure is repeated.For any fixed value of μ n k , the minimum of Eq. ( 6) is given in closed form in terms of the Jacobian matrix G k : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 4 ; 6 3 4 where I is the identity matrix and Êab;k is the estimate of the aberrated electric field.AD-PSM has the same goal, but instead iteratively minimizes the cost function ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 4 ; 5 7 0 where η k is a parameter that encodes the relative importance of minimizing actuator stroke and driving the integrated intensity toward the target value.The minimum of J PSM;k with respect to a k is coincident with the stationary point of L SM;k , 9 meaning that the DM solutions obtained by SM and AD-PSM are, in principle, identical.In this paper, we choose η k in each WFS&C iteration as ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 4 9 6 where η 00 , known as the penalty parameter, is a constant set by the experimenter.The denominator ðk Êab;k k 2 − I target;k Þ 2 scales the cost function to be invariant to the energy in the dark-zone electric field, which helps in practice to obtain good solutions in all WFS&C iterations as Êab;k gradually converges toward zero.

Electric Field Conjugation
EFC attempts to drive the dark zone electric field toward zero, with Tikhonov regularization to mitigate ill-conditioning caused by the presence of actuators that are completely or partially obscured by pupil features.Its cost function for a single correction wavelength is given as ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 3 4 6 where Γ k is the Tikhonov regularization matrix.In the most common case, one chooses Γ k ¼ α k I, making EFC identical to SM with a fixed Lagrange multiplier μ k ¼ 1∕α 2 k .For this case, its solution can be obtained using Eq. ( 7).The general solution for the Jacobian-based formulation of EFC is given as In AD-EFC, one iteratively minimizes a scaled version of the EFC cost function ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 4 ; 2 3 4 Similarly to AD-PSM, the scaling factor 1∕k Êab;k k 2 makes the AD-EFC cost function invariant to the energy in the dark-zone electric field to aid in obtaining numerical solutions in practice.The RMAD adjoint model for J EFC;k is provided in Appendix D.

Experimental Setup
In this section, we provide an overview of the HiCAT testbed and provide details about the implementation of AD-PSM and AD-EFC, the reference Jacobian-based implementations of SM and EFC, and the electric field estimation algorithm used in the estimation step of the WFS&C loop.

HiCAT Testbed
1][12][13][14][15][16][17] These technologies include Lyot coronagraphy, high-order WFS&C for generating and stabilizing dark zones, and low-order wavefront sensing (LOWFS). 18HiCAT operates in a mid-contrast regime (10 −7 to 10 −8 ) which approaches the limit achievable outside of a vacuum environment.The testbed is equipped with two Boston Micromachines Kilo-DMs for high-order sensing and control, with 952 actuators each, making it suitable for our proof-of-concept demonstrations.One DM is placed in a plane conjugate to the HiCAT entrance pupil, while the second DM is placed ∼300 mm farther along the optical axis, corresponding to a Fresnel number N F ≈ 98 at a wavelength of 638 nm.This configuration enables simultaneous control of amplitude and phase aberrations over a dark zone that extends over both halves of the image plane.We conducted our experiments using a Thorlabs MCLS1 laser diode source, which emits monochromatic light with a central wavelength λ 0 ¼ 638 nm.
HiCAT additionally has an IrisAO segmented DM with 37 hexagonal segments, each with controllable piston/tip/tilt, to act as a telescope pupil simulator and to inform experimental efforts devoted to segment-level tolerancing and stabilization. 19,20Figure 5 shows a simplified system layout of HiCAT, including the primary imaging beam path as well as several additional beam paths used by the LOWFS and metrology subsystems.
Our experiments on HiCAT utilized a classical Lyot coronagraph (CLC) design with a hexagonal entrance pupil mask designed to mask the extreme edges of the IrisAO, a circular focal-plane mask, and a circular Lyot stop.Figure 6 shows an overlay of the CLC pupil masks along with a simulated stellar point-spread functions (PSFs).Figure 7 shows example experimental PSFs obtained before and after closed-loop WFS&C using SM, along with the corresponding DM commands.

Algorithm Implementation
We developed a differentiable model of HiCAT using Python, including a hand-derived adjoint model.To facilitate testing, our model was comprised of several sub-modules each with its own individual forward and adjoint models.
1.A model to compute the DM surface resulting from a given set of actuator commands using a fast convolutional representation.This is described in further detail in Appendix C. 2. A model to apply the phase corrections imparted by the in-pupil DM (DM1) and out-ofpupil DM (DM2), including the free-space propagation between the two DMs. 3. A model to propagate the electric field after DM correction through the HiCAT CLC.
We use the semi-analytical Lyot coronagraph model originally described in Ref. 21.
We refer the reader to our earlier work for detailed descriptions of the operations in the forward and adjoint models for the latter two sub-modules. 7A pre-existing numerical model of HiCAT based on the POPPY framework 22,23 served as a reference for calibrating our differentiable model.

Optimization algorithm
We used the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm, 9 as implemented in the SciPy package, 24 as the optimization algorithm for AD-PSM and AD-EFC.L-BFGS is a quasi-Newton algorithm, meaning that it uses the gradient vectors collected during optimization to approximate the inverse Hessian matrix H −1 k of second derivatives, involving substantially less effort than Newton's method with an exact Hessian (see Appendix B).
All numerical optimization algorithms have a termination criterion that determines when the algorithm has reached a location in the parameter space that is sufficiently close to a local optimum.In the L-BFGS implementation used by SciPy, the termination criterion is set by a tolerance parameter ε defined as the magnitude of the largest element of the gradient vector.As ε → 0, L-BFGS carries out a greater number of optimization iterations to terminate closer to the minimum.Choosing a larger value of ε reduces total computation with the tradeoff of a less-optimal DM solution.However, as our experiments showed, in certain instances this can help to regularize the solution by terminating before L-BFGS converges to an overly aggressive DM correction that would otherwise make the WFS&C loop unstable.
The termination criterion presents an additional consideration not present in SM and EFC which exactly minimize their respective cost functions in each WFS&C iteration.The value of the tolerance parameter ε should be chosen to optimally trade-off between accuracy and computational effort.We discuss this in further detail below.

Estimation Algorithm
We used the pairwise probe estimator 3 for the estimation step in all experiments.The pairwise estimator forms a least-squares estimate Êab;k of the focal-plane electric field E ab;k by applying a series of P probing DM commands u p to generate probing electric fields E DM;k ðu p Þ that interfere with E ab;k .The data vector for the least-squares estimate is formed by differencing the images resulting from E DM;k ðu p Þ and E DM;k ð−u p Þ.The estimate of the m'th pixel in the dark zone is then found by finding the least-squares solution of the system E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 7 ; 4 5 6 2 6 6 6 4 where We generated four DM probe functions u p that were optimized to produce probing fields of the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 7 ; 3 5 0 where p ∈ f1;2; 3;4g and sgnfθ x g is the sign of the focal-plane x coordinate.These probing fields, which have anti-Hermitian symmetry, correspond to purely real commands applied to the in-pupil DM, which we can calculate via a regularized least-squares approach for each p, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 7 ; 2 7 7 arg min where α 2 probe is the Tikhonov regularization term for the probe calculation and G 1;k is the Jacobian for the in-pupil DM.This is identical to the EFC problem in Eq. (10) but with the probing field E p DM;k on the right-hand side instead of − Êab;k , and has the solution ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 7 ; 2 0 0 For the experiments with SM and EFC, we applied Eq. ( 16) to compute probe commands for pairwise probe estimation.For experiments with AD-PSM and AD-EFC, we iteratively solved Eq. ( 15) using the AD-EFC framework.In all cases, we set the Tikhonov regularization parameter for probe generation as α 2 probe ¼ 0.7; for the iterative probes, we set the optimizer tolerance to ε ¼ 10 −5 .Figure 8 shows the probe commands u p along with the corresponding magnitude and phase of E DM;k ðu p Þ.
We conducted a series of experiments to compare the contrast performance of AD-PSM and AD-EFC relative to SM and EFC, respectively.All experiments used an annular control region extending from 5.8 λ 0 ∕D LS to 9.8 λ 0 ∕D LS , where λ 0 ¼ 638 nm is the central wavelength of the Thorlabs MCLS1 laser diode and D LS is the Lyot stop diameter.Each experiment consisted of 80 WFS&C iterations.For each experiment, we computed the median and 10th and 90th percentile values of the spatially averaged dark-zone contrast values for the final 50 iterations, which captured the steady-state performance of the WFS&C loop with each algorithm after converging to deepest-possible contrast.Figures 9 and 10 show example contrast vs. iteration time series with AD-PSM/SM and AD-EFC/EFC, respectively, as well as the PSF from the iteration with deepest contrast using AD-PSM and AD-EFC.In both cases, the convergence properties of AD-PSM and AD-EFC are nearly identical to their Jacobian-based counterparts, validating our proposed approach.
For SM and AD-PSM, we chose I target;k ¼ 0.5k Êab;k k 2 as the contrast target, i.e., a factor of two improvement iteration-over-iteration in spatially integrated dark zone contrast.For the Lagrange multiplier line search in SM, we let μ nþ1 k ¼ 1.3μ n k .For AD-PSM, we tested combinations of the penalty parameter η 00 ∈ f10;100;1000g and the nonlinear optimization convergence tolerance ε ∈ f10 −2 ; 10 −3 ; 10 −4 g.We determined these values based on simulations and prior WFS&C experiments on HiCAT.
For EFC and AD-EFC, we selected the Tikhonov regularization matrix as Γ k ¼ α k I, with α 2 k ∈ f10 −2 ; 10 −3 ; 10 −4 g for the first 30 WFS&C iterations and α 2 k ¼ 10 −1 thereafter.We determined through previous experiments with EFC that maintaining an aggressive α k value throughout the experiment was detrimental to the stability of the control loop after reaching the steady-state regime.For AD-EFC, we used the same set of ε values as AD-PSM.
Fig. 8 The probe commands u p for pairwise estimation were optimized so that the resultant dark-zone electric field E DM;k ðu p Þ ¼ sgnfρ x g expfiπðp − 1Þsgnfρ x g∕4g, as described in Sec.3.3.The inner and outer edges of the dark zone are also shown for reference.The probe commands are close to the inverse Fourier transform of the dark zone geometry (an annulus), modulated by a horizontal sinusoid whose phase angle is proportional to the desired piston phase, and projected onto the DM actuator coordinates.Only the pupil-plane DM is modulated.
For AD-EFC and AD-PSM, we compared the value of the cost function for two different starting guesses for the DM correction: a 0 k ¼ a Ã k−1 , i.e., the solution of the previous WFS&C iteration, and a 0 k ¼ 0. The starting guess with the lower of the two cost function values was then selected.In all cases a 0 k ¼ 0 was ultimately chosen as the starting guess.Figures 11 and 12 show the statistics of the post-convergence spatially averaged dark-zone contrast achieved with AD-PSM and AD-EFC for each combination of regularization (η 00 or α k ) and optimization tolerance ε, respectively, compared to reference experiments with SM and EFC.For all parameter combinations, AD-PSM and AD-EFC equaled the contrast performance of SM and EFC, respectively.Moreover, we observed no strongly identifiable trends in achievable contrast as a function of the algorithm parameters, indicating that in all cases, AD-PSM and AD-EFC reached the contrast floor imposed by environmental instabilities.For full contrast versus iteration curves for all experiments, we refer to Fig. 13 in Appendix A.

Discussion
Our experiments were aimed at exploring a relevant subset of the space of free parameters for each algorithm, namely the nonlinear optimization convergence tolerance ε, the Tikhonov regularization α k for AD-EFC, and the penalty parameter η 00 for AD-PSM.In principle, each parameter affects the attainable contrast of the WFS&C loop, but in subtly different ways, which we discuss here.
As discussed in Sec.3.2.1,ε determines the effort that the nonlinear optimization algorithm will expend to find a solution close to the true minimum of the cost function.A smaller value of ε corresponds to a larger number of optimization iterations before termination, and ultimately a larger time interval between DM updates.In principle, on a system, such as HiCAT, where environmental disturbances cause the aberrated electric field to evolve over time scales on the order of seconds or faster, an excessive delay between the estimation step and application of the DM correction can cause a degradation in achievable contrast.However, we observed no meaningful degradation for smaller values of ε.In a real space-borne system, this is unlikely to be a significant consideration because of the much greater electric field stability, and because the total duration of each WFS&C iteration will be dominated by the exposure times needed for the estimation step.
The value ε can also serve as an auxiliary form of regularization, by terminating the optimization algorithm before it reaches an overly aggressive DM correction caused by a noisy electric field estimate, insufficient regularization using α k or η 00 , or both.For instance, in Fig. 12, with α 2 k ¼ 10 −4 (rightmost panel), EFC diverged altogether whereas AD-EFC did not.On the other hand, choosing ε too large can impose an effective contrast floor by limiting the ability of AD-EFC vs. EFC Fig. 12 Median, 10th percentile, and 90th percentile spatially averaged contrast values achieved by AD-EFC (dark orange) as a function of optimizer tolerance, for three different values of the Tikhonov regularization parameter α k .For each α k , we also ran a reference experiment with EFC (light orange); the EFC experiment with α 2 k ¼ 10 −4 diverged and is not shown.The contrast versus iteration time series for the rightmost result in the left pane (α 2 k ¼ 10 −2 , ε ¼ 10 −4 ) is illustrated in Fig. 10.The performance of AD-EFC was statistically equivalent to EFC for α 2 k ¼ 10 −2 and α 2 k ¼ 10 −3 .As discussed in Sec.4.1, the optimizer tolerance prevented AD-EFC from diverging in this case.

AD-PSM vs. SM
Fig. 11 Median, 10th percentile, and 90th percentile spatially averaged contrast values achieved by AD-PSM (dark blue) as a function of optimizer tolerance, for three different values of the penalty parameter η 00 , compared to SM (light blue).The contrast vs. iteration time series for the rightmost result in the left pane (η 00 ¼ 10, ε ¼ 10 −4 ) is illustrated in Fig. 9.In all cases, the contrast performance of AD-PSM was equivalent to that of SM within statistical uncertainty.
the optimization algorithm to converge to appropriately strong corrections.We determined in simulation that this was the case for ε > 10 −2 .
For EFC and AD-EFC, smaller values of the Tikhonov regularization α k correspond to more aggressive correction of the electric field, with the tradeoff of increased sensitivity to small perturbations in the estimated electric field, which reduces the stability of the WFS&C loop.For AD-PSM and SM, the aggressiveness of the WFS&C control loop is set first and foremost by the target contrast level I target;k .For SM, a smaller I target;k corresponds to a larger value of the Lagrange multiplier μ k , and therefore more aggressive correction (recalling from Sec. 2.2 that μ k ¼ 1∕α 2 k ).For AD-PSM, as η 00 tends toward infinity, the goal of achieving the target contrast is enforced more strongly; η 00 is interpretable as tuning the aggressiveness of the control loop up to a maximum level imposed by I target;k .Our experiments indicated that over the range of values considered, the performance of the WFS&C loop was insensitive to the value of η 00 .
In principle, because α k and ε have similar effects on the performance of AD-EFC, there potentially exists a single combination of the two parameters that is optimal in terms of contrast, or perhaps a continuum of combinations with similar performance, with a change in α k compensated by a change to ε in the opposite direction.The same holds true for AD-PSM.

Conclusions and Future Work
In this paper, we reported the first experimental demonstrations of two AD-based wavefront control algorithms, AD-PSM and AD-EFC, using the HiCAT testbed.To within statistical uncertainty, AD-PSM and AD-EFC equaled their Jacobian-based counterparts in dark-zone contrast for all combinations of parameters that we tested.These demonstrations pave the way for future experimental validation at higher contrast.
The analysis in our earlier work indicated that the largest computational gains are realized for DMs with more than 64 × 64 actuators.The DMs currently in use on HiCAT have 34 actuators across the diameter of the active region or 952 actuators per DM in total, which is comparatively low.Therefore, our goal was to validate the fundamental capability of AD-PSM and AD-EFC to reach deep contrast, rather than to demonstrate improved computational efficiency.
In this work and in our earlier work, we utilized the L-BFGS optimization algorithm to minimize the wavefront control cost function because of its desirable convergence properties compared to first-order optimization methods and low storage requirements.Despite this, L-BFGS is known to converge slowly for poorly conditioned problems compared to methods that utilize the exact Hessian matrix. 9However, there exist alternative methods, such as truncated Newton algorithms, with excellent convergence properties for quadratic or nearly-quadratic cost functions, that require only the ability to evaluate Hessian-vector products, rather than the Hessian matrix itself. 9Hessian-vector products can be evaluated using AD in a similar fashion to gradients.Future work will explore such methods as a potentially faster approach than L-BFGS.

Appendix A: Convergence Data for All Experiments
In Sec. 4, we showed the time of spatially averaged contrast vs. WFS&C iteration for an AD-PSM with η 00 ¼ 10 and ε ¼ 10 −4 compared to an experiment with SM (Fig. 9), as well as an AD-EFC experiment with α 2 k ¼ 10 −2 and ε ¼ 10 −4 compared to an EFC experiment with the same value of α 2 k (Fig. 10).We then summarized the steady-state time series statistics (median, 10th percentile, and 90th percentile) for a series of nine AD-PSM and nine AD-EFC runs with different combinations of ðη 00 ; εÞ and ðα 2 k ; εÞ, respectively (Figs. 11 and 12).In this section, Fig. 13 shows the full time series of contrast versus iteration for all 18 AD-PSM and AD-EFC experiments compared to their respective SM and EFC reference experiments.

Appendix B: Equivalence of Jacobian-Based Solutions and Newton's Method
Finding solutions for EFC and SM using the Jacobian matrix is equivalent to minimizing their respective cost functions with respect to a k using Newton's method.Newton's method is a second-order optimization technique that utilizes second-derivative information about the cost function, given by the local Hessian matrix H k ða n k Þ at any point a n k in the DM command parameter space.Given an initial guess for the solution a 0 k , the full Newton update is given as 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 7 ; 7 1 1 For cost functions that are exactly quadratic, including EFC and SM, Newton's method converges in a single iteration.
For general numerical optimization problems, Newton's method is rarely used in practice because forming the Hessian matrix explicitly is expensive.On the other hand, quasi-Newton methods, such as the BFGS algorithm or the limited-memory BFGS (L-BFGS) variant, can approximate H −1 k using changes in ∂J k ∕∂a T k over successive optimization iterations.As a consequence, they are substantially less computationally expensive.Although quasi-Newton algorithms do not converge as rapidly as Newton's method, and in particular can converge slowly for poorly conditioned problems, they are nonetheless superior to purely first-order methods such as steepest descent. 9s we show below, the Hessian matrix for EFC and SM has an analytic expression in terms of the Jacobian G k given by where C is a symmetric, positivedefinite matrix given by Γ T k Γ k for EFC and I∕μ n k for SM.Our approach is to replace this full Newton iteration, requiring computation of the Jacobian, by a series of cheaper quasi-Newton iterations instead, requiring only computation of the gradient ∂J k ∕∂a T k , which we achieve using AD.
The cost function for the EFC algorithm described in Sec. 2 can be written in the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 7 ; 4 6 1 where Γ k is a regularization matrix; for SM, Γ k ¼ I∕ ffiffiffiffi ffi μ k p .For now, we will restrict our attention to the EFC algorithm, but note that the result derived here applies equally as well to minimizing the Lagrangian function for SM with respect to the DM correction a k .
Expanding Eq. ( 18) and recalling that E DZ;k ¼ G k a k þ Êab;k , the cost function has the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 7 ; 3 7 9 where we use the fact that a k is purely real to discard ImfG † k G k g.The Jacobian-based solution is found by finding a k such that ∂J EFC;k ∕∂a k vanishes.We therefore begin by writing down the gradient: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 7 ; 3 0 8 This is a linear system of equations that we can solve for the optimal correction a Ã k : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 7 ; 2 5 6 The Hessian matrix is given as ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 7 ; 2 1 1 Since both terms in H k are positive definite, H k is positive definite as well, confirming that the solution is a minimum of the cost function.
We will now show that the solution a Ã k obtained above is the same as the solution obtained by applying a single iteration of Newton's method to the EFC cost function.Let a 0 k be an initial guess for the solution.Newton's method produces an iterate of the form 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 7 ; 1 1 2 Combining Eqs. ( 20) and ( 22) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 1 1 4 ; 7 2 4 Inserting back into Eq.( 23) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 a ; 1 1 4 ; 6 7 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 b ; 1 1 4 ; 6 3 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 c ; 1 1 4 ; 6 1 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 d ; 1 1 4 ; 5 9 2 Inserting the definition of the Hessian matrix from Eq. ( 22), we see that the Newton iterate a 1 k is identical to the analytical solution in Eq. ( 21) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 a ; 1 1 4 ; 5 5 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 b ; 1 1 4 ; 5 2 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 c ; 1 1 4 ; 5 0 0 As we described earlier, the same result holds if minimizing L SM;k in Sec.2.1 with respect to a k .

Appendix C: Fast Convolutional DM Model
Consider a DM with N A actuators along each side (i.e., N act ¼ N 2 A ) whose surface sðx; yÞ can be modeled as a linear superposition of identical influence functions fðx; yÞ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 1 1 4 ; 4 2 4 sðx; yÞ ¼ For fixed actuator spacing along the horizontal and vertical directions, we can rewrite the above summation as a convolution between a weighted Dirac comb function and the influence function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 1 1 4 ; 3 4 5 sðx; yÞ ¼ fðx; yÞ Ã Fourier transforming both sides transforms the convolution operation into a multiplication: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 a ; 1 1 4 ; 2 8 9 We define sðf x ; f y Þ ≜ F fsðx; yÞg and fðf x ; f y Þ ≜ F ffðx; yÞg: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 1 1 4 ; 1 9 8 We next define the discretized surface and influence function arrays s and f such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 a ; 1 1 4 ; 1 4 0 s½p; q ¼ sðpΔf x ; qΔf y Þ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 b ; 1 1 4 ; 1 0 4 where ∘ denotes element-wise multiplication.Finally, we define the vectors of actuator center coordinates x c and y c such that x c ½m ¼ x m and y c ½n ¼ y n , as well as the array of actuator commands A for which A½m; n ¼ a m;n E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 a ; 1 1 7 ; 6 6 3 s½p; q ¼ f½p; q ∘ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 b ; 1 1 7 ; 6 0 9 A½m; n expf−i2πy c ½nf y ½qg: (33b) We can write this more succinctly as the following sequence of element-wise and matrix products E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 4 ; 1 1 7 ; 5 7 1 where exponentiation is performed element-wise and ab T denotes the outer product of the vectors a and b.
The term in the parentheses is more commonly referred to as the matrix Fourier transform 21 or matrix triple product Fourier transform 25 of A, which we denote as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 5 ; 1 1 7 ; 4 9 7 MFTfA; x c ; y c ; f ) ; t e m p : i n t r a l i n k -; e 0 3 6 ; 1 1 7 ; 4 6 0 The final step is to compute an inverse discrete Fourier transform to obtain the desired discrete DM surface s, which is carried out most efficiently using the inverse fast Fourier transform, yielding the final result E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 7 ; 1 1 7 ; 3 9 8 s ¼ IFFTf f ∘ MFTfA; x c ; y c ; f x ; f y gg: For DMs whose active actuators are a subset of the N A × N A grid modeled above, only the elements of A corresponding to active actuators are set to nonzero values.

C.1 Adjoint Model
The algorithm described in the previous section computes the DM surface resulting from a twodimensional array of actuator commands A, under the assumptions that the influence function is identical across all actuators and that the surface can be approximated as a linear superposition of the actuator influence functions.In the context of gradient-based nonlinear optimization using RMAD, we can derive an adjoint model for this algorithm that computes the derivative A ≜ ∂J∕∂A T for some scalar cost function J, given the derivative s with respect to the surface s.
To begin, we break the forward model into the following sequence: We now apply the RMAD gradient rules 7,26 to each step in reverse order to derive the adjoint model, letting x ≜ ∂J k ∕∂x T for any variable x: E -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 a ; 1 1 4 ; 4 8 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 b ; 1 1 4 ; 4 4 3 Combining and simplifying, we see that the desired gradient is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 4 ; 1 1 4 ; 3 4 5 This gradient is passed to the next block of the adjoint model, which in this case is the adjoint model for propagation through the coronagraph, which, referring to Fig. 4, evaluates the derivatives of J EFC;k with respect to the surfaces s 1;k and s 2;k , respectively.We refer to our earlier work for a derivation of the coronagraph propagation adjoint model. 7

Disclosures
The authors declare no conflicts of interest.

Fig. 5
Fig. 5 Simplified, partially-unfolded layout of the HiCAT testbed.The elements encountered by the primary imaging beam path are highlighted in bold.DM1 and DM2 indicate the in-pupil and out-ofpupil DMs, respectively.Our experiments utilized a coronagraph configuration in which the apodizer is replaced with a flat mirror (see Sec. 3.1).

Fig. 7 Fig. 6
Fig. 7 Experimental on-axis images from HiCAT.(a) Non-coronagraphic image, (b) coronagraphic image prior to the WFS&C loop, (c) coronagraphic image after 80 iterations of SM, and corresponding actuator commands for the in-pupil DM (d), and out-of-pupil DM (e).In panels (b) and (c), the geometrical edge of the focal-plane mask (FPM) with radius 3.34 λ 0 ∕D LS is shown as well as the inner and outer edges of the dark zone at 5.8 λ 0 ∕D LS to 9.8 λ 0 ∕D LS , respectively, where D LS is the Lyot stop diameter and λ 0 ¼ 638 nm.

Fig. 10 (Fig. 9
Fig. 10 (a) Spatially averaged dark-zone contrast versus WFS&C iteration using AD-EFC with α 2 k ¼ 10 −2 and ε ¼ 10 −4 , compared to an experiment with EFC using the same value for α k .The median, 10th, and 90th percentile of the converged datapoints are shown in green.(b) The on-axis PSF corresponding to the iteration with deepest contrast using AD-EFC.

E 7 ; 2 3 8 Ã 1 s 5 Ã 4 A 2 A 6 E 8 J
Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 8 a ; 1 1 ¼ MFTfA; x c ; y c ; f x ; f y g; (38a) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 8 b ; 1 1 7 ; 2 0 0s ¼ f ∘ Ã; (38b) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0This leads to the following adjoint model, following the RMAD adjoint variable rules in Refs.7 and 26:E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 9 a ; 1 1 7 ; 1 5 ¼ FFTfsg; (39a) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 3 9 b ; 1 1 7 ; 1 1 ¼ fÃ ∘ s;(39b)E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 9 c ; 1 1 7 ; 9 ¼ IMFTf Ã; f x ; f y ; x c ; y c g;(39c)where Ã denotes element-wise complex conjugation and IMFT denotes the inverse matrix Fourier transform: Will et al.: High-order coronagraphic wavefront control with algorithmic differentiation. . .J. Astron.Telesc.Instrum.Syst.045004-17 Oct-Dec 2023 • Vol.9(4) Downloaded From: https://www.spiedigitallibrary.org/journals/Journal-of-Astronomical-Telescopes,-Instruments,-and-Systems on 02 Jan 2024 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 0 ; 1 1 4 ; 7 3 6 IMFTf Ã; f x ; f y ; x c ; y c g ≜ expfi2πx c f T x g Ã expfi2πf y y T c g: (40) Combining these expressions, the adjoint model is then E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 1 ; 1 1 4 ; 7 0 ¼ IMFTf fÃ ∘ FFTfsg; f x ; f y ; x c ; y c g: (41) 9 Appendix D: Adjoint Model for EFC Cost FunctionIn Sec.2.2, we describe the cost function for the EFC algorithm for a single correction wavelength.Here, we derive its RMAD adjoint model, which computes the derivative ∂J EFC;k ∕∂E DM;k .We begin by writing the cost function as a series of operations evaluated sequentially:E Q -T A R GE T ; t e m p : i n t r a l i n k -; e 0 4 2 a ; 1 1 4 ; 6 0 DZ;k ¼ E DM;k þ Êab;k ; (42a) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 b ; 1 1 4 ; 5 6 E DZ;k ¼ kE DZ;k k 2 ; (42b) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 c ; 1 1 4 ; 5 4 8c k ¼ Γ k a k ;(42c)E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 d ; 1 1 4 ; 5 3 0J c;k ¼ kc k k 2 ;(42d)E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 e ; 1 1 4 ; 5 1 1

E 5 E
Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 c ; 1 1 4 ; 4 2 3c k ¼ 2c k J c;k ; (43c) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 d ; 1 1 4 ; 4 0 4 a k ¼ Γ T k c k ;(43d)E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 e ; 1 1 4 ; 3 8 5 ĒDZ;k ¼ 2E DZ;k J E DZ;k ; (43e) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 f ; 1 1 4 ; 3 6 DM;k ¼ E DZ;k : RefE DM;k ðu 1 Þg½m ImfE DM;k ðu 1 Þg½m