© ESO, 2014

1. Introduction

Almost twenty years ago Mayor & Queloz (1995) made the first detection of an exoplanet, and more than a thousand are now recognized¹. External occulters will be formidable instruments for performing detailed spectralanalysis of exo-Earths that have previously been detected by other methods.

The first use of an external occulter dates back to Evans (1948), who used an occulting disk set in front of a very small coronagraph of Lyot (1939). Very soon solar astronomers have identified diffraction effects of a circular occulter and considered various kinds of apodized, shaped, or multiple occulters for space missions, as well described in the review paper of Koutchmy (1988). Because of implementation difficulties, toothed or multiple disks are preferred over radially graded ones, which have been described in Purcell & Koomen (1962) and Newkirk Jr & Bohlin (1965). Next to the successful Solar and Heliospheric Observatory mission (Brueckner et al. 1995), future solar coronagraphy envisages formation-flying spacecrafts with an occulter of diameter 1.5 m at about 150 m in front of a small telescope (Vives et al. 2009).

The parameters are completely different for exoplanets,. In rounded numbers, the occulting angle is about 0.2 arcsec against 2000 arcsec for the Sun, and a large 4 m diameter telescope is needed with an occulter of 50 m diameter at a distance of 80 000 km. The positive aspect of exoplanet experiments is that diffraction phenomena are less difficult to simulate numerically than in the solar case. This is due to two reasons. There are fewer than 20 Fresnel zones compared with 8000 in the solar case, and the main source of diffracted light is a point source instead of a huge 1/2° extended object. Nevertheless, the problem remains difficult because of the high expected dynamicrange.

Progression of ideas for exoplanets first followed developments for the Sun, with little interaction, but now the theoretical calculations and instrumental developments have resulted in much more developments than in the solar case. In his review of space astronomy Spitzer (1962) mentioned the technique of apodization envisaged by solar astronomers. Elaborated petal-shaped occulters can already be found in Marchal (1985), and Copi & Starkman (2000) proposed a circular apodizing transmission on a square mask. Studies of shaped and apodized occulters have been reported in many publications, such as Cash (2006), Arenberg et al. (2007), Vanderbei et al. (2007), Soummer et al. (2010), Cash (2011) and Wasylkiwskyj & Shiri (2011), to cite just a few. Although partially transparent petal shaped occulters are envisaged by Shiri & Wasylkiwskyj (2013), the most advanced technological developments are for petaled star shades. The two-step procedure leading to an optimal shaped occulter is well described in Kasdin et al. (2013). These authors explained that they first seek for the optimal variable transmission of an apodized occulter and then chose a sufficient number of petals for the shaped occulter to best approach the theoretical result given by the variable transmission. We here focus only on the first part of this approach, that is, on the search for an ideal transmission, leaving the final design of a shaped occulter to a future work.

Our purpose is to show that the occulter can be optimized advantageously from the telescope focal plane, which will minimize the level of light that is diffracted by the star in a zone of interest of the focal plane for the exoplanet. The resulting illumination on the telescope aperture appears to be apodized from center to edge. The integrated flux on the telescope pupil is no longer minimized. The observation is nevertheless improved because most of the residual flux is trapped behind the geometrical image of the mask, in a blind area for the observation of the exoplanet. Differences in transmission with an occulter that minimizes the flux over the aperture are small, but the gain in the residual background light at the level of the exoplanet range from 1.9 to 50 depending on the spectral bandwidth, which might correspond to a gain in integration time of a factor 3.6 to 2500 for certain observing conditions, inphoton counting mode.

The paper is organized as follows: the fundamental relations that describe the complex amplitudes diffracted by the occulter in the aperture and focal planes are given in Sect. 2. The procedure for optimizing the transmission of a circularly symmetric apodized occulter to minimize the flux on the telescope aperture or for selected regions of the focal plane are described in Sect. 3. Section 4 contains results on the optimization problem and a discussion. Additional information is given in two appendices based on an analytic approach of the Fresnel diffraction of a pure circular disk and of a linearly apodized disk.

2. Expressions of complex amplitudes and intensities in the aperture and focal planes

We denote with 2Ω the overall diameter of the occulter. For an occulting mask, it is convenient to write its transmission as t(r) = 1 − f(r), with the attenuation function f(r) constrained by0 ≤ f(r) ≤ 1 for | r | ≤ Ω, and f(r) = 0 for | r | > Ω. For a wave of unit amplitude, the Fresnel diffraction at the telescope aperture at a distance z from the occulter (see Fig. 1) can be written as $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathit{z}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{i\lambda z}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Omega }}\mathrm{2}\mathit{\pi \xi f}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{\tau }}_{\mathit{z}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{J}}_{\mathrm{0}}\left(\mathrm{2}\mathit{\pi }\frac{\mathit{\xi r}}{\mathit{\lambda z}}\right)\mathrm{d}\mathit{\xi ,}& & \end{array}$(1)where τ_z(r) = exp(iπr²/λz) is a quadratic phase-term corresponding to a diverging lens, and J₀(r) is the Bessel function of the first kind.

