## 1.

## Introduction

The pupils on the Extremely Large Telescope (ELT)-generation telescopes are inevitably segmented and partially shaded by thick support structures due to their large sizes. Due to small subapertures, these secondary mirror support structures, also known as spider legs, lead to fragmentation of the pupil into several disjoint segments and cover partially or even completely up to several subapertures. If a wavefront sensor provides only local information about the wavefront, as the Shack–Hartmann (SH) does, differential pistons between the pupil fragments are not seen by the sensor (are in the null space), and therefore, cannot be controlled.^{1} The pyramid wavefront sensor (PWFS),^{2} due to its increased sensitivity, has been chosen as part of many instruments currently under the development for ELT-sized telescopes. In contrast to the SH sensor, the pyramid sensor does have a footprint of the segmented piston in its data (see Sec. 5.2.2 and Fig. 5). This suggests that the reconstruction of such modes is possible, if the reconstructor is able to extract this information from the sensor data and efficiently use it.

This paper pursues twofold aims. First, various existing reconstructors are analyzed and compared with respect to their ability to provide a stable control of the fragmented piston modes. This analysis exposes a serious lack of tools for the mentioned task, which makes the majority of available reconstruction algorithms helpless in the presence of telescope spiders. Therefore, the second aim of this paper is to provide such a tool.

Let us drop a few words on our framework. First of all, we distinguish between the so-called low-wind effect (LWE) and pupil fragmentation problem. The LWE is a real low-order distortion in the wavefront caused by the heating of the air in the vicinity of spider legs in particular low-wind conditions as observed on the Very Large Telescope.^{3} Pupil fragmentation is, on the contrary, a reconstructor-related phenomenon induced by discontinuity of sensor data. In this paper, we are not considering the LWE, but we focus on the differential pistons related to reconstruction methods only. The latter is caused by the inability of the reconstructor to provide an accurate solution of the wavefront based on fragmented data. Such a failure may have two origins. One reason is the inability of the wavefront sensor to provide the measurements of the jumps between the wavefront parts on pupil fragments, as is in the case with the SH sensor. Or it may as well be a deficiency in the reconstructor itself, which is not able to use in a proper way this kind of information that is though provided by the sensor. Moreover, in this paper, we focus on solutions to the pupil fragmentation problem, which do not require bringing any changes into the mechanical setup of the system, i.e., introducing any additional dedicated focal plane sensors that would measure the fragmented low-order modes during the AO loop, as suggested for the LWE in Refs. 45.6.–7.

Let us consider the pupil fragmentation problem with PWFSs. We demonstrate in this paper that even if this type of sensor is able to provide the information on the jumps between the fragmented wavefronts, its usage does not automatically lead to correct reconstructions of segmented pistons with an arbitrary reconstruction algorithm designed for nonfragmented annular apertures. Based on numerous end-to-end simulations performed with various reconstruction algorithms (summarized in Secs. 3.5 and 4.2), we draw a conclusion that many of the known methods simply fail, meaning that the residuals contain randomly appearing uncontrolled piston modes on the segments. This leads to a significantly reduced correction quality in terms of Strehl ratio, PSF, and contrast. In addition to pupil fragmentation, another challenge that wavefront reconstruction algorithms for the ELTs have to tackle with is the high number of correcting elements that need to be controlled in real time.

Having these two points in mind, in Secs. 3 and 4, we give an overview of the available algorithms and their readiness/ability to operate both fast and stable with high-quality on ELTs. The analysis is performed in the context of the mid-infrared ELT imager and spectograph (METIS) instrument^{8} on the ELT containing a single conjugate adaptive optics (SCAO) system with $74\times 74$ WFS subapertures. The geometry of the planned ELT M4 deformable mirror (DM) is taken into account. The thickness of the six secondary mirror support structures is 50 cm and coincides with the subaperture size.

On the one hand, as shown in Sec. 3.5, the current status of the performed simulations shows that among all possible variants of the interaction-matrix-based matrix–vector multiplication (MVM) methods only one, namely, the minimum variance reconstructor using a zonal interaction matrix and Laplacian regularization, is able to provide a stable control of the segmented low-order modes. However, the well-known drawback of any MVM method is the related computational load. Though for the considered SCAO system, the real-time application of MVM is still doable, it is hardly feasible in case of the planned extreme adaptive optics (XAO) system having tens of thousands of correcting elements to control.

On the other hand, as shown in Sec. 4, there exists a variety of interaction-matrix-free, model-based wavefront reconstructors developed for PWFSs in the recent decade, all of them being computationally much more efficient than interaction-matrix-based approaches. However, since those methods were developed such that they are intrinsically using the approximate forward models of the sensor not including segmented pupils, these methods fail when being applied straightforwardly to segmented sensor data. In Sec. 4.2, we additionally mention several attempts of adapting these methods for wavefront reconstruction on segmented pupils that unfortunately do not yield the expected performance. Therefore, in the presence of spiders, model-based wavefront control algorithms need to be adapted in some sophisticated way in order to handle the differential pistons between the pupil segments.

As a solution allowing for both high-quality and high-speed wavefront control, we suggest in Sec. 5 a hybrid scheme. This approach combines the advantage of the interaction-matrix-based methods of being able to handle pupil segmentation and the advantage of the model-based reconstructors of being fast. The solution, named the split approach, treats the reconstruction of segmented pistons separately from the higher-order modes (or frequencies). Here, the piston-free wavefront reconstruction on segments is provided by some model-based algorithm, e.g., the P-CuReD,^{9}^{–}^{11} as described in Sec. 5.1. In parallel, the segment pistons are reconstructed from the same sensor data with an interaction-matrix-based MVM approach. Since for the direct segment piston reconstruction we are only interested in the modes of order zero, the computational load can be significantly reduced.

In Sec. 5.2, we demonstrate two possibilities toward the formulation of the direct piston reconstruction for segmented pupils. The first one employs the usual setting of the full zonal interaction matrix using a set of dedicated basis functions representing the wavefront. This big matrix is then inverted via standard techniques of regularization and the resulting intermediate control matrix is afterward reduced to a small sized, control matrix relating the sensor data with the vector of segment pistons having only as many entries as pupil segments. Though still requiring the computationally expensive setting up and inversion of a dense matrix, which can be performed offline, the online calculations have linear complexity $O({n}_{a})$ and are very cheap. In the second approach for direct segment piston reconstruction, the initial interaction matrix is formulated in the basis of segment pistons and is, therefore, very small from the start. Resulting in the same number of computations to be performed online, this approach is, in addition, free from the time-consuming (offline) operations involved in the first approach. Sec. 5.4 contains technical details on the implementation of the split approach.

In combination with the P-CuReD, both direct segment piston reconstruction methods have an overall computational effort, which scales linearly. In Sec. 6, we illustrate the performance of the split approach employing the two proposed methods for direct piston reconstruction for segmented pupils with end-to-end closed loop simulation results.

## 2.

## Preliminaries

Let us denote the model of the pyramid sensor by an operator $\mathit{P}:{\mathcal{L}}_{2}({\mathbb{R}}^{2})\to {\mathbb{R}}^{2n}$, which maps real-valued ${\mathcal{L}}_{2}$ functions (wavefronts and residual phases) to a vector of discrete measurements of length $2n$. The measurement process is given by

where $\mathrm{\varphi}$ describes the incoming phase, $\overrightarrow{s}$ is the pyramid sensor measurements, and $\overrightarrow{\eta}$ is the noise in the data. The inverse problem is to reconstruct the wavefront $\mathrm{\varphi}$ from given noisy sensor data $\overrightarrow{s}$. Throughout this paper, we assume the pyramid sensor to operate in closed loop AO to be linear.For a detailed description of the advantages of the PWFS over other sensor types and an overview on the instruments and telescopes, in which this type of sensor is already installed or is planned to be used in the nearest future, we refer the reader to Ref. 12. Here, we only mention briefly that due to the global response, the PWFS is able to sense the differential pistons, which has been successfully demonstrated in the laboratory,^{13}^{,}^{14} supported by numerical simulations, and validated on sky under seeing-limited conditions.^{15} The ability of the PWFS to sense the differential pistons of a segmented mirror and correct for it with an inversion based on the singular value decomposition (SVD) of the measured interaction matrix was first demonstrated in numerical simulations in Ref. 16. Among all WFS types tested in Ref. 17, the PWFS takes the most sensitive measurements of the differential pistons on the segments. Apart from that, compared to the SH WFS, the PWFS provides an increased sensitivity, which leads to higher limiting GS magnitudes and higher sky coverage.^{18}

