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

HAL
open science

# High-order myopic coronagraphic phase diversity (COFFEE) for wave-front control in high-contrast imaging systems

B. Paul, L. M. Mugnier, J. -F. Sauvage, K. Dohlen, M. Ferrari

► To cite this version:

B. Paul, L. M. Mugnier, J. -F. Sauvage, K. Dohlen, M. Ferrari. High-order myopic coronagraphic phase diversity (COFFEE) for wave-front control in high-contrast imaging systems. Optics Express, 2013, 21 (26), pp.31751–31768. 10.1364/OE.21.031751 . hal-01439528

HAL Id: hal-01439528

https://hal.science/hal-01439528

Submitted on 28 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
High-order myopic coronagraphic phase
diversity (COFFEE) for wave-front
control in high-contrast imaging systems.

B. Paul,<sup>1,2,3,*</sup> L. M. Mugnier,<sup>1,3</sup> J.-F. Sauvage<sup>1,3</sup> K. Dohlen<sup>2,3</sup> and M. Ferrari<sup>2,3</sup>

**Abstract:** The estimation and compensation of quasi-static aberrations is mandatory to reach the ultimate performance of high-contrast imaging systems. COFFEE is a focal plane wave-front sensing method that consists in the extension of phase diversity to high-contrast imaging systems. Based on a Bayesian approach, it estimates the quasi-static aberrations from two focal plane images recorded from the scientific camera itself. In this paper, we present COFFEE's extension which allows an estimation of low and high order aberrations with nanometric precision for any coronagraphic device. The performance is evaluated by realistic simulations, performed in the SPHERE instrument framework. We develop a myopic estimation that allows us to take into account an imperfect knowledge on the used diversity phase. Lastly, we evaluate COFFEE's performance in a compensation process, to optimize the contrast on the detector, and show it allows one to reach the 10<sup>-6</sup> contrast required by SPHERE at a few resolution elements from the star. Notably, we present a non-linear energy minimization method which can be used to reach very high contrast levels (better than 10<sup>7</sup> in a SPHERE-like context).

© 2013 Optical Society of America

OCIS codes: (010.7350) Wave-front sensing, (100.5070) Phase retrieval, (110.1080) Active or adaptive optics, (100.3190) Inverse problems, (350.1260) Astronomical optics.

References and links

1. J.-L. Beuzit, M. Feldt, K. Dohlen, D. Mouillet, P. Puget, J. Antici, P. Baudoz, A. Boccaletti, M. Carbillet, J. Char-
ton, R. Claudi, T. Fusco, R. Gratton, T. Henning, N. Hubin, F. Joos, M. Kasper, M. Langlois, C. Moutou, J. Pragt,
P. Rabou, M. Saisse, H. M. Schmid, M. Turatto, S. Udry, F. Vakili, R. Waters, and F. Wildi, "SPHERE: A Planet
Finder Instrument for the VLT," in "Proceedings of the conference In the Spirit of Bernard Lyot: The Direct
Detection of Planets and Circumstellar Disks in the 21st Century.", P. Kalas, ed. (University of California,
Berkeley, CA, USA, 2007).

2. B. A. Macintosh, J. R. Graham, D. W. Palmer, R. Doyon, J. Dunn, D. T. Gavel, J. Larkin, B. Oppenheimer, L. Saddlermyer, A. Sivaramakrishnan, J. K. Wallace, B. Bauman, D. A. Erickson, C. Marois, L. A. Poyneer, and R. Soummer, "the gemini planet imager: from sience to design to construction," in "Adaptive Optics Systems", vol. 7015 (Proc. Soc. Photo-Opt. Instrum. Eng., 2008), vol. 7015.

<footnote n="1"><sup>1</sup>Onera - The French Aerospace Lab,
F-92322 Chatillon, France</footnote>

<footnote n="2"><sup>2</sup>Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille),
UMR 7326, 13388, Marseille, France</footnote>

<footnote n="3"><sup>3</sup>Groupement d'intérêt scientifique PHASE (Partenariat Haute résolution Angulaire Sol et Espace) between Onera, Observatoire de Paris, CNRS, Université Diderot, Laboratoire d'Astrophysique de Marseille and Institut de Planétologie et d'Astrophysique de Grenoble,
France</footnote>

<footnote>baptiste.paul@onera.fr</footnote>
3. O. Guyon, F. Martinache, V. Garrel, F. Vogt, K. Yokochi, and T. Yoshikawa, "The subaru coronagraphic extreme ao (SCExAO) system: wavefront control and detection of exoplanets with coherent light modulation in the focal plane," in "Adaptive Optics Systems II," vol. 7736 (Proc. Soc. Photo-Opt. Instrum. Eng., 2010), vol. 7736.

4. S. Hinkley, B. R. Oppenheimer, N. Zimmerman, D. Brenner, I. R. Parry, J. R. Crepp, G. Vasisht, E. Ligon, D. King, R. Soummer, A. Sivaramakrishnan, C. Beichman, M. Shao, L. C. Roberts Jr., A. Bouchez, R. Dekany, L. Pueyo, J. E. Roberts, T. Lockhart, C. Zhai, C. Shelton, and R. Burruss, "A new high contrast imaging program at palomar observatory," Pub. Astron. Soc. Pacific **123** (2011).

5. R. Soummer, A. Ferrari, C. Aime, and L. Jolissaint, "Speckle noise and dynamic range in coronagraphic images," Astrophys. J. **669** (2007).

6. J.-F. Sauvage, T. Fusco, C. Petit, L. M. Mugnier, B. Paul, and A. Costille, "Focal-plane wave front sensing strategies for high-contrast imaging. experimantal validations on SPHERE," in "Adaptive Optics Systems III," , vol. 8447 (Proc. Soc. Photo-Opt. Instrum. Eng., 2012), vol. 8447.

7. J. K. Wallace, R. S. Burruss, R. D. Bartos, T. Q. Trinh, L. A. Pueyo, S. F. Fregoso, J. R. Angione, and J. C. Shelton, "the gemini planet imager calibration wavefront sensor instrument," in "Adaptive Optics Systems II," , vol. 7736 (Proc. Soc. Photo-Opt. Instrum. Eng., 2010), vol. 7736.

8. P. J. Bordé and W. A. Traub, "High contrast imaging from space : speckle-nulling in a low-aberration regime," Astrophys. J. **638** (2006).

9. A. Give'on, R. Belikov, S. Shaklan, and J. Kasdin, "Closed loop, DM diversity based, wave-front correction algorithm for high contrast imaging systems," Opt. Express **15** (2007).

10. J. Trauger, A. Give'on, B. Gordon, B. Kern, A. Kuhnert, D. Moody, A. Niesser, F. Shi, D. Wilson, and C. Burrows, "Laboratory demonstration of high-contrast imaging for space coronagraphy," in "Tecniques and Instrumentation for Detection of Exoplanets III," , vol. 6693 (Proc. Soc. Photo-Opt. Instrum. Eng., 2010), vol. 6693.

11. S. J. Thomas, A. A. Give'On, D. Dillon, B. Macintosh, D. Gavel, and R. Soummer, "Laboratory test of aplication of electric field conjugation image-sharpening to ground-based adaptive optics," in "Adaptive Optics Systems II," , vol. 7736 (Proc. Soc. Photo-Opt. Instrum. Eng., 2010), vol. 7736.

12. P. Baudoz, J. Mazoyer, M. Mas, R. Galicher, and G. Rousset, "Dark hole and planet detection: laboratory results using the self-coherent camera," in "Ground-based and Airborn Intrusmentation for Astronomy IV," , vol. 8446 (Proc. Soc. Photo-Opt. Instrum. Eng., 2012), vol. 8446.

13. P. Baudoz, A. Boccaletti, J. Baudrand, and D. Rouan, "The self-coherent camera : a new tool for exoplanet detection," in "Proc. IAU Colloquium," (2006).

14. J.-F. Sauvage, L. M. Mugnier, B. Paul, and R. Villecroze, "Coronagraphic phase diversity: a simple focal-plane wavefront sensor," Opt. Lett. **37** (2012).