This relation is equivalent to Eq. (4) of Vanderbei et al. (2007), neglecting here the term of propagation of a plane wave. In general, this integral does not admit an analytical solution, except for f(r) = 1, as shown by Aime (2013) and outlined here in Appendix A. Therefore, Eq. (1), which appears as the Hankel transform of f(r)τ_z(r), must be evaluated numerically. Computations of these transforms are delicate, as discussed by Lemoine (1994) for example, who recommended a sampling of the function based on zeroes of the Bessel function J₀. In practice, after several tests described in Appendix A and B and below, we directly used the function NIntegrate of Mathematica, which performs very well for this kind of calculation, thanks to recent improvements of numerical integration.

This wavefront ψ(r) arrives onto the aperture of the telescope of transmission P(r), and we denote with ϕ(r) = P(r)ψ(r) the complex amplitude of the wave going through the telescope aperture. For the sake of simplicity, we assume in the following a perfectly circular telescope of diameter 2R, centered on the optical axis of the system, so that P(r) = Π(r/R), defining here for convenience that the box distribution Π(r) equals 1 for | r | < 1 and 0 elsewhere. Note that a more realistic telescope with a central obscuration can be easily inserted into the calculations. Using the circular symmetry of the problem, the complex amplitude of the wave in the focal plane can be written as a Hankel transform of ϕ(r) of the form: $\begin{array}{ccc}\mathit{\phi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\tau }}_{\mathit{F}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{i\lambda F}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\mathrm{2}\mathit{\pi \xi \varphi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{J}}_{\mathrm{0}}\left(\mathrm{2}\mathit{\pi }\frac{\mathit{\xi r}}{\mathit{\lambda F}}\right)\mathrm{d}\mathit{\xi }\mathrm{=}\frac{{\mathit{\tau }}_{\mathit{F}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{i\lambda F}}\mathit{\varphi ̂}\left(\frac{\mathit{r}}{\mathit{\lambda F}}\right)\mathit{,}& & \end{array}$(2)where F is the focal length of the telescope. Note that here the phase term τ_F(r) is compensated for by the lens phase function of the telescope, meaning that a simple 2D Fourier transform or Hankel transform denoted by substitutes the more complex Fresnel propagation of Eq. (1).

It is moreover convenient to calibrate the focal plane in terms of angular units θ = r/F on the sky, and the intensity in the focal plane of the telescope can be written $\begin{array}{ccc}\mathrm{\Phi }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\lambda }}^{\mathrm{2}}}{\left|\mathit{\varphi ̂}\left(\frac{\mathit{\theta }}{\mathit{\lambda }}\right)\right|}^{\mathrm{2}}\mathit{,}& & \end{array}$(3)instead of just the Airy pattern $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{0}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\lambda }}^{\mathrm{2}}}{\left|\widehat{{\mathit{\varphi }}_{\mathrm{0}}}\left(\frac{\mathit{\theta }}{\mathit{\lambda }}\right)\right|}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\pi }}^{\mathrm{2}}{\mathit{R}}^{\mathrm{4}}}{{\mathit{\lambda }}^{\mathrm{2}}}{\left(\frac{{\mathit{J}}_{\mathrm{1}}\mathrm{(}\mathrm{2}\mathit{\pi \theta R}\mathit{/}\mathit{\lambda }\mathrm{)}}{\mathit{\pi \theta R}\mathit{/}\mathit{\lambda }}\right)}^{\mathrm{2}}& & \end{array}$(4)for the direct observation without external occulter with a perfect telescope of circular aperture of diameter 2R.

Limiting the analysis to only geometrical aspects, the occulter prohibits the planet observation for an angular radius smaller than Ω /z to the star, and the light coming from the planet passes over the occulter for an angle , a variant of the geometric inner working angle Ω /z of Kasdin et al. (2013), taking into account the telescope size. Note that all the wave and intensity functions defined in this section are illustrated in Fig. 1.

3. Procedure used for optimizing f(r) for integrated intensities