The reconstructors we will be dealing with in Secs. 3 and 4 are using different forward models of the PWFS or its approximations. Here, for the sake of brevity, we intentionally omit describing the forward models of the PWFS since we focus mainly on the achieved performance of various methods under pupil segmentation. For a detailed description of different model approaches and approximations of the WFS data we refer the reader to the devoted works.^{9}^{,}^{10}^{,}^{12}^{,}^{19}^{,}^{20}

Note that the interaction-matrix-based MVM methods are described in Sec. 3 with a much higher level of detail compared to the model-based algorithms in Sec. 4. This is done intentionally by the authors who contributed to the development of model-based algorithms, which turned out to struggle under the circumstances of pupil fragmentation. The focus of this paper is to study, compare, and understand the behavior of different algorithms under pupil fragmentation. As a result of these efforts, this paper serves partially as a review of the currently available wavefront reconstruction algorithms for PWFSs and, in particular, of the present status in the performance those reach on the ELT-era instruments under design.

The connection between incoming wavefronts $\mathrm{\varphi}$, residual wavefronts ${\mathrm{\varphi}}_{\mathrm{res}}$, and the mirror shape $\phi $ in an AO system is given by

Since for ideal compensation the residual wavefront ${\mathrm{\varphi}}_{\mathrm{res}}$ should be equal to zero, the optimal choice for the DM shape is:

Therefore, for the control of DMs, one needs to know either the mirror actuator commands or the shape of the incoming wavefront provided as the solution of the inverse problem (1) for $\mathrm{\varphi}$.

## 3.

## Interaction-Matrix-Based Reconstructors

In this section, we analyze the applicability and performance of the so-called interaction-matrix-based MVM methods for PWFS data and fragmented pupils. These methods involve a registration (or computation) of a WF-to-WFS interaction matrix, its inversion and a subsequent multiplication of the obtained control matrix with a vector of sensor data measurements. In the literature, there have been many variants of the interaction-matrix-based MVM approach presented: statistical estimation or solution in a least-squares sense; zonal or modal control approaches (i.e., the degrees of freedom are actuators/subapertures or modes). The presented overview aims to summarize and compare the performance of the existing interaction-matrix-based MVM methods in case of PWFS data fragmented by spiders. In Sec. 3.1, the generation of a WF-to-WFS interaction matrix is described, and the option of coupling or decoupling this step with the DM is explained in Sec. 3.2. In Sec. 3.3, the simplest least-squares approach and its regularized inversion are specified, whereas Sec. 3.4 deals with more sophisticated statistical approaches. Finally, the quality and speed performance of the interaction-matrix-based MVM algorithms are summarized in Sec. 3.5.

## 3.1.

### Generating the Interaction Matrix

We introduce $({h}_{i})$ as a set of arbitrary basis functions to represent the wavefront, $({h}_{i}^{m})$ as a set of modal/global basis functions, and (IF) denotes the DM influence functions. In order to create the interaction matrix of the system, we need to relate the incoming wavefront with the output (measurements) of the PWFS. We represent the incoming wavefront $\mathrm{\varphi}$ using a set of basis functions $({h}_{i})$. Thus $\mathrm{\varphi}$ can be approximated by

where ${n}_{c}$ indicates the number of used basis functions. The interaction matrix $\mathit{M}\in {\mathbb{R}}^{2n\times {n}_{c}}$ is then given by## Eq. (5)

$$\mathit{M}=(\begin{array}{cccc}{\overrightarrow{s}}_{1}& {\overrightarrow{s}}_{2}& \cdots & {\overrightarrow{s}}_{{n}_{c}}\end{array}),$$## Eq. (6)

$${\overrightarrow{s}}_{i}=\mathit{P}({h}_{i})\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{for}\text{\hspace{0.17em}\hspace{0.17em}}i=\mathrm{1,2},\dots ,{n}_{c},$$The sensor equation, as already mentioned, reads as

To reconstruct the incoming (residual) wavefront $\mathrm{\varphi}$, the matrix $\mathit{M}$ has to be “inverted” and applied to the measurements represented by

with $\overrightarrow{c}={({c}_{i})}_{i=1,\dots ,{n}_{c}}$.After the reconstruction step, one has to derive the actuator commands $\overrightarrow{a}={({a}_{i})}_{i=1,\dots ,{n}_{a}}$ from the reconstruction $\mathrm{\varphi}$, i.e., solve

## Eq. (7)

$$\mathrm{\varphi}(x,y)=\sum _{j=1}^{{n}_{c}}{c}_{j}{h}_{j}(x,y)=\sum _{j=1}^{{n}_{a}}{a}_{j}{\mathrm{IF}}_{j}(x,y).$$If the chosen basis ${h}_{i}$ coincides with the DM influence functions, the vectors $\overrightarrow{c}$ and $\overrightarrow{a}$ are equal.

## 3.2.

### Working with DM Influence Functions

Often, the interaction matrix inversion is coupled with the DM in the sense that for the generation of an interaction matrix one creates a certain (zonal or modal) shape with the DM, which is then sensed by the wavefront sensor. In this approach, one is restricted to wavefront shapes, which can be represented by the DM, i.e., are a linear combination of the DM influence functions:

or using DM modes with with actuator commands $({m}_{l}^{j})$. This results in## Eq. (10)

$$\mathrm{\varphi}(x,y)=\sum _{j=1}^{{n}_{c}}{c}_{j}{h}_{j}^{m}(x,y)=\sum _{j=1}^{{n}_{c}}{c}_{j}\sum _{l=1}^{{n}_{a}}{m}_{l}^{j}{\mathrm{IF}}_{l}(x,y).$$Combining Eq. (1) with Eq. (8) or Eq. (9) produces a DM-to-WFS interaction matrix, which relates the sensor measurements $\overrightarrow{s}$ directly with the command vectors:

## 3.3.

### Deterministic Setting

## 3.3.1.

#### Least-squares pseudoinverse

The least-squares problem

Such a pseudoinversion of the interaction matrix is considered to be the simplest reconstruction algorithm possible. Historically, least-squares with zonal representation was the first approach applied to the wavefront reconstruction problem involving the SH wavefront sensor.^{21} For wavefront reconstruction using PWFS measurements, the least-squares approach allows one to reach high correction accuracy without regularization, at least if the number of degrees of freedom is not very big. With the PWFS, the least-squares reconstructor has provided reasonable results on small scale systems having up to $30\times 30$ subapertures within an 8-m telescope diameter as it was demonstrated in Ref. 22 with modal DM control using Karhunen–Loève (KL) polynomials.

## 3.3.2.

#### Regularized least-squares pseudoinverse

For large-scale AO systems or more sophisticated configurations like multiconjugate adaptive optics (MCAO), the corresponding system matrices have large condition numbers and are difficult to invert. In this case, the conventional least-squares reconstructor performance is not satisfactory, and special treatment is required in the form of a regularization or filtering of unstable modes.

Sensor noise is modeled as a random process obeying zero-mean Gaussian statistics, $\eta \sim \mathcal{N}(0,{C}_{\eta}={\sigma}^{2}I)$, ${C}_{\eta}\in {\mathbb{R}}^{2n\times 2n}$, where ${\sigma}^{2}$ denotes the sensor noise variance. The noise-covariance-weighted least-squares (also known as minimum-norm maximum likelihood^{23}) reconstructor, which minimizes

The pseudoinverse can be computed using the eigendecomposition of ${\mathit{M}}^{\mathrm{T}}{C}_{\eta}^{-1}\mathit{M}$ or the SVD of ${C}_{\eta}^{-1/2}\mathit{M}$. However, in practice, such decompositions are related with expensive computations and difficulties in computing the small singular values with the desired numerical precision. Also, small computation errors in eigen- or singular values lead to instabilities in the reconstruction due to noise amplification. As a regularization method, in the truncated SVD, one filters out the modes corresponding to singular values smaller than a given parameter $\alpha >0$.

Alternatively, a similar effect is achieved by the so-called Tikhonov regularization, well-known in the field of inverse problems.^{24} For solving the inverse problem (1), we consider the least-squares problem

Then the regularized pseudoinverse is derived as

