An analytic expression for coronagraphic imaging through turbulence. Application to on-sky coronagraphic phase diversity

TL;DR

This paper derives an analytic expression for coronagraphic imaging through turbulence, enabling efficient wave-front reconstruction and quasi-static aberration calibration in high contrast exoplanet imaging systems.

Contribution

It extends Roddier's imaging expression to coronagraphic systems, providing a validated, computationally efficient model for turbulence-affected imaging.

Findings

Validated the model with numerical simulations.

Achieved sub-nanometric wave-front reconstruction precision.

Demonstrated potential for real-time aberration calibration during observations.

Abstract

The ultimate performance of coronagraphic high contrast exoplanet imaging systems such as SPHERE or GPI is limited by quasi-static aberrations. These aberrations produce speckles that can be mistaken for planets in the image. In order to design instruments, correct quasi-static aberrations or analyse data, the expression of the point spread function of a coronagraphic instrument in the presence of residual turbulence is most useful. Here we derive an analytic expression for this point spread function that is an extension to coronagraphic imaging of Roddier's expression for imaging through turbulence. We give a physical interpretation of its structure, we validate it by numerical simulations and we show that it is computationally efficient. Finally, we incorporate this imaging model into a coronagraphic phase diversity method (COFFEE) and validate by simulations that it allows wave-front reconstruction in the presence of residual turbulence. The preliminary results, which give a sub-nanometric precision in the case of a SPHERE-like system, strongly suggest that quasi-static aberrations could be calibrated during observations by this method.

Figures

Equations

$$h_{le}=h_a\star h_s,$$

$${\tilde{h}}_a=\exp \left(-\frac{1}{2}{D}_{\varphi }\right),$$

$$h_c\left(\alpha ;{\psi }_u,{\psi }_d\right)={\left|{\mathcal{F}}^{-1}\left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u\right)\right]\right\}(\alpha )\right|}^2,$$

$$h_{sec}(\alpha ,t;{\psi }_u,{\psi }_d)=h_c(\alpha ;{\psi }_u\times {\mathrm{e}}^{\mathrm{i}{\varphi }_t(t)},{\psi }_d).$$

$$h_{lec}(\alpha ;{\psi }_u,{\psi }_d)={⟨h_{sec}(\alpha ,t;{\psi }_u,{\psi }_d)⟩}_t.$$

$$f\otimes g(r)=\iint f^{\ast }(s)g(r+s)\mathrm{d}s.$$

$$\begin{array}{rl}{h}_{lec}(\alpha ;\psi & {}_u,{\psi }_d)=\\ & ⟨{\mathcal{F}}^{-1}(\left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}{\varphi }_t(t)}\right)\right]\right\}\\ & \otimes {\left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}{\varphi }_t(t)}\right)\right]\right\})(\alpha )⟩}_t.\end{array}$$

$$\begin{array}{rl}{\tilde{h}}_{lec}(r;\psi & {}_u,{\psi }_d)=\\ & ⟨\left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}{\varphi }_t(t)}\right)\right]\right\}\\ & {\otimes \left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}{\varphi }_t(t)}\right)\right]\right\}(r)⟩}_t.\end{array}$$

$$\begin{array}{rl}{\tilde{h}}_{lec}(r;{\psi }_u,{\psi }_d)=& \iint {\psi }_d^{\ast }(r_1){\psi }_d(r+r_1)\iint \iint {\mathrm{e}}^{\mathrm{i}2\pi r_1\cdot {\alpha }_1}\\ & \times {\mathrm{e}}^{-\mathrm{i}2\pi (r_1+r)\cdot {\alpha }_2}{\mathcal{M}}^{\ast }({\alpha }_1)\mathcal{M}({\alpha }_2)\\ & \times \iint \iint {\mathrm{e}}^{-\mathrm{i}2\pi r_2\cdot {\alpha }_1}{\mathrm{e}}^{\mathrm{i}2\pi r_3\cdot {\alpha }_2}{\psi }_u^{\ast }(r_2){\psi }_u(r_3)\\ & \times {⟨{\mathrm{e}}^{\mathrm{i}\left[{\varphi }_t(r_3,t)-{\varphi }_t(r_2,t)\right]}⟩}_t\mathrm{d}{r}_3\mathrm{d}{r}_2\mathrm{d}{\alpha }_2\mathrm{d}{\alpha }_1\mathrm{d}{r}_1.\end{array}$$

