# Fresnel diffraction of multiple disks on axis

C. Aime

► To cite this version:

C. Aime. Fresnel diffraction of multiple disks on axis: Application to coronagraphy. Astronomy and Astrophysics - A&A, 2020, 637, pp.A16. 10.1051/0004-6361/201937208 . hal-02903862

HAL Id: hal-02903862

https://hal.science/hal-02903862

Submitted on 21 Jul 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Fresnel diffraction of multiple disks on axis
## Application to coronagraphy

C. Aime

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS UMR 7293 Lagrange, Parc Valrose 06108, Nice, France
e-mail: claude.aime@univ-cotedazur.fr

Received 28 November 2019 / Accepted 4 March 2020

## ABSTRACT

**Aims.** We seek to study the Fresnel diffraction of external occulters that differ from a single mask in a plane. Such occulters have been used in previous space missions and are planned for the future ESA Proba 3 ASPIICS coronagraph.

**Methods.** We studied the shading efficiency of double on-axis disks and generalized results to a 3D occulter. We used standard Fourier optics in an analytical approach. We show that the Fresnel diffraction of two and three disks on axis can be expressed using a Babinet-like approach. Results are obtained in the form of convolution integrals that can be written as Bessel-Hankel integrals; these are difficult to compute numerically for large Fresnel numbers found in solar coronagraphy.

**Results.** We show that the shading efficiency of two disks is well characterized by the intensity of the residual Arago spot, a quantity that is easier to compute and therefore allows an interesting parametric study. Very simple conditions are derived for optimal sizes and positions of two disks to produce the darkest structure around the Arago spot. These conditions are inspired from empirical experiments performed in the sixties. A differential equation is established to give the optimal envelope for a multiple-disk occulter. The solution takes the form of a simple law, the approximation of which is a conical occulter, a shape already used in the SOHO Mission.

**Conclusions.** In addition to quantifying expected results, the present study highlights unfortunate configurations of disks and spurious diffractions that may increase the stray light. Particular attention is paid to the possible issues of the future occulter spacecraft of ASPIICS.

**Key words.** instrumentation: high angular resolution – methods: analytical – Sun: corona – space vehicles: instruments

## 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^{10}$ fainter than that from its star, and, in order to consider exoplanets of the solar system type, an exo-Jupiter is about $10^6$ fainter in the far infrared. The K-corona is a few times $10^{-6}$ of the mean solar brightness close to the solar limb, and decreases to almost $10^{-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 $\Psi_0(r)$ be the complex amplitude of the wave in the plane $z=0$. Its expression $\Psi_z(r)$ after a free-space propagation over a distance $z$ may be written as:

$$ \Psi_z(r) = \Psi_0(r) * \tau_z(r), \quad \text{where} \\ \tau_z(r) = \frac{1}{i\lambda z} \exp\left(\frac{i\pi r^2}{\lambda z}\right), \quad \text{and} \quad \tau_z(r) * 1 = 1, \quad (1) $$

where $i$ is the imaginary unit, $\lambda$ 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\pi z/\lambda)$ has been omitted here. Taking the direct and inverse Fourier transform of Eq. (1), we obtain:

$$ \Psi_z(r) = \mathcal{F}^{-1}[\hat{\Psi}_0(u) \exp(-i\pi\lambda z u^2)], \quad (2) $$