## Eq. (11)

$${\mathit{M}}_{\alpha}^{\u2020}={({\mathit{M}}^{\mathrm{T}}\mathit{M}+\alpha \mathit{I})}^{-1}{\mathit{M}}^{\mathrm{T}},$$For (relatively) small scale systems having up to $30\times 30$ subapertures on an 8-m telescope, high correction quality has been achieved with both an SVD-regularized zonal^{25} and modal least-squares reconstructor using Zernike polynomials^{26} or KL polynomials.^{27}^{–}^{29}

However, due to noise propagation, the least-squares wavefront reconstruction algorithm performs poorly for large-scale or laser-guide star-based AO applications.^{30}^{,}^{31} Hence, there has been a tendency observed in the AO community to prefer the regularized variants of the interaction-matrix-based MVM method taking atmospheric statistics into account.

## 3.4.

### Bayesian Setting

Within a stochastic context, wavefront shapes and sensor noise are independent random processes obeying zero-mean Gaussian statistics, $\mathrm{\varphi}\sim N(0,{C}_{\mathrm{\varphi}})$, $\eta \sim N(0,{C}_{\eta}={\sigma}^{2}I)$, ${C}_{\mathrm{\varphi}}\in {\mathbb{R}}^{{n}_{c}\times {n}_{c}}$, and ${C}_{\eta}\in {\mathbb{R}}^{2n\times 2n}$, with ${\sigma}^{2}$ denoting the sensor noise variance. Such a point of view allows one to use in the reconstruction the prior knowledge of the atmosphere and measurement noise statistics, expressed with the corresponding covariance matrices ${C}_{\varphi}$ and ${C}_{\eta}$, in order to regularize or stabilize the solution.

Two Bayesian statistical approaches, both using a prior probability density assumed on the phase, have been applied to the problem of wavefront reconstruction from sensor data—minimum variance estimation and maximum *a posterior* (MAP) estimation. The minimum variance [or minimum mean-square error (MMSE)] estimator minimizes the variance of the phase estimation error. The MAP estimator identifies the most likely value of $\mathrm{\varphi}$ given the observed data $\overrightarrow{s}$ and prior knowledge on the distribution of $\mathrm{\varphi}$. Since wavefront reconstruction deals with zero-mean Gaussian signal and perturbation, the minimum variance reconstructor providing the minimal MSE coincides with the MAP reconstructor.^{23}^{,}^{32}

In the Bayesian setting, the measurement vector $\overrightarrow{s}$ is a function of the atmospheric turbulence profile. The sensor Eq. (1) using the representation of the wavefront Eq. (4) formulated in terms of arbitrary coefficients ${c}_{i}\in \mathbb{R}$, in the stochastic setting is formulated in terms of wavefront coefficients ${\varphi}_{i}\in \mathbb{R}$ as

The aim in this setting is to compensate the turbulence-induced wavefront error.

The minimum variance/MAP reconstructor minimizes the penalized noise-weighted least-squares functional

## Eq. (12)

$${\mathit{M}}_{\mathrm{MAP}}^{\u2020}={({\mathit{M}}^{\mathrm{T}}{C}_{\eta}^{-1}\mathit{M}+{\mu}_{0}{C}_{\mathrm{\varphi}}^{-1})}^{-1}{\mathit{M}}^{\mathrm{T}}{C}_{\eta}^{-1}.$$Note that besides the pseudoinverse, there exist several other methods based on the above normal equation for solving the inverse problem. The inverse of the phase covariance ${C}_{\mathrm{\varphi}}^{-1}$ must be chosen such that it is physically realistic. Typically, because of singularities (or ill-conditioning) in the turbulence spectra, it is inevitable to assume some discrete approximation on ${C}_{\mathrm{\varphi}}^{-1}$ and an additional regularization, which results in some loss of accuracy but yields stability.^{30}

With the MAP/minimum variance estimators, two kinds of errors are related: the approximation error that tells how well the reconstructor approximates the inverse of the sensing operator $\mathit{P}$ and the noise propagation error related to sensor noise. The sources of the model error are the chosen basis representation of the wavefront and the accuracy of the *a priori* statistical knowledge of the atmosphere.

In contrary to the least-squares approach, a statistical estimation method necessarily needs regularization parameter tuning for an accurate wavefront reconstruction. The numerical simulations indicate that the MAP/minimum variance reconstructor with an optimized parameter ${\mu}_{0}$ performs better (is more stable) than the (noise-weighted) least squares solution.^{33}^{,}^{34}

In the following, we briefly focus on the implementation details of two variants of Bayesian reconstructors that proved to be efficient in astronomical AO.

## 3.4.1.

#### MAP reconstructor with modal control of the DM

In this method, ${C}_{\mathrm{\varphi}}$ is the von Karman wavefront covariance matrix restricted to the modal space of the DM. The parameter ${\mu}_{0}$, which scales as the inverse of the square of the signal-to-noise ratio, allows one to weight the sensor noise and atmospheric priors in a flexible way.

## 3.4.2.

#### Zonal minimum variance estimator with regularized sparse approximation of ${C}_{\mathrm{\varphi}}$

Another variant of MVM that provides high-quality reconstruction is the zonal minimum variance estimator with regularized sparse discrete approximation of ${C}_{\mathrm{\varphi}}^{-1}$, as suggested in Ref. 30, with

where $\mathit{L}$ denotes a discrete Laplacian matrix approximating the Laplacian operator. The constant ${c}_{0}$ is physically interpreted as the strength of the turbulence and additionally normalizes the approximation of ${C}_{\mathrm{\varphi}}$ in order to fit the von Karman turbulence spectrum.^{30}

^{,}

^{37}

Using this covariance approximation corresponds to regularization by the ${l}_{2}$-norm of the Laplacian, i.e., to solve the penalized least-squares functional:

^{37}

This reconstructor using DM influence functions as basis functions is incorporated in YAO, an open-source AO simulation tool written in yorick.^{38} For this study, we have also implemented it in Octopus. In the numerical implementation, stabilization with respect to the wavefront sensor noise is performed by filtering out the columns in the registered interaction matrix corresponding to basis functions whose WFS response was smaller than a certain predefined value (e.g., expressed in a percentage of the maximum registered response). This parameter can be tuned in order to optimize the performance of the reconstructor for different flux settings or spiders’ thickness. Note that in the numerical simulations, the parameter ${c}_{0}$ is heuristically tuned as well for various atmosphere strengths and guide star fluxes.

## 3.5.

### Performance of Interaction-Matrix-Based MVM in the Presence of Spiders

## 3.5.1.

#### Quality

In this section, we summarize and compare the available results of various versions of the interaction-matrix-based MVM methods in the presence of spiders. First, we mention the earlier result achieved in the context of EPICS, the XAO instrument on the ELT in the case of four thick spiders.^{39} It was reported therein that both the zonal and modal reconstructors provide the same quality in the presence of four spider legs. Apart from that, two important points were underlined: the light behind spiders needs to be used in the reconstruction; and the amount of modulation should not be too large, otherwise, the sensor loses its sensitivity to low-order modes.

Next, we report in Table 1 on the performance we obtained with two variants of Bayesian reconstructors in the case of six thick spiders in the context of the METIS instrument on the ELT. The results are presented for a median atmosphere, on-axis correction and high photon flux. Note that the modal MAP result from Octopus and the zonal reconstructor results from YAO were obtained using a regular Fried geometry, whereas the zonal method results from Octopus imply the real ELT M4 geometry.

## Table 1

LE Strehl ratio obtained with three versions of the interaction-matrix-based MVM for pupils with and without spiders.

Method | Zonal MMSE (OCT.) | Zonal MMSE (YAO) | Modal MAP (OCT.) |
---|---|---|---|

No spiders, Fried geometry | 0.89^{40} | 0.89 | |

No spiders, M4 geometry | 0.89 | ||

With spiders, Fried geometry | 0.89^{40} | 0.86 | |

With spiders, M4 geometry | 0.89 |

The zonal minimum variance reconstructor, implemented both in Octopus and YAO, provides the best quality of 0.89 LE Strehl, which barely decreases for fragmented pupil masks. For more YAO simulation results obtained for the METIS case with the zonal minimum variance reconstructor, we refer to Ref. 40.