$${⟨{\mathrm{e}}^{\mathrm{i}\left[{\varphi }_t(r_3,t)-{\varphi }_t(r_2,t)\right]}⟩}_t={\mathrm{e}}^{-\frac{1}{2}{D}_{\varphi }(r_3-r_2)},$$

$$D_{\varphi }(r)={⟨{\mathrm{e}}^{\mathrm{i}\left[{\varphi }_t(r^{\prime },t)-{\varphi }_t(r+r^{\prime },t)\right]}⟩}_t.$$

$$\begin{array}{rl}{\tilde{h}}_{lec}(r;\psi & {}_u,{\psi }_d,D_{\varphi })=\\ & \iint h_a({\alpha }^{\prime };D_{\varphi })\iint {\psi }_d^{\ast }(r_1)\iint {\mathrm{e}}^{\mathrm{i}2\pi r_1\cdot {\alpha }_1}{\mathcal{M}}^{\ast }({\alpha }_1)\\ & \times \iint {\mathrm{e}}^{-\mathrm{i}2\pi r_2\cdot {\alpha }_1}{\psi }_u^{\ast }(r_2){\mathrm{e}}^{-\mathrm{i}2\pi r_2\cdot {\alpha }^{\prime }}\mathrm{d}{r}_2\mathrm{d}{\alpha }_1\\ & \times {\psi }_d(r+r_1)\iint {\mathrm{e}}^{-\mathrm{i}2\pi (r_1+r)\cdot {\alpha }_2}\mathcal{M}({\alpha }_2)\\ & \times \iint {\mathrm{e}}^{\mathrm{i}2\pi r_3\cdot {\alpha }_2}{\psi }_u(r_3){\mathrm{e}}^{\mathrm{i}2\pi r_3\cdot {\alpha }^{\prime }}\mathrm{d}{r}_3\mathrm{d}{\alpha }_2\mathrm{d}{r}_1\mathrm{d}{\alpha }^{\prime }.\end{array}$$

$$\begin{array}{rl}{\tilde{h}}_{lec}(r;\psi & {}_u,{\psi }_d,D_{\varphi })=\\ & \iint h_a({\alpha }^{\prime };D_{\varphi })\left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}2\pi {\alpha }^{\prime }\cdot \mathrm{Id})}\right)\right]\right\}\\ & \otimes \left\{{\psi }_d\times \mathcal{F}\left[\mathcal{M}\times {\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}2\pi {\alpha }^{\prime }\cdot \mathrm{Id})}\right)\right]\right\}(r)\mathrm{d}{\alpha }^{\prime }.\end{array}$$

$$\begin{array}{rl}h& {}_{lec}(\alpha ;{\psi }_u,{\psi }_d,D_{\varphi })=\\ & \iint h_a({\alpha }^{\prime };D_{\varphi }){\left|{\mathcal{F}}^{-1}\left\{{\psi }_d.\mathcal{F}\left[\mathcal{M}.{\mathcal{F}}^{-1}\left({\psi }_u{\mathrm{e}}^{\mathrm{i}2\pi {\alpha }^{\prime }\cdot \mathrm{Id}}\right)\right]\right\}(\alpha )\right|}^2\mathrm{d}{\alpha }^{\prime }.\end{array}$$

$$\boxed{\begin{array}{rl}{h}_{lec}& (\alpha ;{\psi }_u,{\psi }_d,D_{\varphi })=\\ & \iint h_a({\alpha }^{\prime };D_{\varphi })h_c\left(\alpha ;{\psi }_u{\mathrm{e}}^{\mathrm{i}2\pi {\alpha }^{\prime }\cdot \mathrm{Id}},{\psi }_d\right)\mathrm{d}{\alpha }^{\prime },\end{array}}$$

$$h_a\triangleq {\mathcal{F}}^{-1}\left[\exp \left(-\frac{1}{2}{D}_{\varphi }\right)\right]={\mathcal{F}}^{-1}\left[{⟨{\mathrm{e}}^{\mathrm{i}\left[{\varphi }_t(r+r^{\prime },t)-{\varphi }_t(r^{\prime },t)\right]}⟩}_t\right].$$