where $\hat{\Psi}$ stands for the direct Fourier transform $\mathcal{F}[\Psi]$ of $\Psi$, and $\mathcal{F}^{-1}[...]$ 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 $\Omega$ centered on axis. The transmission of this occulter can be written as $1 - \Pi(r/\Omega)$, where $\Pi(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:

$$ \Psi_z(r) = \left[1 - \Pi\left(\frac{r}{\Omega}\right)\right] * \tau_z(r) = 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_z(r) \quad (3) $$

where we use the fact that $\tau_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, $\Psi_z(r)$ is a radial function and can be classically computed as an Hankel integral. We have, for the convolution term in Eq. (3),

$$ \Pi\left(\frac{r}{\Omega}\right) * \tau_z(r) = \frac{2\pi}{i\lambda z} \exp\left(\frac{i\pi r^2}{\lambda z}\right) \int_0^{\Omega/2} \xi \exp\left(\frac{i\pi\xi^2}{\lambda z}\right) J_0\left(\frac{2\pi\xi r}{\lambda z}\right) d\xi. \quad (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 = \Omega^2/\lambda z$. For ASPIICS, $N_F \sim 25000$, 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:

$$ \Psi_z(0) = \exp\left(\frac{i\pi\Omega^2}{4\lambda z}\right). \quad (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 $\Omega$ and $\omega$, 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:

$$ \frac{d}{z} \times \Omega < \omega < \Omega, \quad (6) $$

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:

$$ |AB| = \frac{\omega z - \Omega d}{2(z - d)}. \quad (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 $O_2^-$ immediately before the second disk, at the distance $z-d$, can be written as:

$$ \Psi_{O_2^{-}}(r) = 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_{z-d}(r). \quad (8) $$

The wave in plane $O_2^+$, immediately after the second occulter, is:

$$ \Psi_{O_2^+}(r) = \left[ 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_{z-d}(r) \right] \times \left[ 1 - \Pi\left(\frac{r}{\omega}\right) \right] \quad (9) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic diagram for the Fresnel diffraction of two disks of diameters $\Omega$ and $\omega$ at distances of $z$ and $d$ from observing plane E. AB is the radius of the geometrical dark zone.</figcaption></figure>

and the wave in plane E can be written as:

$$ \Psi_E(r) = \Psi_{O_2^+}(r) * \tau_d(r). \quad (10) $$

Since $\tau_{z-d}(r) * \tau_d(r) = \tau_z(r)$, Eq. (10) simplifies to:

$$ \Psi_E(r) = 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_z(r) - \Pi\left(\frac{r}{\omega}\right) * \tau_d(r) + \left\{ \left[ \Pi\left(\frac{r}{\Omega}\right) * \tau_{z-d}(r) \right] \times \Pi\left(\frac{r}{\omega}\right) \right\} * \tau_d(r). \quad (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:

$$ \Psi_E(r) = 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_z(r) - \left\{ \left[ 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_{z-d}(r) \right] \times \Pi\left(\frac{r}{\omega}\right) \right\} * \tau_d(r). \quad (12) $$

In that form, the complex amplitude $\Psi_E(r)$ can be expressed as the coherent composition, or interference, of two terms, one due to the disk $\Omega$ (propagation over $z$, defined by the convolution $*\tau_z(r)$) and the second due to the disk $\omega$ (propagation over $d$, the convolution $*\tau_d(r)$), the amplitude of the disk $\omega$ being modified by the Fresnel diffraction of $\Omega$ 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 $\Omega$ and $\omega$, 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 $\Omega$ and $\omega$, and its computation is more difficult. Indeed, applying twice the property that the Fourier transform of a 2D circular radial function is a Hankel
<figure><img src="image_3.png" /><figcaption>Fig. 2. Fresnel diffraction patterns of single and double disks. Parameters: $\Omega = 4$ mm, $\omega = 3$ mm, and $z-d = d = 0.1$ m. Top panel: diffraction of disks $\Omega$ (top left quadrant A) and $\omega$ (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 $\pm\Omega/2$ (external black lines), $\pm\omega/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.</figcaption></figure>

transform, we obtain:

$$ \left\{ \left[ \Pi\left(\frac{r}{\Omega}\right) * \tau_{z-d}(r) \right] \times \Pi\left(\frac{r}{\omega}\right) \right\} * \tau_d(r) \\ = \frac{-4\pi^2}{\lambda^2 d(z-d)} \\ \times \exp\left(\frac{i\pi r^2}{\lambda d}\right) \int_{\eta=0}^{\omega/2} \int_{\xi=0}^{\Omega/2} \xi\eta \exp\left(\frac{i\pi\eta^2 z}{\lambda d(z-d)}\right) \exp\left(\frac{i\pi\xi^2}{\lambda(z-d)}\right) \\ \times J_0\left(\frac{2\pi\xi\eta}{\lambda(z-d)}\right) J_0\left(\frac{2\pi r\eta}{\lambda d}\right) d\xi d\eta. \qquad (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.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Red curves: Fresnel diffraction for two disks using Eq. (11) for $\Omega = 3$ mm, $\omega = 2.25$ mm, $z = 1$ m, and $d = 0.5$ m. The black curve shows the diffraction of the disk $\Omega$ alone (Eq. (3)). Top panel: linear scale, bottom: log10 scale. Vertical lines correspond to $\Omega/2$ (black small dashed line), $\omega/2$ (blue dashed line) and the geometrical limit AB (at 0.75 mm, (red line)) of Fig. 1 and Eq. (7).</figcaption></figure>

## 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.

Figure 2 shows the Fresnel diffraction of two disks on-axis compared with diffractions of disks taken alone. The parameters are $\Omega = 4$ mm, $\omega = 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 $\Omega$ has a diameter of 40% of the array. The Fresnel diffraction of the wave after the disk $\Omega$ is calculated using Eq. (2) for the distance $z-d = 0.1$ m. Then the transmission of the second disk $\omega$ 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 $\Omega$ alone, (B) that of the disk $\omega$ 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
C. Aime: Fresnel diffraction of multiple disks

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Fresnel diffraction with two disks for a central region close to the Arago spot, for $\Omega = 1.42$ m, $\omega = 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 $\omega$ alone attenuated by a factor of 321. *Bottom panel:* two curves in log10 scales (black, disk $\Omega$ alone, red, the system of two disks).</figcaption></figure>

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 $\Omega$, the influence of the second disk $\omega$ 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 $\omega$, 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 $|\Psi_E(r)|^2$ 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 $\Omega = 3$ mm, $\omega = 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

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Intensity of the Arago spot computed using Eq. (14). Parameters: $z = 82$ cm, $\Omega = 32$ mm, $\lambda = 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).</figcaption></figure>

that could be obtained if the disk $\Omega$ was alone (Eq. (3)). Outside the central dark zone, the Fresnel diffraction of the first disk $\Omega$ 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 $|\Psi_E(r)|^2$ for the central region very close to the Arago spot for $\Omega = 1.42$ m, $\omega = 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 $\Omega \times \omega$, using Eq. (11). The Fresnel diffraction of the single disk $\omega$ 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 $\omega$ in the umbra of $\Omega$ has indeed the effect to make darker its edges, and the shadow obtained resembles the diffraction of the disk $\omega$ illuminated by a less bright source.

We note that, for the conditions of Fig. 4, the central part of the Fresnel diffraction of $\omega$ can be very well approximated by the function $I_0 \times J_0^2(\pi\omega r/\lambda d)$, where $J_0$ 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(0) = 1$, the intensity at $r = 0$ is the value $I_0$. 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 $\omega$ 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 $\Psi_A(0)$, that is, the Arago spot.

Because $J_0(0) = 1$, we have:

$$ \Psi_A(0) = \exp\left(i\pi\frac{\Omega^2}{4\lambda z}\right) + \exp\left(i\pi\frac{\omega^2}{4\lambda d}\right) - 1 \\ - \frac{4\pi^2}{\lambda^2 d(z-d)} \int_{\eta=0}^{\omega/2} \int_{\xi=0}^{\Omega/2} \xi\eta \exp\left(\frac{i\pi\eta^2 z}{\lambda d(z-d)}\right) \\ \times \exp\left(\frac{i\pi\xi^2}{\lambda(z-d)}\right) J_0\left(\frac{2\pi\xi\eta}{\lambda(z-d)}\right) d\xi d\eta \quad (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 $\Omega = 32$ mm is given as a 2D contour plot in Fig. 5. The second occulter $\omega$ 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.

The results of a much easier experiment are simulated in Fig. 6. In this experiment the disks $\Omega$ and $\omega$ are of fixed size (in the example we took $\Omega = 44$ mm and $\omega = 40$ mm). We assume that the disk $\Omega$ is fixed and that the disk $\omega$ may be moved away

<figure><img src="image_7.png" /><figcaption>Fig. 6. Intensity of the Arago spot observed at $z = 1$ m for a system of two disks of dimensions $\Omega = 4.4$ cm and $\omega = 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.</figcaption></figure>

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 $\Omega \times \omega$ 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 $\Omega$ and $\omega$, the amplitude of the latter term being modified by the Fresnel diffraction of $\Omega$ 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.

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 $\Omega$ being oriented towards the sun and the second disk $\omega$ 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 $\omega = d\Omega/z$ and $\omega = \Omega$, 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 $\Omega$ and $\omega$. For the other higher limit, one can consider that the first disk, $\Omega$, diffracts light towards the second disk, $\omega$, which by interference with the direct light may increase or decrease its illumination.

The best darkening effect appears to be at the mean distance between these two limits, that is, for an optimal $\omega_{opt}$ value:

$$ \omega_{opt} = \Omega \frac{z+d}{2z}. \quad (15) $$
<figure><img src="image_8.png" /><figcaption>Fig. 7. Intensity of the Arago spot as a function of the size of the second disk $\omega$, in a linear (top) and log10 scale (bottom). The first disk of diameter $\Omega = 1.42$ m is at 144.3 m from the telescope, the second disk $\omega$ 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.</figcaption></figure>

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 $\omega_{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 $\omega_{opt}$, increases almost linearly with the distance $z-d$ between disks, as shown in Fig. 8, for $\Omega = 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.

## 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_0$ 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 $\omega_{opt}$ in Eq. (15), but for more generality, let us assume that we choose a direction PQ such as $AM = \alpha MA'$. With $HP = y$, $HA = z_0 - x$,

<figure><img src="image_9.png" /><figcaption>Fig. 8. Rejection factor as a function of the distance $z-d$ for $\Omega = 1.42$ m and $z = 144.3$ m. The straight line corresponds to a linear fit. The corresponding variation of $\omega_{opt}$ is from 1.418 m to 1.413 m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

simple geometric considerations give $HQ = (1 + \alpha)(z_0 - x)$, and we can write:

$$ \frac{HP}{HQ} = y' = -\frac{y}{(1 + \alpha)(z_0 - x)}. \qquad (16) $$

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

$$ y(x) = \frac{\Omega}{2} \times \left(1 - \frac{x}{z_0}\right)^{\frac{1}{1+\alpha}} \qquad (17) $$

where $y(x)$ is the generatrix of the envelope around the rotation axis $x$. If we can assume that $x \ll z_0$, then the solution is a straight line:

$$ y(x) \sim \frac{\Omega}{2} \left( 1 - \frac{x}{(1 + \alpha)z_0} \right) \qquad (18) $$

that generates a conical envelope. An illustration is given in Fig. 10 for $\alpha = 1$. The exact solution is compared with its linear approximation, corresponding to a cone as in LASCO C2, Bout et al. (2000). The value $\alpha = 1$ produces a rather small dark zone, and a larger value of $\alpha$ 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.

## 5. Conclusion

The present work is an analytical study of the diffraction produced by a nonsingle-plane occulter. In particular the Fresnel
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Illustration of the optimal envelope for the solution given in Eq. (17) and $\alpha = 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 $z_0 = 1$ m, and the multiple disks extending up to 0.5 m. These values are for illustration only.</figcaption></figure>

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 $\Omega \times \omega$, one at each end of the satellite. For example the photovoltaic sensors might have been included in a larger first disk acting as $\Omega$, the second disk $\omega$ 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.

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

## References

* Aime, C. 2013, A&A, 558, A138

* Aschwanden, M. 2006, Physics of the Solar Corona: An Introduction with Problems and Solutions (Springer Science & Business Media)

* Baccani, C., Landini, F., Romoli, M., et al. 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 9904, 990450

* Boivin, L. 1978, Appl. Opt., 17, 3323

* Born, M., & Wolf, E. 1999, Principles of Optics, 7th edn. (Cambridge: Cambridge University Press), 461

* Bout, M., Lamy, P., Maucherat, A., Colin, C., & Llebaria, A. 2000, Appl. Opt., 39, 3955

* Cady, E. 2012, Opt. Exp., 20, 15196

* Cash, W. 2006, Nature, 442, 51

* Cox, A. N. 2015, Allen's Astrophysical Quantities (New York: Springer)

* Domingo, V., Fleck, B., & Poland, A. I. 1995, Sol. Phys., 162, 1

* Evans, J. W. 1948, J. Opt. Soc. Am., 38, 1083

* Ferrari, A. 2007, ApJ, 657, 1201

* Ferrari, A., Aime, C., & Soummer, R. 2009, ApJ, 708, 218

* Fort, B., Morel, C., & Spaak, G. 1978, A&A, 63, 243

* Harvey, J. E., & Forgham, J. L. 1984, Am. J. Phys., 52, 243

* Kasdin, N. J., Cady, E. J., Dumont, P. J., et al. 2009, in Techniques and Instrumentation for Detection of Exoplanets IV, Proc. SPIE, 7440, 744005
Kelly, W. R., Shirley, E. L., Migdall, A. L., Polyakov, S. V., & Hendrix, K. 2009, Am. J. Phys., 77, 713

Koutchmy, S. 1988, Space Sci. Rev., 47, 95

Koutchmy, S., & Belmahdi, M. 1987, J. Opt., 18, 265

Landini, F., Baccani, C., Schweitzer, H., et al. 2017, Opt. Lett., 42, 4800

Lenskii, A. 1981, Sov. Astron., 25, 366

Lyot, B. 1932, Z. Astrophys., 5, 73

Mathematica, W. 2009, Wolfram Research (Champaign, Illinois: Wolfram Media Inc.)

Mawet, D., Riaud, P., Absil, O., & Surdej, J. 2005, ApJ, 633, 1191

Mayor, M., & Queloz, D. 1995, Nature, 378, 355

Newkirk, G., Jr., & Bohlin, D. 1963, Appl. Opt., 2, 131

Newkirk, G., Jr., & Bohlin, J. D. 1965, in Astronomical Observations from Space Vehicles, ed. J. L. Steinberg, IAU Symp., 23, 287

Peter, H., Bingert, S., Klimchuk, J., et al. 2013, A&A, 556, A104

Prasad, B. R., Banerjee, D., Singh, J., et al. 2017, Current Science, 113, 613

Purcell, J., & Koomen, M. 1962, J. Opt. Soc. Am., 52, 596

Renotte, E., Baston, E. C., Bemporad, A., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 9143, 91432M

Roddier, F., & Roddier, C. 1997, PASP, 109, 815

Rouan, D., Riaud, P., Boccaletti, A., Clénet, Y., & Labeyrie, A. 2000, PASP, 112, 1479

Rougeot, R., & Aime, C. 2018, A&A, 612, A80

Rougeot, R., Aime, C., Baccani, C., et al. 2018, in Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 10698, 106982T

Rougeot, R., Flamary, R., Galano, D., & Aime, C. 2017, A&A, 599, A2

Rougeot, R., Flamary, R., Mary, D., & Aime, C. 2019, A&A, 626, A1

Shestov, S., & Zhukov, A. 2018, A&A, 612, A82

Shestov, S., Zhukov, A., & Seaton, D. 2019, A&A, 622, A101

Soummer, R., Aime, C., & Falloon, P. 2003, A&A, 397, 1161

Thernisien, A., Colaninno, R., Plunkett, S., et al. 2005, in Solar Physics and Space Weather Instrumentation, Proc. SPIE, 5901, 59011E

Tousey, R. 1965, Ann. Astrophys., 28, 600

Verroi, E., Frassetto, F., & Naletto, G. 2008, J. Opt. Soc. Am. A, 25, 182

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

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

$$ \Psi_A(r) = 1 - \Pi\left(\frac{r}{\Omega}\right) * \tau_z(r) - \Pi\left(\frac{r}{\omega}\right) * \tau_{d+\Delta}(r) - \Pi\left(\frac{r}{\theta}\right) * \tau_d(r) \\ + \left\{ \Pi\left(\frac{r}{\Omega}\right) * \tau_{\Delta}(r) \times \Pi\left(\frac{r}{\omega}\right) \right\} * \tau_{d+\Delta}(r) \\ + \left\{ \Pi\left(\frac{r}{\Omega}\right) * \tau_{2\Delta}(r) \times \Pi\left(\frac{r}{\theta}\right) \right\} * \tau_d(r) \\ + \left\{ \Pi\left(\frac{r}{\omega}\right) * \tau_{\Delta}(r) \times \Pi\left(\frac{r}{\theta}\right) \right\} * \tau_d(r) \\ + \left[ \left\{ \Pi\left(\frac{r}{\Omega}\right) * \tau_{\Delta}(r) \times \Pi\left(\frac{r}{\omega}\right) \right\} * \tau_{\Delta}(r) \times \Pi\left(\frac{r}{\theta}\right) \right] * \tau_d(r). \quad (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 diffrac- tions of disks taken two by two: $\Omega$ with $\omega$, $\Omega$ with $\theta$, and $\omega$ with $\theta$, 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.