The basic idea of an external occulter coronagraph is to avoid letting the stellar light enter the telescope while leaving the light coming from the exoplanet unchanged. In a first approximation, if the angular position of the planet is larger than θ₀, we can expect that the flux coming from the planet is almost unaffected by the occulter, and a fair measurement of the efficiency of the external occulter can be the normalized residual integrated intensity over the telescope aperture, for example the quantity Γ defined by $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathrm{\Delta }\mathit{\lambda }}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{m}}}^{{\mathit{\lambda }}_{\mathrm{M}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\mathit{r}\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\lambda ,}& & \end{array}$(5)where Ψ(r) = | ψ(r) | ² is the wavelength-dependent intensity and Δλ is the spectral bandwidth of the experiment, simply taken equal to λ_M − λ_m here. In this analysis, the same weight is assigned to each wavelength, independently of the brightness of the source, the quantum efficiency of the detector, or the presence of acquisition filters (Soummer et al. 2010). It would be easy to perform a differently weighted average for a more realistic analysis, if necessary. Our numerical computation has been made from λ_m = 380 nm to λ_M = 750 nm and with R = 2 m.

To find an optimal shape for the external occulter, a natural approach would be to minimize Γ with respect to f. Note that this approach has been used in Wasylkiwskyj & Shiri (2011), but this is not the only way to optimize the shape of an occulter. For instance, Vanderbei et al. (2007) and Kasdin et al. (2013) proposed to minimize the upper bound c on the real and imaginary part of ψ(r),∀r ∈ [ 0,R ]. This approach will ensure that the intensity in the aperture plan will be below 2c² and also minimize the flux. We emphasizethat while we decided here to minimize the flux, our approach can be readily adapted to the minimization of the upper bound instead of the whole flux Γ.

Because of the conservation of energy, the flux is conserved from aperture to focal planes. Nevertheless, we are in fact interested in the level of background light produced by the star in the area where we observe the planet. A quantity representative of the residual flux may be given by the residual light γ over a zone of the focal plane, such as $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathrm{\Delta }\mathit{\lambda }}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{m}}}^{{\mathit{\lambda }}_{\mathrm{M}}}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}\mathit{\theta }\mathrm{\Phi }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda ,}& & \end{array}$(6)where θ₁ and θ₂ define the region of interest of the observation (Fig. 1), corresponding for example to the habitable zone around the star, or the possible excursion of a known planet over a star, with probably θ₁ close to θ₀. We therefore have γ ≤ Γ, the equality holding in the limit θ₁ = 0 and θ₂ = ∞.

The two measures of residual light Γ and γ are based on different points of view. If one manages to completelyshut offthe light in the aperture, then the residual light in the observation area will also shut off, but the reverse is false. In other words, minimizing Γ or γ with respect to f will lead to different optimal solutions. We believe that minimizing γ is better when the objective is to observe an exoplanet in a known area of the focal plane.

3.1. Decomposition of f(r) on a basis of functions f_k(r)

As already presented in Jacquinot & Roizen-Dossier (1964) for the systematic search for apodizing properties, a common way to optimize the shape of a function is to force this function to be a weighted sum of basis functions. We here aim to optimize the weights α_k such that $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}{\mathit{f}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{and}\mathrm{0}\mathrm{\le }\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\le }\mathrm{1.}& & \end{array}$(7)Jacquinot & Roizen-Dossier (1964) have shown that, if the expansion contains an infinite number K of terms, “the absolute optimal function according to the criterion chosen” is obtained regardless of the bases used. This result was obtained in the different context of Fraunhofer diffraction, but given the linearity of the equations, the result still holds for our purpose. In practice, however, the number K of elements of the bases must be limited and the choice of basis functions becomes important.

After various tests, we used trapezoidal functions for the f_k(r), equal to 1 for r ≤ Ω_m + (k − 1)Δ and equal to 0 for r ≥ Ω_m + kΔ. The functions linearly decrease from 1 to 0 between Ω_m + (k − 1)Δ and Ω_m + kΔ. These function, illustrated in Fig. 2, are very similar to the binary disk for small Δ. It is the same for their Fresnel diffraction pattern, as shown in Appendix B, but this small difference is sufficient to lead to much better numerical results. Note that these basis functions lead to a transmission that is piece-wise linear as opposed to the use of binary disk that leads to piece-wise constant transmission function. For the numerical analysis, we set Δ equal to 5 cm and used K = 300 functions spreading between Ω_m = 10 m to Ω_M = Ω = 25 m. An example of such functions is given in Fig. 2. The reason for the existence of an innermost opaque section (i.e. t(r) = 0 or f(r) = 1 for r ≤ Ω_m) is of technical origin linked with the spacecraft, as indicated for example by Vanderbei et al. (2007).

Note that expressing function f(r) as a sum of basis function leads to a more general model potentially at the cost of more complex numerical computation. For instance, one might choose to use a more elaborated basis of functions such as prolate functions. This would potentially lead to very good results because those functions are known to provide efficient apodizers, and might also be excellent occulters. Nevertheless, the physical constraints on f(r) as discussed in the following would be more difficult to express. In this particular case, it would be advantageous to compute these constraint on a finite sampling of the function, as in Vanderbei et al. (2007).

3.2. Constraints on f(r) and f_k(r)

The finalattenuation function f(r) has to follow a number of physical constraints. On the one hand, some of those constraints can be easily enforced through a wise choice of the basis functions. For instance the smallest radius of apodization is set by Ω_m, while the largest radius is selected by Ω = Ω_m + KΔ. On the other hand, the attenuation function f(r) has to be equal to 1 for r< Ω_m, and since all f_k(r) functions are equal to 1 for r ≤ Ω_m, we have to enforce the linear constraint ${\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}\mathrm{=}\mathrm{1}$. This constraint can also be expressed as a linear equality of the form 1^⊤α = 1, where 1 is a vector of ones.