$$h_a={\mathcal{F}}^{-1}\left[{⟨{\psi }_t^{\ast }(r^{\prime },t){\psi }_t(r+r^{\prime },t)⟩}_t\right]$$

$$h_{lec}(\alpha ;{\psi }_u,{\psi }_d)=h_c(\alpha ;{\psi }_u,{\psi }_d),$$

$$\begin{array}{rl}{h}_{lec}(\alpha ;\psi & {}_u,{\psi }_d,D_{\varphi })\\ & =\iint h_a({\alpha }^{\prime };D_{\varphi }){\left|{\mathcal{F}}^{-1}\left\{{\psi }_d{\psi }_u{\mathrm{e}}^{\mathrm{i}2\pi {\alpha }^{\prime }\cdot \mathrm{Id}}\right\}(\alpha )\right|}^2\mathrm{d}{\alpha }^{\prime }\\ & =h_a(D_{\varphi })\star {\left|{\mathcal{F}}^{-1}\left\{{\psi }_d{\psi }_u\right\}\right|}^2(\alpha ),\end{array}$$

$${\tilde{h}}_{lec}(r;{\psi }_u,{\psi }_d)={\tilde{h}}_a(r)\times \left({\psi }_d{\psi }_u\otimes {\psi }_d{\psi }_u\right)(r),$$

$$\frac{1}{2}{D}_{\varphi }(r)={\sigma }_{\varphi }^2-A_{\varphi }(r),$$

$$h_a={\mathrm{e}}^{-{\sigma }_{\varphi }^2}{\mathcal{F}}^{-1}\left[{\mathrm{e}}^{A_{\varphi }}\right].$$

$$\begin{array}{rl}{h}_a({\alpha }^{\prime })& \simeq {\mathrm{e}}^{-{\sigma }_{\varphi }^2}{\mathcal{F}}^{-1}\left[1+A_{\varphi }\right]({\alpha }^{\prime })\\ h_a({\alpha }^{\prime })& \simeq {\mathrm{e}}^{-{\sigma }_{\varphi }^2}\delta ({\alpha }^{\prime })+{\mathrm{e}}^{-{\sigma }_{\varphi }^2}{S}_{\varphi }({\alpha }^{\prime })\end{array}$$

$$\begin{array}{rl}h& {}_{lec}(\alpha ;{\psi }_u,{\psi }_d)\simeq \\ & {\mathrm{e}}^{-{\sigma }_{\varphi }^2}\times \left[h_c(\alpha ;{\psi }_u,{\psi }_d)+\iint S_{\varphi }({\alpha }^{\prime })h_c\left(\alpha ;{\psi }_u{e}^{\mathrm{i}2\pi {\alpha }^{\prime }},{\psi }_d\right)\mathrm{d}{\alpha }^{\prime }\right]\end{array}$$

Full text

††pubyear: 2017††pagerange: An analytic expression for coronagraphic imaging through turbulence. Application to on-sky coronagraphic phase diversity.–References

1 Introduction

Direct detection of exoplanets is limited by the presence of speckles on scientific images. These speckles are induced by quasi-static aberrations in the optical system that are not corrected by the adaptive optics loop. On SPHERE on the VLT, the quasi-static aberrations can be calibrated on an internal source thanks to COFFEE (see Paul et al. (2014)), an extension of phase diversity (see Gonsalves (1982) for the principle of phase diversity, and Mugnier et al. (2006) for a review) to coronagraphic imaging. However, the quasi-static aberrations evolve during the observations, so it would be useful to calibrate them during the observations in order to reach the ultimate performance of the instrument. Since COFFEE uses a model of the instrument, an analytic expression of the long exposure point spread function of a coronagraphic instrument in the presence of the residual turbulence of an adaptive optics system is needed.

To the best of our knowledge, no such general analytic expression has been previously published, even though the shape and properties of the coronagraphic point spread function have been explored, e.g. in Perrin et al. (2003); Soummer et al. (2007b), and an analytic approach has previously been developed in Sauvage et al. (2010). However, Sauvage et al. (2010) was based on the hypothesis of a perfect coronagraph, which does not exist, and yields unacceptable errors when used for coronagraphic phase diversity, as shown by Paul et al. (2013b).

