Optimization of starshades: focal plane versus pupil plane
Laboratoire Lagrange, UMR CNRS 7293, Université de Nice SophiaAntipolis,
Observatoire de la Côte d’Azur,
Parc Valrose,
06108
Nice,
France
email: remi.flamary@unice.fr; claude.aime@unice.fr
Received:
20
February
2014
Accepted:
17
July
2014
We search for the best possible transmission for an external occulter coronagraph that is dedicated to the direct observation of terrestrial exoplanets. We show that better observation conditions are obtained when the flux in the focal plane is minimized in the zone in which the exoplanet is observed, instead of for the total flux received by the telescope. We describe the transmission of the occulter as a sum of basis functions. For each element of the basis, we numerically computed the Fresnel diffraction at the aperture of the telescope and the complex amplitude at its focus. The basis functions are circular disks that are linearly apodized over a few centimeters (truncated cones). We complemented the numerical calculation of the Fresnel diffraction for these functions by a comparison with pure circular disks (cylinder) for which an analytical expression, based on a decomposition in Lommel series, is available. The technique of deriving the optimal transmission for a given spectral bandwidth is a classical regularized quadratic minimization of intensities, but linear optimizations can be used as well. Minimizing the integrated intensity on the aperture of the telescope or for selected regions of the focal plane leads to slightly different transmissions for the occulter. For the focal plane optimization, the resulting residual intensity is concentrated behind the geometrical image of the occulter, in a blind region for the observation of an exoplanet, and the level of background residual starlight becomes very low outside this image. Finally, we provide a tolerance analysis for the alignment of the occulter to the telescope, which also favors the focal plane optimization. This means that telescope offsets of a few decimeters do not strongly reduce the efficiency of the occulter.
Key words: methods: analytical / methods: numerical
© 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^{1}. External occulters will be formidable instruments for performing detailed spectralanalysis of exoEarths 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 formationflying 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 petalshaped 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 twostep 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 (1)where τ_{z}(r) = exp(iπr^{2}/λz) is a quadratic phaseterm corresponding to a diverging lens, and J_{0}(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_{0}. 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: (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 (3)instead of just the Airy pattern (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.
Fig. 1
Schematic diagram of an external occulter coronagraph with the principal notations used in the paper. The observation, occulter, aperture, and focal telescope planes are shown. 

Open with DEXTER 
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 θ_{0}, 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 (5)where Ψ(r) =  ψ(r)  ^{2} is the wavelengthdependent 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^{2} 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 (6)where θ_{1} and θ_{2} 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 θ_{1} close to θ_{0}. We therefore have γ ≤ Γ, the equality holding in the limit θ_{1} = 0 and θ_{2} = ∞.
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 & RoizenDossier (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 (7)Jacquinot & RoizenDossier (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.
Fig. 2
Left: illustration of the basis function used to represent the attenuation function f(r) with an example of the combination of three of these functions (right). For clarity of representation, the values of Ω_{m} and Δ are not to scale of those used for the f_{k}(r) basis functions. 

Open with DEXTER 
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 piecewise linear as opposed to the use of binary disk that leads to piecewise 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 . 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 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 abovementioned 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 piecewise linear. Copi & Starkman (2000), Cash (2011) and Wasylkiwskyj & Shiri (2011), on the other hand, used offset highorder polynomial or hyperGaussian 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 (8)and the corresponding intensity becomes (9)with with ψ^{⊤} = [ ψ_{1}(r),...,ψ_{K}(r) ] a rank one matrix of general coefficient , 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 (10)where 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 (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 secondorder 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 closedform 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 (12)and the intergrated density (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 θ_{1} = 0,θ_{2} = ∞, where K_{Γ} = K_{γ}. Our optimization in the focal plane was performed for θ_{1} = 0.1 arcsec and θ_{2} = 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 μ = μ_{0}max_{k,l}  K_{k,l}  to be less dependent on the metric matrix K. μ_{0} was chosen by hand to allow good attenuation performance with a limited dynamic in the solution function f. For the monochromatic study, we set μ_{0} = 10^{8} and for the more complex study on a wide bandwidth, less regularization was necessary with μ_{0} = 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^{2}.
Fig. 3
Radial cut of the attenuation of the functions f(r) for aperture and focal optimization at a wavelength of 562 nm. The dashed curves denoted as ≥0 correspond to an optimization with α ≥0 constraint. 

Open with DEXTER 
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 fcurves 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.
Fig. 4
Top curves: intensity in the telescope aperture plane for an optimization that minimizes the integrated intensity Γ in the aperture plane (blue curve) or γ for optimizing the observation between 0.1 and 0.5 arcsec in the focal plane (red curve), for λ = 562 nm. Note the apodization profile for the γ curve and the inversion of relative position of the curves outside the region of optimization. Bottom curves: focal plane intensity obtained for the two optimizations. An airy pattern corresponding to a 10^{12} fainter exoplanet in the observation area is reported in black. The apodization area defined by the smallest and largest radius of the occulter is also reported in green. Note the inversion of behavior of curves with regard to the overall occulter profile. 

Open with DEXTER 
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 centertolimb 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^{5} 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_{γ}.
Fig. 5
Radial cut of the attenuation of the functions f(r) for aperture and focal optimization for the whole spectral bandwidth between 380 nm and 750 nm. 

Open with DEXTER 
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 θ<θ_{1} in the blind area behind the occulter and that Φ_{γ}(θ) < Φ_{Γ}(θ)∀θ ∈ (θ_{1},θ_{2}) in the observation area of the focal plane (illustrated in green in the figure).
Fig. 6
Same representation as in Fig. 4, but for an optimization made for the whole spectral bandwidth between 380 nm and 750 nm. The behavior of the inversion of the curves is similar to the monochromatic case, although it is somewhat reduced by the effect of the chromatism. 

Open with DEXTER 
Resulting light flux in the aperture (Γ) of focal plane (γ) for all the optimization schemes.
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 onaxis 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 twodimensional 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 zeropadding 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 twodimensional array was filled using an interpolation of the onedimensional 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 twodimensional 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 onaxis 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.
Fig. 7
Parametric plot of the observed intensities in the focal plane: xaxis, aperture plane optimization, yaxis, focal plane optimization. The curve is drawn for the wide spectral bandwidth without positivity constraints. From topright to bottomleft the curve follows the increase of θ, clearly showing the concentration of light behind the occulter for the focal plane optimization. The line y = x clearly shows that the γcurves always give better results than the Γcurves, especially in the observation area (reported in green). 

Open with DEXTER 
Threedimensional 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 apodizationlike 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 offset, 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 onaxis and offaxis of 50 cm and 1 m. The two circles of radius θ_{1} = 0.1 arcsec and θ_{2} = 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.
Fig. 8
Top curves: pupil plane intensities for an offset of 1 m for λ ∈ (380, 750) nm. Density curves below represent the focal plane intensities (log scale) for telescope positions: onaxis (upper curves) and offsets of 50 cm (intermediate curves) and 1 m (lower curves). Left: aperture optimization, right: focal plane optimization. 

Open with DEXTER 
Fig. 9
Average residual light level γ for the focal plane a) and pupil plane b) optimization, as a function of the telescope offset in centimeters for λ ∈ (380, 750) nm. The average is computed over the region from 0.1 to 0.5 arcsec delimited in Fig. 8 by two circles. 

Open with DEXTER 
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 onaxis 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 maximummagnitude minimization as in Kasdin et al. (2013) would also be interesting.
See for example http://exoplanet.eu/
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
 Aime, C. 2013, A&A, 558, A138 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Arenberg, J. W., Lo, A. S., Glassman, T. M., & Cash, W. 2007, Comptes Rendus Physique, 8, 438 [NASA ADS] [CrossRef] (In the text)
 Born, M., & Wolf, E. 1999, Principles of optics: electromagnetic theory of propagation, interference and diffraction of light (CUP Archive) (In the text)
 Boyd, S., & Vandenberghe, L. 2004, Convex optimization (Cambridge university press) (In the text)
 Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357 [NASA ADS] [CrossRef] (In the text)
 Cash, W. 2006, Nature, 442, 51 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Cash, W. 2011, ApJ, 738, 76 [NASA ADS] [CrossRef] (In the text)
 Copi, C. J., & Starkman, G. D. 2000, ApJ, 532, 581 [NASA ADS] [CrossRef] (In the text)
 Evans, J. W. 1948, J. Opt. Soc. Am. (1917–1983), 38, 1083 [NASA ADS] [CrossRef] (In the text)
 Jacquinot, P., & RoizenDossier, B. 1964, Prog. Opt., 3, 29 [CrossRef] (In the text)
 Kasdin, N. J., Vanderbei, R. J., Sirbu, D., et al. 2013, in Proc. IEEE Aerospace Conference (In the text)
 Koutchmy, S. 1988, Space Sci. Rev., 47, 95 [NASA ADS] [CrossRef] (In the text)
 Lemoine, D. 1994, J. Chem. Phys., 101, 3936 [NASA ADS] [CrossRef] (In the text)
 Lyot, B. 1939, MNRAS, 99, 580 [NASA ADS] [CrossRef] (In the text)
 Marchal, C. 1985, Acta Astron., 12, 195 [CrossRef] (In the text)
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [NASA ADS] [CrossRef] (In the text)
 Newkirk Jr, G., & Bohlin, J. 1965, in Annales d’Astrophysique, 28, 234 (In the text)
 Purcell, J., & Koomen, M. 1962, J. Opt. Soc. Am, 52, 596 (In the text)
 Shiri, S., & Wasylkiwskyj, W. 2013, J. Opt., 15, 035705 [NASA ADS] [CrossRef] (In the text)
 Soummer, R., Valenti, J., Brown, R. A., et al. 2010, in SPIE Conf. Ser., 7731 (In the text)
 Spitzer, L. 1962, Amer. Sci., 50, 473 (In the text)
 Vanderbei, R. J., & Shanno, D. F. 1999, Comput. Opt. Appl., 13, 231 [CrossRef] (In the text)
 Vanderbei, R. J., Cady, E., & Kasdin, N. J. 2007, ApJ, 665, 794 [NASA ADS] [CrossRef] (In the text)
 Vives, S., Lamy, P., Koutchmy, S., & Arnaud, J. 2009, Adv. Space Res., 43, 1007 [NASA ADS] [CrossRef] (In the text)
 Wasylkiwskyj, W., & Shiri, S. 2011, J. Opt. Soc. Am. A, 28, 1668 [NASA ADS] [CrossRef] (In the text)
Appendix A: Quality test of the numerical computation using an analytical expression for a sharpedged circular occulter
Fig. A.1
From top to bottom: real, imaginary, squared modulus and unwrapped phase of ψ(r) for a circular occulter of radius Ω = 25 m (the vertical line) at z = 80 000 km for λ = 0.55 μm. The continuous line represents the Lommel series, dots are the result of a direct numerical computation of Eq. (1) using NIntegrate of Mathematica. 

Open with DEXTER 
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: (A.1)At r = Ω, both series converge to the same value given by (A.2)
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)
Fig. B.1
Differences of moduli of Fresnel diffractions for raw and linearly apodized functions for Ω values of 10 m (thick curve) and 10 m (thin curve). 

Open with DEXTER 
Fig. B.2
Representation for r between 0 and 2 m of the moduli of ψ_{k}(r) of linearly apodized occulters of radii 10 m and 10.05 m a). The two curves are almost over imposed, and curve b) represents 100 times their differences. 

Open with DEXTER 
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
Resulting light flux in the aperture (Γ) of focal plane (γ) for all the optimization schemes.
All Figures
Fig. 1
Schematic diagram of an external occulter coronagraph with the principal notations used in the paper. The observation, occulter, aperture, and focal telescope planes are shown. 

Open with DEXTER  
In the text 
Fig. 2
Left: illustration of the basis function used to represent the attenuation function f(r) with an example of the combination of three of these functions (right). For clarity of representation, the values of Ω_{m} and Δ are not to scale of those used for the f_{k}(r) basis functions. 

Open with DEXTER  
In the text 
Fig. 3
Radial cut of the attenuation of the functions f(r) for aperture and focal optimization at a wavelength of 562 nm. The dashed curves denoted as ≥0 correspond to an optimization with α ≥0 constraint. 

Open with DEXTER  
In the text 
Fig. 4
Top curves: intensity in the telescope aperture plane for an optimization that minimizes the integrated intensity Γ in the aperture plane (blue curve) or γ for optimizing the observation between 0.1 and 0.5 arcsec in the focal plane (red curve), for λ = 562 nm. Note the apodization profile for the γ curve and the inversion of relative position of the curves outside the region of optimization. Bottom curves: focal plane intensity obtained for the two optimizations. An airy pattern corresponding to a 10^{12} fainter exoplanet in the observation area is reported in black. The apodization area defined by the smallest and largest radius of the occulter is also reported in green. Note the inversion of behavior of curves with regard to the overall occulter profile. 

Open with DEXTER  
In the text 
Fig. 5
Radial cut of the attenuation of the functions f(r) for aperture and focal optimization for the whole spectral bandwidth between 380 nm and 750 nm. 

Open with DEXTER  
In the text 
Fig. 6
Same representation as in Fig. 4, but for an optimization made for the whole spectral bandwidth between 380 nm and 750 nm. The behavior of the inversion of the curves is similar to the monochromatic case, although it is somewhat reduced by the effect of the chromatism. 

Open with DEXTER  
In the text 
Fig. 7
Parametric plot of the observed intensities in the focal plane: xaxis, aperture plane optimization, yaxis, focal plane optimization. The curve is drawn for the wide spectral bandwidth without positivity constraints. From topright to bottomleft the curve follows the increase of θ, clearly showing the concentration of light behind the occulter for the focal plane optimization. The line y = x clearly shows that the γcurves always give better results than the Γcurves, especially in the observation area (reported in green). 

Open with DEXTER  
In the text 
Fig. 8
Top curves: pupil plane intensities for an offset of 1 m for λ ∈ (380, 750) nm. Density curves below represent the focal plane intensities (log scale) for telescope positions: onaxis (upper curves) and offsets of 50 cm (intermediate curves) and 1 m (lower curves). Left: aperture optimization, right: focal plane optimization. 

Open with DEXTER  
In the text 
Fig. 9
Average residual light level γ for the focal plane a) and pupil plane b) optimization, as a function of the telescope offset in centimeters for λ ∈ (380, 750) nm. The average is computed over the region from 0.1 to 0.5 arcsec delimited in Fig. 8 by two circles. 

Open with DEXTER  
In the text 
Fig. A.1
From top to bottom: real, imaginary, squared modulus and unwrapped phase of ψ(r) for a circular occulter of radius Ω = 25 m (the vertical line) at z = 80 000 km for λ = 0.55 μm. The continuous line represents the Lommel series, dots are the result of a direct numerical computation of Eq. (1) using NIntegrate of Mathematica. 

Open with DEXTER  
In the text 
Fig. B.1
Differences of moduli of Fresnel diffractions for raw and linearly apodized functions for Ω values of 10 m (thick curve) and 10 m (thin curve). 

Open with DEXTER  
In the text 
Fig. B.2
Representation for r between 0 and 2 m of the moduli of ψ_{k}(r) of linearly apodized occulters of radii 10 m and 10.05 m a). The two curves are almost over imposed, and curve b) represents 100 times their differences. 

Open with DEXTER  
In the text 