15. B. Paul, J.-F. Sauvage, and L. M. Mugnier, "Coronagraphic phase diversity: performance study and laboratory demonstration," Astron. Astrophys. **552** (2013).

16. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes : the Art of Scientific Computing* (Cambridge University, 2007).

17. E. Thiébaut, "Optimization issues in blind deconvolution algorithms," in "Astronomical Data Analysis II," , vol. 4847 (Proc. Soc. Photo-Opt. Instrum. Eng., 2002), vol. 4847, pp. 174–183.

18. K. Dohlen, F. Wildi, J. Beuzit, P. Puget, D. Mouillet, A. Baruffolo, A. Boccaletti, J. Charton, R. Claudi, A. Costille, P. Feautrier, M. Feldt, T. Fusco, R. Gratton, M. Kasper, M. Langlois, D. L. Mignant, J. Lizon, N. Hubin. A. Pavlov, C. Petit, J. Pragt, P. Rabou, S. Rochat, R. Roelfsema, J. Sauvage, and H. Schmid, "Sphere system analysis predictions," (AO4ELT 2, 2011).

19. J.-F. Sauvage, L. Mugnier, T. Fusco, and G. Rousset, "Post processing of differential images for direct extrasolar planet detection from the ground," in "Advances in Adaptive Optics II," , vol. 6272, L. Ellerbroek, B. and D. Bonaccini Calia, eds. (Proc. Soc. Photo-Opt. Instrum. Eng., 2006), vol. 6272.

20. S. Bongard, F. Soulez, E. Thiébaut, and E. Pecontal, "3D deconvolution of hyper-spectral astronomical data," Mon. Not. R. Astr. Soc. **418** (2011).

21. R. Soummer, L. Pueyo, A. Sivaramakrishnan, and R. J. Vanderbei, "Fast computation of lyot-style coronagraph propagation," Opt. Express **15** (2007).

22. L. M. Mugnier, A. Blanc, and J. Idier, "Phase diversity: a technique for wave-front sensing and for diffraction-limited imaging," in "Advances in Imaging and Electron Physics," , vol. 141, P. Hawkes, ed. (Elsevier, 2006), chap. 1, pp. 1–76.

23. D. J. Lee, M. C. Roggemann, B. M. Welsh, and E. R. Crosby, "Evaluation of least-squares phase-diversity technique for space telescope wave-front sensing," Astrophys. J. **36** (1997).

24. A. Blanc, T. Fusco, M. Hartung, L. M. Mugnier, and G. Rousset, "Calibration of NAOS and CONICA static aberrations. Application of the phase diversity technique," Astron. Astrophys. **399**, 373–383 (2003).

25. I. Mocœur, L. M. Mugnier, and F. Cassaing, "Analytical solution to the phase-diversity problem for real-time wavefront sensing," Opt. Lett. **34**, 3487–3489 (2009).
# 1. Introduction

Exoplanet imaging is one of the most challenging areas in today's astronomy. The observation of an extremely faint object (the planet) very close to a bright source (the host star) requires the use of an extreme adaptive optics (XAO) system coupled with a high-contrast imaging technique such as coronagraphy. The current generation of instruments dedicated to exoplanets direct imaging (SPHERE on the VLT [1], GPI on Gemini South [2], Subaru SCExAO [3] and Palomar P1640 [4]) aim at detecting massive gaseous planets $10^{-6}$ to $10^{-7}$ times fainter than their host star. In the future, high-contrast imaging instruments on ground based or space based telescopes will perform observation of Earth-like planets, $10^{-9}$ to $10^{-10}$ times fainter than their host star.

The ultimate limitation of a high-contrast imaging instrument lies in its quasi-static aberrations, which originate in imperfections of the optical system such as misalignment or optical surface polishing error. These aberrations, when unseen and thus uncorrected by the AO loop, create long-lived speckles on the detector plane [5], limiting the achievable contrast. Besides, unlike the signal which originates in the residual turbulence (averaged in a long exposure image), these speckles can easily be mistaken for a planet. Thus, to reach the ultimate performance of the imaging system, one must be able to compensate for these aberrations. To perform such a compensation, SPHERE relies on phase diversity [6] to reach a contrast of $10^{-6}$, whereas GPI relies on an interferometric approach [7] for an aimed contrast of $10^{-7}$ on the detector.

Several techniques dedicated to this compensation have been proposed. Closed loop methods, which assume small aberrations ([8–12]), estimate the electric field in the detector plane using at least three images. The technique proposed by Baudoz et al. [13] relies on a modification of the imaging system, but requires only one image. We note that this approach, based on the analysis of fringed speckles, requires a $\sqrt{2}$ oversampling of the coronagraphic images to properly sample the interference fringes. These techniques aim at minimizing the energy in a chosen area (“Dark Hole”), leading to a contrast optimization on the detector in a closed loop process.

The focal plane wave-front sensor we have proposed [14], called COFFEE (for Coronagraphic Focal-plane wave-Front Estimation for Exoplanet detection), requires only two focal-plane images to estimate the aberrations both upstream and downstream of the coronagraph without any modification of the coronagraphic imaging system or assuming small aberrations. In a previous study [15], we presented COFFEE’s early performance and limitations, detailed below, as well as its sensitivity to a realistic experimental environment. In this paper, we present a high order extension of COFFEE and its performance evaluation in a compensation process, in the framework of the quasi-static aberration calibration of a ground-based instrument such as SPHERE. Section 2 presents the modifications that allows COFFEE to overcome its previous limitations [15], allowing an estimation of high order aberrations with nanometric precision for any coronagraphic device. Section 3 presents the noise sensitivity of this extended version of COFFEE using realistic SPHERE-like simulations. The choice of a suitable diversity phase to use is also discussed in this section. Knowing that in a real system, such a diversity phase will not be perfectly introduced, we present, in Section 4, an original approach, hereafter called “myopic”, which significantly improves COFFEE’s robustness to an imprecise knowledge of the diversity phase. Finally, in Section 5, we describe the different compensation processes which can be used with COFFEE. In particular, we present a method of energy minimization in the detector plane that allows the creation of a Dark Hole without any small aberration assumption. Unlike other energy minimization methods [8–12], the one we propose does not rely on the calibration of an interaction matrix, which is sensitive to the position of the coronagraphic image on the detector. This new dark hole method can thus be used on any high contrast instrument without a repetitive dedicated calibration step.
In the SPHERE baseline design, quasi-static aberrations are measured with conventional phase diversity [6] (no coronagraph), which is unable to sense high-order aberrations to a nanometric level. As a result, this high contrast imaging instrument performance will be limited by high-order phase aberrations and not by amplitude aberrations (amplitude variations in a pupil plane). Consequently, the latter are not considered in the simulations presented herein and are currently not estimated by COFFEE.

## 2. Aberrations estimation with COFFEE

COFFEE is based on an inverse problem approach: it estimates the aberrations both upstream and downstream of the coronagraph using two focal-plane images that differ from a known aberration. As described in [15], the two main error sources of the first version of COFFEE have been shown both by simulations and experimentally to be aliasing and modelling error. The former was due to the use of a Zernike basis and prevented COFFEE from estimating high order aberrations. Moreover, the estimation of these high order aberrations is mandatory to optimize the contrast in the detector far from the optical axis (between a few $\lambda/D$ and $20 \lambda/D$ in the case of the SPHERE instrument).

The latter was originating in the image formation model used by COFFEE. The estimations were indeed performed using a perfect coronagraph model, and thus limited by a model error. In practice, COFFEE's use was limited to the apodized Roddier & Roddier coronagraph.

In this Section, we present COFFEE's modifications that allows to get rid of these two limitations. Section 2.1 describes the modification of the maximum a posteriori (MAP) approach on which COFFEE is based, which includes a modification of the basis used for the aberration estimation, now composed of pupil indicator functions (pixels). Such a basis, used with a dedicated regularization metric described in Section 2.2, allows COFFEE to estimate high-order aberrations. Besides, thanks to the modification of the imaging model, described in Section 2.3, COFFEE is now able to perform the estimation for any coronagraphic device.

### 2.1. Criterion expression

