© C. Aime 2020

1. Introduction

Coronagraphy began with the experiments of Lyot (1932). In addition to the study of the solar corona, this technique has become essential for the detection of exoplanets where many theoretical developments have been made since the first detection of an exoplanet by Mayor & Queloz (1995). The goal of these high-contrast techniques is to allow the imaging of a very weak object (solar corona or exoplanet) next to a very bright object (solar disk or star) by getting rid of as much of the stray light diffracted by the bright source at the level of the science object as possible.

The brightness ratio between the object and the neighboring bright source is comparable in the two astronomical domains. The light in the visible domain from an exo-Earth is of the order of 10¹⁰ fainter than that from its star, and, in order to consider exoplanets of the solar system type, an exo-Jupiter is about 10⁶ fainter in the far infrared. The K-corona is a few times 10⁻⁶ of the mean solar brightness close to the solar limb, and decreases to almost 10⁻¹⁰ at an altitude of 5 solar radii (Cox 2015). Moreover, the solar corona is an extremely complex physical environment (Aschwanden 2006) with rapidly evolving structures (Peter et al. 2013).

Requirements in resolution are different, with the detection of exoplanets requiring very large telescopes while observation of the solar corona with apertures of only a few centimeters can already give interesting results. However this is only a difference in scale, and what makes the fundamental difference between the two research fields for instrumental concepts is that the solar disk is a huge extended object while stars are almost unresolved point sources.

This considerably limits the use, for solar observations, of most of the very advanced coronagraphs invented for the detection of exoplanets. Phase-mask coronagraphs, such as those of Roddier & Roddier (1997), Soummer et al. (2003), and Mawet et al. (2005), are inappropriate for observation of the Sun. Indeed these systems, such as for example the four quadrants of Rouan et al. (2000), that can fully cancel the light from a point source on-axis, become totally ineffective for a point source off-axis, and therefore cannot be used for an extended object.

Stray light analysis in solar coronagraphs has been the subject of numerous studies and instrumental developments for more than half a century (Purcell & Koomen 1962; Newkirk & Bohlin 1963, 1965; Tousey 1965). Coupling a Lyot coronagraph with an external occulter goes back to the photometer of Evans (1948), a technique that is found in almost all space-born coronagraphs (Koutchmy 1988), even though the future project VELC Aditya-L1, Prasad et al. (2017) is a pure 15 cm Lyot coronagraph.

The new ESA formation flying Proba-3 ASPIICS project, Renotte et al. (2014), remains in the tradition of an externally occulted Lyot coronagraph. It will consist of two spacecraft flying 144 m apart, the occulter spacecraft, with a 1.42 m disk and the coronagraph spacecraft bearing a Lyot coronagraph with an entrance aperture of 5 cm. The umbra and penumbra solar pattern produced by the external occulter at the level of the coronagraph can be expressed as a convolution, as shown in Aime (2013). Propagating light down to the observing plane requires a much more complicated analytical study (Ferrari 2007; Ferrari et al. 2009).

The current approach involves the coupling of analytic and numerical computation. The procedure consists of two steps. First the response produced by each source point of the solar photosphere in the last observation plane is calculated in a coherent propagation. Rougeot et al. (2017) clearly show differences in response between points on and off axis (their Figs. 3 and 4). The intensities of these point sources coming from the entire solar surface are then summed, because they are mutually incoherent. A series of recent papers made use of similar procedures (Rougeot et al. 2018; Shestov & Zhukov 2018; Shestov et al. 2019). Because the goal is to observe the corona very close to the limb, the effects of scattered light are far from negligible. In particular, Rougeot et al. (2019) showed that the quality of the optics is crucial for the Lyot coronagraph of ASPIICS, and that the rejection is mainly due to the effect of the external occulter.

Optimization of the shape or the structure of a simple or multiple external occulters is a long-standing goal of solar coronography in space (Purcell & Koomen 1962; Tousey 1965; Fort et al. 1978; Lenskii 1981) and remains active today (Thernisien et al. 2005; Verroi et al. 2008).

Tousey (1965) used a serrated-edge external occulter and obtained the first observations of the outer solar corona in the absence of an eclipse. These shaped occulters, unlike petal occulters for exoplanet studies (Cash 2006; Kasdin et al. 2009) did not have the expected success for solar astronomy. The reason sometimes invoked is the poor quality of their realization, but it might instead be the fact that the number of teeth used was not large enough, as shown by Rougeot & Aime (2018). Indeed, for the case of ASPIICS, a 64-tooth serrated occulter was found to give poorer results than a sharp-edged disk, while occulters with one thousand teeth are expected to give excellent results.

Multiple disks were considered as the optimal external occulting system, Purcell & Koomen (1962). Newkirk & Bohlin (1963) wrote that “the maximum reduction occurs just when the second and smaller disk completely obscures the objective lens from the edge of the first disk”. More recently, Thernisien et al. (2005) presented experimental and numerical optimization of an external occulter made of three disks on axis. Their numerical simulation is based on Fresnel filtering techniques.

In the present study, our approach is mainly analytic, and numerical results are obtained computing Bessel-Hankel integrals in Sect. 3. This is because of the difficulty of a purely numerical calculation at very large Fresnel numbers corresponding to solar experiments. The question of sampling was analyzed in detail in Rougeot & Aime (2018), with illustrations in their Appendix A. The problem worsens in the case of double or multiple diffractions. The analytical approach is also more general, and allows a more detailed analysis; for example, the use of the Arago spot in Sect. 4. Nevertheless, at very low Fresnel numbers, a 2D numerical illustration is possible and useful for a simple illustration of the Fresnel diffraction pattern produced by two disks, and we use it in an example.