Another more complicated constraint is the fact that 0 ≤ f(r) ≤ 1, ∀r. These block constraints, thanks to the simple structure of the f_k functions, can be expressed as a linear inequality constraints Aα ≥ b with $\begin{array}{ccc}{A}\mathrm{=}\left[\begin{array}{c}\\ \\ \end{array}\right] {C}\mathrm{=}\left[\begin{array}{c}\\ & & & & \\ & & & & \\ & & & & \\ & & & & \\ & & & & \end{array}\right]\mathit{,} \mathrm{and} {b}\mathrm{=}\left[\begin{array}{c}\\ \\ \\ \\ \\ \\ \end{array}\right]\mathit{.}& & \end{array}$Another constraint that has been used in previous works is the smoothness of the transmission function. This is typically enforced using derivatives of the function, but in our case smoothness can be readily enforced using a classical quadratic regularization term α^⊤α. Indeed, minimizing the euclidean norm of the α vector under the constraint 1^⊤α = 1 tends to promote similar values on all α_k, leading to a smooth transmission function (linear apodization for an infinite regularization). This kind of regularization is equivalent to the constraint on the curvature of f described in Kasdin et al. (2013).

An additional constraint that has been enforced in previous literature is to force f to be monotonic, which in this case is a decreasing function w.r.t. the radius r. This constraint, also enforced in Kasdin et al. (2013), has an extremely simple form in our formulation, in effect it is the positivity constraint α ≥ 0. To have a grasp of the performance loss due to this additional constraint, we performed an optimization with and without this constraint in the numerical experiments.

Note that in this work, we made the design choice of using linear apodization functions as basis functions for the transmission. These choices lead to the above-mentioned physical constraints on the weights α. An equivalent formulation has been proposed in Vanderbei et al. (2007) and Kasdin et al. (2013), where the basis function are trapezoidal functions overlapping such that the final transmission is also piece-wise linear. Copi & Starkman (2000), Cash (2011) and Wasylkiwskyj & Shiri (2011), on the other hand, used offset high-order polynomial or hyper-Gaussian functions.

Below, we present the optimization problems derived from minimizing the intensity Γ in the telescope aperture and minimizing the intensity γ in a region of the focal plane. The optimal functions are denoted as f_Γ(r) and f_γ(r) for the optimal attenuation w.r.t. the aperture plane and w.r.t. a selected region of the focal plane, respectively.

3.3. Minimizing the integrated starlight in the telescope aperture plane: f_Γ(r)

Substituting the decomposition of f(r) of Eqs. (7) into (1), we obtain $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}\left(\mathrm{1}\mathrm{-}{\mathit{\tau }}_{\mathit{z}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Omega }}{\mathit{f}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\tau }}_{\mathit{z}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{J}}_{\mathrm{0}}\left(\mathrm{2}\mathit{\pi }\frac{\mathit{\xi r}}{\mathit{\lambda z}}\right)\mathrm{d}\mathit{\xi }\right)\\ & \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}{\mathit{\psi }}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}\end{array}$(8)and the corresponding intensity becomes $\begin{array}{ccc}\mathrm{\Psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& {\left|\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\right|}^{\mathrm{2}}\mathrm{=}{\left|\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}{\mathit{\psi }}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right|}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{k,l}\mathrm{=}\mathrm{1}}^{\mathit{K,K}}{\mathit{\alpha }}_{\mathit{k}}{\mathit{\alpha }}_{\mathit{l}}{\mathit{\psi }}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\psi }}_{\mathit{l}}\mathrm{(}\mathit{r}{\mathrm{)}}^{}\\ & \mathrm{=}& {{\alpha }}^{\mathrm{\top }}{{K}}_{\mathit{\psi }}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\alpha }\mathit{,}\end{array}$(9)with with ψ^⊤ = [ ψ₁(r),...,ψ_K(r) ] a rank one matrix of general coefficient ${\mathit{K}}_{\mathit{l,k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\psi }}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\psi }}_{\mathit{l}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, where ψ_k(r) is the Fresnel diffraction of a wave of unit amplitude for the basis function occulter f_k(r). The flux Γ can then be expressed as $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathrm{\Delta }\mathit{\lambda }}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{m}}}^{{\mathit{\lambda }}_{\mathrm{M}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\mathit{r}{{\alpha }}^{\mathrm{\top }}{{K}}_{\mathit{\psi }}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\alpha }\mathrm{d}\mathit{r}\mathrm{d}\mathit{\lambda }\mathrm{=}{{\alpha }}^{\mathrm{\top }}{{K}}_{\mathrm{\Gamma }}{\alpha }\mathit{,}& & \end{array}$(10)where ${{K}}_{\mathrm{\Gamma }}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathrm{\Delta }\mathit{\lambda }}{{}^{\mathrm{\int }}\mathit{\lambda }_{\mathrm{m}}}_{{\mathit{\lambda }}_{\mathrm{M}}}^{}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{R}}\mathit{r}{{K}}_{\mathit{\psi }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\lambda }$ is a K × K matrix integrating all the surface of the aperture and all the wavelength in (λ_m,λ_M). The optimization problem is finally defined as $\begin{array}{ccc}& & \underset{{\alpha }}{\min } {{\alpha }}^{\mathrm{\top }}\mathrm{(}{{K}}_{\mathrm{\Gamma }}\mathrm{+}\mathit{\mu }{I}\mathrm{)}{\alpha }\\ & & \mathrm{s}\mathrm{.}\mathrm{t}\mathrm{.} {\mathrm{1}}^{\mathrm{\top }}{\alpha }\mathrm{=}\mathrm{1} \mathrm{and} {A}{\alpha }\mathrm{\ge }{b}\mathit{,}\end{array}$(11)where I is the identity matrix and μ is a positive regularization parameter. As discussed in Wasylkiwskyj & Shiri (2011), regularization is important in this case because the matrix K_Γ is poorly conditioned and would lead to difficulties when solving the problem. Moreover, as discussed in the previous section, the μ parameter weight the quadratic regularization on the α vector, promoting smoothness in the transmission function. This smoothness has also been shown to be important in practice in Kasdin et al. (2013), where an equivalent smoothness constraint was inserted in the optimization problem. Also note that in practice one can easily express the smoothness constraint of Kasdin et al. (2013) in our optimization problem instead of using a regularization term. Since with the trapezoidal basis, the second-order derivatives are exactly computed using the finite difference between adjacent α components, the corresponding constraints can be obtained by adding lines to the A matrix that contain finite difference operators, and lines to b that contain the constraint parameter σ of Kasdin et al. (2013).