Most of the notations of this article are coherent with [14, 15]. We consider a coronagraphic imaging system made of four successive planes denoted by A (circular entrance pupil of diameter $D_u$), B (coronagraphic focal plane), C (Lyot Stop), and D (detector plane). The optical aberrations are considered as static and introduced in pupil planes A and C. The coronagraphic device is composed of a focal plane mask located in plane B and a Lyot Stop in plane C. No particular assumption is made on the pupil shape or intensity, which can be calibrated using data recorded from the instrument. We note that this model does not consider out of plane aberrations and the corresponding amplitude aberrations (which originate in Fresnel effect), but as said previously, we consider in this paper the case of a ground based instrument (such as SPHERE) limited by phase aberrations.

COFFEE requires only two images $i_c^{\text{foc}}$ and $i_c^{\text{div}}$ recorded on the detector (plane D) that, as in phase diversity, differ from a known aberration $\phi_{\text{div}}$, to estimate aberrations both upstream ($\phi_u$) and downstream ($\phi_d$) of the coronagraph.

In this paper, we consider the case of the instrument calibration, assumed to be performed at high signal to noise ratio (SNR) value, with a monochromatic source, emitted from a single-mode laser fiber. Since the impact of the source finite size on the estimation is not significant on the aberrations estimation [15], we consider here that this calibration is performed with an unresolved object, and use the following imaging model:

$$i_c^{\text{foc}} = \alpha_{\text{foc}}h_{\text{det}} * h_c(\phi_u, \phi_d) + n_{\text{foc}} + \beta_{\text{foc}} \quad (1)$$

$$i_c^{\text{div}} = \alpha_{\text{div}}h_{\text{det}} * h_c(\phi_u + \phi_{\text{div}}, \phi_d) + n_{\text{div}} + \beta_{\text{div}}$$
where $\alpha_p$ is the incoming flux (p is for “foc” or “div”), $h_c$ the coronagraphic on axis “point spread function” (PSF) of the instrument (which is the response of a coronagraphic imaging system to a point source), $h_{det}$ the known detector PSF, $n_{foc}$ and $n_{div}$ are the measurement noises and comprise both detector and photon noises, $\beta_p$ is a unknown uniform background (offset), and $\star$ denotes the discrete convolution operation.

COFFEE is based on a maximum a posteriori (MAP) approach: it estimates the aberrations $\phi_u$ and $\phi_d$ as well as the fluxes $\alpha = [\alpha_{foc}, \alpha_{div}]$, and the backgrounds $\beta = [\beta_{foc}, \beta_{div}]$ that maximize the posterior likelihood $p(\alpha, \beta, \phi_u, \phi_d | i_c^{foc}, i_c^{div})$ of the data. For practical issues, it is more convenient (and equivalent) to minimize the opposite of the logarithm of the posterior likelihood, or neg-log-likelihood $J(\alpha, \beta, \phi_u, \phi_d) = -\ln[p(\alpha, \beta, \phi_u, \phi_d | i_c^{foc}, i_c^{div})]$ which includes regularization terms $\mathcal{R}(\phi_u)$ and $\mathcal{R}(\phi_d)$ designed to enforce smoothness of the sought phases:

$$ (\hat{\alpha}, \hat{\beta}, \hat{\phi}_u, \hat{\phi}_d) = \arg\min_{\alpha, \beta, \phi_u, \phi_d} [J(\alpha, \beta, \phi_u, \phi_d)] \quad (2) $$

where

$$ J(\alpha, \beta, \phi_u, \phi_d) = \frac{1}{2} \left\| \frac{i_c^{\text{foc}} - (\alpha_{\text{foc}}h_{\text{det}} \star h_c(\phi_u, \phi_d) + \beta_{\text{foc}})}{\sigma_n^{\text{foc}}} \right\|^2 \\ + \frac{1}{2} \left\| \frac{i_c^{\text{div}} - (\alpha_{\text{div}}h_{\text{det}} \star h_c(\phi_u + \phi_{\text{div}}, \phi_d) + \beta_{\text{div}})}{\sigma_n^{\text{div}}} \right\|^2 \\ + \mathcal{R}(\phi_u) + \mathcal{R}(\phi_d) \quad (3) $$

$\|x\|^2$ denotes the sum of squared pixel values of map $x$, $\sigma_n^{\text{foc}}$, and $\sigma_n^{\text{div}}$ are the noise standard deviation maps of each image. The corresponding variances can be computed as a sum of the photon and detector noise variances. The former can be estimated as the image itself thresholded to positive values, and the latter can be calibrated prior to the estimation.

Any aberration $\phi$ is expanded on a basis $\{b_m\}$. In [15], we showed that the use of a truncated Zernike basis for the reconstruction led to a strong aliasing error, let alone the inability to estimate high frequency aberrations. In this paper, the phase is expanded on pixel indicator functions in the pupil plane: $\phi = \sum_m \phi^m b_m$ (with $\phi^m$ the value of the m-th pixel in the pupil). Such a basis, used with the proper regularization metrics, will allow COFFEE to estimate high order aberrations and strongly reduce the aliasing error, as shown in the following.

The minimization of metric $J(\alpha, \beta, \phi_u, \phi_d)$ of Eq. (3) is performed by means of a limited memory variable metric (BFGS) method ([16, 17]), which is a fast quasi-Newton type minimization method. It uses the analytical expression of gradients $\frac{\partial J}{\partial \phi_u}$, $\frac{\partial J}{\partial \phi_d}$, $\frac{\partial J}{\partial \alpha}$ and $\frac{\partial J}{\partial \beta}$, which we have calculated, to estimate $\phi_u$, $\phi_d$, $\alpha$ and $\beta$ (Implementation details can be found in Appendix A).

Sauvage et al. [14] established that a suitable diversity phase $\phi_{\text{div}}$ for COFFEE was a mix of defocus and astigmatism: $\phi_{\text{div}} = a_4^{\text{div}}Z_4 + a_5^{\text{div}}Z_5$ with $a_4^{\text{div}} = a_5^{\text{div}} = 0.8$ rad RMS, introduced upstream of the coronagraph. In this paper, such a diversity phase will be used for a start; the optimal phase diversity to be used with COFFEE will be discussed later.

## 2.2. Regularization metric

The use of a pixel basis for the phase reconstruction is required for COFFEE to estimate high order aberrations. However, this leads to a large number of unknowns, which in turn calls for a regularization metric in order to reduce the noise sensitivity. We chose a regularization metric that is based on the available a priori knowledge on the quasi-static aberrations. Indeed, they can be reasonably assumed to be Gaussian, homogeneous and thus endowed with a power spectral density (PSD) $S_{\phi_k}$ (where k stands for u (upstream) or d (downstream)), which is usually
assumed to follow a power law:

$$ \begin{cases} S_{\phi_k} \propto \frac{1}{v^n} \\ \langle \phi_k \rangle = 0 \end{cases} \qquad (4) $$

with v the spatial frequency and $\langle \phi_k \rangle$ the mean of $\phi_k$. The regularization term $\mathcal{R}(\phi_k)$ can thus be written as:

$$ \mathcal{R}(\phi_k) = \frac{1}{2} \sum_{v} \frac{|\mathcal{F}[\phi_k](v)|^2}{S_{\phi_k}(v)}, \qquad (5) $$

where $\mathcal{F}$ represent the Fourier transform operation. In order to be able to cope with any pupil shape, we implement this metric in direct space rather than Fourier space as follows; for n=2 in Eq. (4), we obtain:

$$ \mathcal{R}(\phi_k) = \frac{\mu_k}{2} \| \nabla \phi_k(r) \|^2. \qquad (6) $$

And for n=4 in Eq. (4),

$$ \mathcal{R}(\phi_k) = \frac{\mu_k}{2} \| \Delta \phi_k(r) \|^2. \qquad (7) $$