A decisive experimental breakthrough was made with the use of a three-dimensional occulter for the coronagraph LASCO-C2 of the SOHO Mission, Domingo et al. (1995). In their experimental study, Bout et al. (2000) found that three-dimensional conical external occulters gave performances superior to triple- disk occulters, a conclusion previously reached by Koutchmy & Belmahdi (1987).

However, a numerical study of the Fresnel diffraction of a three-dimensional occulter is very difficult, because it is no longer in the field of classical Fourier optics, which must assume plane-to-plane propagations. The alternative may be the approach of Baccani et al. (2016) and Landini et al. (2017), which consisted of substituting planar cuts of the occulter to its true 3D shape, and computing the successive Fresnel diffractions numerically from one cut to another. This constitutes a first possible approach, and we adopt it here to derive an optimal shape for a 3D occulter.

The paper is organized as follows. Section 1 contains this introduction. In Sect. 2 we give the formal expressions for the Fresnel diffraction for two disks on axis, using a Babinet approach. We show the interest of the Arago spot to quantify the darkness of the central shadowed zone. Section 3 is devoted to a parametric analysis of the spot of Arago depending on the configuration in size and position of the two disks on axis. A simple law of optimization is obtained from this study. In Sect. 4 we use this latter to derive the shape of the envelope of multiple occulters that optimally reduce the level of the Arago spot. Finally, a discussion is provided in Sect. 5. An appendix gives the analytic expression for three disks.

2. Fresnel diffraction of one and two disks on axis

As already indicated above, Fresnel diffraction is an effect of the coherent nature of the light. Here we reiterate a few basic relations of Fourier optics.

Let Ψ₀(r) be the complex amplitude of the wave in the plane z = 0. Its expression Ψ_z(r) after a free-space propagation over a distance z may be written as:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_z(r)& ={\mathrm{\Psi }}_0(r)\ast {\tau }_z(r),\mathrm{where}\hfill \\ \hfill {\tau }_z(r)& =\frac{1}{\mathrm{i}\lambda z}exp\left(\frac{\mathrm{i}\pi r^2}{\lambda z}\right),\mathrm{and}{\tau }_z(r)\ast 1=1,\hfill \end{array}$$(1)

where i is the imaginary unit, λ the wavelength of the light, and the asterisk represents two-dimensional convolution, to be performed on the components x and y of r(x, y). A constant term exp(2iπz/λ) has been omitted here. Taking the direct and inverse Fourier transform of Eq. (1), we obtain:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_z(r)={\mathcal{F}}^{-1}[{\widehat{\mathrm{\Psi }}}_0(u)exp(-\mathrm{i}\pi \lambda z{u}^2)],\end{array}$$(2)

where $\widehat{\mathrm{\Psi }}$ stands for the direct Fourier transform ℱ[Ψ] of Ψ, and ℱ⁻¹[…] the inverse Fourier transform of the quantity inside the brackets. This approach is commonly used in direct numerical computations using Fast Fourier Transforms.

2.1. Fresnel diffraction of a single disk

Let us consider an incident wave of unit amplitude impinging on an opaque disk of diameter Ω centered on axis. The transmission of this occulter can be written as 1 − Π(r/Ω), where Π(r) is the window function equal to 1 for |r| < 1/2 and 0 elsewhere. According to Eq. (1), the wave at the distance z can be written as:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_z(r)=[1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)]\ast {\tau }_z(r)=1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_z(r)\end{array}$$(3)

where we use the fact that τ_z(r)*1 = 1. This equation shows that the Fresnel diffraction of a disk can be expressed as 1 minus the diffraction of the complementary hole occulter, which is an illustration of Babinet theorem for Fresnel diffraction.

Moreover, Ψ_z(r) is a radial function and can be classically computed as an Hankel integral. We have, for the convolution term in Eq. (3),

$$\begin{array}{cc}\hfill \mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_z(r)=\frac{2\pi }{\mathrm{i}\lambda z}& exp\left(\frac{\mathrm{i}\pi r^2}{\lambda z}\right){\int }_0^{\mathrm{\Omega }/2}\xi exp\left(\frac{\mathrm{i}\pi {\xi }^2}{\lambda z}\right)J_0\left(\frac{2\pi \xi r}{\lambda z}\right)\mathrm{d}\xi .\hfill \end{array}$$(4)

One interest of the Babinet approach is that the integral is to be calculated from zero to the value of the radius of the occulter, instead of being over an infinite support (from the radius to infinity). The numerical computation of this integral is delicate for large Fresnel numbers N_F = Ω²/λz. For ASPIICS, N_F ∼ 25 000, meaning that the complex exponential quadratic term makes about 8000 rounds in the complex plane. The Lommel series can be used to compute Eq. (4), as well described in the reference book of Born & Wolf (1999). Aime (2013) showed that the present numerical computation of integrals, for example the function NIntegrate of Mathematica (2009), can give excellent results. The value for r = 0 of Eq. (3) gives the amplitude of the famous Arago spot:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_z(0)=exp\left(\frac{i\pi {\mathrm{\Omega }}^2}{4\lambda z}\right)·\end{array}$$(5)