The purpose of this letter is firstly to present a general analytic expression for coronagraphic imaging through turbulence – which is an extension to the coronagraphic case of the well-known expression by Roddier (1981) for imaging through turbulence – and secondly to show that it can be integrated into COFFEE to measure quasi-static aberrations with nanometric precision during observations.

Moreover, the expression that we derive in this letter could be of great use for at least two other applications than wave-front sensing. One is the design, simulation and optimisation of coronagraphic systems. The other is the processing of high contrast images based on an imaging model, as proposed by Ygouf et al. (2013), which is essential for the detection and characterisation of exoplanets.

2 Classical imaging through turbulence and instantaneous coronagraphic point spread function

2.1 Classical imaging through turbulence

Let us start with a reminder of the non-coronagraphic case of imaging through turbulence. The Fourier transform of Equation (4.15) of the classical text by Roddier (1981) shows that the long exposure point spread function $h_{le}$ of a telescope through a (possibly residual) turbulent atmosphere is

[TABLE]

where $h_{s}$ is the point spread function of the optical system (which accounts for quasi-static aberrations), $\star$ is the convolution operator, and $h_{a}$ is the so-called atmospheric point spread function, which is the inverse Fourier transform of the atmospheric transfer function $\tilde{h}_{a}$:

[TABLE]

$D_{\phi}$ being the phase structure function of the (possibly residual) atmospheric turbulence.

This expression is appealing because it untangles the influence of the atmospheric turbulence from the influence of the quasi-static aberrations on the global point spread function.

We show below that an analogous formula can be derived for coronagraphic systems.

2.2 Instantaneous point spread function of a telescope equipped with a coronagraph

Let us consider an imaging system on a telescope equipped with a coronagraph as described in Figure 1. We will denote by $P_{u}$ the amplitude transmission in a pupil plane before the coronagraph, or “upstream pupil”, $\phi_{u}$ the phase aberrations in this plane, which will be called “upstream aberrations” hereafter, and $\psi_{u}=P_{u}\mathrm{e}^{\mathrm{i}\phi_{u}}$ the complex amplitude in this plane. Similarly, we will denote by $P_{d}$ the downstream transmission in a pupil plane after the coronagraph, or “downstream pupil” (also called the Lyot stop plane), $\phi_{d}$ the phase aberrations in this plane, which will be called “downstream aberrations” hereafter, and $\psi_{d}=P_{d}\mathrm{e}^{\mathrm{i}\phi_{d}}$ the multiplicative term – due to the downstream aberrations – that affects the complex amplitude of the wave that has propagated from the entrance aperture to this plane. For the sake of simplicity, we want to keep the same orientation of the images from one focal plane to the following, and the same pupil image orientation from one pupil plane to another. For this reason, we rotate the axis of the downstream plane by 180° with respect to the upstream plane, and similarly we rotate the orientation of the focal mask plane by 180° with respect to the final focal plane. This axis change is taken into account mathematically by choosing an inverse Fourier transform ($\mathcal{F}^{-1}$) instead of a forward one ($\mathcal{F}$) to model the propagation of light from a pupil to a focal plane. This convention allows for the image orientation to always remain the same. Of course, our main result (Equation (11)) is independent of this orientation convention.

Let us denote by $h_{c}$ the coronagraphic point spread function of the instrument without turbulence. It writes

[TABLE]

where $\alpha$ is the vector of angular coordinates of a position on the sky, and $\mathcal{M}$ the coronagraphic focal plane mask. This equation simply expresses the propagation from successive planes to the next ones and the quadratic detection, illustrated in Figure 1. For notation simplicity we take $\lambda=1$ in the equations, that is to say we use the same coordinate system for the pupil phases and for the transfer functions.

In the presence of turbulence, the instantaneous point spread function of this instrument writes

[TABLE]

Here $h_{sec}(\alpha,t,\psi_{u},\psi_{d})$ ($sec$ stands for “short exposure coronagraphic”) is the instantaneous point spread function taken at coordinate $\alpha$ in the detector focal plane at time $t$ and $\phi_{t}(t)$ is the instantaneous residual turbulence phase at time $t$ if the telescope is ground-based and equipped with adaptive optics.

3 Long exposure coronagraphic point spread function

In current high contrast systems such as SPHERE or GPI, the exposure time is large with respect to the typical evolution time of corrected turbulence, so that the signal that is actually recorded by the detector is the long exposure point spread function, which is an average over time of the short exposure point spread function:

[TABLE]

Here $h_{lec}$ is the long exposure coronagraphic point spread function, and $\left\langle\right\rangle_{t}$ is the averaging over time. Let us denote by $f^{*}$ the complex conjugate of any function $f$, and by $\otimes$ the bi-dimensional correlation product:

[TABLE]

In developed form, the long exposure coronagraphic point spread function writes, thanks to Wiener-Khintchine’s theorem:

[TABLE]

To somewhat simplify calculations, we will consider the optical transfer function instead of the point spread function. Let $\tilde{h}_{lec}~{}=~{}\mathcal{F}(h_{lec})$ be the Fourier transform of the point spread function. The integral over space commutes with the averaging over time, so

[TABLE]

In order to express the average over time, we must fully develop this expression. We introduce the pupil-plane variables $r_{1}$, $r_{2}$ and $r_{3}$, the focal-plane variables $\alpha_{1}$ and $\alpha_{2}$. After some manipulations, the fully developed expression reads

[TABLE]

Following Roddier (1981), we assume that turbulence is an ergodic stationary process (which is a very reasonable assumption in the case of residual turbulence after an extreme adaptive optics system), so we write

[TABLE]

where $D_{\phi}$ is the turbulent phase structure function, defined by

[TABLE]

In order to be able to separate variables $r_{2}$ and $r_{3}$ in Equation (8), we take the inverse Fourier transform of $\mathrm{e}^{-\frac{1}{2}D_{\phi}}$. We denote $\alpha^{\prime}$ the conjugate variable of $r_{3}-r_{2}$, and $h_{a}$ is the atmospheric point spread function defined by Equation (2). We will come back to the meaning of $h_{a}$ in Section 4. We obtain:

[TABLE]

In a more compact form, using $\mathrm{Id}$ to denote the identity function, this equation can be written:

[TABLE]

To obtain the point spread function back from this optical transfer function, we just have to take the inverse Fourier transform, then apply Wiener-Khintchine’s theorem again. Finally, the long exposure coronagraphic point spread function reads

[TABLE]

This expression gives the long exposure coronagraphic point spread function as a function of three deterministic parameters, namely: upstream aberrations, downstream aberrations, and residual turbulence-induced phase structure function.

4 Physical interpretation

In order to interpret Equation (10) physically, we can re-write it in the more compact following form:

[TABLE]

where $h_{c}\!\left(\alpha;\psi_{u}\mathrm{e}^{\mathrm{i}2\pi\alpha^{\prime}\cdot\mathrm{Id}},\psi_{d}\right)$ is the coronagraphic point spread function of the instrument in the absence of turbulence, but with a tilt $\alpha^{\prime}$ added to the upstream aberrations, that is to say with the light coming from the star being tilted by an angle $\alpha^{\prime}$.

A physical interpretation of the atmospheric point spread function defined in Eq. (2) can shed some light on the expression of Eq. (11). Recall that

[TABLE]

If we note $\psi_{t}(r,t)=\exp\left(\mathrm{i}\phi_{t}(r,t)\right)$ the contribution to the electric field at position $r$ and time $t$ by the atmospheric turbulence, we can re-write $h_{a}$ as :

[TABLE]

Assuming stationariness and ergodicity, we recognize the inverse Fourier transform of the auto-correlation of $\psi_{t}$. Hence we can identify $h_{a}$ as the energy spectral density of the turbulence-induced complex field. In other words, $h_{a}(\alpha^{\prime})$ gives the fraction of the energy of light that is diffracted in each direction $\alpha^{\prime}$.

Finally, Equation (11) can be interpreted as an incoherent plane wave decomposition: the long exposure coronagraphic point spread function $h_{lec}$ is a weighted sum of coronagraphic point spread functions $h_{c}$, without any turbulence, but with an upstream tilt $\alpha^{\prime}$. The weight on any of these tilted point spread functions is the fraction of the light energy $h_{a}(\alpha^{\prime})$ that the atmosphere scatters in the direction of the corresponding tilt $\alpha^{\prime}$.

Moreover, the expression of Equation (11) separates the turbulent part and the coronagraphic part of the problem: $h_{a}$ codes for the characteristics of the turbulent atmosphere, while $h_{c}$ codes for the characteristics of the instrument. This separation of the point spread function into an atmospheric point spread function and a point spread function of the instrument is strikingly similar to the one found in the Roddier expression for non-coronagraphic imaging.