Here, r denotes the pupil plane position which will be omitted in the following for the sake of simplicity. $\nabla$ and $\Delta$ represent the gradient and the Laplacian operators, respectively. The balance parameter $\mu_k$ will be called “hyperparameter” hereafter. For both cases, derivatives $\nabla\phi_k$ and $\Delta\phi_k$ are computed as finite differences between neighboring points, and summations are limited to points whose computation requires only pixels inside the pupil. In this paper, we consider a PSD decrease as $1/v^2$ ($n=1$ in Eq. (4)), which corresponds to a classical assumption for optical surface polishing errors, according to K. Dohlen et al. [18], who measured and characterized the PSD of the SPHERE optical system. Additionally, identification between Eqs. (5) and (6) yields the analytic value of the hyperparameter $\mu_k$:

$$ \mu_k = \frac{1}{\sigma_{\nabla\phi_k}^2}, \qquad (8) $$

where $\sigma_{\nabla\phi_k}^2$ is defined as $\sigma_{\nabla\phi_k}^2 = \sigma_{\nabla_x\phi_k}^2 + \sigma_{\nabla_y\phi_k}^2$, with $\sigma_{\nabla_x\phi_k}^2$ and $\sigma_{\nabla_y\phi_k}^2$ the variances of $\nabla(\phi_k)$ in directions x and y, respectively. The fact that the hyperparameter is given by Eq. (8) stems from the assumption that the phase $\phi_k$ is statistically homogeneous, and is whitened by the differentiation in Eq. (6) [19, 20]. One can notice that $\sigma_{\nabla\phi_k}^2$ can be analytically computed from $S_{\phi_k}$ and the phase variance $\sigma_{\phi_k}^2$. Thus, this regularization does not require any manual tuning.

## 2.3. Coronagraphic image formation model

To perform the minimization of criterion J in Eq. (3), the image formation model used by COFFEE (Eq. (1)) requires the expression of a coronagraphic PSF $h_c$. Let r be the pupil plane position vector and $\gamma$ the focal plane position vector. The entrance pupil function $P_u$ is such that:

$$ P_u(r) = \Pi\left(\frac{2r}{D_u}\right) \Phi(r), \qquad (9) $$

with $\Pi$ the disk of unit radius, $D_u$ the entrance pupil diameter, and $\Phi$ a known apodization function. The electric field in the entrance pupil can be written as:

$$ \Psi_A(r) = P_u(r) e^{j\phi_u(r)}. \qquad (10) $$

The electric field in the detector plane $\Psi_D$ is obtained by propagating $\Psi_A$ through each plane of the coronagraphic imaging system: the signal is first focused on the coronagraphic focal
plane mask $\mathcal{M}$; then, the electric field is propagated through the Lyot Stop pupil $P_d(r)$ ($P_d(r) = \Pi(2r/D_d)$ with $D_d$ the Lyot Stop pupil diameter). The electric field in the detector plane $\Psi_D$ can thus be written as:

$$ \Psi_d(\gamma) = \mathcal{F}^{-1} \left\{ \mathcal{F} \left[ \mathcal{F}^{-1}(\Psi_A(r)) \cdot \mathcal{M} \right] P_d(r) e^{j\phi_d(r)} \right\}, \quad (11) $$

where $\mathcal{F}^{-1}$ is the inverse Fourier transform operation. For the sake of simplicity, spatial variables $r$ and $\gamma$ will be omitted in the following.

The coronagraphic PSF $h_c$ is the square modulus of $\Psi_D$:

$$ h_c = \left| \mathcal{F}^{-1} \left\{ \mathcal{F} \left[ \mathcal{F}^{-1}(\Psi_A) \cdot \mathcal{M} \right] P_d e^{j\phi_d} \right\} \right|^2 \quad (12) $$

In Eq. (12), $\mathcal{M}$ can easily be adapted to represent any coronagraphic device, allowing COF-FEE to be used with a broad class of high contrast imaging instruments.

### 3. Performance evaluation

This Section presents the performance of the new extension of COFFEE presented in Section 2. Section 3.1 gathers the different parameters used for these simulations. In Section 3.2, the impact of the noise on COFFEE's estimation is quantified, while the optimal phase diversity $\phi_{div}$ to be used with COFFEE is studied in Section 3.3. Section 3.4 presents COFFEE's sensitivity to a difference between the prior assumed in the phase reconstruction and the true phase.

#### 3.1. Parameters and criteria

<table><caption>Table 1 gathers the parameters used for the simulations presented in this section:</caption><tr><td>image size</td><td>64 × 64 $\lambda/D$ (128 × 128 pixels, Shannon-sampled)</td></tr><tr><td>Light spectrum</td><td>Monochromatic, wavelength $\lambda = 1589$ nm</td></tr><tr><td>Entrance pupil</td><td>$D_u = 64$ pixels</td></tr><tr><td>Lyot stop pupil</td><td>$D_d = D_u$</td></tr><tr><td>Aberration upstream of the coronagraph ($\phi_u$)</td><td>WFE<sub>u</sub> = 50 nm RMS</td></tr><tr><td>Aberration downstream of the coronagraph ($\phi_d$)</td><td>WFE<sub>d</sub> = 20 nm RMS</td></tr><tr><td>Coronagraph</td><td>Apodized Lyot Coronagraph (ALC), focal plane mask angular diameter $d = 4.52\lambda/D$</td></tr></table>

<caption>Table 1. COFFEE: simulation parameters</caption>

These parameters have been chosen so that the following simulations are representative of the SPHERE instrument. The chosen coronagraph (ALC) is the one designed for the considered wavelength on SPHERE, and the apodization function used in the image formation model (Figure 1) is the one designed for this coronagraph. In this paper, we consider the case of a high-contrast imaging instrument calibration prior to the scientific observation, so a monochromatic source is considered. It is worth mentioning that COFFEE could easily be adapted to polychromatic images, although such a study is beyond the scope of this paper. Such an adaptation would require a modification of the image formation model (Eq. (1)), in which a polychromatic coronagraphic PSF would be computed from several monochromatic PSF for different wavelengths. In order to properly model the ALC coronagraph, the coronagraphic PSF $h_c$ (Eq. (12)) is computed using the method developed by R. Soummer et al. [21]. This approach allows an accurate numerical representation of Lyot-style coronagraphs by accurately sampling the coronagraphic
<figure><img src="image_3.png" /><figcaption>Fig. 1. Apodized Lyot coronagraph: apodization function used in this paper, computed from an experimental image recorded on SPHERE</figcaption></figure>

focal plane mask, which can hardly be done using the common Fast Fourier Transform (FFT) algorithm. Such an operation would indeed require the manipulation of very large arrays.

For each simulation, coronagraphic images are computed from randomly generated aberrations $\phi_u$ and $\phi_d$ using the imaging model presented in Eqs (1) and (12). Using these two images, COFFEE performs the phase estimation by minimizing criterion $J$ of Eq. (3).

In order to quantify the reconstruction accuracy, we define the reconstruction error $\epsilon_k$ as the RMS value of $\phi_k - \hat{\phi}_k$, where $k$ is either $u$ (upstream) or $d$ (downstream), $\phi_k$ is the simulated aberration and $\hat{\phi}_k$ its estimation made by COFFEE. In this Section, every reconstruction error value is an average value, computed from ten independent randomly generated phases to make sure that the result obtained is independent of the phase realization. The PSD $S_\phi$ of each generated phase is such that $S_\phi \propto 1/v^2$.

## 3.2. Noise propagation

The ultimate limitation of the estimation performed by COFFEE lies in the propagation of noise present in the images. As mentioned in Section 2.2, the use of a suitable regularization metric ensures the smoothness of the phase, limiting the propagation of the noise from the images to the estimated aberrations. In this Section, we first demonstrate that the analytic value of the regularization metric hyperparameter (see Section 2.2) is the one that gives the smallest reconstruction error. We consider here an incoming flux $\alpha = 10^7$ photons and a read-out noise (RON) of standard deviation $\sigma_{det} = 5$ e$. This value, slightly higher than the expected RON on a SPHERE-like system, is chosen to strengthen the impact of the hyperparameter value. Besides, photon noise is added in the simulated coronagraphic images. The hyperparameter is here such as $\mu_k = \chi/\sigma_{\nabla\phi_k}^2$, with $\chi \in [10^{-3}; 10^3]$. A reconstruction error value is computed for each value of $\chi$.