The corresponding intensity remains that of the incident wave. This calculation, first made by Poisson, is at the origin of the confirmation of Fresnel theory. Many publications can be found on this Arago spot; see for example Harvey & Forgham (1984) and Kelly et al. (2009).

2.2. Generalization to the Fresnel diffraction of two disks

Let us now consider the system of two disks on axis schematized in Fig. 1, of diameters Ω and ω, and set at distances z and d from the observing plane, respectively. From a geometrical point of view, the blue triangle corresponds to the region where the second occulter is in the shadow of the first one and masks it from the on-axis point A. The blue triangle corresponds to the geometrical condition:

$$\begin{array}{c}\hfill \frac{d}{z}\times \mathrm{\Omega }<\omega <\mathrm{\Omega },\end{array}$$(6)

Fig. 1. Schematic diagram for the Fresnel diffraction of two disks of diameters Ω and ω at distances of z and d from observing plane E. AB is the radius of the geometrical dark zone.

in agreement with the finding of Tousey (1965) and Newkirk & Bohlin (1963). Point B is the position in the field beyond which the first occulter is no longer shielded by the second one. A simple calculus gives:

$$\begin{array}{c}\hfill |\mathrm{AB}|=\frac{\omega z-\mathrm{\Omega }d}{2(z-d)}·\end{array}$$(7)

AB gives the geometrical radius of a central dark circular zone and presents similarities with the Boivin (1978) radius for serrated occulters, although these two dark zones have completely different physical origins. We see below that there are conditions for this zone to really be a dark zone.

Let us come back to diffraction effects. The wave in the plane ${\text{O}}_2^{-}$ immediately before the second disk, at the distance z − d, can be written as:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{{\mathrm{O}}_2^{-}}(r)=1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{z-d}(r).\end{array}$$(8)

The wave in plane ${\text{O}}_2^{+}$, immediately after the second occulter, is:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{{\mathrm{O}}_2^{+}}(r)=[1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{z-d}(r)]\times [1-\mathrm{\Pi }\left(\frac{r}{\omega }\right)]\end{array}$$(9)

and the wave in plane E can be written as:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_E(r)={\mathrm{\Psi }}_{{\mathrm{O}}_2^{+}}(r)\ast {\tau }_{\mathrm{d}}(r).\end{array}$$(10)

Since τ_z − d(r)*τ_d(r) = τ_z(r), Eq. (10) simplifies to:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_E(r)=& 1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_z(r)-\mathrm{\Pi }\left(\frac{r}{\omega }\right)\ast {\tau }_{\mathrm{d}}(r)\hfill \\ \hfill & +\{[\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{z-d}(r)]\times \mathrm{\Pi }\left(\frac{r}{\omega }\right)\}\ast {\tau }_{\mathrm{d}}(r).\hfill \end{array}$$(11)

As for Babinet’s theorem, the Fresnel diffraction of two disks can be expressed using only complementary hole occulters. The interest of such a formulation for a numerical calculation is that the integrals are limited to the radii of the disks.

A slight reordering of this equation allows for a more physical presentation of the result. Assembling the last two terms of Eq. (11), we obtain:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_E(r)=& 1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_z(r)\hfill \\ \hfill & -\{[1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{z-d}(r)]\times \mathrm{\Pi }\left(\frac{r}{\omega }\right)\}\ast {\tau }_{\mathrm{d}}(r).\hfill \end{array}$$(12)

In that form, the complex amplitude Ψ_E(r) can be expressed as the coherent composition, or interference, of two terms, one due to the disk Ω (propagation over z, defined by the convolution τ_z(r)) and the second due to the disk ω (propagation over d, the convolution τ_d(r)), the amplitude of the disk ω being modified by the Fresnel diffraction of Ω onto it, the term inside square parentheses.

The calculation for two disks leading to Eq. (11) can be generalized to three disks centered on axis, in a system similar to the one described by Newkirk & Bohlin (1963) or Thernisien et al. (2005). The result is given in Appendix A. Unfortunately there is not a simple law to express the diffraction of three disks from that of two, and no recurrent relationship can be established from N disks towards N + 1 disks, and therefore here we limit our analysis to Eq. (11) and the Fresnel diffraction of two disks.

The two first convolution terms of Eq. (11) correspond to the Fresnel diffractions of holes of diameters Ω and ω, and can be easily computed using Mathematica, as already indicated. The last term corresponds to the diffraction of the wave that would travel through the two holes Ω and ω, and its computation is more difficult. Indeed, applying twice the property that the Fourier transform of a 2D circular radial function is a Hankel transform, we obtain:

$$\begin{array}{cc}& \{[\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{z-d}(r)]\times \mathrm{\Pi }\left(\frac{r}{\omega }\right)\}\ast {\tau }_{\mathrm{d}}(r)\hfill \\ \hfill & =\frac{-4{\pi }^2}{{\lambda }^{2}d(z-d)}\hfill \\ \hfill & \times exp\left(\frac{\mathrm{i}\pi r^2}{\lambda d}\right){\displaystyle {\int }_{\eta =0}^{\omega /2}{\int }_{\xi =0}^{\mathrm{\Omega }/2}}\xi \eta exp\left(\frac{\mathrm{i}\pi {\eta }^{2}z}{\lambda d(z-d)}\right)exp\left(\frac{\mathrm{i}\pi {\xi }^2}{\lambda (z-d)}\right)\hfill \\ \hfill & \times J_0\left(\frac{2\pi \xi \eta }{\lambda (z-d)}\right)J_0\left(\frac{2\pi r\eta }{\lambda d}\right)\mathrm{d}\xi \mathrm{d}\eta .\hfill \end{array}$$(13)

