# Simulation of inverse Fredholm reconstruction in a vignetting zone: application to ASPIICS

C. Theys, C. Aime, R. Rougeot, H. Lantéri

► To cite this version:

C. Theys, C. Aime, R. Rougeot, H. Lantéri. Simulation of inverse Fredholm reconstruction in a vignetting zone: application to ASPIICS. Astronomy and Astrophysics - A&A, 2022, 665, 10.1051/0004-6361/202141217. insu-04004082

HAL Id: insu-04004082

https://insu.hal.science/insu-04004082

Submitted on 24 Feb 2023

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
# Simulation of inverse Fredholm reconstruction in a vignetting zone: application to ASPIICS

C. Theys<fnref n="1" />, C. Aime<fnref n="1" />, R. Rougeot<fnref n="2" />, and H. Lantéri<fnref n="1" />

Received 30 April 2021 / Accepted 4 July 2022

## ABSTRACT

**Aims.** This work deals with image reconstruction in a vignetting zone where the point spread function becomes evanescent and the image undergoes a Fredholm transformation. The application of this method is aimed at the reconstruction of the solar corona in the vignetting zone of the ASPIICS coronagraph. It extends on a previous paper in several aspects.

**Methods.** We used a matrix formalism for the exact inversion of the Fredholm integral. The stray light appears there as a bias. We performed two procedures: either the direct processing of the biased data or their processing following the subtraction of the bias. In the first case, the statistics follow a Poisson distribution and a Kullback-Leibler divergence was used; in the second case, we were led to use a simplifying Gaussian statistic. In both cases, a physical regularization using a Strehl criterion was implemented and this improved the results. Image reconstruction in the vignetting area is done in the case of a perfect coronagraph for two diameters of the internal occulter, but also in the case of formation flight errors and optical defects that are present but ignored by the inversion procedure.

**Results.** Poisson and Gauss models both give much better results than simple flux compensation. For the Poisson model, unexpected pseudo-fringes are present in the reconstructed raw image but are greatly reduced using regularization. The Gaussian model (using de-biased data) is found to give better results, no matter whether it is the regularized or non-regularized version of the algorithm that is used. Despite a high level of stray light, the internal occulter of a smaller dimension allows us to approach much closer to the solar edge without too great a loss in terms of quality in the outer regions. This conclusion remains true in the case of optical micro-defects leading to speckles in the PSF because that has only a slight impact on the images in the vignetting area. In the case of formation flying errors, the Fredholm inversion is more affected by these for the small internal occulter than for the larger one.

**Conclusions.** The method proposed for the Fredholm inversion is general and can be transposed to other systems using external occulters. An application of this method to the imaging of exoplanets is generally envisaged.

**Key words.** techniques: high angular resolution – methods: numerical – techniques: image processing – Sun: corona

## 1. Introduction

In this paper, we study the reconstruction of an object based on an observation that does not result from a relationship of mere convolution between the object and the response of the instrument but, rather, requires a complete Fredholm integral of the first kind for its description, Groetsch (2007). In the present application, the object is the solar corona and the instrument is the coronagraph ASPIICS, using an external occulter in front of a Lyot (1933) coronagraph, as per the configuration first proposed by Evans (1948).

The external occulter creates an artificial eclipse for the coronagraph set in its shadow and allows for a better reduction of the straylight, as analyzed by Newkirk & Bohlin (1963) and Koomen et al. (1975), for example. The drawback of such system is the difficulty in carrying out an observation very close to the solar limb due to vignetting. The aperture of the telescope becomes progressively reduced to a crescent shape by the umbra of the external occulter located in front of it. This was well illustrated in Llebaria et al. (2004, 2006) for the LASCO-C2 coronagraph of the SOHO mission, Brueckner et al. (1995). In spite of these vignetting effects, which do not exist in a pure Lyot (1933) coronagraph, the strong reduction of the stray light prevails; and an Evans configuration has been adopted in almost all space experiments, as noted by Koutchmy (1988).

As a first approximation, the loss of the inner observing angle close to the limb is the ratio of the diameter of the telescope to the diameter of the occulter, in units of solar radius. A small vignetting zone requires a very large occulter that is set far away from the telescope. Such a result cannot be obtained with a single spacecraft and the formation flying project ESA's Proba-3 ASPIICS is an excellent way to solve this problem, Lamy et al. (2010); Renotte et al. (2014); Galano et al. (2018).

The first satellite holds a 1.42-m diameter occulting disk and the second satellite, positioned at 144-m away in its umbra, carries a 5-cm Lyot coronagraph, where the final observation is made, Galano et al. (2019). ASPIICS aims to observe the solar corona very close to the solar limb, as has never been done before beyond eclipse conditions. The first proposals for ASPIICS (Lamy et al. 2007, 2010, 2017) carried the ambition of an observation of the solar corona very close to the limb, as low as 1.01 $R_⊙$ (solar radius), and even the goal of including chromospheric and prominence physics. The current project, for reasons of reliability in the event of positioning errors, will place the telescope pupil in the shadow cone at a position where the shadow disk is significantly larger than the pupil. The drawback is that this increases the blind area near the solar limb. Recent studies have envisaged observations starting at no sooner than 1.05 $R_⊙$ or even 1.1 $R_⊙$, depending on the experimental parameters used, Shestov et al. (2021) and Rougeot (2020).

<footnote n="1">¹ CNRS UMR 7293 Lagrange Observatoire de la Côte d'Azur – Université Côte d'Azur 06108 Nice Cedex 2, France
e-mail: celine.theys@univ-cotedazur.fr</footnote>

<footnote n="2">² European Space Research and Technology Center, European Space Agency, Keplerlaan 1, 2201 Noordwijk, The Netherlands</footnote>
In the region of vignetting, the point spread function (PSF) degrades rapidly before disappearing when the shadow of the occulter covers the entire entrance pupil. While the main effect is geometrical, diffraction makes the response more complex and the frontier between light and darkness less well defined. The spatial variation of the PSF in the region of vignetting is very fast, and solutions considered in the field of adaptive optics, such as in La Camera et al. (2015), Denis et al. (2015), Thiébaut et al. (2016), cannot apply to the present work.

The present paper follows and complements the study of Aime et al. (2019), based on the assumption of a Richardson-Lucy type algorithm, (Richardson 1972; Lucy 1974), applied and adapted to the case of a non-constant PSF. This keeps the complete matrix formalism and upgrades the analysis with regard to several fundamental aspects. A realistic solar object is used instead of simple point sources. The stray light present in the observation is taken into account, and not only the direct Fredholm transformation as in the former study. This stray light increases since the internal occulter used is smaller and a compromise between resolution near the solar limb and reduced stray light is to be considered.

This work has a signal processing aspect and a physical aspect that are considered in that order in this paper. In Sect. 2, we describe the Fredholm inversion in the case of a space variant PSF. The matrix formalism and the reconstruction algorithm are given there. Different statistical models adapted to the raw data (Poisson model) and to the de-biased data (Gauss-Poisson and Gaussian models) are presented in Sect. 3. The PSFs and biases calculated for ASPIICS with two values of the internal occulter are given in Sect. 4. Numerical results performed on an eclipse image are shown in Sect. 5. A robustness study in the case of instrumental defects is presented in Sect. 6 for formation flying errors that are present but ignored in the inversion procedure and for the surface roughness of the optics. Prospects for the extension of the matrix formalism to other instruments are proposed in Sect 7. Our conclusions are given in Sect. 8. Peculiarities of the PSFs in the region dominated by diffraction are detailed in the Appendix.

## 2. Fredholm inversion principle

As stated in the introduction, the objective of this work is to reconstruct the object $X(\alpha, \beta)$ in the vignetting area of an externally occulted coronagraph. The intensity $Y(\alpha, \beta)$ that is formed in the final observing plane of the coronagraph can be written as:

$$Y(\alpha, \beta) = \iint X(\xi, \eta)h(\alpha, \xi, \beta, \eta)d\xi d\eta + B(\alpha, \beta), \quad (1)$$

where $h(\alpha, \xi, \beta, \eta)$ is the space-variant PSF in the region of vignetting, the integral a two-dimensional (2D) Fredholm integral of the first kind and $B(\alpha, \beta)$ is the stray light that appears here as a bias. This equation corresponds to Eq. (3) of Aime et al. (2019), with the addition of stray light. We note that dummy variables $(\xi, \eta)$ and coordinates $(\alpha, \beta)$ are associated with the object and image planes, respectively. In the following, we give the elements of the matrix approach to discretize Eq. (1).

### 2.1. Matrix formalism

We note that in the case of an invariant PSF, the Fredholm relation in Eq. (1) can be simplified to a convolution relation, allowing for a fast resolution in Fourier space via FFTs. This

is no longer the case for the variant PSF considered here and matrix manipulations must be carried out. The discretized form of Eq. (1) can be written as:

$$y_{m,n} = \sum_{i=1}^{M} \sum_{j=1}^{N} x_{i,j} h_{m,i,n,j} + b_{m,n} \quad (2)$$