Figure 2 shows that the analytic value of the hyperparameter ($\chi = 1$) is the one that gives the smallest reconstruction error for both upstream and downstream aberrations. When the regularization metrics are under balanced ($\chi < 1$), the prior information is not “strong” enough in the minimization to prevent noise propagation in the estimated aberrations. One can notice here that without a regularization metric ($\chi = 0$), the reconstruction error would have been unacceptable. On the other hand, when the regularization metrics are over balanced ($\chi > 1$), their impact is too strong, and prevents the estimation of the high frequency components of $\phi_u$ and $\phi_d$. In the following, considering the result of this simulation, all the estimations performed by COFFEE will be done using the analytic value of the hyperparameter ($\chi = 1$).

In Fig. 3, we present the evolution of the reconstruction errors with respect to the total incoming flux. As previously, photon noise and detector noise ($\sigma_{det} = 1$ e$^{-}$) are added in the simulated images used by COFFEE to perform the estimation.

In Fig. 3, one can see that in a low flux regime, both reconstruction errors upstream and
<figure><img src="image_4.png" /><figcaption>Fig. 2. Regularization metrics: reconstruction error upstream (red line) and downstream (blue line) of the coronagraph as functions of the hyperparameter value</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Reconstruction error (solid red line) upstream (left) and downstream (right) as a function of the incoming flux $\alpha$. For comparison, $1/\sqrt{\alpha}$ (magenta dashed line) theoretical behaviour is plotted for photon noise only. The dotted-dashed black line represents the WFE value upstream and downstream of the coronagraph.</figcaption></figure>
downstream reach a saturation level, which correspond to the WFE value, rather than becoming arbitrarily large. Such a behaviour originates in the regularization metric: in a low flux regime, the speckles that originate in aberrations upstream of the coronagraph, which are used by COFFEE to estimate the aberrations, can hardly be distinguished from the noise. The less these speckles are visible in the images, the less the corresponding aberrations can be estimated using COFFEE. The limit is reached when no speckles are visible: in such a case, the estimated phase tend to zero and the reconstruction error is equal to the WFE value. Notice that without regularization metric, the reconstruction errors would have been much stronger in a low flux regime, due to noise propagation (as presented in Figure 2, in the case of an under-balanced hyperparameter).

When the flux increases, the reconstruction error evolution is proportional to $1/\alpha$, which correspond to a photon noise limited regime. For very high flux values ($\alpha > 10^9$ photons), the upstream reconstruction error seems to reach another saturation level. We have shown that this saturation, slightly above 1 nm RMS can be attributed to numerical difficulties in the minimization due to the very high dynamic range of the noise variance in the criterion. If it were of practical interest, these difficulties could be solved to reach even better accuracies.

In this paper, we consider the case of a high-contrast imaging instrument calibration, performed off-line prior to the observation with a high SNR value. Thus, considering the noise propagation behaviour presented above, the following simulations will be performed with an incoming flux value $\alpha = 10^9$ photons and a read-out noise (RON) of standard deviation $\sigma_{det} = 1 e^{-}$.

### 3.3. Choice of a diversity phase

This section aims at studying the sensitivity of COFFEE to the diversity phase $\phi_{div}$, which was until now a mix of defocus and astigmatism: $\phi_{div} = a^{div}(Z_4 + Z_5)$ with $a^{div} = 0.8$ rad RMS (202 nm RMS at $\lambda = 1589$ nm), introduced upstream of the coronagraph. This choice has been made following Sauvage et al. [14], who demonstrated that for a perfect coronagraph model and low order aberrations, such a diversity phase allowed a suitable criterion shape for the minimization. Indeed, the use of this diversity phase instead of defocus alone enlarges the global minimum, leading to an easier criterion minimization. In this Section, we study the influence of the diversity phase on the reconstruction accuracy more thoroughly and for realistic high order aberrations and coronagraph.

In classical phase diversity (no coronagraph), the optimal diversity phase depends on several parameters such as the signal to noise ratio (SNR), the level of the aberrations and their PSD [22]. A theoretical work, based on the computation of the Cramer-Rao lower bound (following Lee et al. [23]) could be performed to determine an optimal diversity phase; however, such a study would assume that there are no local minima in the criterion. Since we know that such minima appears in the criterion when the diversity phase amplitude is small, we adopt, in this section, a more practical approach to determine a suitable diversity phase for an aberration estimation with nanometric precision.

We will consider different diversity phases: a diversity phase composed of defocus alone, $\phi_{div} = a_{div}Z_4$, and a diversity phase composed of a mix of defocus and astigmatism, $\phi_{div} = a^{div}(Z_4 + Z_5)$. For each diversity phase, the evolution of the reconstruction errors with the diversity phase amplitude $a_{div}$ value will be plotted for 3 different WFE$_u$ value upstream of the coronagraph. The parameters used in this simulation are gathered in Table 1.

Figure 4 shows the evolution of both reconstruction errors upstream and downstream of the coronagraph with respect to the amplitude $a^{div}$. Here, the reconstruction error is due to noise and to local minima, which are gradually removed when the diversity phase amplitude increases, leading to an improvement of the estimation accuracy. When the diversity phase amplitude
<figure><img src="image_6.png" /><figcaption>Fig. 4. Reconstruction errors upstream (left) and downstream (right) of the coronagraph as functions of the amplitude $a^{\text{div}}$ of a diversity phase composed of defocus only ($Z_4$, dashed lines), or composed of a mix of defocus and astigmatism ($Z_4 + Z_5$, solid lines).</figcaption></figure>

$a^{\text{div}}$ is high enough, all local minima are removed. Then, the reconstruction error reaches a saturation level which correspond to the level of noise in the images.

One can notice that with a diversity composed of defocus and astigmatism, the reconstruction error decreases faster than in the case of a diversity composed of a defocus alone. In order to have $\epsilon_u < 5$ nm RMS, $a_{\text{div}}$ must be greater than $2 \times \text{WFE}_u$ when the diversity is composed of defocus and astigmatism, and greater than $2.5 \times \text{WFE}_u$ when it is composed of defocus alone.

This result confirms what has been shown in a simple case with a perfect coronagraph model [14]: a diversity phase composed of a mix of defocus and astigmatism enlarges the global minimum and pushes away local minima, making criterion *J* (Eq. (3)) easier to minimize, and thus allowing a more accurate estimation of the aberrations both upstream and downstream of the coronagraph.

### 3.4. Sensitivity to the a priori assumption accuracy

As mentioned in Section 2.2, the *a priori* assumption about the aberration PSD considered by COFFEE is derived from [18]. However, this assumption will not be exact for all optical surface; for instance, the PSD that correspond to the polishing error of the Very Large Telescope primary mirrors follow a $1/v^3$ power law, as demonstrated by Bordé and Traub [8].

In this section, the sensitivity of COFFEE to the validity of the assumed PSD is evaluated: using parameters listed in Table 1, coronagraphic images are computed with aberrations generated with three different power laws: $1/v$, $1/v^2$ and $1/v^3$. Then, using these simulated images, COFFEE perform the aberration estimation assuming a PSD following a $1/v^2$ power law (as described in Section 2.2).

<table><caption>Table 2. Comparison of COFFEE's estimation accuracy when the *a priori* knowledge on the aberration's PSD is not perfectly accurate. As previously, *k* stands for *u* (upstream) or *d* (downstream).</caption><thead><tr><td>Aberration's PSD</td><td>ε<sub>u</sub></td><td>ε<sub>d</sub></td></tr></thead><tbody><tr><td>S<sub>φk</sub> ∝ 1/v</td><td>3.08 nm RMS</td><td>5.85 nm RMS</td></tr><tr><td>S<sub>φk</sub> ∝ 1/v<sup>2</sup></td><td>2.28 nm RMS</td><td>4.09 nm RMS</td></tr><tr><td>S<sub>φk</sub> ∝ 1/v<sup>3</sup></td><td>1.54 nm RMS</td><td>3.40 nm RMS</td></tr></tbody></table>
Table 2 shows the results of this evaluation. One can see that COFFEE is not very sensitive to the a priori accuracy, since the reconstruction errors remain indeed small when the PSD assumed by COFFEE does not perfectly match the actual PSD.