This term is difficult to evaluate because of the double integral and the very fast oscillating complex exponential terms for small values of (z − d) compared with disk diameters. Indeed, for parameters similar to ASPIICS, a propagation of 7 centimeters (the width of the external toroidal occulter), the Fresnel numbers N_F would be larger than 50 millions.

2.3. Numerical illustrations

Figures 2–4 give illustrations of Fresnel diffractions that can be observed in plane (E) of Fig. 1. A global 2D numerical illustration of the whole pattern is given in Fig. 2 and radials cuts in Figs. 3 and 4 for different configurations. Detailed explanations follow.

Fig. 2. Fresnel diffraction patterns of single and double disks. Parameters: Ω = 4 mm, ω = 3 mm, and z − d = d = 0.1 m. Top panel: diffraction of disks Ω (top left quadrant A) and ω (bottom left quadrant B) alone. The right two quadrants (C) show the diffraction of the system of two disks. Vertical and horizontal lines correspond to ±Ω/2 (external black lines), ±ω/2 (blue dashed lines) and the geometrical limit of the dark zone (red lines) AB of Fig. 1 and Eq. (7). Bottom panel: zoom on the central region around the Arago spot. In this example, the gain in darkness is about 8.

Fig. 3. Red curves: Fresnel diffraction for two disks using Eq. (11) for Ω = 3 mm, ω = 2.25 mm, z = 1 m, and d = 0.5 m. The black curve shows the diffraction of the disk Ω alone (Eq. (3)). Top panel: linear scale, bottom: log10 scale. Vertical lines correspond to Ω/2 (black small dashed line), ω/2 (blue dashed line) and the geometrical limit AB (at 0.75 mm, (red line)) of Fig. 1 and Eq. (7).

Fig. 4. Fresnel diffraction with two disks for a central region close to the Arago spot, for Ω = 1.42 m, ω = 1.40524 m, z = 144.3 m, and z − d = 3 m. The top panel is in linear scale; the red dots correspond to Eq. (11) and the continuous curve to the Fresnel diffraction of the disk ω alone attenuated by a factor of 321. Bottom panel: two curves in log10 scales (black, disk Ω alone, red, the system of two disks).

Figure 2 shows the Fresnel diffraction of two disks on-axis compared with diffractions of disks taken alone. The parameters are Ω = 4 mm, ω = 3 mm, z = 0.2 m and d = 0.1 m. This system corresponds to a Fresnel numbers of 290 for the diffraction of the first disk onto the second. A dark central zone of radius AB of Fig. 1 and Eq. (7) is clearly visible in top Fig. 2. For this example, a fully numerically computation was made using arrays of 32 768 × 32 768 points. It is the only numerical calculation in this paper. The programming language used is MATLAB and the computer is a machine with two Intel Xeon(R) @ 2.3 GHz × 56 cores and 503.8 G of RAM operated under Ubuntu.

The transmission of the disks in the computation is either 0 or 1. Disk Ω has a diameter of 40% of the array. The Fresnel diffraction of the wave after the disk Ω is calculated using Eq. (2) for the distance z − d = 0.1m. Then the transmission of the second disk ω is applied to the wave, and the final result is obtained computing the propagation over a same distance d = 0.1 m using again Eq. (2). The representation is divided in 3 parts giving (A) the diffraction of the disk Ω alone, (B) that of the disk ω alone, and (C) the diffraction for the two disks on-axis. A dark central zone appears in the right part (C) of the top figure of Fig. 2. This circular zone is well inside the red lines giving the geometrical limits AB of Fig. 1 and Eq. (7). For the outer part of the diffraction, near the external black lines corresponding to Ω, the influence of the second disk ω becomes only marginal, the diffraction being dominated by the first disk. A zoom of the central part of the diffraction pattern is shown in Fig. 2, bottom. There, part (C) appears as a darkened version (about 8 times fainter) of the pattern due to the second disk ω, drawn in the lower left part (B) of the figure. Diffraction rings in (C) are noisier than in (B). This is due to numerical computational problems, mainly because of the limited number of points in arrays which limits both the field and the sampling. The noise increases also because of the two successive propagations for (C), instead of just one for (A) and (B).

Figure 3 gives an example of a radial cut |Ψ_E(r)|² using the analytic expression of Eq. (11). Here the Bessel-Hankel integrals are computed using the function NIntegrate, that is a general numerical integrator of Mathematica. The same Xeon based computer is used as for getting Fig. 2 with MATLAB but the Mathematica license we had did not allowed the machine to run in parallel on all cores (only 16 out of 56 available). Results are much more precise than the direct numerical calculation using Fourier transforms. Parameters for this system are Ω = 3 mm, ω = 2.25 mm, z = 1 m and d = 0.5 m. The diffraction of the first disk onto the second corresponds to a low Fresnel number of 33. Geometrical limits are given in the figure. This representation clearly illustrates the Fresnel diffraction in the external region where the distance to the center is greater than the radius AB of the dark zone given in Eq. (7). The result of the Fresnel diffraction for the two disks is compared with the Fresnel diffraction that could be obtained if the disk Ω was alone (Eq. (3)). Outside the central dark zone, the Fresnel diffraction of the first disk Ω tends to dominate the pattern progressively. The ratio of the two Arago spots is 7.4 here, but this is not a focal point of this illustration.