where $y_{m,n}$, $x_{i,j}$ and $b_{m,n}$ are the $M \times N$ discretized values of $Y(\alpha, \beta)$, $X(\xi, \eta)$, and $B(\xi, \eta)$, respectively, while $h_{m,i,n,j}$ are the $M^2 \times N^2$ discretized values of $h(\alpha, \xi, \beta, \eta)$.

To use a matrix formalism, we arrange the arrays in vectors. The order chosen is lexicographic, that are the values of each array of dimension $M \times N$ (read from left to right and from top to bottom) are arranged in this order in a vector. An image of size $M \times N$ then yields a vector $\mathbf{y}$ of length $M \times N$:

$$\mathbf{y} = [y_{1,1}, y_{1,2}, \dots, y_{2,1}, \dots, y_{M,N}]^T. \quad (3)$$

In doing so, Eq. (2) can be written in the following matrix form:

$$\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{b}, \quad (4)$$

where $\mathbf{y}, \mathbf{x}$ and $\mathbf{b}$ are vectors of size $M \times N$ containing the discretized values of the image and the object. $\mathbf{H}$ is the matrix containing the elements $h_{m,i,n,j}$, of dimension $MN \times MN$:

$$\mathbf{H} = \begin{pmatrix} h_{1,1,1,1} & h_{1,1,1,2} & \cdots & h_{1,N,1,N} \\ h_{1,1,2,1} & h_{1,1,2,2} & \cdots & h_{1,N,2,N} \\ \vdots & \vdots & \ddots & \vdots \\ h_{M,1,M,1} & h_{M,1,M,2} & \cdots & h_{M,N,M,N} \end{pmatrix}. \quad (5)$$

Each column of $\mathbf{H}$ corresponds to one PSF ordered lexicographically. The first column corresponds to the PSF of the point source at the upper left corner of the object ($i = j = 1$), while the last column corresponds to the PSF of the point source at the bottom-right corner of the object ($i = M, j = N$). It should be noted that inevitable boundary problems exist in the discrete representation of the Fredholm integral, as they exist for the simple convolution relation.

### 2.2. Reconstruction algorithm

The estimation of $\mathbf{x}$ given the data, $\mathbf{y}$, matrix, $\mathbf{H}$, and bias, $\mathbf{b}$, namely, Eq. (4), is a so-called ill-posed problem. In fact, the direct inversion of this equation provides unacceptable solution due to very bad conditioning of $\mathbf{H}$. The classical method to solve such an inverse problem is to use an iterative algorithm to minimize a criterion which measure a discrepancy between the observation $\mathbf{y}$ and the linear model $\mathbf{H}\mathbf{x}$. Iterations have to be stopped before instabilities in the solution appear unless an explicit regularization is used, (Demoment 1989; Titterington 1985). Since we are dealing with intensities, a non-negativity constraint has to be imposed on $\mathbf{x}$. The problem is then to solve:

$$\min_{\mathbf{x}} J(\mathbf{x}), \quad \text{s.t} \quad x_i \ge 0, \quad \forall i, \quad (6)$$

with $J$ a criterion depending on $\mathbf{x}$, and, of course, on $\mathbf{y}$ as well. As is classically done, we decompose $J$ into two terms:

$$J(\mathbf{x}) = J_1(\mathbf{x}) + \gamma J_2(\mathbf{x}), \quad (7)$$

where $J_1(\mathbf{x})$ and $J_2(\mathbf{x})$ are, respectively, the data consistency and the regularization terms. The coefficient $\gamma$ makes it possible
to adjust the weight of $J_1$ with respect to $J_2$, (see Titterington 1985). Depending on the statistical model for the noise, the data consistency term $J_1(\mathbf{x})$ takes different forms, as we show in Sect. 3.

The regularization term $J_2$ is independent from the data consistency term, that is, it is a regularization term and must, in the context of an ill-posed problem, impose a constraint of smoothness on the solution $\mathbf{x}$, Titterington (1985). The choice of $J_2$ is further discussed in Sect. 5.2, once the particular form of the ASPIICS PSF has been presented (in Sect. 4).

$J(\mathbf{x})$ is assumed to be convex relative to $\mathbf{x}$ and, therefore, the simplest way to minimize it is to use a gradient descent algorithm with non-negativity constraint. The procedure is described in Lantéri et al. (2002) and summarized here.

The Karush Kuhn Tucker (KKT) conditions, (Karush 1939; Kuhn & Tucker 1951), once resolved, can be reduced to an unique equation that the optimal value of the solution, $\mathbf{x}^*$, must satisfy, that is, for a component $r$:

$$[\nabla_{\mathbf{x}} J(\mathbf{x}^*)]_r \mathbf{x}_r^* = 0, \quad r = 1, \dots, MN. \qquad (8)$$

This condition will then serve as a descent direction in the iterative algorithm:

$$\mathbf{x}_r^{k+1} = \mathbf{x}_r^k + \alpha_r^k \mathbf{x}_r^k [-\nabla_{\mathbf{x}} J(\mathbf{x}^k)]_r, \qquad (9)$$

with $\alpha_r^k$, the descent step-size chosen to have $J(\mathbf{x}^{(k+1)}) < J(\mathbf{x}^k)$ at each iteration $k$. The search procedure is carried out in an interval computed at each $k$ to ensure $\mathbf{x}^{k+1} \ge 0$, starting from a positive initial value, making it an interior point algorithm, (see Armijo 1966; Lantéri et al. 2002).

Making use of Eq. (7), the gradient $\nabla_{\mathbf{x}} J$ necessary for the implementation of the algorithm of Eq. (9), can be split into the two terms:

$$\nabla_{\mathbf{x}} J = \nabla_{\mathbf{x}} J_1 + \gamma \nabla_{\mathbf{x}} J_2. \qquad (10)$$

The gradient of the data consistency term $\nabla_{\mathbf{x}} J_1$ depends on the given model and is given in Sect. 3, $\nabla_{\mathbf{x}} J_2$ will be given in the Sect. 5.2.

## 3. Data modeling

This section is dedicated to the development of statistical models for raw data and for debiased data further on. We show in what follows that the drawback of subtracting the bias is that it makes it more difficult to set up of the statistical model needed for the Fredholm inversion.

### 3.1. Raw data (raw) and the Kullback-Leibler divergence

We may consider that the discretized values of $y_{m,n}$ in Eq. (2) correspond to a classical integrated intensity in the pixels of the detector. These intensities, converted into units of photoelectrons, undergo a photodetection process, which is well represented, for an ideal detector by a Poisson distribution $P$. We remark that in the case of this model, the $y_{raw}$ data are non-negative.

The corresponding noisy data are denoted as $\mathbf{y}_{raw}$:

$$\mathbf{y}_{\text{raw}} \sim \mathcal{P}(\mathbf{H}\mathbf{x} + \mathbf{b}). \qquad (11)$$

Under the assumption of Eq. (11), the raw data follow a Poisson statistics. Assuming moreover that all pixels of the image $y_i$ are

independent, the likelihood or conditional probability $p(\mathbf{y}_{\text{raw}}/\mathbf{x})$ may be written as:

$$p(\mathbf{y}_{\text{raw}}/\mathbf{x}) = \prod_{i=1}^{MN} \frac{(\mathbf{H}\mathbf{x} + \mathbf{b})_{i}^{(\mathbf{y}_{\text{raw}})_i}}{(\mathbf{y}_{\text{raw}})_i!} \exp(-(\mathbf{H}\mathbf{x} + \mathbf{b})_i). \qquad (12)$$

In order to achieve a minimization algorithm, we take the opposite of the log likelihood to obtain $J_1$, in Eq. (7):

$$J_1(\mathbf{x}) = -\log(p(\mathbf{y}_{\text{raw}}/\mathbf{x})), \qquad (13)$$

$$\propto \sum_i -(\mathbf{y}_{\text{raw}})_i \log(\mathbf{H}\mathbf{x} + \mathbf{b})_i + (\mathbf{H}\mathbf{x})_i, \qquad (14)$$

and the gradient w.r.t $\mathbf{x}$, to get $\nabla_{\mathbf{x}} J_1$, Eq. (10) is:

$$\nabla_{\mathbf{x}} J_1 = \mathbf{H}^T (-\mathbf{y}_{\text{raw.}}/(\mathbf{H}\mathbf{x} + \mathbf{b}) + 1_{MN}), \qquad (15)$$

where ./ indicates element-wise division.

This expression of the gradient of $J_1$, $\nabla_{\mathbf{x}} J_1$ associated with the gradient of $J_2$, Eq. (29), will be used in the iterative reconstruction algorithm, Eq. (9), to process raw data.

We note that in the non-regularized case, namely, $\gamma = 0$ and with a descent step-size $\alpha_r^k = 1$, the algorithm obtained from Eq. (9) in the case of a Poisson statistic is equivalent – although not strictly equal – to the so-called Richardson Lucy algorithm.

### 3.2. Debiased data (deb)

Since we assume a knowledge of the bias, that is a numerical computation of the straylight, it seems an a priori natural choice to subtract it as part of the pre-processing of the data. Let us denote the new bias-removed data, $\mathbf{y}_{deb}$:

$$\mathbf{y}_{\text{deb}} \sim \mathcal{P}(\mathbf{H}\mathbf{x} + \mathbf{b}) - \tilde{\mathbf{b}} \sim \mathcal{P}(\mathbf{H}\mathbf{x}) + (\mathcal{P}(\mathbf{b}) - \tilde{\mathbf{b}}), \qquad (16)$$

where $\tilde{\mathbf{b}}$ is the estimate of the straylight. Contrary to the Poisson case, the $\mathbf{y}_{deb}$ data can be negative, that is the case once the bias is removed from the data. We can approximate $\mathcal{P}(\mathbf{b}) - \tilde{\mathbf{b}}$ with a Gaussian law of zero mean and mean $\tilde{\mathbf{b}}$, because of the large number of photons in this solar experiment. Two statistics are derived from there, presented below.

#### 3.2.1. Gauss-Poisson model

We considered the data after removing the bias, Eq. (16). Due to the subtraction of an estimate of the bias, some pixels became negative and the observed quantity can no longer be written as a pure Poisson process. The model of Eq. (16) then becomes a Gauss-Poisson model, such that:

$$\mathbf{y}_{\text{deb}} \sim \mathcal{P}(\mathbf{H}\mathbf{x}) + \mathbf{g}, \qquad (17)$$

with $\mathbf{g} \sim N(0, \text{diag}(\mathbf{b}))$, where $\text{diag}(\mathbf{b})$ is the diagonal covariance matrix with diagonal $\mathbf{b}$.

This Gauss-Poisson model has been the subject of several astrophysical publications in e.g Lantéri & Theys (2005); Benvenuto et al. (2008). The probability function in this case is the convolution of the Poisson and Gauss probability laws for a pixel, as per Eq. (3) of Benvenuto et al. (2008):

$$p(\mathbf{y}_i/\mathbf{x}) = \sum_{n_i} \frac{1}{\sqrt{\mathbf{b}_i} \sqrt{2\pi}} \exp\left(-\frac{(\mathbf{y}_i - n_i)^2}{2\mathbf{b}_i}\right) \frac{(\mathbf{H}\mathbf{x})_{i}^{n_i}}{n_i!} \exp(-\mathbf{H}\mathbf{x})_{i}, \qquad (18)$$
and the gradient for **x** is, for a pixel, *i*:

$$[\nabla J_1(\mathbf{x})]_i = \sum_j \left(h_{ji} - h_{ji} \frac{r_j}{(H\mathbf{x})_j}\right) = 0, \qquad (19)$$

with

$$r_j = \frac{p_j}{q_j}, \qquad (20)$$

$$p_j = \sum_{n_i} \exp\left(-\frac{(y_j - n_i)^2}{2\mathbf{b}_j}\right) \frac{n_i}{n_i!} (H\mathbf{x})_j^{n_i} \qquad (21)$$

and

$$q_j = \sum_{n_i} \exp\left(-\frac{(y_j - n_i)^2}{2\mathbf{b}_j}\right) \frac{1}{n_i!} (H\mathbf{x})_j^{n_i}. \qquad (22)$$

It is important to note that the terms to be calculated in the gradient expression above can lead to numerical instability for large values of $n_i$. This leads us to consider a simplified model.

### 3.2.2. Gaussian model

A simplified and more tractable model for debiased data is a full Gaussian model with a covariance matrix depending on the data:

$$\mathbf{y}_{\text{deb}} \sim \mathcal{N}(\mathbf{H}\mathbf{x}, \text{diag}(\mathbf{y}_{\text{deb}})), \qquad (23)$$

where diag($\mathbf{y}_{\text{deb}}$) is the covariance matrix with $\mathbf{y}_{\text{deb}}$ on the diagonal. Then, the conditional probability, assuming that all pixels of the image $y_i$ are independent, is:

$$p(\mathbf{y}_{\text{deb}}/\mathbf{x}) = \prod_{i=1}^{MN} \frac{1}{\sqrt{2\pi(\mathbf{y}_{\text{deb}})_i}} \exp\left(-\frac{(\mathbf{y}_{\text{deb}} - H\mathbf{x})_i^2}{2(\mathbf{y}_{\text{deb}})_i}\right), \qquad (24)$$

In order to achieve the minimization algorithm, we take the opposite of the log likelihood:

$$J_1(\mathbf{x}) = -\log(p(\mathbf{y}_{\text{deb}})) \qquad (25)$$

$$\propto \sum_i \frac{(\mathbf{y}_{\text{deb}} - \mathbf{H}\mathbf{x})_i^2}{(\mathbf{y}_{\text{deb}})_i}, \qquad (26)$$

and the gradient for **x**:

$$\nabla_{\mathbf{x}} J_1 = -\mathbf{H}^T \text{diag}(1./\mathbf{y}_{\text{deb}}) \mathbf{y}_{\text{deb}} + \mathbf{H}^T \text{diag}(1./\mathbf{y}_{\text{deb}}) \mathbf{H}\mathbf{x}. \qquad (27)$$

In the same way as for the KL divergence, this expression of $\nabla_{\mathbf{x}} J_1$, associated with the gradient of $J_2$, will be used in the iterative reconstruction algorithm, namely, Eq. (9).

## 4. Application to ASPIICS

In order to derive the values of the matrix **H** and the straylight **b** for ASPIICS, we made use of the numerical model developed by Rougeot et al. (2017). This approach was further used as a starting point for studies including optical defects of the coronagraph (Rougeot et al. 2017, 2019; Shestov & Zhukov 2018; Shestov et al. 2021).

The program, written in MATLAB, computes the propagation of a wave originating from a point source at the level of the Sun through the whole ASPIICS system until the final observing

<figure><img src="image_2.png" /><figcaption>Fig. 1. Expected level of the corona (model of Cox 2015) and straylights for the small (SIO) and the large (LIO) internal occulter</figcaption></figure>

plane, with a spatial sampling of 0.7497 × 0.7497 arcsec. This model takes into account the Fresnel diffraction of the external occulter mounted on the Occulter Spacecraft and its propagation across the Lyot-style coronagraph carried by the Coronagraph Spacecraft. The numerical model makes it possible to compute the resulting PSF for each point source, inside the solar disc or in the corona, and by integrating the stray light attributed to the entire solar disc. In this section, we assume that the complete coronagraph system is perfect. The effects for an imperfect system is analyzed in Sect. 6.

In the present study, we perform the calculations for two values for the radius of the IO, corresponding to 1.0206 and 1.0759 times the size of the EO image (physical values of 1.662 mm and 1.752 mm inside the experiment). The larger the diameter of the IO, the better the straylight rejection is, with the drawback of preventing the observation of the corona very close to the solar limb. In the following, we refer to these IOs as the Small IO (SIO) and the Large IO (LIO). The inversions of the Fredholm transform are made on rectangular arrays of 128 (transverse) by 236 (radial) points (96 × 177) arcsec, both for the object and the evanescent PSFs. The radial dimension is so to cover the two different vignetting areas of LIO and SIO. The matrix **H** are of size 30 208 × 30 208, since an image of dimension $N \times M$ requires the use of a matrix of dimension $(N \times M)^2$. The results for SIO given in the present paper have never been published before.

### 4.1. Stray light, **b**

Stray lights corresponding to the SIO or LIO internal occulter for the perfect coronagraph are given in Fig. 1. Stray light curves exhibit the characteristic two-bump shape that is observed in all studies for ASPIICS. As expected, the stray light is lower for LIO than for SIO. The maximal value is $1.3210^{-6}$ at $\rho = 1.0542$ (in units of solar radius) for SIO, and $1.112710^{-7}$ at $\rho = 1.1112$ for LIO. We note that in a small region, the second bump of straylight of LIO overcomes the stray light of SIO.

To compare the stray light with respect to the corona, the expected theoretical level of the K-corona in white light is plotted in the same figure. It is drawn following the model of Cox (2015).

### 4.2. Evanescent PSF and the matrix **H**

The combined effect of the external and internal occulters produces a vignetting on the telescope aperture of the Lyot coronagraph. Since the internal occulter (IO) is placed on the image
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Radial cuts of a few PSFs computed for SIO (blue) and LIO (red). Point sources are distant one another of $\Delta\rho = 7.810^{-3}$ Rsun, or about 7.5 arcsec. Envelopes (black curves) correspond to Strehl ratios shown in Fig. 3.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Integrated flux $\Phi_G$ (geometrical model) and $\Phi_D$ (end to end diffraction computation) for the small (SIO) and the large (LIO) internal occulter. Corresponding Strehl ratios $S_D$ and $S_G$ are also given. At this scale, geometric and diffraction curves are almost identical.</figcaption></figure>

of the external occulter (EO) in the coronagraph and is larger than it, we can consider, as a first geometric approximation, that the IO defines the vignetting domain. In that region, the PSF evolves both in shape and integrated intensity, ultimately disappearing close to the limb. Figure 2 represents radial cuts of a few PSFs for the two cases SIO and LIO. For SIO, the PSF vanishes below roughly 1.04 solar radius, while the same phenomenon occurs already at 1.1 solar radius for LIO. In this representation, only a few PSF are drawn, for the sake of clarity. A 2D image of a PSF will be shown in Appendix A for comparison between perfect and imperfect optics.