When $S_{\phi_k}$ follows a $1/v$ power law, we note a slight increase in the reconstruction error. We have checked that, as expected, this increase is due to very high frequencies in the phase which are over-regularized and thus not reconstructed. In the conditions of our simulations, these very high frequencies are actually beyond the spatial frequencies that the DM is able to control, so that in closed loop this HF error on the reconstructed phase would have no impact.

We note that when $S_{\phi_k}$ follow a $1/v^3$ power law, reconstruction errors are lower than when the correct $1/v^2$ power law is considered. Again, the reconstruction errors originate mostly in the estimation of the high frequency aberrations, which give birth to lower energy speckle than the low frequency aberrations. When the PSD follows a $1/v^3$ power law, the quantity of high frequencies decreases in the aberrations to estimate, leading to an improved reconstruction error.

### 4. Circumventing calibration errors of the diversity phase: the myopic approach

The value of the diversity phase $\phi_{\text{div}}$ is one of the few inputs COFFEE needs to perform the phase estimation. Thus, an imprecise calibration of $\phi_{\text{div}}$ will lead to an error on the estimated aberration. In the case of classical phase diversity, this error is the one that drives the total error budget [24]. In [15], we demonstrated that an error $\epsilon_{\text{div}}$ on the knowledge of $\phi_{\text{div}}$ was leading to a reconstruction error of about $\epsilon_{\text{div}}/2$ on both upstream and downstream estimated aberrations. Now that both aliasing and model errors have been tackled in COFFEE, the diversity calibration error would be the most important one in the error budget.

The most convenient way to introduce the diversity phase $\phi_{\text{div}}$ on the instrument is to modify the reference slopes of the AO loop to introduce a calibrated aberration, as described in [15]. The accuracy of such a process is thus limited by the DM's ability to achieve a given shape, leading to an error $\epsilon_{\text{div}}^{\text{DM}}$ on the phase diversity $\phi_{\text{div}}^{\text{DM}}$ actually introduced. As this error will always be present on an AO system (thermal evolution of the DM, ageing of influence functions, inability to re-calibrate them regularly), we have adapted COFFEE to make it able to perform a joint estimation of an error on $\phi_{\text{div}}$. This approach, called hereafter "myopic estimation", consist in a slight modification of the criterion $J$ to be minimized (Eq. (3)), in which an additional unknown parameter $\phi_{\epsilon}$ (called hereafter diversity error phase) is introduced:

$$ J(\alpha, \beta, \phi_u, \phi_d, \phi_{\epsilon}) = \frac{1}{2} \left\| \frac{i_c^{\text{foc}} - (\alpha_{\text{foc}}h_{\text{det}} * h_c(\phi_u, \phi_d) + \beta_{\text{foc}})}{\sigma_n^{\text{foc}}} \right\|^2 + \frac{1}{2} \left\| \frac{i_c^{\text{div}} - (\alpha_{\text{div}}h_{\text{det}} * h_c(\phi_u + \phi_{\text{div}}^{\text{cal}} + \phi_{\epsilon}, \phi_d) + \beta_{\text{div}})}{\sigma_n^{\text{div}}} \right\|^2 + \mathcal{R}(\phi_u) + \mathcal{R}(\phi_d) + \mathcal{R}(\phi_{\epsilon}), \quad (13) $$

where $\phi_{\text{div}}^{\text{cal}}$ is the calibrated diversity phase: $\phi_{\text{div}} = \phi_{\text{div}}^{\text{cal}} + \phi_{\epsilon}$. $\mathcal{R}(\phi_{\epsilon})$ is an optional regularization metric designed to enforce our knowledge that $\phi_{\epsilon}$ should be small and smooth. Using the gradient $\frac{\partial J}{\partial \phi_{\epsilon}}$ (whose computation is trivial knowing $\frac{\partial J}{\partial \phi_u}$), COFFEE is able to perform a joint estimation of $\phi_{\epsilon}$ along with the previously estimated parameters $\alpha, \beta, \phi_u$ and $\phi_d$.

In this paper, the estimated phases $\phi_u$ and $\phi_d$ are expanded on a pixel basis, which allows the estimation of high-order aberrations. However, since the diversity phase is composed of low order aberration, one can expect that the error $\phi_{\epsilon}$ will be mainly composed of low order aberrations. Thus, in order to optimize COFFEE's performance, the estimation of $\phi_{\epsilon}$ has been implemented in two different ways:
* $\phi_{\epsilon}$ can be expanded on a pixel basis if high-order patterns (such as dead actuators) are expected in the diversity error phase. In this case, the regularization metric $\mathcal{R}(\phi_{\epsilon})$ used will have the same expression as $\mathcal{R}(\phi_{u})$, presented in Section 2.2.

* If only low order aberrations are expected in $\phi_{\epsilon}$, it can be expanded on a truncated Zernike basis, composed of a few Zernike modes (typically up to 5 if $\phi_{\text{div}}$ is composed of defocus ($Z_4$) and astigmatism ($Z_5$)). With this basis, which reduces the number of parameters to be estimated, no regularization metric dedicated to $\phi_{\epsilon}$ is needed.

We validate this myopic approach by a realistic simulation: using parameters gathered in Table 1, we simulate coronagraphic images, and consider that the diversity phase used for the simulation is not perfectly known in the estimation stage. The coronagraphic simulated diversity image is computed with a diversity $\phi_{\text{div}} = \phi_{\text{div}}^{\text{cal}} + \phi_{\epsilon}$, where $\phi_{\epsilon} = a_{\epsilon}(Z_4 + Z_5)$, considering that the amplitude of the phase diversity is not perfectly known. In this simulation, $a_{\epsilon} = 0.04$ rad RMS (10 nm RMS).

COFFEE's phase reconstruction is then performed using both simulated images considering that the diversity phase is equal to $\phi_{\text{div}}^{\text{cal}}$, first without the myopic estimation, then with an estimation of $\phi_{\epsilon}$ expanded on pixel basis and on a truncated Zernike basis composed of 4 modes: tip, tilt (which allows an estimation of a differential tip-tilt between the two images), defocus and astigmatism.

<table><caption>Table 3. Comparison of COFFEE estimation accuracy with and without myopic estimation of $\phi_{\epsilon}$ when the diversity phase is not perfectly known.</caption><thead><tr><th></th><th>ε<sub>u</sub></th><th>ε<sub>d</sub></th></tr></thead><tbody><tr><td>φ<sub>ε</sub>: no estimation</td><td>6.92 nm RMS</td><td>8.51 nm RMS</td></tr><tr><td>φ<sub>ε</sub>: estimation on a pixel-wise map</td><td>2.55 nm RMS</td><td>5.45 nm RMS</td></tr><tr><td>φ<sub>ε</sub>: estimation on 4 Zernike modes</td><td>2.64 nm RMS</td><td>4.16 nm RMS</td></tr></tbody></table>

The reconstruction errors corresponding to each reconstruction are gathered in Table 3. One can clearly see here that the myopic estimation approach significantly improves the reconstruction accuracy, roughly by a factor 2, allowing an optimal use of COFFEE even when the diversity phase $\phi_{\text{div}}$ is not perfectly known.

The usefulness of this myopic approach will be further illustrated in Section 5.2, where we show that the estimation of an error on the diversity phase allows us to improve the contrast on the detector plane in a compensation process.

## 5. Closed loop quasi-static aberration compensation using COFFEE

In this Section, we present COFFEE's performance in a closed loop compensation process. We consider here the case of the calibration of a SPHERE-like instrument. Two coronagraphic images are simulated using randomly generated aberrations, whose PSD follows a $1/v^2$ law. Then, using these two simulated coronagraphic images, the aberrations upstream $\phi_u$ and $\phi_d$ are estimated using COFFEE.

Once the phase reconstruction is performed, we modify the DM actuator voltages to compensate for the estimated aberrations and thus optimize the contrast in the detector plane. This compensation is performed using two different techniques: in Section 5.1, we use the conventional phase conjugation. In Section 5.2, we minimize the energy in a chosen area in order to optimize the contrast in the selected region of the detector plane. To perform such a compensation, we have developed a method dedicated to high-contrast imaging instruments that does
not rely on any small aberration approximation. These simulations are performed using the parameters gathered in Tables 1 and 4.