Figure 4 gives |Ψ_E(r)|² for the central region very close to the Arago spot for Ω = 1.42 m, ω = 1.40524 m, z = 144.3 m and z − d = 3 m. These parameters correspond to a diffraction at very high Fresnel number (1.2 million), and took several days of computing time with Mathematica. The top figure of Fig. 4 shows in red dots the central regions of Fresnel diffraction for the 2 disks Ω × ω, using Eq. (11). The Fresnel diffraction of the single disk ω at the distance d, computed with Eq. (3), is drawn as a continuous black line. To make curves superimposable, the single disk diffraction pattern was divided by a constant factor (here 321). After that, the two curves are remarkably similar. Figure 4, bottom figure, shows the two raw diffraction patterns in a log10 scale and the constant factor appears there as an offset between curves.

The effect of stray light reduction expected by the use of two disks is well obtained. No wonder with that because it is the idea developed by the first inventors of the multiple disks system. However this result comes out here from the theoretical expression of Eq. (11). Reasoning in terms of the Maggi-Rubinowich approach, well described in Born & Wolf (1999) and used by Cady (2012) and Rougeot & Aime (2018) for example, provides a physical foundation for the former empirical approach. It states that, for an opaque occulter, the diffraction comes mainly from the edges. Positioning the second occulter ω in the umbra of Ω has indeed the effect to make darker its edges, and the shadow obtained resembles the diffraction of the disk ω illuminated by a less bright source.

We note that, for the conditions of Fig. 4, the central part of the Fresnel diffraction of ω can be very well approximated by the function $I_0\times J_0^2(\pi \omega r/\lambda d)$, where J₀ is the zero-order Bessel function of the first kind. This expression was first obtained by Harvey & Forgham (1984) and retrieved by Aime (2013), as the first term of the Lommel series giving the Fresnel diffraction of a disk, as described in section 8.8 “The three-dimensional light distribution near focus” of Born & Wolf (1999). In that formula, since J₀(0) = 1, the intensity at r = 0 is the value I₀. The Arago spot gives the value of the intensity of the whole central pattern. In Fig. 4 this approximation exactly matches values for the single disk ω computed with Eq. (3) and have not been plotted there.

These numerical results show that the intensity of the spot of Arago is a convenient measure of the level of the dark zone that can be obtained using two disks. Of course this level depends on the disks configuration (size and distance), which is the subject of study in the following section.

3. Spot of Arago for two disks on axis

The integral expression of Eq. (11) somewhat simplifies for the value of the central point Ψ_A(0), that is, the Arago spot.

Because J₀(0) = 1, we have:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{\mathrm{A}}(0)=& exp\left(\mathrm{i}\pi \frac{{\mathrm{\Omega }}^2}{4\lambda z}\right)+exp\left(\mathrm{i}\pi \frac{{\omega }^2}{4\lambda d}\right)-1\hfill \\ \hfill & -\frac{4{\pi }^2}{{\lambda }^{2}d(z-d)}{\displaystyle {\int }_{\eta =0}^{\omega /2}{\int }_{\xi =0}^{\mathrm{\Omega }/2}}\xi \eta exp\left(\frac{\mathrm{i}\pi {\eta }^{2}z}{\lambda d(z-d)}\right)\hfill \\ \hfill & \times exp\left(\frac{\mathrm{i}\pi {\xi }^2}{\lambda (z-d)}\right)J_0\left(\frac{2\pi \xi \eta }{\lambda (z-d)}\right)\mathrm{d}\xi \mathrm{d}\eta \hfill \end{array}$$(14)

which remains a difficult integration because it still contains the very fast oscillating quadratic complex term.

An example of the result of the computation of Eq. (14) for z = 82 cm and Ω = 32 mm is given as a 2D contour plot in Fig. 5. The second occulter ω varies between 22 and 36 mm, and z − d varies between 0.8 and 16.4 cm. The figure highlights a triangular region where the position and size of the second occulter can efficiently reduce the level of the Arago spot. This region corresponds to the blue triangle in Fig. 1 and closely matches the findings of the early inventors of the multiple disk systems in the 1960’s. The log10 scale allows better visualization of the variation of the gain. It would be very interesting to verify these results experimentally, however this is a difficult task because of the number of disks of different sizes required.

Fig. 5. Intensity of the Arago spot computed using Eq. (14). Parameters: z = 82 cm, Ω = 32 mm, λ = 0.55 μm. Linear (top) and log10 (bottom) representations are shown. The triangular region corresponds to the one in Fig. 1 and conditions given in Eq. (6).

The results of a much easier experiment are simulated in Fig. 6. In this experiment the disks Ω and ω are of fixed size (in the example we took Ω = 44 mm and ω = 40 mm). We assume that the disk Ω is fixed and that the disk ω may be moved away so as to vary the distance z − d, using a simple translation mount. The Arago spot is observed at z = 1 m and its intensity is given as a function of z − d. The dashed vertical line corresponds to the position where the second disk enters the blue triangle. For distances z − d that are not large enough, the Ω × ω system, instead of making the Arago spot darker, actually makes it brighter by up to a factor of 1.33. This is a surprising result that comes from the equations, but can be understood in a heuristic way. As highlighted by Eq. (12), the resulting Fresnel diffraction can be seen as the coherent composition of the complex amplitudes coming from the disks Ω and ω, the amplitude of the latter term being modified by the Fresnel diffraction of Ω onto it. Depending on whether the amplitudes arrive in phase or in opposition, this produces constructive or destructive interference at the origin of the oscillatory behavior of the curve. The unwanted effect becomes important when the second disk approaches the blue triangle of Fig. 1 while remaining outside.