Figure 3 gives integrated quantities for the small (SIO) and the large (LIO) internal occulter. The integrated flux correspond to the integral of the PSF, variable from zero to 1 in the vignetting zone, this integral being normalized to 1 outside this region. In the same graph pseudo Strehl ratios are given, also variable from zero to 1 in the vignetting zone. This Strehl ratio was already drawn in Fig. 2 as the maximal value of the PSF, being normalized to 1 outside the vignetting zone. We note that these two representations correspond to different

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Expected level of the corona (denoted by C) and observed intensities obtained as the product $C \times \Phi_D + SL$, with $\Phi_D$ the integrated flux (end-to-end diffraction computation) and $SL$ is the solar straylight, for the small (SIO) and the large (LIO) internal occulter.</figcaption></figure>

normalizations of the PSF (integral or maximal point normalized to 1). The exact curves obtained after propagation of waves through the entire coronagraph system are denoted with an index *D*, $\Phi_D$, and $S_D$, with *D* for Diffraction. Also, $\Phi_G$ and $S_G$ correspond to the simple geometrical model. At this scale, geometric and diffraction curves are almost identical. A more detailed analysis of the particular behavior of PSFs in the vignetting domain is given in Appendix A.

Figure 4 gives the level of the corona from the model of Cox (2015), namely, the level of this signal transmitted by the coronagraph after the Fredholm transformation and the final observed intensity after the addition of straylight, which appears as a bias in the observation. These curves are drawn in a linear scale. Denoting by *C* the true level of the solar corona, this operation is summarized as $C \times \Phi + SL$, where *Φ* and *SL* are the values for SIO and LIO (shown in Figs. 1 and 3).

These curves highlight the loss of signal near the solar limb caused by the use of LIO in comparison with SIO. On the contrary, it also highlight the much greater bias due to straylight when using the SIO occulter. This figure clearly demonstrates the fact that a compromise has to be found for the size of the internal occulter between straylight and closeness of observation near the solar limb. A zoom on the central part of matrix **H** is given for LIO and SIO (in Fig. 5, top and bottom panels, respectively).

### 4.3. Object *x*

The inversion of the Fredholm transformation has to be verified on a realistic representation of the solar corona near the limb. In their recent study, Shestov et al. (2021) based their simulation on a direct image obtained during a real solar eclipse and synthetic object that present finer details than the eclipse image, Shestov et al. (2019). The object used in the present study lies halfway between these two objects.

We used a solar corona picture obtained by D. Sabers, R. Royer and M. Druckmüller during the total solar eclipse of March 9, 2016 in Indonesia, where the fine details of the coronal structure have been enhanced by image processing. We impose the synthetic object to follow the dynamic range and radial decrease of intensity of the true corona, as done by Shestov et al. (2021) on their synthetic image. An illustration of this object is given in Fig. 6 represented directly and in polar coordinates,
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Region of 1000 × 1000 points of the 30 208 × 30 208 matrix for LIO (**top panel**) and SIO (**bottom panel**) internal occulters.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Solar eclipse observation of March 9, 2016 (thanks to Pr Miloslav Druckmüller) from which the data material is chosen. **Top panel:** direct image. **Bottom panel:** representation in polar coordinates of the highlighted region around the limb. The object **x** is a small fragment drawn from this image.</figcaption></figure>

as in DeForest et al. (2018), for example. Because of the size of the matrix **H**, a small part of this image can be used as a reference object in the present study. An example is given in Fig. 7. We took this region because it is radially homogeneous and allows a fair comparison between both internal occulters. It corresponds to a small region of 96 arcsec × 177 arcsec of the solar corona near the solar limb (from 1.0206 to 1.2041 in solar radius). Another part of the image is used further on in the paper. For the purposes of our study, these images are re-sampled with the sampling of 0.7497 arcsec per pixel used for the PSFs. The intensity levels have been modified to recover the dynamic of the corona given in Fig. 1. To do so, we forced the radial mean image to follow the model of Cox (2015) by multiplying the image by a factor dependent on ρ. Of course, other images of the corona could be used, such the ones shown in Habbal et al. (2010), Boe et al. (2020) for example, or other pictures from Druckmüller's eclipse home page.

## 4.4. Simulated observation of coronal features y

Figures 7 and 8 show the data material used for the simulations, in accordance with Eqs. (4) and (11). Figure 7 shows the images corresponding to the noiseless classical intensity for the SIO and LIO cases. Top panel of the figure corresponds to the effects of the vanishing PSF alone, while bottom panel of the figure includes the straylight, which appears very strong in the SIO case.

Figure 8 (top panel) gives an example of the observed images of Fig. 7 after photodetection, for an average of $10^4$ photons per pixel. In the bottom panel, the figure shows the same image with a subtraction of the bias. It should be noted that in this procedure, the calculation corresponds to performing a Poisson transformation of the bias (photonisation) and subtraction of its mean value. The top panel of Fig. 8 is the biased observation, bottom panel is the de-biased one. We note that in these two figures, the upper third of the images reproduce the same object in order to facilitate the comparison of the successive effects.

The subtraction of an estimate of the bias get rid of the vertical structures due to the scattered light, at the expense of a few negative values. At this high level of light, the debiased data gives an image much closer to the reality of the solar corona than the raw data.

# 5. Numerical results

This section is dedicated to the results of numerical simulations. The data material used, **y**, **H**, and **b** are those mentioned more precisely in the previous section.

## 5.1. Correction of data by simple flux normalization

A very simple compensation of the variation of the flux of PSFs due to vignetting can be applied to the debiased noisy observation of Fig. 8 (bottom panel). The result is shown in Fig. 9. Here, again, the upper third of the image corresponds to the object. The observed images are simply divided by the flux $\Phi_D$ of Fig. 3 for SIO and LIO. To prevent divisions by zero or by very small terms, a constant value (here: 0.01) is added to $\Phi_D$.

The result is that, in the vignetting region, the variation of the integrated flux of the PSF is compensated to a large extent, but
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Triple representation of the same region of the solar corona: the object, its Fredholm transformation for SIO and LIO. The solar limb is on the left. *Top panel:* without the stray light. *Bottom panel:* with the stray light (from Fig 1). Color bar is in solar units (1 at disk center).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Triple representation of the same region of the solar corona: the object and the observations for SIO and LIO. *Top panel:* raw photodetected image for an average of $10^4$ photons per pixel. *Bottom panel:* De-biased data (subtraction of an estimate of the straylight, inducing a few negative values). Color bar is given as the number of photons.</figcaption></figure>

not its variable blurring effect. The image quality sharply deteriorates with increasing vignetting effects. Moreover, the bias subtraction has made some of the data to become negative, as it can be seen in the range of the color bar (about 2000 pixels of the 30 000 pixels of the image are negative).

What we seek to obtain in the present study is to also compensate for this lack of image quality, or, in other words, to perform a reconstruction of the observed image inside the vignetting zone and recover a positive image. This fully justifies the choice of using a more sophisticated reconstruction algorithm, as presented in Sect. 2.2. The application to ASPIICS makes it possible to choose the term of regularization adapted to the object that is to be reconstructed.

## 5.2. Choice of the regularization term

The regularization term $J_2$, as its name suggests, must regularize, namely, it must prevent the appearance of high frequencies in the reconstruction during the iterative process. In addition, it must take into account the specificities of the object to be reconstructed. The solar corona is an object with smooth variations and this justifies the choice of the distance, L2. This choice of distance also guarantees the convexity of the functional.

More specifically, we chose for $J_2(x)$ a discrepancy function between the object, $x$, to a smoothed version of it called $A_{Rx}$, which corresponds to the effect of a perfect circular aperture of

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Correction of data by simple flux normalization: debiased observations of Fig. 8 are divided by the flux $\Phi_D$ for SIO and LIO (given in Fig. 3). Correction factors are limited to 100 to prevent division by zero. Some values are negative. *From top to bottom:* object, the results for SIO and LIO.</figcaption></figure>
radius, R. At each iteration, the matrix, $\mathbf{A}_R$, operates the con-
volution of $\mathbf{x}$ with an Airy pattern corresponding to the PSF of
a telescope of radius R. Here, R becomes a parameter that has
to be defined. After a few tests, we found that the value of R =
2.5 cm was suitable. This corresponds to the entrance aperture
of the ASPIICS coronagraph. However, we did not carry out an
exhaustive study on the value of this parameter. This regulariza-
tion term prevents the amplification of high spatial frequencies
beyond the cut-off of the perfect instrument.

We further give a particular and complex structure to $J_2(\mathbf{x})$ by including the Strehl, $S_i$, in the regularization function, as follows:

$$
J_2(\mathbf{x}) = \sum_i S_i (\mathbf{x}_i - [\mathbf{A}_R \mathbf{x}]_i)^2, \quad (28)
$$

and the gradient is given by:

$$
[\nabla_x J_2]_r = -S_r [(\mathbf{A}_R + \mathbf{A}_R^T)\mathbf{x}]_r + S_r [(I + \mathbf{A}_R^T \mathbf{A}_R)\mathbf{x}]_r. \quad (29)
$$

Here, $S_i$ plays a complementary role to the regularization fac-
tor $\gamma$ of Eq. (7). The product $\gamma S_i$ weights the strength of the
regularization according to the image quality. The regularization
coefficient then decreases from $\gamma S_i = \gamma$ in the external region
where $S_i = 1$ toward $\gamma S_i = 0$ in the shadowed region close
to the solar limb where $S_i = 0$. There, the algorithm becomes
unregularized.

In practice, the maximum value of each column of **H** (hence,
of each PSF) gives the quantity associated with the generalized
Strehl ratio described in Sect. 4.4, which is represented in Fig. 3.
The Strehl ratio, as with the PSF, is the same for a given radial
distance to the solar edge. Therefore, in the matrix formulation,
we denote using $S_i$ the value of the Strehl ratio in the column $i$.
We note that due to the lexicographic approach, $S_i$ is periodic of
the period, $N$.

We can see that in a region characterized by very strong
vignetting, the PSFs are dominated by diffraction effects and
present the surprising behavior shown in Fig. A.3 (as already
discussed in Sect. 4.2). These PSFs are no longer centered on
their corresponding point sources and their responses are shifted
towards areas of less vignetting (see Fig. A.3). This aspect
and related consequences are further discussed in Sects. 5.4
and 5.5. Simulations that are not presented here indeed con-
firm the efficiency of such a regularization by including the
variable term $S_i$ in comparison with a fixed term $S_i = 1$, for
example.

5.3. Implementation of the algorithmic reconstruction

The reconstruction is then performed according to the algo-
rithm of Eq. (9). The regularization used is that described in
the Eq. (28). Thus, to speed up the calculation at each iteration,
$k$, the convolution of $\mathbf{x}$ by an Airy function (the term $\mathbf{A}_R\mathbf{x}$) is
performed by spatial filtering using matlab 2D discrete cosine
transform dct2.

The quality criterion chosen is the reconstruction error, $\epsilon^k$, calculated by:

$$
\epsilon^k = \frac{\|\mathbf{x}_t - \mathbf{x}^k\|}{\|\mathbf{x}_t\|}, \quad (30)
$$

where $\mathbf{x}_t$ is the true object. The calculation of the error is
restricted to areas where $\Phi_D$ is not too small, leading to two
different areas for LIO (from 1.1050 to 1.2034 solar radius) and
SIO (from 1.0472 to 1.2034 solar radius).

In the case of early stopping, that is, without regularization,
$\gamma = 0$, the chosen reconstruction object $\hat{\mathbf{x}}$ is such that:

$$
\hat{\mathbf{x}} = \mathbf{x}_{\text{opt}}, \tag{31}
$$

$$
k_{\text{opt}} = \arg(\min(\epsilon^k)). \quad (32)
$$

In the case of a regularized procedure, $\gamma \neq 0$, the reconstructed
object is given by:

$$
\hat{\mathbf{x}} = \mathbf{x}_{\gamma \to \infty}, \tag{33}
$$

$$
\gamma_{\text{opt}} = \arg(\min_{\gamma} (\epsilon^{k \to \infty})). \quad (34)
$$

We draw attention to the fact that early stopping requires some a
priori knowledge of $\mathbf{x}_r$, and is therefore difficult to implement in
real cases. Conversely, the regularization procedure in its princi-
ple stabilizes the error and the algorithm can be stopped at any
time after a large number of iterations is made, without allowing
the errors to dominate.

A key problem is thus a correct choice of $\gamma$. In the simula-
tions, a numerical set of regularization coefficient is used in the
algorithm and the value chosen is the one giving the minimum
of the error as a function of $\gamma$ for large $k$, as shown in Eq. (34).

The regularization coefficient will depend on both the data
and the functions $J_1$ and $J_2$. In the present case, the value $\gamma =
10^{-6}$ is found to be the best for the Kullback functional, $\mathbf{y}_{\text{raw}}$
and the value $\gamma = 10$ is found to be the best for the Gaussian
functional, $\mathbf{y}_{\text{deb}}$. The order of magnitude difference in $\gamma$ for the
two models comes from the different expressions of $J_1$, Eqs. (14)
and (26).

5.4. Object reconstruction using raw data

The results obtained in this section are those obtained with the
raw noisy data, that is, without the subtraction of the estimate of
the bias. Furthermore, Eq. (14) and the corresponding Poisson
model described in Sect. 3.1 are used for processing the data.

Figure 10 gives the result of the direct inversion by min-
imisation of $J_1(\mathbf{x})$, that is, using Eq. (7) with $\gamma = 0$, and
stopping the iterations when the minimum of the criterion $\epsilon^k$ of
Eq. (30) is obtained. Early stopping allows to obtain what is the
best reconstruction that can be achieved in this un-regularized
inversion.

In addition to a noisy granular structure, there is a surpris-
ing quasi-deterministic phenomenon of pseudo-fringes, which is
especially clearly visible for the LIO case. The interpretation of
the presence of these pseudo-fringes was not immediate. After
several trials and errors based on simulations, we came to the
conclusion that these pseudo fringes arose mainly due to the
very particular shape of the PSF in the area of strong vignetting,
especially, with the variation of flux of the PSF. Indeed, when a
simulation is carried out by imposing a normalization of the PSF,
that is to say a constant integral also in the vignetting zone, this
phenomenon disappears. The phenomenon of parasitic fringes
also disappears if the images are processed without the presence
of bias. These simulation tests are not transcribed here for the
sake of brevity. It can therefore be deduced from these observa-
tions that the phenomenon of parasitic fringes is linked both to
the evanescent nature of the PSFs and to the presence of additive
bias. A more in-depth study of these aspects would be very inter-
esting from a data processing point of view, but this obviously
goes beyond the objectives of this work.

The resolution for most of these disturbances was to imple-
ment the regularization, namely, to use the full Eq. (7) with
$\gamma = 10^{-6}$. The result is shown in Fig. 11. For SIO, the fringes
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Reconstructed images from the raw data given in Fig. 8 (top panel), Poisson model (KL divergence of Sect. 3.1, $\gamma = 0$, early stopping). *From top to bottom:* object and the best reconstructions for SIO and LIO.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Reconstructed images from debiased data of Fig. 8 (bottom panel) based on the Gaussian model of Sect. 3.2.2, with $\gamma = 0$, early stopping. *From top to bottom:* object and the best reconstructions for SIO and LIO.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Regularized reconstructed images from the raw data given in Fig. 8 (top panel), Poisson model (KL divergence of Sect. 3.1, $\gamma = 10^{-6}$). *From top to bottom:* object and the best reconstructions for SIO and LIO.</figcaption></figure>

have completely disappeared, while for LIO they remain partially present, but for a zone outside the domain of exploitation of these observations, namely, below the geometric cutoff due to vignetting.

## 5.5. Object reconstruction using debiased data

We have proposed two models for the debiased data: a Gaussian Poissonian model (Sect. 3.2.1) and a pure Gaussian model (Sect. 3.2.2). In practice, the implementation of the first model leads to a digital dead end, due to the very large number of photons; in particular, $(H\mathbf{x})_i$ and $n_i$ can take values of several thousand photons, so the terms $(H\mathbf{x})_i^{n_i}/n_i!$ cannot be evaluated numerically. We therefore had to abandon this for the purposes of reconstructing the solar corona. Consequently, the results

obtained in this section are those obtained with the Gaussian model described in Sect. 3.2.2.

As for the raw data, $y_{\text{raw}}$, Fig. 12 gives the result of the direct inversion by minimization of $J_1(\mathbf{x})$, that is using Eq. (7) with $\gamma = 0$, and stopping the iterations when the minimum of the criterion $\epsilon^k$ of Eq. (30) is obtained. As already indicated for the Poisson case, this simulation assumes the knowledge of the object $\mathbf{x}_f$, which makes it possible to obtain the best reconstruction in this non regularized inversion. Unlike what is obtained in the previous processing based on raw data, no parasitic fringes are present here and we find that these non-regularized results are already quite promising.

Figure 13 gives the result obtained by applying a regularization with a weight $\gamma = 10$. As noted previously, the optimal numeric values of $\gamma$ are those which give the smallest error, as the number of iterations increases, ($k = 1000$ here). These results are also excellent. We can note a surprising appearance of very weak fringes for LIO in the completely obscured zone where no measurement can be expected.

In Fig. 14, another part of the solar corona is processed with the debiased data (see Sect. 3.2.2). Reconstructions are given for LIO (Fig. 14, top panel) and SIO (Fig. 14 bottom). Again, the use of SIO allows for a reconstruction of the corona closer to the solar limb.

## 5.6. Error analysis as a function of $\rho$

A quantitative estimate of the reconstruction error is given in Fig. 15. For each point of the reconstructed image, an error taking the form of Eq. (30) is computed and then its average across the image is drawn as a function of the radial position, $\rho$.