<table><caption>Table 4. COFFEE: parameters used for the compensation simulations of Section 5</caption><tr><td>Lyot stop pupil</td><td>$D_d = 0.96D_u$</td></tr><tr><td>Incoming flux</td><td>$10^9$ photons</td></tr><tr><td>Detector noise</td><td>$\sigma_{det} = 1 e^-$</td></tr><tr><td>Deformable mirror (DM)</td><td>$41 \times 41$ actuators, Gaussian-shaped influence functions</td></tr></table>

In order to perform a realistic simulation, we consider here that an error $\phi_{\epsilon}$ is made on the diversity phase: the coronagraphic diversity image is computed with a diversity phase $\phi_{\text{div}}^{\text{sim}} = \phi_{\text{div}} + \phi_{\epsilon}$. In this section, we consider that $\phi_{\epsilon} = a_{\epsilon}(Z_4 + Z_5)$, with $a_{\epsilon} = 0.04$ rad RMS (10 nm RMS).

COFFEE estimates the aberrations considering that the calibrated diversity phase is equal to $\phi_{\text{div}}$, and jointly searches for the diversity phase error $\phi_{\epsilon}$, as described in Section 4.

## 5.1. Phase conjugation

Conventional phase conjugation aims at compensating for the aberrations upstream of the coronagraph $\phi_u$ in an iterative process. After its criterion minimization, COFFEE gives an estimation $\hat{\phi}_u$. The aberrations upstream of the coronagraph at the iteration $i+1$ are thus given by:

$$ \phi_u^{i+1} = \phi_u^i - g\hat{\phi}_u^{\text{DM}}, \quad (14) $$

where $g$ is the gain of the iterative process ($g = 0.5$ in this simulation) and $\hat{\phi}_u^{\text{DM}}$ is the aberration introduced by the DM in the entrance pupil plane. Such an aberration corresponds to the best representation of $\hat{\phi}_u$ achievable by the deformable mirror. Let $F$ be the DM's influence matrix and $T$ its generalized inverse. The aberration introduced by the DM can be computed as follows:

$$ \hat{\phi}_u^{\text{DM}} = F T \hat{\phi}_u. \quad (15) $$

Figure 5 presents the result of the simulation of a compensation performed by phase conjugation after 5 iterations of the loop (the average computation time for one loop iteration is 2 minutes). On this figure, one can see that the aberration compensation performed using COFFEE estimation allows a significant improvement in the coronagraphic images (figure 5, left). In particular, the average contrast plot (figure 5, right) shows that after compensation, the performance in the area controlled by the DM ($\pm 20\lambda/D$) is very close to the one that would be obtained from a coronagraphic image computed without aberrations.

We note that the level of contrast reached after compensation presented in Fig. 5 exceeds the SPHERE instrument specification for off-line calibration, which relies on a phase estimation based on classical phase diversity (no coronagraph). In particular, classical phase diversity is not able to estimate high-order frequencies, and thus will not compensate for speckles located beyond $8\lambda/D$. The use of COFFEE, which allows a compensation in the whole area controlled by the DM in the detector plane (as shown in Fig. 5), could thus improve the SPHERE instrument performance.

## 5.2. Creation of a Dark Hole on the detector

Speckle nulling iterative techniques aim at minimizing the energy in a chosen area of the detector in order to facilitate exoplanet detection in this area, called a “Dark Hole” (DH). To
<figure><img src="image_7.png" /><figcaption>Fig. 5. COFFEE: simulation of aberration compensation using classical phase conjugation. Left: coronagraphic images before (top) and after five iterations of the closed loop process (bottom, logarithmic scale for both images). Right: normalized average raw contrast computed from the images without coronagraph (solid red line), with coronagraph before (solid green line) and after compensation (solid blue line). For comparison, the contrast computed from a coronagraphic image computed without any aberrations is plotted as well (dashed black line).</figcaption></figure>

create this DH, the methods developed until now, which rely on a small aberration approximation, minimize the energy during an iterative process [8, 9]. This process requires several iterations and is based on the knowledge of an interaction matrix between the detector plane and the DM pupil plane to create the Dark Hole on the detector.