Fig. 6. Intensity of the Arago spot observed at z = 1 m for a system of two disks of dimensions Ω = 4.4 cm and ω = 4 cm, as a function of the distance z − d between the first and second disk. The vertical red dashed line corresponds to the distance when the second disk enters the blue triangle of Fig. 1.

Figure 7 corresponds to an experiment which takes the key components of ASPIICS, that is, two satellites 144.3 m apart, but where the current occulting spacecraft is replaced by a system bearing two sharp-edged disks, one at each end of the satellite; the first disk Ω being oriented towards the sun and the second disk ω towards the second spacecraft. In the example of Fig. 7, the distance between the disks is 50 cm. The curve gives the intensity of the spot of Arago as a function of the diameter of this second disk. The two vertical red dashed lines correspond to ω = dΩ/z and ω = Ω, the geometrical limits given in Eq. (6). Here again, if the values are not correctly chosen, the system of two disks may increase the value of the Arago spot, which can become brighter than with a single disk. Oscillatory increase of the Arago spot is obtained near the edges of the blue triangle in Fig. 1. For the lower limit, the heuristic interpretation has already been given for Fig. 6, that is, the combination of waves coming from the edges of Ω and ω. For the other higher limit, one can consider that the first disk, Ω, diffracts light towards the second disk, ω, which by interference with the direct light may increase or decrease its illumination.

Fig. 7. Intensity of the Arago spot as a function of the size of the second disk ω, in a linear (top) and log10 scale (bottom). The first disk of diameter Ω = 1.42 m is at 144.3 m from the telescope, the second disk ω is at d − z = 50 cm from the first one. Vertical red dashed lines correspond to limits given in Eq. (6), and the blue dashed line corresponds to the optimal value of Eq. (15). There is an increase in brightness outside these limits.

The best darkening effect appears to be at the mean distance between these two limits, that is, for an optimal ω_opt value:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{opt}}=\mathrm{\Omega }\frac{z+d}{2z}·\end{array}$$(15)

This position is indicated by a blue dashed line in Fig. 7. The curve is rather flat, and provided that the size of the second disk does not differ too greatly from the optimal value ω_opt, the rejection factor remains very good.

The rejection factor, which we define as the inverse of the value of the Arago spot at the optimal value ω_opt, increases almost linearly with the distance z − d between disks, as shown in Fig. 8, for Ω = 1.42 m and z = 144.3 m. The linear fit is 2.156 + 106.3 (z − d). This expression gives the correct value of 321 found in Fig. 4 for z − d = 3 m.

Fig. 8. Rejection factor as a function of the distance z − d for Ω = 1.42 m and z = 144.3 m. The straight line corresponds to a linear fit. The corresponding variation of ωopt is from 1.418 m to 1.413 m.

4. The optimal envelope for multiple disks

In this section we seek to derive an optimal envelope for an ensemble of disks. To this aim, we follow the approach of Landini et al. (2017), which is to replace the 3D structure by a series of 2D cuts. We assume that this remains valid within the approach of Fourier optics, an assumption that would require further investigation. Let us consider the schematic representation of Fig. 9. Here, z = z₀ is the distance from the first disk occulter in O to the point of observation on axis A, and HP is the Nth disk at the distance x from O. According to the condition given in Eq. (15), the top P’ of the Nth + 1 disk (such as H’P’) should be in the direction PM in the same blue triangle of Fig. 1. The point M should be preferably the middle of AA’, according to ω_opt in Eq. (15), but for more generality, let us assume that we choose a direction PQ such as AM = αMA’. With HP = y, HA = z₀ − x, simple geometric considerations give HQ = (1 + α)(z₀ − x), and we can write:

$$\begin{array}{c}\hfill -\frac{\mathrm{HP}}{\mathrm{HQ}}=y^{\prime }=-\frac{y}{(1+\alpha )(z_0-x)}·\end{array}$$(16)

Fig. 9. Geometric elements for the derivation of the differential equation of Eq. (16) making use of the condition given in Eq. (15). The blue triangle has the same meaning as in Fig. 1.

With y(0) = Ω/2, the differential equation gives the simple solution:

$$\begin{array}{c}\hfill y(x)=\frac{\mathrm{\Omega }}{2}\times {(1-\frac{x}{z_0})}^{\frac{-1}{1+\alpha }}\end{array}$$(17)

where y(x) is the generatrix of the envelope around the rotation axis x. If we can assume that x ≪ z₀, then the solution is a straight line:

$$\begin{array}{c}\hfill y(x)\sim \frac{\mathrm{\Omega }}{2}(1-\frac{x}{(1+\alpha )z_0})\end{array}$$(18)

that generates a conical envelope. An illustration is given in Fig. 10 for α  =  1. The exact solution is compared with its linear approximation, corresponding to a cone as in LASCO C2, Bout et al. (2000). The value α  =  1 produces a rather small dark zone, and a larger value of α may be preferably used to increase it, unless observation of only the outer corona is the goal of the experiment. Importantly, only the last disk is seen from the aperture.