The optimization problem defined Eq. (11) is a constrained quadratic optimization Boyd & Vandenberghe (2004). The linear constraints forbid the use of the closed-form solution as proposed in Wasylkiwskyj & Shiri (2011), but there exist several efficient optimization strategies in the literature, among which the active set approach presented in Vanderbei & Shanno (1999).

3.4. Minimizing the starlight in a region of the telescope focal plane: f_γ(r)

To express, the intensity in the focal plane, we used the same decomposition of the function f(r) as a sum of functions f_k(r) as above, the complex amplitude in the focal plane can be written as a sum of elementary function φ_k(r), Fourier transforms of P(r)ψ_k(r) according to Eqs. (7) and (2). The intensity in the focal plane of Eq. (3) becomes $\begin{array}{ccc}\mathrm{\Phi }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\lambda }}^{\mathrm{2}}}{\left|\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{\alpha }}_{\mathit{k}}{\mathit{\phi }}_{\mathit{k}}\left(\frac{\mathit{\theta }}{\mathit{\lambda }}\right)\right|}^{\mathrm{2}}\mathrm{=}{{\alpha }}^{\mathrm{\top }}{{K}}_{\mathit{\phi }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\alpha }\mathit{,}& & \end{array}$(12)and the intergrated density $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathrm{\Delta }\mathit{\lambda }}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{m}}}^{{\mathit{\lambda }}_{\mathrm{M}}}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}}^{{\mathit{\theta }}_{\mathrm{2}}}\mathit{\theta }{{\alpha }}^{\mathrm{\top }}{{K}}_{\mathit{\phi }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\alpha }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }\mathrm{=}{{\alpha }}^{\mathrm{\top }}{{K}}_{\mathit{\gamma }}{\alpha }\mathit{.}& & \end{array}$(13)The resulting optimization problem is exactly the same that of Eq. (11), but with a different metric K_γ.

The optimal functions f_Γ(r) and f_γ(r) are different because the solutions α are different. Because of conservation of energy between aperture and focal planes, these functions become identical in the special case θ₁ = 0,θ₂ = ∞, where K_Γ = K_γ. Our optimization in the focal plane was performed for θ₁ = 0.1 arcsec and θ₂ = 0.5 arcsec.

4. Results and discussion

4.1. Numerical computations

In the following, we denote Ψ_Γ(r) and Φ_Γ(θ) the intensities obtained in the aperture and focal plane of the telescope using an external occulter of transmission 1 − f_Γ(r) minimizing the integrated intensity Γ and described in Sect. 3.3. Similarly, we denote Ψ_γ(r) and Φ_γ(θ) the results obtained using an occulter transmission 1 − f_γ(r), as described in Sect. 3.4.