Here, we propose a new method to minimize the energy that does not rely on any small aberration assumption, allowing us to deal with high amplitude phase aberrations such as the ones created by the DM dead actuators (which are currently part of the SPHERE instrument's limitations). Besides, this method takes into account both upstream and downstream aberrations. This compensation method, coupled with COFFEE would be particularly adapted to the calibration of a high-contrast imaging instrument.

Let us define the energy in the focal plane $\mathcal{E}_{DH}$ in the DH as:

$$ \mathcal{E}_{DH} = \alpha \sum_{m,n \in \text{DH}} |\Psi_{DH}(m,n)|^2 \quad (16) $$

With $m, n$ the pixel position in the DH (for the sake of simplicity, these variables will be omitted in the following). $\alpha$ is the incoming flux and $\Psi_{DH}$ the electric field in the DH area which, using notations of Section 2.3, is given by:

$$ \Psi_{DH}(v) = \mathcal{F}^{-1} \left\{ \mathcal{F}^{-1} \left( P_u e^{j(\phi_u + \psi(v))} \right) \mathcal{M} \right\} P_d e^{j\phi_d} \quad (17) $$

$\psi$ represents the aberration introduced by the DM: $\psi(v) = Fv$, with $v$ the set of voltages applied to the DM actuators. Thus, creating a DH on the detector means finding the set of voltages $v_{DH}$ that minimize the energy $\mathcal{E}_{DH}(v)$, knowing the relation between the entrance pupil plane and the detector plane of the high contrast imaging system.

Here, we minimize $\mathcal{E}_{DH}(v)$ by means of the same optimization method as the one which was used by COFFEE to perform its estimation, which is a limited memory variable metric (BFGS) method [16, 17]. Such an operation gives us the voltage $v_{DH}$, and thus the aberration $\psi(v_{DH})$ to
introduce using the DM to create the DH. Unlike others energy minimization methods, the one we propose here does not require an interaction matrix calibration between the detector plane and the DM pupil plane. In particular, using COFFEE, we are able to estimate the aberrations downstream of the coronagraph, and among them the coronagraphic PSF position on the detector, i.e. downstream tip-tilt. This estimation is then taken into account in our compensation method (Eqs (16) and (17)). Thus, variation in the aberrations downstream of the coronagraph will not require any particular action during the compensation process, since these variations will be estimated by COFFEE along with the aberrations upstream of the coronagraph.

As mentioned previously, COFFEE does not estimate amplitude aberrations. However, the compensation method described in this section can easily be adapted to take into account amplitude aberrations upstream of the coronagraph $\xi$, by modifying Eq. (17):

$$ \Psi_{\text{DH}}(v) = \mathcal{F}^{-1} \left\{ \mathcal{F} \left[ \mathcal{F}^{-1} \left( P_{it} e^{j(\phi_u + \psi(v)) + \xi} \right) \mathcal{M} \right] P_d e^{j\phi_d} \right\} \quad (18) $$

With $\xi$ the amplitude aberrations. Thus, provided these aberrations are known, amplitude aberrations does not limit the performance of the compensation method presented in this Section.

We now validate this new energy minimization method by simulation. It requires the knowledge of the aberrations upstream $\phi_u$ and downstream $\phi_d$ of the coronagraph, and of the incoming flux $\alpha$. In this section, we use COFFEE to perform the estimation of these parameters, and then minimize the energy in order to create a DH in the detector plane. The energy $\mathcal{E}_{\text{DH}}$ is minimized between $5\frac{\lambda}{D}$ and $20\frac{\lambda}{D}$ to create the Dark Hole in the right part of the focal plane. As in the phase conjugation case (Section 5.1), such a compensation is performed in an iterative process, where the aberrations upstream of the coronagraph at the iteration $i+1$ are given by:

$$ \phi_{it}^{i+1} = \phi_{it}^{i} + gF_{\text{vDH}}. \quad (19) $$

As previously, we consider $g = 0, 5$ in this simulation.

The result of this compensation after 5 iterations of the loop is presented in Fig. 6 (the average computation time for one loop iteration is 3 minutes). In the targeted area (between $5\frac{\lambda}{D}$ and $20\frac{\lambda}{D}$), this method allows a significant improvement: when the compensation is performed by conventional phase conjugation (Figure 5), the average contrast in the same area is $7.2 \times 10^{-7}$. The use of our new compensation method allow a contrast improvement by a factor 10 after compensation (Figure 6, bottom left). Besides, the interest of the myopic approach is illustrated in Fig. 6: when the closed loop process is performed without estimation of an error on the diversity phase, the performance decreases by a factor 2 (solid cyan line).

If a very high contrast level is required, this compensation technique can thus be used instead of conventional phase conjugation to calibrate the instrument. Besides, in order to push down the dark hole floor, it is possible to narrow the energy minimization area, as mentioned by Bordé et al. [8].

## 6. Conclusion

In this paper, an extended version of our coronagraphic phase diversity, nicknamed COFFEE, has been presented. The use of a regularized pixel basis in the estimation allows COFFEE to estimate high order aberrations with nanometric precision (Section 2). Besides, thanks to a modification of the coronagraphic PSF used in the imaging model, COFFEE is no longer limited to a particular coronagraphic device. COFFEE's performance has been studied and discussed in Section 3, while Section 4 has described a so-called myopic extension of COFFEE, which consists in a joint estimation of an error on the diversity phase in order
<figure><img src="image_8.png" /><figcaption>Fig. 6. COFFEE: simulation of aberration compensation by minimizing the energy in a chosen area (between $5/2$ and $20/2$ in the right part of the focal plane). Left: coronagraphic images before (top) and after five iterations of the closed loop process (bottom, logarithmic scale for both images). Right: normalized average raw contrast computed in the dark hole area from the images without coronagraph (solid red line), with coronagraph before (solid green line) and after compensation. After minimization, the average contrast in the Dark Hole is $3.1 \ 10^{-8}$ when the myopic approach is used (solid blue line), and $6.6 \ 10^{-8}$ otherwise (solid cyan line). For comparison, the contrast computed from a coronagraphic image computed without any aberrations is plotted as well (dashed black line).</figcaption></figure>

to improve COFFEE's accuracy in a real system, where the diversity phase is not perfectly known. Lastly, in Section 5, the achievable contrast optimization on a SPHERE-like system using COFFEE in a compensation process has been studied using realistic simulations. In the latter section, we have presented a new compensation method which minimizes the energy in a chosen area of the detector through a non-linear minimization, in order to reach higher level of contrast than those that can be obtained using phase conjugation.

The experimental validation of this high-order and myopic version of COFFEE is ongoing, and aims at demonstrating the ability of COFFEE to estimate both low and high order aberrations, and to compensate for them. Several perspectives are currently considered for this work. With an adaptation of the coronagraphic imaging model, COFFEE can be extended to work on ground-based, long exposure images with residual turbulence induced aberrations. Another perspective lies in optimization of the computation time required for the aberration estimation, possibly following I. Mocœur et al. [25]. These two improvements will allow COFFEE to work on-line, in closed loop during the scientific exposure. A further perspective is to extend COFFEE to the estimation of amplitude aberrations, mandatory to reach the very high levels of contrast required for exo-earth imaging, in order to create a dark hole area on the detector using the method presented in this paper.

## Acknowledgments

The authors would like to thank several key players of the SPHERE instrument, namely Thierry Fusco, David Mouillet, Jean-Luc Beuzit and Marc Ferrari for simulating discussions and support, and the Région Provence-Alpes-Côte d'Azur for partial financial support of B. Paul's scholarship. This work was partly funded by the European Commission under FP7 Grant
A. Gradients expression

The numerical minimization of criterion J (Eq. (3)) requires the analytic expression of gradients $\frac{\partial J}{\partial \phi_u}$, $\frac{\partial J}{\partial \phi_d}$, $\frac{\partial J}{\partial \alpha}$ and $\frac{\partial J}{\partial \beta}$ to estimate the aberrations upstream $\phi_u$ and downstream $\phi_d$ of the coronagraph, as well as the incoming flux $\alpha$ and the residual background $\beta$. Let us rewrite here the expression of criterion J:

$$ J(\alpha, \beta, \phi_u, \phi_d) = \frac{1}{2} \left\| \frac{i_c^{\text{foc}} - (\alpha_{\text{foc}}h_{\text{det}} \star h_c^{\text{foc}} + \beta_{\text{foc}})}{\sigma_n^{\text{foc}}} \right\|^2 + \frac{1}{2} \left\| \frac{i_c^{\text{div}} - (\alpha_{\text{div}}h_{\text{det}} \star h_c^{\text{div}} + \beta_{\text{div}})}{\sigma_n^{\text{div}}} \right\|^2 \\ + \mathcal{R}(\phi_u) + \mathcal{R}(\phi_d) \\ = J^{\text{foc}} + J^{\text{div}} + \mathcal{R}(\phi_u) + \mathcal{R}(\phi_d) \tag{20} $$

The expressions of $\frac{\partial J}{\partial \alpha}$ and $\frac{\partial J}{\partial \beta}$ can be found in [15]. The calculation of gradients $\frac{\partial J}{\partial \phi_u}$ and $\frac{\partial J}{\partial \phi_d}$ is performed following what have been done in [15]: we derive $J^{\text{foc}}$, and then deduce the gradients' expressions of $J^{\text{div}}$ using a trivial substitution. The notations used here are the ones introduced in Section 2:

$$ \frac{\partial J^{\text{foc}}}{\partial \phi_d} = 2\mathfrak{S} \left\{ \psi_0^* - \varepsilon \psi_d \mathcal{F} [\mathcal{M} \mathcal{F}^{-1} (\psi_u)]^* \times \mathcal{F} \left[ \frac{\partial J^{\text{foc}}}{\partial h_c^{\text{foc}}} (\Psi_0 - \varepsilon \Psi_c) \right] \right\} \\ \frac{\partial J^{\text{foc}}}{\partial \phi_u} = 2\mathfrak{S} \left\{ \psi_0^* \mathcal{F} \left[ \frac{\partial J^{\text{foc}}}{\partial h_c^{\text{foc}}} (\Psi_0 - \varepsilon \Psi_c) \right] \right\} \\ - \varepsilon \psi_u^* \mathcal{F} \left[ \mathcal{M}^* \mathcal{F}^{-1} \left( \Psi_d^* \mathcal{F} \left\{ \frac{\partial J^{\text{foc}}}{\partial h_c^{\text{foc}}} [\Psi_0 - \varepsilon \Psi_c] \right\} \right) \right] $$

with:

$$ \frac{\partial J^{\text{foc}}}{\partial h_c^{\text{foc}}} = \frac{1}{\sigma_n^{\text{foc}} 2} [\alpha h_{\text{det}} (\alpha h_{\text{det}} \star h_c^{\text{foc}} - i_c^{\text{foc}})] \tag{22} $$

and:

$$ \begin{align} \psi_u &= P_u e^{j\phi_u} \notag \\ \psi_d &= P_d e^{j\phi_d} \qquad \Psi_d = \mathcal{F}^{-1}(\psi_d) \tag{23} \\ \psi_0 &= P_u e^{j(\phi_u + \phi_d)} \qquad \Psi_0 = \mathcal{F}^{-1}(\psi_0) \notag \\ \Psi_c &= \mathcal{F}^{-1}\{\psi_d \mathcal{F}[\mathcal{M} \mathcal{F}^{-1}(\psi_u)]\} \notag \end{align} $$

The regularization metric expression $\mathcal{R}(\phi_k)$ (k is for u (upstream) or d (downstream)) is given by Eq. (6). Its gradient $\frac{\partial \mathcal{R}}{\partial \phi_k}$ can be written as:

$$ \frac{\partial \mathcal{R}}{\partial \phi_k} = \mu_k \| \Delta \phi_k(r) \|. \tag{24} $$

where $\Delta$ represent the Laplacian operator.