Fig. 10. Illustration of the optimal envelope for the solution given in Eq. (17) and α = 1. The exact curve (1) and linear approximation (2) of Eq. (18) are shown as the blue and red dashed lines, respectively. The drawing is made for z0 = 1 m, and the multiple disks extending up to 0.5 m. These values are for illustration only.

5. Conclusion

The present work is an analytical study of the diffraction produced by a nonsingle-plane occulter. In particular the Fresnel diffraction produced by an incident wave plane impinging on two disks on axis is described analytically. A few numerical examples are given. The problem encountered is that evaluation of the Bessel-Hankel integrals expressing the complex amplitude of the wave is not straightforward. Furthermore, these integrals take a lot of time to be obtained, even with a powerful computer. An optimization of the calculation of these integrals should be done to make the calculation faster if possible, but this is outside the scope of the present work.

Here, we show that the Arago spot, that is, a single central point, is an excellent indicator of the level of darkness obtained by such a system. This allows a simple parametric study of the level of darkening as a function of the geometry of the system (dimensions of the disks, separation between them, distance to the observation plane).

Our results match those intuitively obtained by the early inventors of the multiple disk system in the 1960’s. The difference is that the results shown here are from a purely numerical calculation based on a theoretical expression obtained using basic theoretical Fourier optics. Moreover, quantitative results are obtained for the rejection performance, also highlighting unwanted light amplification for misconfigured disks.

An analytical expression was also obtained for three disks on axis, but there was no numerical calculation because this expression is really too difficult to evaluate. Nevertheless, using simplifying assumptions, the envelope for a multiple occulting system is obtained, but without giving the value of the rejection.

The most perplexing result of this work may be the finding of disk configurations that have the opposite effect to what is expected, that is, a brighter zone instead of a dark zone. The torus-shaped external occulter chosen for ASPIICS may fall into these unwanted configurations, as can be deduced from Fig. 11 of Baccani et al. (2016). From the telescope aperture, a non-null part of the edge of the mask will be seen from the telescope aperture, and not only the last disk as in Fig. 10. The present work cannot precisely quantify the effect this may induce, but we cannot exclude an amplification of the brightness of the central area at the pupil of the Lyot coronagraph.

Another point of interest is that this phenomenon of light amplification exists also when the second disk is larger than the first one. As a consequence, it cannot be excluded that the satellite itself may produce inappropriate enhancement of the Fresnel diffraction pattern of its occulting disk. This is very difficult to estimate numerically. Laboratory experiments would be very useful, starting with the one described in Fig. 6. To take into account the effects of the full 3D structure, a small-scale model could be realized and its Fresnel diffraction studied in a montage that conserves the optical etendue, as described by Landini et al. (2017).

A very efficient setting would have been to include the occulter spacecraft of ASPIICS inside a system of two disks Ω × ω, one at each end of the satellite. For example the photovoltaic sensors might have been included in a larger first disk acting as Ω, the second disk ω being set in the place of the actual toroidal occulter. The gain in rejection factor, for a distance between disks of 1 m (the width of the occulter spacecraft is 90 cm) would have been of the order of 100, as shown in Fig. 8. Whether this would have been feasible or not is a technical question that goes beyond the objectives of this work.

As already indicated, this paper gives results for a single incident plane wave. The calculation of umbra and penumbra patterns on the pupil, and the residual stray light in the final focal plane by the incoherent addition of all point sources of the whole Sun remains to be investigated.

This could be done for the two-disk occulter proposal envisaged above, for example, but the calculation of the diffraction using Eq. (11) would take an enormous amount of time. The speed of calculation of the Bessel Hankel integrals should be increased. Indeed, this could be done using a faster computer, but another solution would be to obtain a gain in speed by optimizing the parameters of the NIntegrate function, or by replacing some parts of the results with acceptable approximations, as we have seen for the area around the Arago spot in Fig. 4.

For the calculation of the umbra and penumbra, it might be sufficient to calculate the Fresnel diffraction pattern on a relatively small number of points, because the pattern is smoothed by the solar photospheric image. This involves the assumption that the two disks are close enough in comparison with the distance to the second spacecraft, so that the problem can be treated as a simple convolution, as in Aime (2013). Even then, the complete calculation down to the last focal plane of the coronagraph, as described in Rougeot et al. (2017), would be more intricate. These developments are left for a further study.

Acknowledgments

Thanks are due to Raphaël Rougeot for interesting discussions and to the anonymous referee for very constructive comments.

References

Appendix A: Analytic expression for the Fresnel diffraction of three disks