5 Special cases

5.1 Non-turbulent point spread function

As a first particular case, let us assume that there is no turbulence, so that $h_{a}$ is reduced to a Dirac distribution. Indeed, in the case where there is no turbulence, $\phi_{t}$ is a constant, so $D_{\phi}=\left\langle\left[\phi_{t}(r^{\prime})-\phi_{t}(r+r^{\prime})\right]^{2}\right\rangle_{r^{\prime}}=0$. Then, according to Equation (2), $h_{a}=\delta.$ Hence, by use of Equation (11),

[TABLE]

which is precisely the classical expression for a coronagraphic point spread function in the absence of turbulence.

5.2 Non-coronagraphic optical transfer function

As a second particular case, if we now consider a non-coronagraphic instrument, $\mathcal{M}=1$, so Equation (10) reads

[TABLE]

which is identical to Equation (1), with $h_{s}=\left|\mathcal{F}^{-1}\left\{\psi_{d}\psi_{u}\right\}\right|^{2}$, that is to say that the aberration is the sum of the upstream and downstream aberrations. The long exposure optical transfer function is thus

[TABLE]

which is the classical formula by Roddier (1981) of the long exposure transfer function for imaging through turbulence; indeed the autocorrelation of $\psi_{d}\psi_{u}$ is the optical transfer function of the instrument due to diffraction and quasi-static aberrations.

5.3 Approximation in the case of weak turbulence

Additional insight on Equation (11) can be gained by assuming that the residual turbulence phase after the adaptive optics is a small perturbation. It is straightforward to see that

[TABLE]

where $A_{\phi}(r)=\left\langle\phi_{t}(r^{\prime})\phi_{t}(r+r^{\prime})\right\rangle_{r^{\prime}}$ is the autocorrelation of the phase, hence

[TABLE]

Then, if the turbulence is small, $A_{\phi}$ is small compared to $1$, and we can perform a first order Taylor expansion:

[TABLE]

where $S_{\phi}$ is the power spectrum density of the turbulent phase $\phi_{t}$, defined as the (inverse) Fourier transform of $A_{\phi}$. Thus we can express $h_{lec}$ :

[TABLE]

This means that $h_{lec}$ is approximately the non-turbulent coronagraphic point spread function, with a corrective additive term that takes into account the power spectrum density of the turbulent phase, all this dampened by the coherent energy $\mathrm{e}^{-\sigma_{\phi}^{2}}$.

6 Numerical validation and efficiency considerations

6.1 Numerical validation

To validate our point spread function model, we test it against an average of short exposure point spread functions, each one of them with a different outcome of the residual turbulence-induced phase. Each short exposure point spread function is computed using a matrix Fourier transform as proposed by Soummer et al. (2007a) in order to perform an accurate computation of the field on the mask focal plane. A few averages of point spread functions are displayed on Figure 2. They are simulated with a Lyot coronagraph of diameter $3\lambda/D$, the residual turbulence of a SPHERE-like adaptive optics and upstream and downstream white noise phase aberrations of variance $0.1\;\textrm{rad}^{2}$.

Let us define the convergence error as follows. We denote by $h_{lec}$ the long exposure point spread function using the analytic formula of Equation (11). We denote $h_{sseN}$ the average of $N$ short exposures. We define the error as $\mathrm{err}_{N}=\mathrm{max}\left(\left|h_{sseN}-h_{lec}\right|/h_{lec}\right),$ where the division is taken pixel by pixel. We have performed a simulation of the evolution of error versus $N$, the number of short exposures. The error tends to zero as $N^{-1/2}$, which is consistent with the central limit theorem. The error goes below $10^{-3}$ just after $N=10^{7}$ short exposures.

6.2 Computing cost

The evolution of the convergence error gives us an easy criterion to quantify the comparative computing costs of the long exposure coronagraphic point spread function and the empirical average of short exposures. Let us take the computing cost of a short exposure point spread function as the unit computing cost. Then, if we want an error of less than $10^{-3}$ on the point spread function, approximately $10^{7}$ short exposure point spread functions must be averaged, for a cost of $10^{7}$. The analytic formula for the long exposure point spread function has a total cost of the number of points on which the phase structure function is known. This implies that, for square images of $512\times 512$ pixels, our exact formula is about 38 times less costly to evaluate than an average of short exposures.