The modal MAP reconstructor (as implemented in Octopus) achieves the same LE Strehl ratio of 0.89 in the case without spiders. However, in contrary to the zonal approach, its performance in the presence of spiders is not as good. After running multiple tests over a set of regularization parameters, the best LE Strehl ratio that we were able to achieve is 0.859. In the residual screens, we always observed uncompensated random segmented pistons. Note that in order to obtain reasonable results one has to set the illumination parameter to some value lower than 0.5, i.e., to use data from the subapertures partially covered by the spiders. If one does not use this information, the reconstruction gets unstable because of the poorly corrected segmented pistons.

Similar results with the MAP reconstructor have been recently reported in Ref. 41 for an 8-m telescope. First, the authors therein underline the importance of using the light from the subapertures partially covered by the spiders, in the same way that we observed in our tests. Second, the quality of the spider-free case could not be achieved in the presence of spiders. The authors in Ref. 41 conclude that the modal MAP reconstructor as it is known is not able to overcome the fragmentation problem completely. It is presumed that a dedicated change of basis functions, which would be defined on each segment instead of the entire pupil, could potentially be helpful, but no attempts have been performed so far in this direction.

An intermediate conclusion is that only one version of the interaction-matrix-based methods, namely the zonal minimum variance reconstructor, handles the pupil fragmentation problem without almost no loss in quality.

## 3.5.2.

#### Speed

Interaction matrix inversion is related to a high computational load, which can, however, to a large extent be performed offline. Assuming that an interaction-matrix-based MVM algorithm is implemented using conventional matrix inversions and MVMs, the method requires approximately $O({n}_{a}^{3})$ operations to compute the control matrix and $O(n{n}_{a})$ operations to apply the control matrix.^{30} The computational complexity becomes a limitation of these approaches if the values of ${n}_{a}$ are of order 10,000 as currently under consideration for XAO systems on ELTs.

## 4.

## Advanced Model-Based Reconstructors

Apart from the interaction-matrix-based methods, there exists plenty of recently developed algorithms based on the mathematical analysis of the forward models of the PWFS. The main feature of these algorithms is a low computational complexity resulting in an adequate handling of wavefront reconstruction on AO systems in real time and still guaranteeing high-quality and robustness of the methods. A detailed overview and comparison of various algorithms can be found in Ref. 42. Until lately, these methods have been extensively tested with nonfragmented pupils. The question is: are they also able to provide reasonable reconstructions with fragmented pupil masks?

Keeping the structure similar to the one in the previous section, let us first start with a brief overview of the methods. Note that, for the sake of brevity, these are presented without any mathematical details, which can be found in the provided references. Then we will have a look at quality and speed these methods provide for nonfragmented pupils. Finally, we discuss the performance of the model-based methods in the presence of spiders.

We consider the following algorithms:

• preprocessed cumulative reconstructor with domain decomposition (P-CuReD)

^{9}^{–}^{11}• convolution with the linearized inverse filter (CLIF)

^{9}^{,}^{12}^{,}^{43}• pyramid Fourier transform reconstructor (PFTR)

^{12}^{,}^{43}• finite Hilbert transform reconstructor (FHTR)

^{9}• singular value type reconstructor (SVTR)

^{44}• a number of linear iterative methods applied to pyramid sensors like conjugate gradient for the normal equation (CGNE), steepest descent (SD), Landweber iteration, and Kaczmarz versions of the previous mentioned algorithms

^{45}^{–}^{47}• nonlinear Landweber and Landweber–Kaczmarz iteration.

^{48}

The first three methods, the P-CuReD, the CLIF, and the PFTR are all using Fourier domain analysis of the underlying (modulated and nonmodulated) pyramid sensor operators and their approximations. The CLIF and PFTR suggest different numerical implementations of the same idea of reconstructing the wavefront spectrum from the data spectrum using the known Fourier domain relation between both. The PFTR applies the corresponding inverse Fourier domain filter to the data spectrum, the CLIF algorithm works in the spatial domain. The P-CuReD correlates the pyramid with the SH sensor and by employing the Fourier domain relation between both, suggests a simple convolution-type conversion of pyramid sensor data to SH-like data. Subsequently, the CuReD is applied to the preprocessed pyramid data, a method that was originally developed for the SH-WFS.^{49}^{–}^{51}

The FHTR and the SVTR are based on a simplification of the nonmodulated pyramid sensor model represented as the finite Hilbert transform of the incoming wavefront. The FHTR applies a direct inversion formula of the finite Hilbert transform, the SVTR uses a version of the involved operators’ SVD for wavefront estimation. In contrast to MVM approaches, the algorithms consist of analytical derivations of the mentioned tools for inversion.

In the iterative methods, we use well-known mathematical algorithms for solving the inverse problem of wavefront reconstruction from pyramid sensor data, namely CGNE, SD method, SD-Kaczmarz method, Landweber iteration, and Landweber–Kaczmarz iteration. All algorithms are applicable to pyramid sensors with and without modulation, have comparable numerical complexity, and provide similar reconstruction quality. Outstanding are the Kaczmarz versions of the algorithms, which enable a sophisticated combination of the two data sets ${\overrightarrow{s}}_{x}$ and ${\overrightarrow{s}}_{y}$ for reconstruction. Although all mentioned model-based reconstructors first were introduced based on a linearization of the pyramid sensor model, for Landweber and Landweber–Kaczmarz iteration, we additionally considered their nonlinear versions in order to develop nonlinear reconstructors.

## 4.1.

### Quality and Speed Performance Without Spiders

For a detailed comparison of the methods with respect to the reconstruction quality, we refer to the corresponding publications. Here we only briefly mention that in closed loop simulations all the model-based methods achieve a quality close to that of interaction-matrix-based approaches with some of them even slightly outperforming the latter. In Table 2, we provide a comparison of wavefront reconstruction algorithms for pyramid sensors in terms of their computational complexities. From this table, one can see that all recently developed model-based reconstruction algorithms require much fewer computations to be performed than an MVM implementation. The P-CuReD algorithm reaches the highest quality results in a variety of tested conditions, and at the same time has the lowest, namely linear, computational complexity. In Refs. 52 and 53, AO simulation tools’ users compared the performances of the modal MAP method with the P-CuReD for XAO settings. In Ref. 53, the MVM and P-CuReD algorithm give the same reconstruction quality with only very slight discrepancies in some of the test cases. Moreover, in Refs. 9–11, it was shown that the P-CuReD provides a significantly improved quality in the low-flux cases compared to an interaction-matrix-based reconstructor, and also that it converges to high Strehl ratios faster than the tested MVM approach.

## Table 2

Comparison of the currently existing algorithms for wavefront reconstruction from pyramid sensor data in terms of their flexibility and computational complexity.

Algorithm | Modulation | Complexity | ||
---|---|---|---|---|

No | Small | Large | ||

Interaction-matrix-based MVM | ✓ | ✓ | ✓ | $O({n}_{a}^{2})$ |

FHTR and SVTR | ✓ | $O({n}_{a}^{3/2})$ | ||

Iterative methods | ✓ | ✓ | ✓ | $O({n}_{a}^{3/2})$ |

CLIF | ✓ | ✓ | ✓ | $O({n}_{a}^{3/2})$ |

PFTR | ✓ | ✓ | ✓ | $O({n}_{a}\mathrm{log}\text{\hspace{0.17em}}{n}_{a})$ |

P-CuReD | ✓ | ✓ | ✓ | $O({n}_{a})$ |

Another point we would like to stress here is that all our model-based reconstructors are free from the time-consuming precomputation of matrices and fine tuning of the regularization parameters associated with MVM approaches. Since there are no intrinsic regularization parameters necessary to consider during the reconstruction process, no optimization is needed if atmospheric conditions change. This facilitates the usage of the methods by external users, as confirmed by Refs. 52 and 53 with P-CuReD applications and Ref. 54, where the CuReD method, which was originally developed for the SH sensor^{49}^{,}^{50} and nowadays also constitutes a part of the P-CuReD method for pyramid sensor, has been successfully tested on-sky.

## 4.2.

### Adapting Advanced Reconstructors to Segmented Pupils