Figure 15 gives the errors corresponding to the LIO (top panel) and SIO (medium panel) internal occulters. In each figure, the PSF geometrical limits of vignetting are drawn as vertical dashed lines (for LIO from 1.10093 to 1.17312 solar radius, and for SIO from 1.04249 to 1.11468 solar radius). The two regimes, below and above the lower geometrical limits, clearly appear.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Regularized reconstructed images from debiased data of Fig. 8 (bottom panel) based on the Gaussian model of Sect. 3.2.2, with $\gamma = 10$. *From top to bottom:* object and the best reconstructions for SIO and LIO.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Result obtained in another part of Fig. 6. Regularized reconstructed images from debiased data (Gaussian model of Sect. 3.2.2, $\gamma = 10$). *Top panel:* for LIO; *top to bottom:* object, data, and reconstructed object. *Bottom panel:* for SIO, with the same representations.</figcaption></figure>

Both algorithms, based on Poisson and Gauss statistics, give better results than the simple correction, which is only a satisfactory minimum. The Gaussian model is an approximate model, but it performs a bias subtraction before processing, giving,

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Error as a function of $\rho$, for LIO (top panel) and SIO (medium panel). Here, SC: simple correction; REG and NR: regularized and non-regularized algorithms; raw (red curves) and deb (blue curves): raw data and data after bias subtraction, respectively. Vertical dashed lines represent the geometrical limits of vignetting. *Bottom panel:* zoom on the error of the above figures. Raw data: red and magenta curves. Bias subtracted data: Blue and cyan curves. See more details in the text.</figcaption></figure>

better results than the exact Poisson statistic. In particular, it does not present the problem of fringes which seems, at least in part, to be linked to the presence of the bias. Despite a strong stray light, the SIO gives a much better result than the LIO, in terms of the proximity of observation to the solar limb. However, if we analyze the error very closely (as in Fig. 15, bottom panel) the LIO shows an improvement over the SIO outside $\rho > 1.14$ solar units, but it is very slight.

## 6. Inverse Fredholm reconstruction for the imperfect coronagraph

In the previous section, the ASPIICS coronagraph was assumed to be perfect for the Fredholm inversion. To discuss further, we
now consider some “defaults” from the optics or the satellite fly-
ing formation. These aspects have been analysed by Shestov et al.
(2021) and Rougeot (2020). We limit here our analysis to forma-
tion flying errors and the effect of inevitable residual roughness
of the optical surface of the telescope of the coronagraph.

6.1. Effects of formation flying errors

The novelty of the Proba-3 formation flying mission is that the external occulter is not mechanically attached to the coronagraph. While it offers of a clear advantage in terms of the stray light reduction and vignetting (Rougeot et al. 2017), it also comes with small and varying misalignment and off-pointing errors between the occulter and the instrument, due to the imperfect in-flight formation of the two satellites. The detailed description of the different errors and their respective impact on the diffraction is presented in Shestov & Zhukov (2018). As a reminder, the ideal configuration is when the entrance aperture is well centered in the umbra and the coronagraph boresight perfectly points toward the center of the sun and, thus, the center of the external occulter. Here, we consider the lateral displacements and off-pointing of the coronagraph.

Such formation errors make both the response of the corona-
graph in the vignetting region and the straylight be asymmetric.
On one side, the occulter masks more the inner field of view. The
stray light is thus weaker and the vignetting more pronounced
with respect to the ideal configuration. On the other side, it is
the opposite, as the occulter edge is located closer to the rim of
the solar disk. The direct consequence is that the bias **b** and the
matrix **H** of the PSFs responses are naturally different from the
ideal case that was previously considered. In the following, we
denote, with **H**<sub>FF</sub> and **b**<sub>FF</sub>, respectively, the new matrix and stray
light bias in the case of formation flying errors. Thus, Eq. (4)
becomes:

$$
y = \mathbf{H}_{\mathrm{FF}}x + \mathbf{b}_{\mathrm{FF}}. \tag{35}
$$

Moreover, the misalignment and coronagraph pointing attitude are time-varying and not steady, due to the periodic actuation of both satellites to maintain and control the formation (with a characteristic time of ten seconds). It is thus incorrect to assume that the new matrix, **H**<sub>FF</sub>, and the estimated bias **b**̃<sub>FF</sub> are known a priori at any given time.

To simulate such a situation, which is likely to occur in real-
ity, we consider the case where the observations are made with
positioning errors (i.e., with **H**<sub>FF</sub> and **b**<sub>FF</sub>) and the Fredholm
inversion is performed as if the observations were perfect (i.e.,
with **H** and **b**). This approach corresponds to refraining from the
“inverse crime”, as it is called in the field of signal processing
(see, e.g., Colton et al. 1998 and Wirgin 2004). Therefore, the
de-biased data of Eq. (16) now appear as follows:

$$
y_{\text{deb}} \sim \mathcal{P}(\mathbf{H}_{\text{FF}}\mathbf{x} + \mathbf{b}_{\text{FF}}) - \mathbf{b}, \quad (36)
$$

where we synthesize the data $y_{\text{deb}}$ with $\mathbf{H}_{\text{FF}}$ and $\mathbf{b}_{\text{FF}}$ which contain formation flying errors effects, from which we subtract the stray light $\tilde{\mathbf{b}}$ expected for the perfect coronagraph.

In order to perform the Fredholm inversion in the case of
such formation flying errors, we need to compute **H**<sub>FF</sub> and **b**<sub>FF</sub>.
Here, we considered an extreme but realistic case of a 5mm lat-
eral displacement between the entrance aperture and the occulter,
combined to an off-pointing of the coronagraph boresight by
5 arcsec (Galano et al. 2019). The directions of both errors have
been chosen such that their effects are cumulative. The data

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Expected level of straylight for SIO in blue, and LIO in red. The plain curves show the nominal case without formation flying errors, as in Fig.1. The dotted curves correspond to the case with the formation flying error, in the region where the straylight is weaker due to the displacement of the external occulter.</figcaption></figure>

material used here were computed with the MATLAB code of
Rougeot (2020; see Chapter 5 in particular). We note that the
straylight is no more axis-symmetric in the case of a formation
flying error: it is brighter in one part of the field of view and
weaker in the other part. For the purposes of illustration, we
only considered the latter case. Figure 16 shows cuts of the stray-
light, namely, **b**<sub>FF</sub>, which is compared to the ideal case, for both
occulters SIO and LIO. The impact is more pronounced in the
LIO case, as the relative difference reaches 17% at the stray light
peak.

Figure 17 shows the best results that can be obtained when inverting the Fredholm transform using **H** and **b** instead of the unknown **H**<sub>FF</sub> and **b**<sub>FF</sub> true values. Obviously, the results are not as good as in the case of the perfect coronagraph, but we do not see an exponential increase of errors. The effect appears more visible for SIO, in particular in the region around 1.06 solar radius. This is probably due to the more important error in the bias subtraction there. Pseudo fringes are slightly visible in y<sub>deb</sub>, and are amplified in the final reconstruction despite the use of a strong regularization term. The larger internal occulter appears more favourable in that case.

6.2. Effects of surface roughness of the telescope

The surface roughness of optical surfaces produces micro phase
defaults whose diffraction degrades the PSF structures. The
modeling of such random micro-structures is non-trivial, as
described in Rougeot et al. (2019). The same approach is used
in this paper.

The surface roughness is modeled by a 2D random func-
tion which represents micro-deviations of glass thickness with
respect to the ideal surface. Such defects are well known in optics
to produce speckles. A rough approximation is to consider that
the final degraded PSF is the sum of the ideal PSF plus a speckle
figure, the two figures interfering in terms of complex amplitude.

For the almost perfectly polished surfaces of the corona-
graph, the considered defaults are very small. The root mean
square (rms) amplitude is on the order to λ/300, but with an
extremely small correlation length of λ/2. The resulting speckle
pattern is of a very small amplitude, while extending very far
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Processing under unknown formation flying errors. Top panel: for LIO (from top to bottom), the object, the image $y_{deb}$ of Eq. (36), and the reconstructed image using **H** and **b** instead of the unknown **H**$_{FF}$ and **b**$_{FF}$ true values. Bottom panel: same details but for SIO. The wrong subtraction of the bias slightly visible in $y_{deb}$ appears to be largely responsible for the rise in errors in the final result.</figcaption></figure>

away. The speckle structure, which is very weak compared to the Airy spot, is only revealed very far from the center of the response, as was first clearly illustrated in Fig. A1 in the appendix of Rougeot et al. (2019). Indeed, the degraded PSF due to surface roughness was found to be almost identical to the perfect PSF up to an angular distance of 2 arcmin and the speckle regime becomes clearly dominant outside of 4 arcmin.