In addition, it should be noted that, since the long exposure is an integral, and thus, in practice, a sum, it is very easy to compute it in parallel on several processors. We made all calculations in parallel on 16 cores.

Finally, it should also be noted that a useful approximation can accelerate the computing of the long exposure point spread function by a great deal. Indeed, when the tilt $\alpha^{\prime}$ is sufficiently greater than the radius of the focal-plane mask of the coronagraph, the point spread function can be well approximated by a shifted non-coronagraphic point spread function that can be computed once and for all, so the sum in Equation (11) for the computation of the long exposure point spread function must actually only be computed on a square of typical width $8\lambda/D$. For a Nyquist-rate sampled image, that is only $16\times 16$ short exposure point spread functions, with a maximum relative error of less than $10^{-3}$ on each pixel for SPHERE-like parameters.

7 Application: estimation of quasi-static aberrations with COFFEE in the presence of residual turbulence

7.1 Motivation

In current systems such as SPHERE, quasi-static aberrations are only corrected during day-time. For an upgrade, or for future high contrast instruments on extremely large telescopes, it would be useful to correct them during the scientific acquisition. To achieve this, we combine COFFEE (see Paul et al. (2013a), Paul et al. (2014)) with the long exposure phase diversity proposed by Mugnier et al. (2008).

7.2 Simulation results

We have incorporated the long exposure coronagraphic imaging model developed above into the coronagraphic phase diversity method COFFEE. In the following, we perform a preliminary validation of the quasi-static aberration estimation with COFFEE in the presence of residual turbulence by means of simulation. The parameters of the simulations are the following: the phase map size is $64\times 64$ pixels, the incoming flux is $10^{9}$ photoelectrons, the standard deviation of the readout noise is 1 electron, the (supposedly monochromatic) wavelength is 1589 nm. The phase structure function is typical of the residual turbulence of the extreme adaptive optics system of the SPHERE instrument of the VLT, the coronagraph is a Roddier & Roddier phase mask, and there are no downstream aberrations. The upstream aberration is 50 nm RMS. Figure 3 shows the true upstream aberrations, the estimated upstream aberrations, and the difference between the upstream aberrations and the estimated ones, magnified a hundred times. The error in this case is much less than 1 nm RMS, which is a very encouraging result. Additionally, the usefulness of this reconstruction is even better than what the estimated error suggests. Indeed, if we examine the Fourier transform of the difference on Figure 3, we notice that the error is mainly located on the spatial frequencies that are on the border of the array, that is to say, those that are not corrected by the deformable mirror. Quantitatively, the error in the corrected zone accounts for less than 10 % of the total estimation error variance.

8 Conclusion and perspectives

To summarize this work, we have derived an analytic expression to model long exposure image formation for a coronagraphic telescope through the residual turbulence of an extreme adaptive optics system. This model is general and valid for any coronagraphic mask. It only depends on three deterministic quantities: quasi-static aberrations before and after the coronagraph, and the residual phase structure function. Its expression shows, quite nicely, that the long exposure coronagraphic point spread function is the weighted sum of coronagraphic point spread functions without turbulence with a weight distribution given only by the turbulence-induced residual phase structure function. It is computationally effective, allowing a time gain of at least an order of magnitude on standard sized images. Lastly, we have obtained numerical results that suggest that this long exposure model can be applied to the on-sky measurement and correction of aberrations in the context of ground-based direct exoplanet detection.

These encouraging results suggest that we carry on working in two directions. On the one hand, we will explore the robustness of the COFFEE estimations to an error on the phase structure function of the residual atmospheric turbulence. On the other hand, we will demonstrate experimentally that COFFEE can estimate aberrations through turbulence, first on a laboratory bench and then on sky.

Acknowledgements

The PhD work of O. Herscovici-Schiller is co-funded by CNES and ONERA. We thank J.-M. Le Duigou for his support. This work received funding from the E.U. under FP7 Grant Agreement No. 312430 OPTICON, from the CNRS (Défi Imag’In), and from ONERA in the framework of the VASCO research project. We thank Dr Faustine Cantalloube for useful discussions. We thank the reviewer for his careful review and constructive comments.