In this section, we describe our first approaches to overcome the difficulties caused by pupil fragmentation or spider effects when the width of the spiders is no longer negligible compared to the size of the subapertures. Numerical investigations were performed for the ELT having a primary mirror divided into six segments as shown in Fig. 1. As reconstruction method, we use the P-CuReD algorithm. Almost as important as the width of the spider legs is their placement. If the spiders are in parallel with the $x$ and $y$ axes, some basic attempts described in the following succeed. In contrast, for arbitrarily located spiders, we examine poor wavefront reconstruction suffering from differential piston influences. Note that the following described attempts for wavefront reconstruction in the presence of six telescope spiders were only investigated with application of the P-CuReD algorithm. Different reconstruction algorithms will behave differently.

One idea is to make the illumination factor necessary for the usage of subapertures for wavefront reconstruction very low, i.e., utilizing (almost) all available subapertures for reconstruction, also those being less illuminated. This method combined with the reconstruction algorithm P-CuReD does not lead to success since the light suffers from obstruction effects especially at the boundaries of the pupil segments. The wavefront reconstruction is adversely influenced on the whole pupil by differential piston effects as seen in Fig. 2 for different illumination factors.

Another approach to overcome the effects of pupil segmentation is data interpolation or interpolation of the reconstructed phase under the spider legs. Instead of using defective measurements under the legs provided by the pyramid sensor, we generate data or reconstructions artificially in these areas. For that purpose, we considered bilinear and spline interpolation. Unfortunately, these approaches do not eliminate differential piston effects.

A more sophisticated attempt incorporates the pyramid sensor model in the measurement continuation process by applying an iterative measurement extension method already presented in Refs. 44 and 55. The idea in Ref. 44 was originally developed for wavefront reconstruction on annular, nonsegmented telescope pupils. The basic concept is to generate artificial but pyramid related data under the spider legs by application of the finite Hilbert transform, which is a simplification of the Fourier optics-based pyramid sensor model. Thus the provided data correspond better to pyramid measurements as it is the case for a general interpolation procedure described above. However, the simulation results are quite similar to the approach using bilinear or spline interpolation, and differential piston effects are developing within time and make the reconstruction poor.

An additional experiment is to replace the obstructed data by zeros, i.e., it is assumed that the wavefront is planar in the areas obstructed under the spider legs. Note that zero padding is successfully used for the central obstruction induced by a secondary mirror for some of the reconstructors mentioned in this section, e.g., the SVTR. In contrast to obstruction induced by spiders, a central obstruction does not cause segmented mirrors and hence does not induce differential pistons. Several numerical simulations show that the results of this approach differ significantly for different spider leg locations. In particular, if one considers four spider legs parallel to $x$ and $y$ axes (as those first investigated for SH sensors in Ref. 1), the zero padding approach gives satisfying reconstruction quality with the pyramid sensor. Considering the six ELT legs, the approach again suffers from differential piston development as shown in Fig. 3.

Approaches like jump minimization between the segments by boundary integral coupling do bring quality improvements but not as high as we have hoped for. For SH sensors, this method will precisely be described in an upcoming paper of our group. Actuator coupling/slaving of those actuators that are situated at the boundary of spider legs is a hot topic for pyramid sensors on ELTs, especially for wavefront sensing at shorter wavelengths than K-band.

Altogether, the above-described methods do not satisfactorily handle the impact of 1 or 2 subaperture thick spiders on the wavefront correction performance. The negative influences of randomly appearing piston modes on the wavefront reconstruction gave us the motivation to develop ideas to separately examine the segmented pistons. After trying all the described methods, a way for successfully eliminating spider obstruction effects was found. The idea is to reconstruct the segmented piston modes separately from other frequencies in the wavefront. Within these attempts, which we recap as split approach, one still can use the fast interaction-matrix-free wavefront reconstruction algorithms presented in this section and at the same time obtain stable wavefront correction in the presence of spiders.

In the next section, we introduce the basic concept of the split approach and two different algorithms tested for direct segment piston reconstruction. Due to the high-end performance of the P-CuReD, both in terms of quality and speed, we choose this method as one of the components in our split approach aimed at segment-piston-free reconstruction of the wavefront for segmented pupils.

## 5.

## Split Approach

In this section, we describe robust algorithms that allow to compute optimal mirror configurations from signals containing telescope spider obstruction by dividing the wavefront reconstruction into two parts:^{56}

1. Piston-free wavefront reconstruction on each segment: This wavefront reconstruction method handles all modes seen by the PWFS except modes of order zero (piston modes) on each segment.

2. Direct segment piston reconstruction: Here, we focus on modes of order zero on each segment solely.

The reason why we suggest to split piston reconstruction from the full phase reconstruction is twofold. On the one hand, we want to make stable wavefront reconstruction in the presence of spiders feasible with the computationally efficient model-based algorithms presented in Sec. 4. Those algorithms were developed using the forward model of the sensor for annular apertures that do not include segmented pupils and spider effects. Straightforward attempts as described in Sec. 4.2 of applying these algorithms to sensor data “spoiled” by spider obstruction failed so far. On the other hand, it was recognized that the interaction-matrix-based approaches, as described in Sec. 3, are able to correct for differential pistons but are very time-consuming. The related computational effort may be affordable for the METIS system but hardly feasible for XAO systems. Therefore, our goal was to combine the P-CuReD (or any other model-based wavefront reconstruction method) and the advantages of interaction-matrix-based reconstruction to obtain a fast and robust reconstruction approach for segmented pupils having less computational complexity than an MVM.

We recall that $\mathit{P}$ describes the pyramid sensor operator, $\mathrm{\varphi}$ is the incoming phase, and $\overrightarrow{s}$ is the pyramid sensor measurements. Usually, all wavefront modes (frequencies) that are seen by the PWFS and afterward corrected are treated within the same wavefront reconstruction process. In the split approach, we separate the incoming wavefront $\mathrm{\varphi}$ into the parts:

where ${\mathrm{\varphi}}_{i}$ indicates the piston-free wavefront reconstruction on segment ${\mathrm{\Omega}}_{i}$. The corresponding reconstruction procedure denoted by ${\tilde{\mathit{P}}}^{\u2020}$ is provided by any of the existing model-based algorithms (mentioned in Sec. 4), which provide high-quality reconstruction on each segment. The term ${p}_{i}$ describes the corresponding piston information on every segment for $i=\mathrm{1,2},\dots ,k$ calculated independently using direct segment piston reconstruction methods that will be described in Sec. 5.2 and are in the following denoted by $\mathbf{\Pi}$.Hence, we separate the whole wavefront reconstruction into two parts:

## Eq. (14)

$$\mathrm{\varphi}={\mathrm{\varphi}}_{\mathrm{pistonfree}}+{\mathrm{\varphi}}_{\mathrm{piston}}={\tilde{\mathit{P}}}^{\u2020}\overrightarrow{s}+\mathbf{\Pi}\overrightarrow{s},$$^{19}

^{,}

^{57}

## 5.1.

### Piston-Free Wavefront Reconstruction on Segments

In this section, we focus on piston-free wavefront reconstruction on segments. For the estimation of the wavefront ${\mathrm{\varphi}}_{\mathrm{pistonfree}}$, we can use any of the existing fast model-based wavefront reconstruction algorithms mentioned in Sec. 4. These methods have shown exceptional wavefront correction quality on nonsegmented pupils, i.e., on the full annular telescope aperture in the following denoted by $\mathrm{\Omega}$. The spider legs divide the aperture into segments meaning that for $k$ spider legs we obtain $k$ disjoint segments indicated by ${\mathrm{\Omega}}_{i},i=\mathrm{1,2},\dots ,k$ (for instance the six spider legs of the ELT shown in Fig. 1). For segmented pupils, the wavefront reconstruction method ${\tilde{\mathit{P}}}^{\u2020}$—in addition to standard wavefront reconstruction requirements—fulfills two conditions:

1. The method is implemented on segments ${\mathrm{\Omega}}_{i}$, $i=\mathrm{1,2},\dots ,k$ instead of the full mask $\mathrm{\Omega}$.

2. The reconstruction ${\mathrm{\varphi}}_{i}$ on every segment ${\mathrm{\Omega}}_{i}$, $i=\mathrm{1,2},\dots ,k$ needs to be piston-free.

Note that condition 2 does not constitute a restriction since one can always compute the local piston information of a full segment from the reconstructed wavefront and subtract it afterward. Condition 1 may possibly be attenuated to an algorithm being implemented on the full aperture but dividing the reconstruction into segments thereafter. In this case also, the elimination of zero-order modes on the segments can be performed separately from the reconstruction process. However, we clearly want to point out that we did not investigate this idea in detail. Until now, we only used a segment-piston-free reconstruction method implemented on the segments. An analysis of these considerations will be part of a subsequent study.