Note that while obtaining the matrices K_Γ and K_γ is a computationally intensive problem, we can use heavily parallel computing to estimate as a first step the ψ_k(r) and φ_k(θ) functions through numerical integration. This was performed using Mathematica, with K computational jobs computed by ≈100 parallel processes on eight core Intel processors. In the numerical experiments, both ψ_k(r) and φ_k(θ) were regularly sampled on 3001 samples from respectively r ∈ (0,R + 1) and θ ∈ (0,0.644) arcsec. Note that ψ_k(r) was computed outside of the telescope so that we can used it for the tolerance analysis in Sect. 4.4. Moreover, the behavior of the different optimization strategies outside of the telescope aperture is also interesting. The wavelength λ has also been sampled regularly with 21 samples in the visible light interval (380,750) nm. After computing ψ_k(r) and φ_k(θ), a numeric 2D integration can be performed to obtain the matrices K_Γ and K_γ. The regularization parameter μ in the optimization problem was chosen adaptively with μ = μ₀max_k,l | K_k,l | to be less dependent on the metric matrix K. μ₀ was chosen by hand to allow good attenuation performance with a limited dynamic in the solution function f. For the monochromatic study, we set μ₀ = 10^-8 and for the more complex study on a wide bandwidth, less regularization was necessary with μ₀ = 10^-10 because of a more complex problem and a better conditioning of the K matrix.

Finally, the code for all the numerical experiments will be freely available on the authors website to promote reproducible research².

4.2. Optimization for monochromatic light

For a single monochromatic wave, here λ = 562 nm, we have represented the apodizing curves f_Γ(r) and f_γ(r) in Fig. 3 and the resulting intensities Ψ_Γ(r) and Ψ_γ(r) (top curves), and Φ_Γ(θ) and Φ_γ(θ) (bottom curves) in Fig. 4. We emphasize that these curves were obtained assuming a perfect circular telescope of radius R. The f-curves are only slightly different and are really similar to a linear apodization. Moreover, in the monochromatic case, enforcing a monotonic f leads to the same results for the thickness of the trace, which will not be the case anymore for a wide spectral bandwidth, as we show below.

The results for the Ψ and Φ-curves are quite different. The residual starlight Ψ_Γ(r) is very well attenuated across the telescope pupil, as found by Vanderbei et al. (2007). In contrast, Ψ_γ(r) presents a center-to-limb decrease which suggests an apodized structure.It is interesting to note the behavior of the curves beyond the region of optimization. There, the relative position of the curves is reversed very rapidly, the focal plane optimization becoming lower than the pupil plane optimization. This is an aspect that deserves further studies.

The intensity in the focal plane is also very interesting. While f_Γ tends to minimize the intensity at all angles θ, f_γ only focuses on the region of interest, and a strong residual light is observed in the blind area behind the occulter.

The resulting flux Γ,γ for all the different optimizations is given in the upper part of Table 1. Interestingly, f_Γ leads to a flux Γ in the aperture plane 10⁵ times fainter than f_γ. Nevertheless, in the focal plane the flux γ in the observation area is ≈50 times fainter when using f_γ. This clearly shows the importance of using focal plane optimization for optimal apodization.

4.3. Optimization for a wide spectral bandwidth

Next we investigated the two optimization strategies on a wide bandwidth. For this we compute a regular sampling in the interval λ ∈ (380, 750) nm and computing the average flux across λ. Note that both a sampling of 20 and 100 values were investigated, and while we used in the experiments a sampling of 100 for a better precision, the final results were very similar. In practice the variation of the wavefront is extremely smooth w.r.t.λ and the selected sampling in sufficient.

For multiple wavelengths, λ ∈ (380, 750) nm, the results are very different. First the optimal functions are much more similar, as illustrated in Fig. 5. The functions also tend to oscillate and the monotonic variants are this time different from the free variant of the optimal functions.

The resulting intensities in the aperture and focal plan are plotted in Fig. 6. In the aperture plane, the intensity is still apodized for f_Γ, but the flux is much more attenuated than in the monochromatic case. In the focal plane, the intensities retain the same tendencies, but the gain due to the focal plane optimization is lower. Interestingly, enforcing a monotonic function through α ≥ 0 leads to a clear loss in terms of attenuation for both f_Γ and f_γ.

Intermsof quantitative performances, Table 1 (lower part) shows that the flux in the observation area of the focal plane is still more attenuated with f_γ with a gain of 1.9 in the chromatic case (1.6 with positivity constraint). While the gain in performance is smaller than in the monochromatic case, it is far from negligible when observing faint objects around a star. Interestingly, in the chromatic case, the α ≥ 0 constraint leads to a loss of 2 in attenuation (for both f_Γ and f_γ), proving that relaxing the constraint can lead to additional gain as long as the apodization is physically realizable.

An additional visualization comparing the performances of f_Γ versus f_γ in the focal plane is given in Fig. 7. The intensity Φ_γ(θ) is plotted as a function of Φ_Γ(θ) for θ ∈ (0,0.7) arcsec. The red line allows a comparison between the methods; if a curve is below this line, it means that Φ_γ(θ) < Φ_Γ(θ) and vice versa. We clearly see that Φ_γ(θ) > Φ_Γ(θ) for θ<θ₁ in the blind area behind the occulter and that Φ_γ(θ) < Φ_Γ(θ)∀θ ∈ (θ₁,θ₂) in the observation area of the focal plane (illustrated in green in the figure).