The same behavior is observed for PSFs in the vignetting region. We give an example in Fig. 18 for the LIO case and a point source situated at 1.1011 in units of solar disk. Pinned speckles are visible at the location of the airy rings, which appear to be broken; see Bloemhof et al. (2001) and Soummer et al. (2007) for more details. These speckles appear very far from the center of the response in Fig. 18 (bottom panel). It is necessary to make a comparison with the perfect case shown in the same Fig. 18 (top panel) to see them with some difficulty. The region shown in the figure is much larger (384 × 290 arcsecs) than the region used for the computation of PSFs to establish the matrix **H**.

For the order of magnitude of the microstructures considered in this analysis, these speckle effects do not modify the matrix, **H**, which is quite the same as for the perfect instrument. There

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** PSF from the point source at 1.1011 Rsun obtained for perfect optics (top panel) and optics with surface roughness expected for ASPIICS (bottom panel). The pinned speckles are only visible with difficulty far away from the pic of the response, in the region $x > 1.35$ in the image. The field is of about 384 × 290 arcsec.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** In black, left y-axis, stray lights for the perfect and surface roughness (speckle) cases. Right y-axis, integrated flux $\Phi$ and Strehl values, $S$, for the two experiments are found to be identical.</figcaption></figure>

is a noticeable effect for the straylight, **b**, which is very sensitive to the long-range behavior of diffraction originating from all point sources from the solar disk, as shown in Rougeot (2020). Here, this effect also remains quite weak in the region concerned
by the Fredholm inversion, as shown in Fig. 19. It can be seen
there that the straylight curve is slightly above that of the per-
fect instrument while integrated fluxes and Strehl for the perfect
or speckle case remain strictly identical. Because of these small
differences, we have not offered the same study here as for the
formation flying errors.

**7. Perspectives**

This section considers the application of Fredholm inversion to
other solar coronagraphs and then to exoplanet imaging using
starshades.

7.1. *Applicability to other externally occulted solar coronagraphs*

In the present work, we deal with the case of imagery in the vignetting area of ASPIICS. The C2 and C3 coronagraphs from SOHO-LASCO, which are externally occulted coronagraph with strong vignetting, Llebaria et al. (2006), could be treated in a similar way. The problem will be to obtain a satisfactory estimate of **H** and **b**. Their computation will be more delicate than for ASPIICS for two reasons. One is the presence of the support of the external occulter in the field. The other is the short distance between the occulters and the telescope aperture, leading to very large Fresnel numbers, a difficulty analysed from a theoretical and numerical point of view in Aime (2020) and in Rougeot & Aime (2018). A study could also be done for the new concept of coronagraph METIS aboard the Solar Orbiter ESA-NASA observatory, Antonucci et al. (2020), but a complete study of the propagation would need to be carried out.

It is also possible to envisage a similar study for a coro-
nagraph using an external serrated occulter. Such an occulter
have long been considered for its qualities of straylight rejec-
tion, Newkirk & Bohlin (1965). According to Rougeot & Aime
(2018), to be efficient, an external serrated occulter must have
a very large number of teeth; typically, a 1000-tooth occulter is
expected to provide a gain of 100 in the darkness of the shadow.
Such a free-flying occulter would produce an angular period-
icity of the impulse response of the order of a thousandth of
a circumference, that is a periodicity of barely 6 arcseconds
around the solar limb, a structure that the matrix, **H**, con-
sidered in the present work could handle numerically without
difficulty.

7.2. Inverse Fredholm reconstruction for exoplanet imaging
using starshades

Another totally different field of application may be the detec-
tion of exoplanets and the study of circumstellar environments.
Indeed, if we consider that circumstellar fields may be very com-
plex and structured, the problem of imaging of exoplanets does
not simplify to imaging an ensemble of point sources, and the
inversion of the Fredholm transformation may be of interest there
for observations using starshades, Cash (2006). The same prob-
lem of an evanescent PSF arises, in a rather complex manner.
In this case, the shaped occulters have a low number of petals
(for example 20 in the computation of Cady 2012) with a par-
ticularly studied profile to obtain a very dark shadow, Vanderbei
et al. (2007), Flamary & Aime (2014). These occulters are typi-
cally of tens of meters to tens of thousands of kilometers in front
of a telescope of a few meters in diameters (see Arenberg et al.
2021 for a very pedagogic discussion of this issue).

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

**Fig. 20.** Rough illustration of vignetting for exoplanet imaging using starshade. *Top-left panel:* example of an external occulter with 12 petals and a shaped profile following a Sonine curve of the form $(1 - \rho^2)^2$. The background image gives its Fresnel diffraction pattern at the level of the pupil of the telescope. The discs give the illumination on the pupil for 3 directions of point sources in the sky (in fact the pupil is fixed and the shadow moves according to the direction of observation, but the relative effect is the same). *Bottom panel:* PSFs corresponding to the three above directions on the sky.

Although the physical characteristics are very far from the
problem of ASPIICS studied in this paper, it can be very
similar for the numerical aspects. The use of occulter with a
small number of petals will lead to a very structured vignetting.
The calculus has to be extremely precise for the straylight (the
term **b**), but probably less demanding for the vanishing PSF and
the matrix, **H**. If the external occulter rotates, as in the Starshade
Rendezvous Mission project, Seager et al. (2019), and if the
exposure time is so to include an integer number of rotations,
the observation recovers a radial symmetry similar to that of
ASPIICS. For a non-spinning occulter, there will be an angu-
lar periodicity depending on the number of petals. In both cases
the main problem will remain a precise estimation of **H** and **b**.

A rough illustration of the principle is made in Fig. 20,
showing the results obtained by a simple Fresnel filtering. We
present a shaped external occulter, with its Fresnel diffraction
at the level of the telescope aperture, the light that this aper-
ture would receive for three directions of source points around
the star, and corresponding PSFs. Reversing the Fredholm trans-
formation as we simulated here for ASPIICS could presumably
improve near-star imagery close to the inner working angle of
the system. It is likely that more elaborate techniques using the
Maggi-Rubinowicz approach (Born & Wolf 1999, Cady 2012,
Rougeot & Aime 2018) will need to be used to obtain precise
results in a realistic case.
## 8. Conclusion

In this paper, we develop a general matrix method of image reconstruction for data that have undergone a Fredholm transformation. The application is made for ASPIICS, a test solar object is drawn from an eclipse observation and the effect of vignetting is made using the same procedure as in Rougeot et al. (2017). The stray light is taken into account and simulations are carried out for two sizes of the internal occulter of Lyot's coronagraph.

Signal processing and modeling aspects are deeply analyzed, in particular, the very important question of whether it is better to reconstruct the observation from raw data or from data from which an estimate of the straylight has been suppressed. We used specific algorithms for these two types of data, with a subtraction (or without) of an estimate of the stray light. A physical regularization was used in both cases. According to the study of errors we can derive a notable advantage from the procedure where a subtraction of an estimate of the straylight is done before processing the data; and this is regardless of the fact that we had to use a mere Gaussian model instead of the fundamental Poisson model. In summary, the results are excellent. Despite the fact that it produces a greater amount of stray light, the preference goes to the small internal occulter. However, this is a solution that will not be chosen in the final project, for reasons of experimental safety and tolerance to positioning errors.

The main study was carried out for the purposes of a perfect experiment, for which we know exactly the PSFs and the level of straylight. We have also processed the case were the observation and inversion were not realized using the same PSF and bias. This was done in Sect 6.1 for formation flying errors, using $H_{FF}$ and bias, $b_{FF}$, to synthesize the data and just $H$ and bias, $b$, used to reconstruct the object. Of course the results are not as good as if every thing is known, but the image reconstruction remains interesting. A study made for optical quality defects, in the form of microstructures producing speckles in the PSF, showed that these effects were rather negligible in the region of vignetting. This can be explained by the fact that the main effect, in the vignetting zone, is the geometric variation in the shape of the pupil while the wave defects intervene as second-order effects.

This work led us to define a procedure for processing future ASPIICS data. However, this study remains theoretical, as in practice, such a fine knowledge of the PSFs may not be achievable. The next treatments will have to be done based on actual observations.

We considered generalizations of this work for other solar coronagraphs using external occulters, such as SOHO-LASCO C2 or even the new coronagraph METIS. Because of the difficulties of calculation for high Fresnel numbers involved in these experiments, as well as because of a more complex geometry, the calculation of $H$ and $b$ will be difficult to obtain there. It would be simpler to consider a system similar to ASPIICS equipped with a serrated occulter as computed by Rougeot & Aime (2018), or even with two circular occulters, as studied by Aime (2020), but such projects are not envisaged at the present time.

A very promising perspective on this work might be to apply the inverse Fredholm technique to imaging circumstellar fields and exoplanets, Arenberg et al. (2021). This point is presented briefly in Sect 7.2 and this is a study that we intend to develop in a future work.

**Acknowledgements.** We thank Prof. RNDr. Miloslav Druckmüller, CSc. of the Institute of Mathematics, Faculty of Mechanical Engineering Brno University of Technology, Czech Republic for the use of eclipse images.

See: www.zam.fme.vutbr.cz/~druck/. Thanks are due to Dr. Mamadou N'Diaye for interesting discussions.

## References

Aime, C. 2020, A&A, 637, A16

Aime, C., Theys, C., Rougeot, R., & Lantéri, H. 2019, A&A, 622, A212