For the segment-piston-free wavefront reconstruction, we use the P-CuReD^{9}^{–}^{11} applied to each segment. With linear computational complexity, this algorithm is the fastest method available for wavefront reconstruction from pyramid sensor data and, at the same time, it provides a reconstruction quality close to the theoretical limits.

## 5.2.

### Direct Segment Piston Reconstructors

Now we introduce two methods for direct segment piston reconstruction denoted by $\mathbf{\Pi}$ providing the piston information on the disjoint segments ${({p}_{i})}_{i=\mathrm{1,2},\dots k}$ divided by $k$ spider legs. All of them follow the idea of an interaction-matrix-based reconstruction using MVM. In contrast to the conventional approach, we now do not focus on the reconstruction of the complete wavefront, but only on the reconstruction of piston modes on segments. The first method (see Sec. 5.2.1 for details) uses a zonal interaction matrix containing the sensor response to every single zonal basis function. The second one (see Sec. 5.2.2 for details) is rather modal with a very small number of modes used as basis. Namely, we use a segment basis consisting of $k$ modes, where each mode depicts a piston on a given segment.

## 5.2.1.

#### Direct segment piston reconstructor I: single-poke-approach

We start with wavefront reconstruction using a full zonal interaction matrix as already described in Sec. 3 and transform this algorithm to the reconstruction of segment pistons solely.

For the so-called single-poke-approach, we measure the response of the pyramid sensor for every basis function. More precisely, we use the full interaction matrix $\mathit{M}$ of the system as described in Eqs. (5) and (6) computed for a set of zonal basis functions. The amplitude corresponding to the mentioned basis functions has to be small to ensure a linear response of the pyramid sensor. The control matrix ${\mathit{M}}_{\alpha}^{\u2020}\in {\mathbb{R}}^{{n}_{c}\times 2n}$ is created as in Sec. 3.4.2 using the squared Laplacian as a regularization term.

Since we are only interested in the reconstruction of segment pistons and omit the reconstruction of other modes, the dimension of the problem is drastically reduced. Assume we have the complete wavefront $\overrightarrow{\mathrm{\varphi}}={\mathit{M}}_{\alpha}^{\u2020}\overrightarrow{s}$ reconstructed. Extraction of the segment piston information $\overrightarrow{p}={({p}_{i})}_{1\le i\le k}$ from the known wavefront $\overrightarrow{\mathrm{\varphi}}$ is obtained by averaging the phase values within each segment. This step is modeled as a multiplication of $\overrightarrow{\mathrm{\varphi}}\in {\mathbb{R}}^{{n}_{c}}$ with a matrix $\mathit{Q}\in {\mathbb{R}}^{k\times {n}_{c}}$, where the $i$’th row of $\mathit{Q}$ contains a vector representation of the segment ${\mathrm{\Omega}}_{i}$ divided by the number of active subapertures on ${\mathrm{\Omega}}_{i}$ (for the averaging). The application of the matrix $\mathit{Q}$ leads to

## Eq. (15)

$$\overrightarrow{p}=\mathit{Q}\overrightarrow{\mathrm{\varphi}}=\mathit{Q}{\mathit{M}}_{\alpha}^{\u2020}\overrightarrow{s}\u2255{\mathbf{\Pi}}_{1}\overrightarrow{s},$$In Fig. 4, we provide an illustration to an Octopus simulation using the split approach with the P-CuReD and the first direct segment piston reconstructor (DSPR). Here we focus on the reconstructions in the initial loop steps, when the sensor is yet far from its linear regime and the loop is only getting closed. One can see that the reconstruction is reasonable and provides a stable convergence. The corresponding closed loop results for a longer simulation are presented in Sec. 6.

## 5.2.2.

#### Direct segment piston reconstructor II: segment-poke-approach

Inspired by the simple relation Eq. (15) connecting segment piston values $\overrightarrow{p}$ with the sensor data $\overrightarrow{s}$ through a small-size matrix ${\mathbf{\Pi}}_{1}$, we want to formulate another direct segment piston reconstruction approach, which will allow us to skip the computationally expensive and time-consuming step of setting up the full interaction matrix of the system. Rewriting Eq. (15) as

we see that, formally, it is possible to define a small interaction matrix of the system in the basis consisting of only a few segment pistons ${p}_{i}$, $i=1,\dots ,k$.We again consider arbitrarily located telescope spiders having $k$ spider legs, i.e., dividing the aperture into $k$ disjoint segments ${\mathrm{\Omega}}_{i}$, $i=\mathrm{1,2},\dots ,k$. The effects of a piston offset with amplitude $c$ on a single segment are described by

With the measurement vectors ${\overrightarrow{s}}_{i}\in {\mathbb{R}}^{2n}$ for every segment ${\mathrm{\Omega}}_{i}$, $i=\mathrm{1,2},\dots ,k$, we obtain a piston interaction matrix ${\mathit{M}}^{p}\in {\mathbb{R}}^{2n\times k}$ consisting of

The pyramid response to a segment poke is shown in Fig. 5. While in the previous method a very time-consuming precomputation of the full interaction matrix is required, the matrix necessary for this approach only has dimension $2n\times k$, and hence its computation is much cheaper. Therefore, setting up the matrix containing the WFS response to segment piston modes can be recomputed quickly for changing seeing conditions and readily performed online.

Now the piston reconstruction on every segment is described as minimization of

## 5.3.

### Fast and Robust Wavefront Reconstruction Under Pupil Segmentation Using the Split Approach

The general scheme of the split approach for wavefront reconstruction using one of the above-introduced direct segment piston reconstruction methods is described by Algorithm 1.

## Algorithm 1

Split approach.

choose segment-piston-free wavefront reconstruction method ${\tilde{\mathit{P}}}^{\u2020}$ |

choose direct segment piston reconstruction method $\mathbf{\Pi}$ |

for all time steps do |

get measurements $\overrightarrow{s}$ |

for segments $i=\mathrm{1,2},\dots ,k$ |

substract global piston for every ${\overrightarrow{\mathrm{\varphi}}}_{i}$ if necessary |

end for |

end for |

The segment piston control matrices ${\mathbf{\Pi}}_{i}$, $i=\mathrm{1,2}$ are dense but only of dimension $k\times 2n$, which leads to an optimization of the proposed approaches with respect to computational complexity. The expensive steps can be precomputed offline, thus, the algorithms scale linearly. This is a clear advantage to the computationally expensive wavefront reconstruction using full interaction-matrix-based MVM approaches as indicated in Fig. 6. For a moderately large-scale SCAO system such as METIS, the gain in the computational efficiency provided by the split approach is of order ${10}^{4}$. For an extremely large-scale SCAO system as the planned XAO system, the corresponding gain is of order ${10}^{5}$.

## 5.4.

### Details on the Realization

At the end, we will focus on a few numerical implementation details of the split approach, namely the incorporated pyramid sensor model, illumination factor, loop gain, and phase ambiguity.

For the proposed direct segment piston reconstruction methods, the nonlinear PWFS model including interference effects can be taken into account for the calculation of the piston control matrices ${\mathrm{\Pi}}_{i}$, $i=\mathrm{1,2}$.

Our numerical simulations and those performed in Ref. 41 show that for wavefront estimation under pupil fragmentation it is important to choose an appropriate illumination factor determining that subapertures are active and therefore used for wavefront reconstruction. As already discussed, we perform the segment-piston-free reconstruction on disjoint segments, hence we only use subapertures that are illuminated at least 75%. In contrast, for the DSPR, usage of measurements on less illuminated subapertures is crucial.

The loop gains for the two parts of the split approach have to be considered separately. We identify an optimal loop gain for the segment-piston-free wavefront reconstruction using the P-CuReD and another one for the direct segment piston reconstruction. At the beginning of a closed loop simulation, it is crucial to avoid phase ambiguity caused by the sinusoidal part of the pyramid data, meaning that piston offsets of size $2\pi $ radians in the phase cannot be distinguished by the pyramid sensor but heavily influence the image quality.^{58} Due to the nonlinearity of the PWFS (“optical gain”) and the need to correct low frequencies fast, we use a higher integrator control loop gain for the DSPR in the first iterations. This results in a stronger emphasis on the correction of low-order modes and provides an adequate control of piston offsets for data corresponding to larger phases.