4.4. Tolerance analysis to positioning errors of the telescope

A study of the tolerance sensitivity to positioning errors for a 4 m circular telescope out of its on-axis position was carried out up to an offset of 1 m.

Because of this offset, the resulting pupil plane and focal plane intensities are no longer circularly symmetric, and the computation of the final focal plane image can no longer use the Hankel transformation of Eq. (2). The focal plane image must be computed taking the modulus squared of the two-dimensional Fourier transform of the complex amplitude of the wave on the telescope aperture.

In practice, the study was made using a discrete Fast Fourier Transform of an array of 1024 × 1024 points corresponding to an overall zone of 20 × 20 m, inside which the telescope of diameter 2 m was defined by 1 or 0 transmitting pixels over about 205 points in diameter. As a result, altogether 33 081 points were set equal to 1. This percentage of zero-padding is enough to obtain a satisfactorily sampled focal plane image.

Positioning errors were simulated by moving the complex amplitude obtained in the shadow of the occulter over the telescope aperture in one direction, by steps of 1 pixel or about 2 cm. The two-dimensional array was filled using an interpolation of the one-dimensional complex amplitude of the wavefront computed with Eq. (1).

The focal plane image was obtained taking the modulus squared of the Fourier transform of the wave on the telescope aperture. The effect of wavelength was taken into account afterward. The diffraction pattern increases in size with the wavelength. For that, a two-dimensional interpolation of the result was used to resample the image. Moreover, a scaling factor in intensity inversely proportional to the square of the wavelength was applied, as in Eq. (4). As a result, the integrated flux in the focal plane becomes wavelength independent and equal to that crossing through the telescope aperture. As before, for the sake of simplicity, we assumed that the received flux is constant over the spectral bandwidth.

Figures 8 and 9 present results obtained for transverse displacements of the telescope up to 1 m from its nominal position. We recall that the occulter shapes optimize the measurement when the telescope is at its nominal on-axis position, either for the integrated intensity on the telescope aperture, or for a region in the focal plane between 0.1 and 0.5 arcsec. In both cases the computation was made for the whole bandwidth of λ ∈ (380, 750) nm.

Three-dimensional representations of the normalized intensities in the telescope aperture were used to enhance differences of intensities in the results of the two optimizations, and the apodization-like pattern obtained for the focal plane optimization. Both curves were clipped to make the central parts of the images well visible. As expected, there is much less light for the aperture plane minimization when the telescope is at its nominal position. Since the optimization was calculated for a 4 m diameter telescope, there is a rise in the intensity at the edges for this large offset. This increases becomes somewhat surprisingly in favors of the focal plane optimization for a large off-set, as already discussed in Fig. 6.

In Fig. 8 focal plane intensities are also represented in color levels, with a common color scale for the telescope on-axis and off-axis of 50 cm and 1 m. The two circles of radius θ₁ = 0.1 arcsec and θ₂ = 0.5 arcsec define the angular area inside which the residual intensity is integrated. The focal plane optimization is found to always be better than the aperture plane optimization. From comparing these curves, it is clear that the residual diffracted light remains concentrated longer inside the first circle for the focal plane optimization, which is especially clear for the 50 cm offset images. This behavior is confirmed by Fig. 9 which represents the integrated intensity in this region with a logarithmic scale.

5. Conclusion

We addressed the problem of optimizing of the shape of a circularly symmetric apodized occulter. To this end, we proposed a general expression for the attenuation function, as a weighted sum of basis functions. We expressed the optimization problem as the minimization of the flux, first in the pupil plane, as is commonly done in the literature, and then in the focal plane of a 4 m diameter telescope, for a region between 0.1 to 0.5 arcsec in which the exoplanet is supposed to be observed.

Numerical experiments show that the focal plane optimization leads to a better attenuation of the starlight than the pupil plane optimization. One interesting result is that while the focal plane optimization leads to a strong flux in the pupil plane, most of this light is concentrated in the area behind the occulter in the focal plane, which is not a problem for a perfect telescope. Finally, the robustness of our approach to lateral misalignment of the occulter was investigated and showed that the focal plane optimization also yields a better robustness.

As already indicated, the optimization study was made for a 4 m telescope assumed to be at its on-axis position. A simple procedure to make the experience less sensitive to defect position would be to perform the optimization for a larger aperture than is actually used. A better procedure would be to estimate the statistics of pointing errors and include them in the optimization procedure. Then a better tolerance to pointing errors would be obtained, probably at the expenses of a tolerable loss in the starlight rejection. Nevertheless, our analysis showed that telescope offsets of a few decimeters will not strongly reduce the efficiency of the occulter.

This conclusion was obtained for the condition that the telescope has a perfectly circular aperture and another study is required to model a telescope with central obscuration. Moreover, it would be interesting to study the effect of a a more realistic telescope with small defaults of phase and amplitude. Deviations from this perfection could lead to somewhat different results, but the exigency of quality is limited to a maximum of 0.5 arcsec, or about twenty Airy rings for a 4 m diameter telescope. Finally, with the knowledge of the exact PSF of the telescope, one can readily adapt our optimization approach. A comparison of pupil vs focal optimization for the maximum-magnitude minimization as in Kasdin et al. (2013) would also be interesting.