Let Ω, ω, and θ be the diameters of the three disk occulters. Let us assume for simplicity that there is an identical separation Δ between successive disks, and that d is still the distance between the last disk (θ) and the plane A. It can be shown after a few calculations that the wave on the telescope aperture can be written as:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{\mathrm{A}}(r)=& 1-\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_z(r)-\mathrm{\Pi }\left(\frac{r}{\omega }\right)\ast {\tau }_{\mathrm{d}+\mathrm{\Delta }}(r)-\mathrm{\Pi }\left(\frac{r}{\theta }\right)\ast {\tau }_{\mathrm{d}}(r)\hfill \\ \hfill & +\{[\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{\mathrm{\Delta }}(r)]\times \mathrm{\Pi }\left(\frac{r}{\omega }\right)\}\ast {\tau }_{\mathrm{d}+\mathrm{\Delta }}(r)\hfill \\ \hfill & +\{[\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{2\mathrm{\Delta }}(r)]\times \mathrm{\Pi }\left(\frac{r}{\theta }\right)\}\ast {\tau }_{\mathrm{d}}(r)\hfill \\ \hfill & +\{[\mathrm{\Pi }\left(\frac{r}{\omega }\right)\ast {\tau }_{\mathrm{\Delta }}(r)]\times \mathrm{\Pi }\left(\frac{r}{\theta }\right)\}\ast {\tau }_{\mathrm{d}}(r)\hfill \\ \hfill & +[(\{[\mathrm{\Pi }\left(\frac{r}{\mathrm{\Omega }}\right)\ast {\tau }_{\mathrm{\Delta }}(r)]\times \mathrm{\Pi }\left(\frac{r}{\omega }\right)\}\ast {\tau }_{\mathrm{\Delta }}(r))\times \mathrm{\Pi }\left(\frac{r}{\theta }\right)]\ast {\tau }_{\mathrm{d}}(r).\hfill \end{array}$$(A.1)

Here also, as for two disks, the diffraction can be described by complementary screens, with reference to Babinet’s theorem. Three terms correspond to the Fresnel diffractions of the disks taken independently; three terms correspond to Fresnel diffractions of disks taken two by two: Ω with ω, Ω with θ, and ω with θ, and are similar to the last term in Eq. (11). The last term of Eq. (A.1) corresponds to the diffraction of the wave passed through the three holes. It would be very difficult to evaluate this term numerically, which is one of the reasons for limiting this study to Fresnel diffractions with two disks.

All Figures

	Fig. 1. Schematic diagram for the Fresnel diffraction of two disks of diameters Ω and ω at distances of z and d from observing plane E. AB is the radius of the geometrical dark zone.
In the text

Fig. 2. Fresnel diffraction patterns of single and double disks. Parameters: Ω = 4 mm, ω = 3 mm, and z − d = d = 0.1 m. Top panel: diffraction of disks Ω (top left quadrant A) and ω (bottom left quadrant B) alone. The right two quadrants (C) show the diffraction of the system of two disks. Vertical and horizontal lines correspond to ±Ω/2 (external black lines), ±ω/2 (blue dashed lines) and the geometrical limit of the dark zone (red lines) AB of Fig. 1 and Eq. (7). Bottom panel: zoom on the central region around the Arago spot. In this example, the gain in darkness is about 8.

In the text

	Fig. 3. Red curves: Fresnel diffraction for two disks using Eq. (11) for Ω = 3 mm, ω = 2.25 mm, z = 1 m, and d = 0.5 m. The black curve shows the diffraction of the disk Ω alone (Eq. (3)). Top panel: linear scale, bottom: log10 scale. Vertical lines correspond to Ω/2 (black small dashed line), ω/2 (blue dashed line) and the geometrical limit AB (at 0.75 mm, (red line)) of Fig. 1 and Eq. (7).
In the text

Fig. 4. Fresnel diffraction with two disks for a central region close to the Arago spot, for Ω = 1.42 m, ω = 1.40524 m, z = 144.3 m, and z − d = 3 m. The top panel is in linear scale; the red dots correspond to Eq. (11) and the continuous curve to the Fresnel diffraction of the disk ω alone attenuated by a factor of 321. Bottom panel: two curves in log10 scales (black, disk Ω alone, red, the system of two disks).

In the text

	Fig. 5. Intensity of the Arago spot computed using Eq. (14). Parameters: z = 82 cm, Ω = 32 mm, λ = 0.55 μm. Linear (top) and log10 (bottom) representations are shown. The triangular region corresponds to the one in Fig. 1 and conditions given in Eq. (6).
In the text

	Fig. 6. Intensity of the Arago spot observed at z = 1 m for a system of two disks of dimensions Ω = 4.4 cm and ω = 4 cm, as a function of the distance z − d between the first and second disk. The vertical red dashed line corresponds to the distance when the second disk enters the blue triangle of Fig. 1.
In the text

Fig. 7. Intensity of the Arago spot as a function of the size of the second disk ω, in a linear (top) and log10 scale (bottom). The first disk of diameter Ω = 1.42 m is at 144.3 m from the telescope, the second disk ω is at d − z = 50 cm from the first one. Vertical red dashed lines correspond to limits given in Eq. (6), and the blue dashed line corresponds to the optimal value of Eq. (15). There is an increase in brightness outside these limits.

In the text

	Fig. 8. Rejection factor as a function of the distance z − d for Ω = 1.42 m and z = 144.3 m. The straight line corresponds to a linear fit. The corresponding variation of ωopt is from 1.418 m to 1.413 m.
In the text

	Fig. 9. Geometric elements for the derivation of the differential equation of Eq. (16) making use of the condition given in Eq. (15). The blue triangle has the same meaning as in Fig. 1.
In the text

	Fig. 10. Illustration of the optimal envelope for the solution given in Eq. (17) and α = 1. The exact curve (1) and linear approximation (2) of Eq. (18) are shown as the blue and red dashed lines, respectively. The drawing is made for z0 = 1 m, and the multiple disks extending up to 0.5 m. These values are for illustration only.
In the text