The last step of the algorithm contains the projection of the reconstructed wavefront on the DM. One can, therefore, either solve Eq. (7) or simply evaluate $\mathrm{\varphi}$ at the actuator positions. We did choose the second approach. To be able to control also actuators outside the reconstruction area, we smoothly extend the reconstruction $\mathrm{\varphi}$ to a larger domain covering all used DM actuator positions.

## 6.

## Numerical Validation

## 6.1.

### Simulation Environment and Parameters

We use the simulation tool Octopus^{35}^{,}^{36} provided by the European Southern Observatory and study the METIS instrument having an SCAO system incorporated. As one of the three first light instruments of the ELT, METIS will allow investigations of exoplanets concerning physical and chemical properties like weather, temperature, seasons, or the composition of their atmospheres. It will, among others, focus on proto-planetary disks, the formation of planets, and the Solar System as well as the growth of supermassive black holes. The science wavelengths of the instrument range from 3 to $19\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the wavefront sensing is performed at $2.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. For the METIS instrument, an annular mask is currently considered with an outer diameter of 37 m and an inner diameter of 11.1 m (30%), i.e., all edges of the real ELT primary mirror having a diameter of 39 m are cropped such that we obtain a circular area. As a sensing device, we simulate a PWFS having $74\times 74$ subapertures corresponding to a subaperture size of 0.5 m. The simulation grid size is selected as $740\times 740$ pixels on the aperture resulting in a resolution of 0.05 m per pixel for a 37-m telescope.

The phase screens of our simulation are based on the von Karman atmospheric model having 35 layers and an outer scale of ${L}_{0}=25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The layer model of the atmosphere ranging from 30 m to 26.5 km above the observatory’s platform was provided by the European Southern Observatory. Two seeing conditions were simulated: median with the Fried parameter ${r}_{0}$ at $0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}=0.157\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and very bad with ${r}_{0}$ at $0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}=0.097\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. We consider two different photon fluxes of 600 and 100 incident photons per subaperture per frame corresponding to a guide star magnitude of 8 or 10, respectively.

A key element of every AO system is the DM. It provides fast steering capabilities to compensate for wavefront aberrations caused by atmospheric turbulence and telescope perturbations in real time and hence allows one to optimize the telescope performance. The actuator positions of the M4 DM used in our Octopus simulations are shown in Fig. 7. We use the M4 influence functions internally incorporated in Octopus. Altogether, we have 3874 active subapertures and 5190 active actuators in use. All test case parameters are summarized in Table 3.

## Table 3

Simulation parameters used in the tests.

Simulation parameters | |
---|---|

Telescope diameter | 37 m |

Central obstruction | 11.1 m, i.e., 30% |

Pupil mask | ELT pupil with 50-cm wide spiders |

Pupil segments | 6 |

Science target | On-axis (SCAO) |

WFS | PWFS |

Sensing band | K ($2.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) |

Evaluation bands | K ($2.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) |

L ($3.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) | |

N ($10.0\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) | |

Modulation | $4\lambda /D$ |

Controller | Integrator |

Atmospheric model | von Karman |

Number of simulated layers | 35 |

Outer scale ${L}_{0}$ | 25 m |

Atmosphere | Median |

Fried radius ${r}_{0}$ at $\lambda =500\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ | 0.157 m, 0.097 m |

Coherence time ${\tau}_{0}$ at $\lambda =500\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ | 5.35 ms |

Number of subapertures | $74\times 74$ |

Minimum subaperture illumination | 45% |

Number of active subapertures | 3874 out of 5476 |

Linear size of simulation grid | 740 pixels |

DM geometry | ELT M4 model |

DM delay | 1 frame |

Number of active actuators | 5190 |

Detector read-out noise | 1 electron/pixel |

Background flux | 0.431 photons/pixel/frame or $4.19\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ADU}/{\mathrm{m}}^{2}/\mathrm{ms}/{\mathrm{arcsec}}^{2}$ |

Frame rate | 500 Hz |

Photon flux or guide star brightness | [100, 600] photons/subaperture/frame or [8, 10] mag |

Zenith angle | [0 deg,30 deg,60 deg] |

Simulation time | 2 s (1000 iterations) |

## 6.2.

### Numerical Results

A numerical performance analysis and a comparison of the approaches in terms of reconstruction quality for the above-described simulation environment are provided in Tables 4–Table 56 and Fig. 8. The three tables present the long-exposure Strehl ratios evaluated at three different wavelengths $\lambda =2.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, $\lambda =3.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, and $\lambda =10.0\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ correspondingly. In each of these three tables, the test cases cover two photon fluxes of 600 and 100 photons/subaperture/frame and three zenith angles of 0 deg, 30 deg, and 60 deg. Each table demonstrates the results obtained with the P-CuReD for a spider-free simulation as benchmark and the spider effected results obtained within the split approach using each of the DSPR.

## Table 4

Long-exposure Strehl ratios in K-band (2.2 μm) after 1000 closed loop simulation steps. As reconstruction method, we use the P-CuReD. We compare the results for a benchmark simulation without spiders using the P-CuReD only and employing the split approach combined with the proposed DSPRs in the presence of telescope spiders.

Photon flux | Zenith angle (deg) | No spiders | DSPR I | DSPR II | Theoretical limit |
---|---|---|---|---|---|

600 | 0 | 0.8851 | 0.8775 | 0.8654 | 0.8882 |

30 | 0.8690 | 0.8597 | 0.8462 | 0.8820 | |

60 | 0.7799 | 0.7660 | 0.7412 | 0.7956 | |

100 | 0 | 0.8741 | 0.8650 | 0.8523 | |

30 | 0.8578 | 0.8473 | 0.8327 | ||

60 | 0.7712 | 0.7551 | 0.7307 |

## Table 5

Long-exposure Strehl ratios in N-band (10 μm) after 1000 closed loop simulation steps using the P-CuReD in all simulations.

Flux | Zenith angle (deg) | No spiders | DSPR I | DSPR II | ESO goal | ESO requirement | Theoretical limit |
---|---|---|---|---|---|---|---|

600 | 0 | 0.9941 | 0.9937 | 0.9930 | 0.9946 | ||

30 | 0.9932 | 0.9927 | 0.9919 | 0.9943 | |||

60 | 0.9880 | 0.9871 | 0.9855 | 0.60 | 0.9901 | ||

100 | 0 | 0.9935 | 0.9930 | 0.9923 | |||

30 | 0.9926 | 0.9920 | 0.9911 | 0.95 | 0.93 | ||

60 | 0.9875 | 0.9864 | 0.9848 |

## Table 6

Long-exposure Strehl ratios in L-band (3.7 μm) after 1000 closed loop simulation steps using the P-CuReD in all simulations.

Flux | Zenith angle (deg) | No spiders | DSPR I | DSPR II | ESO goal | ESO requirement | Theoretical limit |
---|---|---|---|---|---|---|---|

600 | 0 | 0.9577 | 0.9547 | 0.9500 | 0.9605 | ||

30 | 0.9515 | 0.9478 | 0.9425 | 0.9583 | |||

60 | 0.9157 | 0.9095 | 0.8990 | 0.57 | 0.9277 | ||

100 | 0 | 0.9535 | 0.9499 | 0.9449 | |||

30 | 0.9472 | 0.9429 | 0.9371 | 0.80 | 0.60 | ||

60 | 0.9121 | 0.9052 | 0.8944 |

From the tables, we see that, compared to the case without spiders, some small loss in quality is present when pupil segmentation due to spiders is taken into account and solved with the split approach. Moreover, one can see that the loss increases slightly for larger zenith angles and that the DSPR I provides slightly higher Strehl ratios compared to the DSPR II. Figure 8 shows the corresponding short-exposure Strehl ratios over the simulation time demonstrating that the DSPR I is indeed more stable compared to the DSPR II. More precisely, the performance of the two methods is not identical because phase averaging over the segments and matrix inversion are interchanged in the DSPR II compared to the DSPR I. As confirmed in Fig. 9 showing the residual pistons on each segment for the DSPR I, the reconstructions obtained with the split approach are free from large differential piston effects.