---

Acknowledgments

We want to thank the referee Jeremy Kasdin for his very constructive comments, and in particular for his suggestion to study the effect of telescope misalignment. This work was granted access to the HPC and visualization resources of Centre de Calcul Interactif hosted by Université Nice Sophia Antipolis (http://calculs.unice.fr/en).

References

Appendix A: Quality test of the numerical computation using an analytical expression for a sharp-edged circular occulter

Aime (2013) has shown, using an approach similar to Born & Wolf (1999), that an analytic expression using Lommel series could be derived for the Fresnel diffraction of a sharp circular occulter of the form f(r) = f_•(r,Ω) = Π(r/ Ω), defining here Π(r) as the function equal to 1 for | r | ≤ 1 and 0 otherwise. At the wavelength λ, the complex amplitude ψ_•(r,Ω) of the wave at the distance z for

• (i)

  inside the geometrical umbra (r< Ω) and,

• (ii)

  outside (r> Ω), can be written using two series: $\begin{array}{ccc}{\mathit{\psi }}_{\mathrm{•}}\mathrm{\left(}\mathit{r,}\mathrm{\Omega }\mathrm{\right)}& \mathrm{=}& \\ \mathrm{\left(}\mathit{i}\mathrm{\right)}& & \mathit{\tau }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\tau }\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{n}}{\left(\frac{\mathit{r}}{\mathrm{\Omega }}\right)}^{\mathit{n}}{\mathit{J}}_{\mathit{n}}\left(\frac{\mathrm{2}\mathit{\pi }\mathrm{\Omega }\mathit{r}}{\mathit{\lambda z}}\right)\\ \mathrm{\left(}\mathit{ii}\mathrm{\right)}& & \mathrm{1}\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\tau }\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{n}}{\left(\frac{\mathrm{\Omega }}{\mathit{r}}\right)}^{\mathit{n}}{\mathit{J}}_{\mathit{n}}\left(\frac{\mathrm{2}\mathit{\pi }\mathrm{\Omega }\mathit{r}}{\mathit{\lambda z}}\right)\mathrm{·}\end{array}$(A.1)At r = Ω, both series converge to the same value given by $\begin{array}{ccc}{\mathit{\psi }}_{\mathrm{•}}\mathrm{\left(}\mathrm{\Omega }\mathit{,}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\exp \left(\mathit{i}\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\lambda z}}\right){\mathit{J}}_{\mathrm{0}}\left(\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\lambda z}}\right)\right)\mathit{.}& & \end{array}$(A.2)

These series are alternative series for the real and imaginary terms, and according to the Leibniz estimate, an upper bound of the error for the sum limited to n = N terms is given by the absolute value of the N + 1 term. The convergence is easy when r/ Ω is smaller than 1, for the first sum for r< Ω inside the geometrical umbra, or the inverse outside. The main error occurs near the transition zone r = Ω.

In Fig. A.1 we compare the analytical and numerical results for an example of a circular occulter of 50 m diameter set at 80 000 km, and for a wavelength of 550 nm. The Lommel series were computed for 100 terms, and there is no difference between the two methods. Results were drawn for the real, imaginary, modulus, and unwrapped phase of the Fresnel diffraction. Note that this is just an example, in particular, the real and imaginary part of the diffraction pattern are highly sensitive to a small variation of any of the parameters. The modulus of the intensity pattern clearly shows the visible Poisson spot at the center of the pattern. The unwrapped phase presents a strong phase variation in the center of the umbra and a quieter structure outside the geometric umbra, that is for a region utilized by the telescope for the planet observation. The consistency between the two analytic and numerical results is recognized here as a proof of the quality of the numerical calculation.

Appendix B: Numerical calculation of linearly apodized basis functions f_k(r)

As indicated in the body of the paper, we tried several forms for the set of apodized basis functions f_k(r), and decided to use the simple trapezoidal functions described in the body of the paper because they affect the final solution the least. We did notsucceed in deriving an analytic form for the Fresnel diffraction of these functions, but comforted by the results obtained on the Fresnel diffraction of Π(r/ Ω), we used the direct numerical calculation of Mathematica.

We recall that these functions are linearly apodized disks of radii varying from 10 m to 25 m in steps of 5 cm, the apodization applies for 5 cm at the edge. The resulting ψ_k(r) are very similar to those obtained for the functions Π(r/ Ω). Differences between moduli are lower than 1%, as shown in Fig. B.1 for two extreme cases, but essential to obtain proper results on diffraction patterns. In Fig. B.2 we focus on the central part of the diffraction pattern, for r ≤ 2 m. Curves are drawn for a raw and linearly apodized occulter of 10 m and 10.05 m.

All Tables

All Figures