Antonucci, E., Romoli, M., Andretta, V., et al. 2020, A&A, 642, A10

Arenberg, J. W., Harness, A. D., & Jensen-Clem, R. M. 2021, J. Astron. Teles. Instrum. Syst., 7, 021201

Armijo, L. 1966, Pacific J. Math., 16, 1

Bayanna, A. R., Mathew, S. K., Sankarasubramanian, K., et al. 2011, Exp. Astron., 29, 145

Benvenuto, F., Camera, A. L., Theys, C., et al. 2008, Inverse Prob., 24, 035016

Bloemhof, E. E., Dekany, R. G., Troy, M., & Oppenheimer, B. R. 2001, ApJ, 558, L71

Boe, B., Habbal, S., & Druckmüller, M. 2020, ApJ, 895, 123

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

Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357

Cady, E. 2012, Opt. Expr., 20, 15196

Cash, W. 2006, Nature, 442, 51

Colton, D. L., Kress, R., & Kress, R. 1998, Inverse Acoustic and Electromagnetic Scattering Theory (Berlin: Springer)

Cox, A. N. 2015, Allen's Astrophysical Quantities (Springer)

DeForest, C., Howard, R. A., Velli, M., Viall, N., & Vourlidas, A. 2018, ApJ, 862, 18

Demoment, G. 1989, IEEE Trans. Acoustics Speech Signal Process., 37, 2024

Denis, L., Thiébaut, E., Soulez, F., Becker, J.-M., & Mourya, R. 2015, Int. J. Comput. Vision, 115, 253

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

Flamary, R., & Aime, C. 2014, A&A, 569, A28

Galano, D., Bemporad, A., Buckley, S., et al. 2018, SPIE, 10698, 106982Y

Galano, D., Jollet, D., Mellab, K., et al. 2019, in 10th International Workshop on Satellite Constellations and Formation Flying

Groetsch, C. W. 2007, J. Phys. Conf. Ser., 73, 012001

Habbal, S. R., Druckmüller, M., Morgan, H., et al. 2010, ApJ, 719, 1362

Karush, W. 1939, M. Sc. Dissertation, University of Chicago, USA

Koomen, M. J., Detwiler, C. R., Brueckner, G. E., Cooper, H. W., & Tousey, R. 1975, Appl. Opt., 14, 743

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

Kuhn, H. W., & Tucker, A. W. 1951, in Second Berkeley Symposium on Mathematical Statistics and Probability, 481

La Camera, A., Schreiber, L., Diolaiti, E., et al. 2015, A&A, 579, A1

Lamy, P., Vivès, S., Koutchmy, S., & Arnaud, J. 2007, ASP Conf. Ser., 368, 271

Lamy, P., Damé, L., Vivès, S., & Zhukov, A. 2010, SPIE Conf. Ser., 7731, 773118

Lamy, P. L., Vivès, S., Curdt, W., et al. 2017, SPIE Conf. Ser., 10565, 105650T

Lantéri, H., & Theys, C. 2005, J. Adv. Signal Process., 2005, 643143

Lantéri, H., Roche, M., & Aime, C. 2002, Inverse Prob., 18, 1397

Llebaria, A., Lamy, P. L., & Bout, M. V. 2004, SPIE Conf. Ser., 5171, 26

Llebaria, A., Lamy, P., & Danjard, J. F. 2006, Icarus, 182, 281

Lucy, L. B. 1974, AJ, 79, 745

Lyot, B. 1933, JRASC, 27, 265

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

Newkirk, Jr., G., & Bohlin, J. D. 1965, IAU Symp., 23, 287

Renotte, E., Baston, E. C., Bemporad, A., et al. 2014, SPIE Conf. Ser., 9143, 91432M

Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55

Rougeot, R. 2020, PhD thesis, Université Côte d'Azur, France

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

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

Seager, S., Kasdin, N. J., Booth, J., et al. 2019, BAAS, 51, 106

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

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

Shestov, S. V., Zhukov, A. N., Inhester, B., Dolla, L., & Mierla, M. 2021, A&A, 652, A4

Soummer, R., Ferrari, A., Aime, C., & Jolissaint, L. 2007, ApJ, 669, 642

Thiébaut, É., Denis, L., Soulez, F., & Mourya, R. 2016, SPIE Conf. Ser., 9909, 99097N

Titterington, D. M. 1985, A&A, 144, 381

Vanderbei, R. J., Cady, E., & Kasdin, N. J. 2007, ApJ, 665, 794

Weisstein, E. W. 2003, https://mathworld.wolfram.com/

Wirgin, A. 2004, ArXiv e-prints [arXiv:math-ph/0401050]
<figure><img src="image_22.png" /><figcaption>Fig. A.1. Zoom on the very low values of the flux for the small internal occulter (SIO), shown in Fig. 3. Below $\rho = 1.0425$ Rsun (dark point) the light is geometrically blocked by the occulter (red points) but faint diffracted light enters the entrance aperture (blue points).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. A.2. Point source at the zero geometrical transmission position ($\rho = 1.0425$ of Fig. A.1) and corresponding shifted response.</figcaption></figure>

## Appendix A: PSFs in the extreme zone of vignetting

In this appendix, we describe the behavior of PSF in the extreme vignetting zone and, in particular, the existence of a response where simple geometry would give total darkness. Figure 3 of Sect. 4.2 already shows that the flux $\Phi$ and the Strehl for diffraction computations extended over a larger zone than expected from a pure geometrical point of view.

The term "geometric" corresponds to the area of the pupil free of occultation, as shown in Bayanna et al. (2011) and Shestov et al. (2021). It is indeed a problem of intersecting circles (the telescope aperture and the external occulter). According to Weisstein (2003), we have:

$$ \Phi_G(R, r, d) = 1 + \frac{d}{\pi r^2} \sqrt{r^2 - \frac{(d^2 + r^2 - R^2)^2}{4d^2}} + \frac{1}{\pi} \cos^{-1}\left(\frac{d^2 + r^2 - R^2}{2dr}\right) - \frac{R^2}{\pi r^2} \cos^{-1}\left(\frac{d^2 - r^2 + R^2}{2dR}\right), \quad (A.1) $$

where R is the radius of the external occulter, increased by internal occulters to 1.02060 (SIO) or 1.0759 (LIO) times its size, r is the radius of the Lyot coronagraph aperture, reduced to 0.97 times its size by the Lyot stop and d as the angular radial position, $\rho$, converted in meters by a factor which corresponds, for a point at the aperture center, to the radius of the projected solar image onto the external occulter. It depends, of course, on the (variable) distance to the Sun. For the parameters of our study, this radius was 0.6718-m.

The difference between $\Phi_G$ and $\Phi_D$ is clearly visible in Fig. A.1 where PSFs vanish. The example given in the figure is for SIO, but a similar behavior is observed for LIO. As a result, diffraction makes a PSF response to exist where the geometrical calculation gives a total darkness. For example, Fig. A.1 shows that the geometric transmission becomes equal to zero below $\rho = 1.0425$ Rsun. Fig. A.2 shows the response to a point source at that exact position. In the figure, the position of the point source is marked with a white circle, and the response is surprisingly shifted to the right of that point.

This behavior is general for the whole region of strong vignetting. An illustration of that is given in Fig. A.3 that shows five positions of point sources (vertical lines) and cuts of corresponding PSFs for SIO (x-axis between 1.03 to 1.08) and five for LIO (x-axis between 1.08 and 1.14). To easily match point sources with corresponding PSFs, vertical lines are drawn with heights corresponding to the maximum value of corresponding responses and drawn in the same color. For both SIO and LIO cases, PSFs are shifted to the right (toward larger radial

<figure><img src="image_24.png" /><figcaption>Fig. A.3. Positions of point sources (vertical lines) and corresponding PSFs drawn with the same color. Responses are shifted to almost the same position (vertical green lines) for SIO and LIO, respectively.</figcaption></figure>

positions) and appear at the same angular position (green dashed lines). For example, for LIO, the five point sources are situated at positions: 1.0933, 1.0948, 1.0964, 1.0979, and 1.0995, while the responses all have their maximum at about 1.11, in units of solar radius. For SIO, the same behavior is present, shifted by about 0.0156 solar unit.

These PSFs, that exist only because of diffraction effects of the external occulter, are very weak, thus they have little impact on the final image. However these weak responses have a noticeable effect for the inversion of the Fredholm transform. Indeed, since there is no more bijection between responses and point sources, it is difficult to recover the direction from which an intensity has come.

A probable explanation for this behavior may be that The Fresnel diffraction by the occulter induces effects on both the modulus and the phase of the diffracted wave at the entrance aperture of the coronagraph. In the bright region of the Fresnel diffraction, the phase variation is very small, while it rotates many times $2\pi$ in the shadowed region. This phase effect produces a tilt of the wave front. The light is directed toward the centre of the pattern. This effect is at the origin of the well known spot of Arago behind the occulting disk. As the pupil of the coronagraph enters into the shadowed region, the diffracted wave gets more and more affected by this phase variation, which amplifies the effect of the tilt.