Based on the provided reconstruction quality and stability features, the conclusion is that the DSPR I method is preferable. That is why we focused on the DSPR I in further numerical studies with Figs. 10 and 11 presenting acceptable METIS results obtained at a zenith angle 60 deg under very bad seeing conditions (nominal ${r}_{0}$ at $0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}=0.064\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$), which corresponds to ${r}_{0}$ at $2.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}=0.379\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$.

However, one should keep in mind that in case the atmospheric or observational conditions change, registration and inversion of an interaction matrix will be required. This process is very time-consuming and may lead to a loss of precious operational time on the telescope. On the contrary, in the DSPR II method, the underlying interaction matrix is much smaller, therefore, significantly less operations are needed for its registration and inversion. Since the DSPR II method provides acceptable reconstruction quality as well, it can be used at least as a substitute of the DSPR I method when the latter requires some updates to be performed.

In addition, Tables 5 and 6 contain the ESO goals and ESO requirements. Note that both are defined as performance evaluations over at least 15 min of telescope operations under nominal conditions. These include some additional error sources, e.g., wind induced vibrations or noncommon path aberrations (NCPAs), which were not taken into account in the presented simulations. Therefore, we cannot straightforwardly compare our results that evaluate only 2 s of operation, with the ESO goals and requirements, but we can still see that there is a big safety gap allowed for vibrations and other error terms reducing the quality. We would like to point out that we are aware of those additional error sources and consider them in separate dedicated tests. For instance, the influence of wind shake and NCPAs on the METIS performance has been previously analyzed in Ref. 40. In this paper, we focus on eliminating the impact of spiders in the first turn. Analyzing and reducing different error terms one-by-one allows one to achieve confidence in the final performance of the instrument.

Moreover, in each of the tables, we show the roughly estimated theoretical limit of the achievable long-exposure Strehl ratio. These rough estimates are obtained using the assumption that in the high-flux case the reachable quality is limited from above by the two main error sources, the fitting, and the temporal delay error. Following the approach in Ref. 34, we evaluate the fitting error by

where $d$ denotes the average actuator’s distance, ${r}_{0}$ is the Fried radius, and the delay error by where $\tau $ is the delay and ${\tau}_{0}$ is the coherence time of the atmospheric turbulence.Furthermore, we want to remark that the results presented in the tables were obtained with the same loop gains for different guide star magnitudes and zenith angles, which clearly underline the stability of the algorithms. The DSPR II was experienced to be more sensitive to the choice of the loop gain than the DSPR I. The quality results of several numerical simulations can even be slightly improved by applying a loop gain being optimized with respect to the special parameter choices of the individual test cases.

In order to quantify the achieved raw contrast, we plot in Fig. 12 the simulated radially averaged PSF obtained with the DSPR I method. The area of interest, where the correction actually happens, is defined by the number of actuators and the telescope diameter as ${n}_{a}(\lambda /D)$. For $\lambda =3.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the best corrected central area of the PSF has a square shape with a side length of 1546 mas. The maximum achieved contrast is of order ${10}^{-6}$ within the control area of the DM. Figure 13 shows the difference in the simulated radially averaged point spread function (PSF) at a wavelength of $\lambda =3.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ obtained with the DSPR II method compared to the DSPR I in the case of a bright guide star and zenith angle 0 deg.

## 7.

## Conclusions

The first part of this paper provided a brief overview on the ability of existing algorithms to reconstruct wavefronts from PWFS data in the presence of support structures for the secondary mirror, the so-called spiders. It is known that a partial or complete shading of subapertures, depending on the thickness of the spider legs, leads to difficulties in controlling piston modes on disjoint pupil segments. We have shown that on the one hand, at least some variants of the interaction-matrix-based MVM approaches are able to handle spider obstruction effects successfully, though they are known to be computationally demanding. On the other hand, there exist several fast interaction-matrix-free model-based algorithms, which provide high-quality reconstruction in case of annular apertures without spiders. Unfortunately, they run into problems when dealing with sensor data (partially) shaded by spider legs. In end-to-end simulations, we observe in this case random uncontrolled segment pistons in the residuals. For the high-contrast large-scale AO systems in design, a combination of the advantages of both reconstructor types, i.e., the ability to reconstruct wavefronts with high precision and speed in the presence of spiders, is highly desirable.

In the second part of this paper, we presented a solution, the so-called split approach, which combines the advantages of the interaction-matrix-based and free methods and provides a high-quality wavefront reconstruction in the presence of spiders with little computational demand. The solution is based on the idea to split the reconstruction of pistons on the pupil segments from the reconstruction of higher-order modes (frequencies). For the latter part, one can use any of the available fast model-based algorithms. Especially, beneficial here is the P-CuReD, which has the smallest computational load and also provides the best reconstruction quality in case no spider legs shade the pupil.

For segment piston reconstruction, we presented two approaches. One of them requires the registration and inversion of a full interaction matrix. However, the control matrix is afterward projected to the space of segment pistons only, hence reduced in size. Another one uses a small size interaction matrix instead, which is registered in the basis of segment pistons. The resulting control matrices are in both methods of small size and the direct segment piston reconstruction step has a linear complexity.

Combined with the P-CuReD used for reconstruction of higher-order modes, the number of computations required for the complete wavefront reconstruction scales linearly with the number of controlled actuators. This represents a big advantage of the split approach compared to the usual interaction-matrix-based MVM, whose complexity scales quadratically with the number of actuators. Although MVM approaches may still be computationally doable for relatively large-scale systems like METIS having a $74\times 74$ pyramid sensor in 2026, they are hardly feasible for extremely large systems like XAO having $200\times 200$ subapertures and a corresponding number of actuators to be controlled in real time. The presented split approach causes, in contrast, no difficulties for the real-time implementation even on the largest AO systems of ELT-sized telescopes.

Moreover, the split approach makes existing model-based phase reconstruction algorithms developed for nonsegmented pupils suitable for wavefront control in the presence of telescope spiders. Alternatively, another idea to overcome the effects of pupil fragmentation with the model-based algorithms consists in an appropriate adjustment of the underlying forward models. A corresponding extension of the algorithms will be necessary and may allow for a stable control of segmented piston offsets without using interaction-matrix-based attempts. The question whether such an extension is possible and for which methods it is applicable needs further investigations.

The analysis of the split approach was done in the context of the instrument METIS for sensing in the K-band. Further investigations about DSPRs using pyramid sensors in different sensing wavelengths are of high interest for several instruments in development for ELTs and therefore planned future work. Although, the theory of the split approach is nonspecific to the instrument and therefore applicable to any AO system, the wavefront reconstruction may possibly suffer from phase ambiguity or other wavelength dependent effects of the pyramid sensor.

## Acknowledgments

The authors are grateful to Miska Le Louarn from the European Southern Observatory (ESO) and the Austrian Adaptive Optics team for support. This work was partially funded by the Austrian Federal Ministry of Science and Research (HRSM) and the Austrian Science Fund (FWF) F68-N36, project 5.

## References

## Biography

**Victoria Hutterer** received her BS and MS degrees in mathematics from Johannes Kepler University Linz, Austria, in 2014 and 2015, respectively. She is a PhD student at the Johannes Kepler University Linz. Her current research interests include inverse problems, adaptive optics, pyramid wavefront sensors, wavefront reconstruction algorithms, and iterative and nonlinear algorithms.

**Iuliia Shatokhina** received her BS and MS degrees in physics from the National University of Kyiv-Mohyla Academy in 2006 and 2007, respectively, and her PhD in mathematics from Johannes Kepler University Linz, Austria, in 2014. She is a research assistant at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences. Her current research interests include inverse problems, adaptive optics, pyramid wavefront sensors, and wavefront reconstruction algorithms.

**Andreas Obereder** received his MS degree in mathematics in 1998 and his PhD in mathematics from Johannes Kepler University Linz, Austria, in 2003. He is a research assistant at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences. His current research interests include inverse problems, astronomical adaptive optics, and wavefront reconstruction algorithms.

**Ronny Ramlau** received his MS degree in mathematics in 1994 and his PhD in mathematics from the University of Potsdam, Germany, in 1997. He is a university professor at Johannes Kepler University Linz, Austria, where he leads the Industrial Mathematics Institute. Also, he is a scientific director at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences, leading the Transfer Group. He is the author of more than 50 journal papers and has written 3 book chapters. His current research interests include nonlinear inverse problems, medical imaging, regularization methods, and iterative methods.