# How to describe Fourier filtering wavefront sensor nonlinearity with interaction matrices

Mahawa Cisse, Olivier Fauvarque, Vincent Chambouleyron, Nicolas Levraud, Charlotte Z. Bond, Benoit Neichel, Jean-François Sauvage, Thierry Fusco

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

## To cite this version:

Mahawa Cisse, Olivier Fauvarque, Vincent Chambouleyron, Nicolas Levraud, Charlotte Z. Bond, et al.. How to describe Fourier filtering wavefront sensor nonlinearity with interaction matrices. Journal of Astronomical Telescopes Instruments and Systems, 2025, 11 (02), pp.029001. 10.1117/1.jatis.11.2.029001. hal-05040075

HAL Id: hal-05040075

https://hal.science/hal-05040075v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# How to describe Fourier filtering wavefront sensor nonlinearity with interaction matrices

Mahawa Cisse<sup>a,b,*</sup> Olivier Fauvarque,<sup>c,*</sup> Vincent Chambouleyron,<sup>d</sup>

Nicolas Levraud,<sup>e</sup> Charlotte Z. Bond,<sup>f</sup> Benoit Neichel,<sup>b</sup>

Jean-François Sauvage,<sup>a,b</sup> and Thierry Fusco<sup>g</sup>

**ABSTRACT.** Fourier filtering wavefront sensors (WFSs) are a class of highly sensitive sensors that can significantly enhance adaptive optics (AO) performance, particularly in low-flux regimes. However, their nonlinear behavior limits their effectiveness for measuring high-amplitude phases. This paper presents a method for characterizing the nonlinearity of Fourier filtering WFSs within the widely used matrix formalism. The nonlinearity arises from an over-modulation effect depending on the phase being measured. Consequently, the matrix describing the sensor outside its linearity range must also depend on this phase. We first propose a theoretical framework, derived from light propagation equations, to construct a reconstructor capable of accounting for the WFS's nonlinear responses. This analytical approach yields an exact expression for the reconstructor within the matrix formalism, referred to as the specific matrix, as it depends on the phase to be measured, which makes it impractical to use. Therefore, a portion of the paper is dedicated to deriving an approximation of the specific matrix. A method for approximating the specific matrix using data fusion with a focal plane camera is introduced. Simulation results demonstrate the efficacy of this approach when applied to a nonlinear WFS, such as the unmodulated pyramid WFS, under challenging seeing conditions (up to 1.3").

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JATIS.11.2.029001]

**Keywords:** adaptive optics; wavefront sensor; extremely large telescope

Paper 24123G received Aug. 1, 2024; revised Feb. 23, 2025; accepted Mar. 10, 2025; published Apr. 8, 2025.

## 1 Introduction

The angular resolution of ground-based telescopes depends, among other factors, on the performance of their adaptive optics (AO) systems. Specifically, the capabilities of the deformable mirror (DM) and wavefront sensor (WFS) are essential in measuring and correcting phase aberrations caused by atmospheric turbulence. We focus on a particular class of WFSs known as Fourier filtering WFSs. These include, for example, the pyramid WFS (PWFS)¹ and its variants,²⁻⁵ whether modulated or not, as well as Zernike WFS.⁶⁻⁸ The interest in these sensors arises

<footnote><sup>a</sup>DOTA, ONERA, Salon Cedex Air, France</footnote>

<footnote><sup>b</sup>Aix Marseille University, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote><sup>c</sup>Ifremer, RDT Research and Technological Development, Plouzané, France</footnote>

<footnote><sup>d</sup>University of California Santa Cruz, Santa Cruz, California, United States</footnote>

<footnote><sup>e</sup>LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, Meudon, France</footnote>

<footnote><sup>f</sup>UK Astronomy Technology Centre, Edinburgh, United Kingdom</footnote>

<footnote><sup>g</sup>DOTA, ONERA, Université Paris Saclay, Palaiseau, France</footnote>

<footnote>*Address all correspondence to Mahawa Cisse, mahawacisse13@gmail.com; Olivier Fauvarque, olivier.fauvarque@ifremer.fr</footnote>
from their excellent noise propagation properties and their designation as the WFSs of choice for the AO systems in future Extremely Large Telescopes (ELTs).

In AO, the measurement process can be described as follows: the incident phase enters the optical system of the WFS. The electromagnetic field then propagates through the system before being detected by a 2D sensor. During detection, noise (such as readout noise or photon noise) may be added to the signal. The resulting “signal plus noise,” now a digital quantity, is passed through a numerical operator known as the reconstructor, **R**, which inverts it to estimate the phase. To meet the real-time constraints of an AO system, the reconstructor **R** is typically assumed to be linear. Mathematically, this process can be expressed as

$$s = h(\phi) + \text{Noise}, \qquad (1)$$

$$\hat{\phi} = \mathbf{R}[h(\phi) + \text{Noise}], \qquad (2)$$

where *s* is the signal, *h* is a function of the phase *φ*, and $\hat{\phi}$ is the estimated phase.

Assuming a linear operator, the estimation error can be written as

$$\epsilon = \hat{\phi} - \phi, \qquad (3)$$

$$=\mathbf{R}[\text{Noise}] + (\mathbf{R}[h(\phi)] - \mathbf{Identity}[\phi]). \qquad (4)$$

If *h* is invertible, the reconstructor **R** becomes $h^{-1}$, yielding

$$\mathbf{R} = h^{-1} \rightarrow \mathbf{R}h = \mathbf{Identity}. \qquad (5)$$

In this case, the only source of error arises from the propagation of noise through the reconstructor. This noise propagation has been widely studied⁹,¹⁰ and is closely linked to the sensitivity of the WFS, which should be maximized to improve the signal-to-noise ratio.⁷,¹¹,¹²

However, another significant source of error arises when the input-output model of the WFS, given by Eq. (1), is not correct. If the model does not accurately represent the system, the relationship can be expressed as

$$s = h'(\phi) + \text{Noise} \neq h(\phi) + \text{Noise},$$

where $h'$ is an unknown function describing the actual system behavior, whereas *h* remains the only available model. This mismatch introduces modeling errors that add to the estimation error. Addressing these errors involves a dual challenge: constructing an accurate input-output model of the WFS for on-sky operation and developing mathematical methods to invert the model effectively.

Recent research into wavefront reconstruction has explored the use of artificial intelligence, particularly those leveraging deep learning, to provide a promising middle ground between purely analytical methods and traditional AO system calibration techniques. Neural networks trained on extensive datasets of wavefront propagation simulations have shown the capability to address the inherent nonlinear behavior of Fourier-based WFSs. Notable studies¹³⁻¹⁶ have successfully employed convolutional and transformer networks to mitigate nonlinearity and enhance reconstruction accuracy, even in open- and closed-loop systems. Specifically, these studies demonstrated the feasibility of using a neural network to close the AO loop using an unmodulated PWFS on an experimental bench. These advancements highlight the potential for deep learning to enhance real-time AO system performance. Although the deep-learning approach has shown promising results, it does not provide insights into the behavior of nonlinear WFSs. Moreover, these approaches have yet to be practically implemented in on-sky applications. The limited understanding of WFS behavior, combined with the challenges of on-sky implementation for deep-learning methods, underscores the robustness and reliability of analytical and calibration-based approaches for operational use.

An alternative approach uses the physical laws describing light propagation within the WFS to model its behavior. Although this can yield sufficiently accurate and invertible models in some cases, such as the PWFS¹ and its iterative inversion algorithms¹⁷⁻¹⁹ or Zernike WFS,⁶ it does not entirely eliminate modeling errors because they do not describe the effective path of light in the WFS.
To address this problem, most AO systems rely on a calibration approach, constructing an “interaction matrix” that accounts for the effective light propagation from the DM to the WFS. This method assumes linearity therefore, under this hypothesis, the function $h$ from Eq. (1), becomes the interaction matrix **H**. Two cases can be considered

1. The matrix **H** is partially invertible. In that case, the reconstructor will be the pseudo-inverse of the matrix **H**: $R = H^\dagger = (H'H)^{-1}H'$. Two cases arise: (i) no eigenmodes are filtered in the process, leading to the correction of all controlled modes. The only source of noise is the noise or (ii) eigenmodes are filtered in the inversion, which adds an estimation error.

2. The matrix **H** only describes the input/output relationship in a specific phase regime.

Unfortunately, these points will lead to modeling and reconstruction errors, as representing the behavior of a WFS with a matrix assumes that the sensor responds linearly to the incident phase. Although this assumption is valid for some sensors over a large range of phase amplitude,²⁰ it appears that most WFSs, especially the Fourier filtering WFSs,³ have signals inherently nonlinear with respect to the phase.

## 1.1 Outline

This paper aims to unify the matrix formalism with the inherent nonlinearity of Fourier filtering WFSs within a single coherent framework. Specifically, we demonstrate how the behavior of these WFSs can be described using a matrix, regardless of the amplitude of the incoming phase. Furthermore, we show how this formalism can enhance the reconstruction process by providing an accurate estimate of the optimal reconstructor, **R**.

Section 2 introduces the nonlinearity of Fourier filtering WFSs and its effects on phase estimation. This discussion establishes the fundamental characteristics of the target reconstructor to be computed. The reconstructor is defined as the matrix that best represents the WFS response within the given phase regime.

In Sec. 3, the analytical expression for the matrix that accurately captures the nonlinearity of Fourier filtering WFSs is presented. This matrix, called the “specific matrix”, is shown to depend on the phase being measured. Its dependence on the phase (which is the unknown) makes it impractical for implementation.

Therefore, in Sec. 4, we propose a method to approximate the specific matrix. The proposed approach is based on data fusion and enables the derivation of a reconstructor—the hybrid reconstructor—that closely approximates the specific matrix. This section provides a simulation-based demonstration of this method using an AO system akin to that of the Very Large Telescope (VLT).

# 2 Handling the Nonlinearity of the Fourier Filtering Wavefront Sensors

## 2.1 General Framework of Fourier Filtering Wavefront Sensing

Fourier filtering WFSs use an object placed in the Fourier (or focal) plane, referred to as the filtering mask, to convert phase fluctuations of the incoming light into intensity variations in the subsequent pupil plane (see Fig. 1).

For instance, in the case of a PWFS, the filtering process is supplemented by a modulation stage.¹ This consists of a mirror in the pupil plane that cyclically applies a tip/tilt aberration to the incoming light. This additional degree of freedom, represented by a weighting function $w$,²¹ was introduced to expand the sensor's linear range. The weighting function encodes the path of the Point Spread Function (PSF) on the focal plane due to tip/tilt modulation.

The following notations are adopted in this framework:

* $\mathbb{I}_p$: The telescope's pupil indicator function, which equals 1 within the pupil and 0 outside.

* $w$: The modulation weighting function, normalized to 1 for energy conservation.

* $m$: The transparency function of the filtering mask, constrained by $|m| < 1$.

* $\phi_r$: The reference phase, around which the system operates. For simplicity, we will assume that the reference phase is zero.
<figure><img src="image_3.png" /><figcaption>Fig. 1 Schematic view (in 1D) of a Fourier filtering optical system.</figcaption></figure>

* $\phi$: The incident phase, expressed as $\phi = \delta\phi + \phi_r$, where $\delta\phi$ denotes the difference between the incident and reference phase.

* $I(\phi)$: The intensity detected at the WFS output for a given phase $\phi$.

Using these notations, the intensity at the detector can be formulated as a function of the incident phase and WFS parameters using the diffraction laws as in²¹

$$I(\phi)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'|} \mathbb{I}_P|_{\vec{r}} e^{i(\phi|_{\vec{r}'} - \phi|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \hat{w}|_{\vec{r}'-\vec{r}}, \quad (6)$$

where $\vec{r}', \vec{r}$ and $\vec{R}$ are the 2D spatial variables. $\hat{}$ is the Fourier transform, and $\vec{}$ the complex conjugate.

To work around the operating point, a tare operation is performed on the intensity data. This involves subtracting the reference intensity, yielding the “reduced intensity,” defined as

$$\Delta I(\delta\phi; \phi_r) \equiv I(\phi) - I(\phi_r), \quad (7)$$

$$=I(\delta\phi + \phi_r) - I(\phi_r). \quad (8)$$

This reduced intensity may be expressed as a function of the WFS parameters as in²¹

$$\Delta I(\delta\phi; \phi_r)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'|} \mathbb{I}_P|_{\vec{r}} e^{i(\phi_r|_{\vec{r}'} - \phi_r|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \hat{w}|_{\vec{r}'-\vec{r}} (e^{i(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} - 1). \quad (9)$$

This equation serves as the input-output relationship for the WFS. Estimating the phase depends on the ability to invert this relation, i.e., to determine the phase being measured $\delta\phi$ as a function of the reduced intensity, $\Delta I(\delta\phi; \phi_r)$.

## 2.2 One-Dimension Toy Model

To simplify the analysis, we use a 1D model to illustrate the matrix formalism. In this model (see Fig. 2), the incident phase $\phi$, the phase being measured $\delta\phi$, and the reference phase $\phi_r$ are treated as real scalars, denoted as $x$, $\delta x$, and $x_r$, respectively. Similarly, the intensity $I$ is simplified to a scalar function $f$. The correspondences are as follows:

$$\begin{aligned} \phi &\leftrightarrow x, \\ \delta\phi &\leftrightarrow \delta x, \\ \phi_r &\leftrightarrow x_r, \\ I(\phi) &\leftrightarrow f(x). \end{aligned}$$

In this 1D model, the WFS output, i.e., the reduced intensity, becomes

$$\Delta I(\delta\phi; \phi_r) \leftrightarrow f(x_r + \delta x) - f(x_r). \quad (10)$$

To mimic the behavior of Fourier filtering WFSs, the function $f$ is assumed to be positive, nonlinear, and saturating, as shown by the black curve in Fig. 2.
<figure><img src="image_4.png" /><figcaption>Fig. 2 One-dimensional representation of the input-output WFS relation (black curve). Classic linear modeling of the WFS around its operating point (blue curve). Linear modeling of the WFS for the phase-to-be-measured δx (red curve).</figcaption></figure>

### 2.2.1 Linear modelling of Fourier filtering WFS

Returning to the equation of propagation, the input-output relationship can be approximated within a linear framework, as demonstrated in Ref. 22. This approach is valid only around the reference phase $\phi_r$, where small phase variations ($\delta\phi \ll 1$) allow for the linearization of Eq. (6). The linearized expression for the reduced intensity is

$$ \Delta I(\delta\phi; \phi_r)_{\text{lin}|\vec{R}} = \int_{\mathbb{R}^2} d^2\vec{r}' (\mathbb{I}_p|_{\vec{r}}\delta\phi|_{\vec{r}}) 2 \mathbf{Im} \left[ \tilde{m}|_{\vec{R}-\vec{r}} \int_{\mathbb{R}^2} d^2\vec{r}' \mathbb{I}_p|_{\vec{r}'}\hat{m}|_{\vec{R}-\vec{r}'} \hat{w}|_{\vec{r}'-\vec{r}} e^{i(\phi_r|_{\vec{r}'}-\phi_r|_{\vec{r}})} \right]. \quad (11) $$

To invert this relationship, the usual technique is to perform a very low amplitude push-pull for each phase, {$\Phi_i$}, around the reference phase and to collect the differential intensities into the matrix known as the interaction matrix, $\mathbf{M}$

$$ \begin{aligned} \mathbf{M}_i &\equiv \frac{I(\phi_r + \varepsilon\Phi_i) - I(\phi_r - \varepsilon\Phi_i)}{2\varepsilon} & \text{with } \varepsilon \ll 1 \\ &= \frac{\Delta I(\varepsilon\Phi_i; \phi_r) - \Delta I(-\varepsilon\Phi_i; \phi_r)}{2\varepsilon} \\ &= \frac{\Delta I_{\text{lin}}(\varepsilon\Phi_i; \phi_r) - \Delta I_{\text{lin}}(-\varepsilon\Phi_i; \phi_r)}{2\varepsilon}, \end{aligned} \quad (12) $$

This yields

$$ \Delta I(\delta\phi; \phi_r)_{\text{lin}} = \mathbf{M}\mathbf{a}(\delta\phi), \quad (13) $$

where $\mathbf{a}(\delta\phi)$ is a vector corresponding to the projection of the phase to be measured on the {$\Phi_i$} basis.

The reconstructor $\mathbf{R}$ is then the pseudo-inverse of $\mathbf{M}$. The following relation:

$$ \mathbf{a}(\delta\phi) = \mathbf{M}^{\dagger} \Delta I(\delta\phi; \phi_r)_{\text{lin}}, \quad (14) $$

constitutes an estimation of the phase from the reduced intensities when $\delta\phi$ is in the linear regime of the WFS.

For the 1D model, $\mathbf{M}$ corresponds to the derivative of $f$ at $x_r$, expressed as

$$ \mathbf{M}(0; \phi_r) \leftrightarrow f'(x_r). \quad (15) $$

In Fig. 2, the blue line represents the linearization of $f$ around $x_r$, whereas the scalar pseudo-inverse of $\mathbf{M}$ becomes the inverse of the function $f$ as follows:
$$ \mathbf{M}(0; \phi_r)^{\dagger} \leftrightarrow \frac{1}{f'(x_r)}. \qquad (16) $$

Applying this reconstructor to estimate $\delta x$, we find

$$ \tilde{\delta}x = \frac{1}{f'(x_r)} [f(x_r + \delta x) - f(x_r)]. \qquad (17) $$

However, the linear model leads to errors when $\delta x$ moves away from $x_r$. These errors typically result in an underestimation of $\delta x$, consistent with the underestimation observed in practice (see Sec. 4.1.1 for the optical gain (OG) discussion).

## 2.3 Modelling/Reconstruction Error Outside the Linear Regime

The reconstructor $\mathbf{M}^{\dagger}$ provides accurate phase estimation when the phase amplitude is small relative to the reference phase ($\delta\phi \ll 1$). However, as the phase amplitude increases, the linear approximation of the reduced intensity is no longer valid, and $\mathbf{M}^{\dagger}$ fails to estimate $\delta\phi$ accurately. This limitation arises from the intrinsic nonlinearity of Fourier filtering WFSs, specifically due to the term: $(\delta\phi|_{\hat{\mathbf{r}}'} - \delta\phi|_{\hat{\mathbf{r}}}) - 1$ from Eq. (6). This nonlinearity introduces errors in the modeling and reconstruction processes, manifesting in two distinct ways as $\delta\phi$ increases

* Sensitivity loss and underestimation: nonlinearity causes a reduction in the WFS's sensitivity. This phenomenon is characterized by modal coefficients forming a diagonal matrix known as the OG matrix.²³ Although various methods exist to correct this sensitivity loss,²⁴⁻²⁶ they tend to lose accuracy as the phase amplitude increases.

* Mode coupling: when the phase amplitude is large, nonlinearity induces cross-talk between modes. This means each mode {$\Phi_i$} is no longer estimated independently; instead, its estimation is distorted by the presence of other modes in $\hat{\delta}\phi$.

## 2.4 Specific Matrix

This work aims to show that it is possible to address the modeling and reconstruction errors caused by nonlinearity while remaining within the matrix formalism. The key challenge is to determine the matrix that accurately represents the behavior of the WFS when measuring $\delta\phi$, regardless of its amplitude. Because nonlinearity can be viewed as self-modulation (cf Sec. 2.6), which alters the sensitivity of the WFS, the matrix must depend on the phase being measured. For this reason, the matrix is called the “specific matrix” and is a function of both the reference phase $\phi_r$ and the measured phase $\delta\phi$.

The specific matrix, $\mathbf{M}(\delta\phi; \phi_r)$, satisfies the following relationship:

$$ \Delta I(\delta\phi; \phi_r) = \mathbf{M}(\delta\phi; \phi_r)\mathbf{a}(\delta\phi), \qquad (18) $$

where $\mathbf{M}(\delta\phi; \phi_r)$ is the matrix that best describes the WFS response for the given phase regime. In the linear case, $\mathbf{M}(\delta\phi; \phi_r)$ reduces to the interaction matrix $\mathbf{M}$, and the reconstructor $\mathbf{R}$ becomes its pseudo-inverse. Outside the linear regime, the following relation:

$$ \mathbf{a}(\delta\phi) = \mathbf{M}(\delta\phi; \phi_r)^{\dagger}\Delta I(\delta\phi; \phi_r), \qquad (19) $$

which ensures an accurate estimation of $\delta\phi$ from the reduced intensity, regardless of the phase amplitude.

Under this framework, the pseudo-inverse of the specific matrix provides an accurate reconstruction of the controlled modes of the phase $\delta\phi$, with the remaining error attributed solely to aliasing, which is negligible when using Fourier filtering WFSs.²⁷,²⁸

## 2.5 One-Dimension Toy Model to Obtain the Expression of the Specific Matrix

To understand the concept of the specific matrix, consider the 1D model. For exact phase estimation, the linearization around the reference phase (blue curve in Fig. 2) must be replaced by a linearization at the measurement point $\delta x + x_r$ (red curve in Fig. 2). This curve is locally linear and, therefore, invertible. Its slope corresponds to the 1D equivalent of the specific matrix
$$ \mathbf{M}(\delta\phi; \phi_r) \leftrightarrow \frac{f(x_r + \delta x) - f(x_r)}{\delta x}. \qquad (20) $$

This rate of change, which depends on $\delta x$, captures the behavior of the WFS for the specific phase being measured. In mathematical terms, this can be expressed as

$$ \frac{f(x_r + \delta x) - f(x_r)}{\delta x} = \int_0^1 dt f'(x_r + t\delta x). \qquad (21) $$

Extending this to the full WFS framework, where derivatives in 1D correspond to push-pull interaction matrices [see Eq. (15)], the specific matrix becomes

$$ \mathbf{M}(\delta\phi; \phi_r) = \int_0^1 dt \mathbf{M}(0; t\delta\phi + \phi_r). \qquad (22) $$

For the full mathematical demonstration of the specific matrix, the reader may refer to Appendix.

In practice, this integral can be approximated using a discrete sum over $n$ evenly spaced steps

$$ \mathbf{M}_n(\delta\phi; \phi_r) \equiv \frac{1}{n} \sum_{i=1}^{n} \mathbf{M}\left(0; \frac{i}{n+1} \delta\phi + \phi_r\right). \qquad (23) $$

## 2.6 Physical Interpretation

The specific matrix provides a linear representation of the WFS response for any phase amplitude. It is not an interaction matrix but an average of push-pull matrices computed at various operating points along the path from the reference phase $\phi_r$ to the measured phase $\phi = \phi_r + \delta\phi$.

This averaging process can be interpreted as a form of modulation or self-modulation. In traditional tip/tilt modulation, the interaction matrix is averaged over variable operating points induced by tip/tilt phases. Thus, the resulting interaction matrix is the sum of incoherent propagation. Similarly, by averaging push-pull matrices computed at different operating points, the specific matrix is also the result of an incoherent sum. Thus, the specific matrix represents the effect of self-modulation, capturing the changes in sensitivity caused by the wavefront itself.

## 3 First Results Using the Specific Matrix

Although the specific matrix formalism can be applied to any Fourier filtering WFSs, our focus in this study is on the PWFS. This choice is motivated by the PWFS being the most widely used Fourier filtering WFSs in current AO systems. Furthermore, the PWFS is expected to play a crucial role in the next generation of AO systems, making it a relevant and practical choice for our simulations and analysis.

By employing the numerical approximation of the specific matrix, as described in Eq. (23), it becomes feasible to evaluate its ability to characterize the behavior of the WFS within the intensity space. To achieve this, Fig. 3 compares the end-to-end (E2E) reduced intensity of an unmodulated PWFS considering a high-amplitude phase $\delta\phi$, following the Kolmogorov statistics, to intensity patterns obtained using different linear models based on various matrices:

* Interaction matrix: $\mathbf{M}\delta\phi = \Delta I(\delta\phi)$, representing the classic approach using the interaction matrix around the reference phase $\phi_r$.

* Push-pull matrix around the phase: $\mathbf{M}(0; \delta\phi + \phi_r)\delta\phi = \Delta I(\delta\phi)$, which employs the push-pull matrix around the phase being measured, $\mathbf{M}(0; \delta\phi + \phi_r)$. This matrix corresponds to an intuitive linear approximation of the input-output relationship and corresponds to the green dashed curve in Fig. 2.

* Average matrix: $\mathbf{M}_{av}\delta\phi = \Delta I(\delta\phi)$, defined as the average of the interaction and push-pull around the phase being measured matrices $\mathbf{M}_{av} = \frac{1}{2}(\mathbf{M} + \mathbf{M}(0; \delta\phi + \phi_r))$. This matrix serves as a compromise between the two aforementioned reconstructors.

* Specific matrix approximation: $\mathbf{M}_n(\delta\phi; \phi_r)\delta\phi = \Delta I(\delta\phi)$, based on linear models constructed from approximations of the specific matrix $\mathbf{M}_n(\delta\phi; \phi_r)$, where $n$ indicates the number of push-pull matrices used in the numerical computation (e.g., $n \in [1,5,10,20]$).
<figure><img src="image_5.png" /><figcaption>**Fig. 3** (a) Input phase $\delta\phi$: amplitude 1.51 rad RMS and E2E intensity; (b) intensity cut obtained with $M$ and $M(0; \delta\phi)$; (c) intensity cut obtained with and $M_{av}$, $M_1$, $M_5$, and $M_{20}$. Top: Input phase and corresponding E2E intensity. Middle: Intensity cuts obtained with the interaction matrix $M$ in rad dashed and with the push-pull matrix around the phase $M(0; \delta\phi)$ in blue. Bottom: Intensity cuts obtained with the average matrix $M_{av}$ in yellow dashed and with the different approximation of the specific matrix $M_n$ with $n \in \{1; 5; 20\}$.</figcaption></figure>

<table><caption>**Table 1** Simulation parameters.</caption><tr><td>Telescope</td><td>8 m</td></tr><tr><td>Pupil resolution</td><td>160 x 160 pixels</td></tr><tr><td>Phase basis</td><td>1000 Karhunen-Loève (KL) modes</td></tr><tr><td>Sensing path</td><td>$\lambda = 850$ nm - 40 x 40 subpupils in D</td></tr></table>

For these simulations, the reference phase $\phi_r$ is set to zero, and its contribution is omitted from the matrix expressions. Each reconstructor is computed using the Moore-Penrose algorithm<sup>29,30</sup> to derive the pseudo-inverse of the respective matrices. The WaveFront Error (WFE) is then evaluated for the controlled modes.

The simulation parameter can be found in Table 1.

These simulations demonstrate the following key insights:

* The classic linear interaction matrix, $M$, fails to accurately reproduce the WFS intensity [see Fig. 3(b), red dashed line] due to the WFS operating outside its linear regime. Consequently, its pseudo-inverse leads to substantial phase reconstruction errors, as shown in Table 2.

<table><caption>**Table 2** WFE in radian RMS using different reconstructors. The reference phase equals zero. The incident phase is shown in Fig. 3(a), and its RMS norm is 1.51 rad RMS.</caption><thead><tr><th>Reconstructer</th><th>$M^\dagger$</th><th>$M(0; \delta\phi)^\dagger$</th><th>$M_{av}^\dagger$</th><th>$M_1^\dagger$</th><th>$M_5^\dagger$</th><th>$M_{10}^\dagger$</th><th>$M_{20}^\dagger$</th></tr></thead><tbody><tr><td>WFE (rad RMS)</td><td>0.98</td><td>0.84</td><td>0.54</td><td>0.42</td><td>0.13</td><td>0.07</td><td>0.04</td></tr></tbody></table>
* The push-pull matrix around the phase, $M(0; \delta\phi)$ [see Fig. 3(b), blue dashed line], also fails to accurately replicate the E2E intensity, resulting in poor phase estimation performance (see Table 2).

* The average matrix, $M_{av}$, offers a marked improvement [see Fig. 3(c), yellow dashed line], achieving a better approximation of the WFS intensity. Its pseudo-inverse reduces the WFE by almost half compared with the interaction matrix (see Table 2), yielding results comparable to those achieved by the matrix $M_1$, which is computed around the median operating point $M_1 = M(0; \delta\phi/2)$ [see Fig. 3(c), green line].

* Increasing the number of steps, $n$, in the numerical computation of the discrete specific matrix $M_n$ further enhances the model's accuracy. For example, using $n = 5$, the WFE is reduced by a factor of nearly 10 compared with the interaction matrix (see Table 2), demonstrating the effectiveness of the specific matrix in representing the WFS's behavior. With $n = 20$, the reconstructor achieves unwrapped phase estimates [see Fig. 3(c), cyan dotted line].

It has been shown that it is possible to describe the nonlinearity of Fourier filtering WFSs in the matrix formalism. To do this, the matrix describing the WFS outside the low-phase regime has to depend on the phase being measured, which is the reason why this matrix is called the specific matrix. Its analytical expression has been given and shows that the specific matrix can be obtained as an average of push-pull matrices computed around different operating points. Such a fact confirms that nonlinearity may be understood as self-modulation. In the next section, several applications of this theoretical tool will be developed.

# 4 Practical Implementation of the Specific Matrix

This section focuses on the practical implementation of the specific matrix. Because the specific matrix depends on knowledge of the wavefront being estimated, its practical implementation is not feasible. The goal of this section is to develop a reconstructor that closely approximates the specific matrix, enabling more accurate and precise phase estimation.

The first example involves approximating the specific matrix by combining the OG formalism with information from a focal plane camera and a statistical approach to obtain the coupling terms. This method aims to construct a hybrid reconstructor that integrates frame-by-frame and statistical approaches to approximate the specific matrix at each iteration of the AO loop. Results for a SPHERE-type AO system will be presented. The final part discusses how to implement this tool in laboratory experiments.

## 4.1 Simulation of a Hybrid Phase Reconstructer

The dependency on the phase to be measured makes the practical use of the specific matrix difficult. Our approach consists of considering the specific matrix as the target reconstructor and exploring potential approximations that allow practical implementation. In particular, we will start by adapting the OG formalism described below. The classical OG approach, which generally only considers the diagonal of the gain matrix, will be completed by the knowledge of the cross-talk terms to build a fully dense gain matrix. This allows the coupling terms to be compensated. The challenge now consists of finding practical ways to obtain this full gain matrix on a frame-by-frame basis as the turbulence evolves dynamically.

The OG formalism was introduced to compensate for the loss of sensitivity in operation. To account for the change of regime between calibration and on-sky operation, the idea was to use another push-pull matrix that best describes the regime where the reconstructor will be used. In this manner, the new matrix is the one around the phase being measured $M(0; \delta\phi + \phi_r)$. Yet, instead of using this full matrix, an approximation of this matrix is used. Indeed, the information about the loss of sensitivity can be found in the diagonal of the matrix resulting from the projection of the matrix around the phase to the command matrix, as given in the following equations:

$$ T(0; \delta\phi + \phi_r) = M^{\dagger}M(0; \delta\phi + \phi_r), \qquad (24) $$

$$ OG(0; \delta\phi + \phi_r) = \text{diag}(T(0; \delta\phi + \phi_r)). \qquad (25) $$
This OG vector is then used to update the command matrix at each iteration. To make it feasible during operation, previous studies have suggested the use of an abacus containing the appropriate OG values depending on the observation conditions.²⁴,²⁵ It demonstrates the possibility of extracting the OG from the PSF image, also called a gain scheduling camera (GSC), thereby eliminating the need for interaction matrix computation during operation. Following this idea, we propose the use of a hybrid method, which consists of using a GSC for frame-by-frame estimation of the modal gains, complemented by a statistical approach to estimate the coupling terms using the specific matrix formalism.

### 4.1.1 Approximation of the specific matrix through the gain formalism

The initial step in constructing our hybrid reconstructor involves determining the modal gains. The formalism used here follows the approach using the matrix formalism for obtaining the OG²⁵ [cf Eq. (24)]. However, rather than employing the matrix computed around the operating point, the specific matrix is used. To compute the OG from the PSF image (hereafter referred to as OG<sub>GSC</sub>), we use the convolutional model¹²,²⁵ and, in particular, Eq. (14) from Ref. 25 and reminded here:

$$t^i = \frac{\langle \mathbf{IR}_\phi \star \phi_i | \mathbf{IR}_{\text{calib}} \star \phi_i \rangle}{\langle \mathbf{IR}_{\text{calib}} \star \phi_i | \mathbf{IR}_{\text{calib}} \star \phi_i \rangle},$$

where $t^i$ is the diagonal term of the gain matrix, $IR_\phi$ the impulse response of the sensor around the phase $\phi$, and $IR_{\text{calib}}$ the impulse response at calibration. The $\star$ stands for the convolutive product.

A key distinction in this study is that, instead of directly using the gains derived from the PSF, we adjust the gains. These are derived by combining the gains obtained from the GSC (OG<sub>GSC</sub>) with those from the calibration process. Previous work has shown that OG<sub>GSC</sub> provides an excellent estimation of the gains arising from a push-pull operation around the phase itself. However, we have demonstrated that the best estimation of these gains is achieved using the specific matrix.

The OG values should be understood as representing a loss of sensitivity induced by the wavefront itself. Consequently, accurate gains result from projecting the optimal reconstructor of that wavefront onto the interaction matrix. Within the matrix formalism, the specific matrix serves as the best reconstructor, as demonstrated in the previous section. Thus, the true OG must be derived from this matrix. Furthermore, we have shown in Ref. 31 that a better estimation of the gains is obtained from the specific matrix. An approximation for these gains can be obtained by averaging the interaction matrix and the push-pull matrix calculated around the operating point. This process is described as follows:

$$\begin{align}
\mathbf{OG}_{av} &= \frac{1}{2}(\mathbf{diag}(\mathbf{M}^\dagger\mathbf{M}) + \mathbf{diag}(\mathbf{M}^\dagger\mathbf{M}(0; \delta\phi + \phi_r))) \nonumber \\
&= \frac{1}{2}(\mathbf{diag}(\mathbf{I}_d) + \mathbf{OG}(0; \delta\phi + \phi_r)). \tag{26}
\end{align}$$

This represents a significant deviation from previous methods, as prior OG estimations were not optimal. Although the impact of these differences is minor when working with small phases ($\delta\phi \ll 1$), it becomes significant for systems operating around large phase amplitudes.

The simulations are conducted around a flat reference wavefront, $\phi_r = 0$, which simplifies the equations by removing the reference phase. The simulated wavefront aberrations are constructed using a power spectral density (PSD) model with a power-law exponent of -2.5 to emulate AO residuals. The resulting equivalent Strehl ratio is 55% in the H band.

Figure 4(a) shows the OG obtained through various methods for a modulated PWFS with a modulation radius of $r_{\text{mod}} = 3\lambda/D$. The OG values derived from the specific matrix (black curve) are higher than those obtained using the matrix $M(0; \delta\phi)$ around the residual phase (blue curve). In addition, the OG<sub>GSC</sub> values (green curve) align well with those estimated around the current operating point, consistent with previous results.²⁵ Finally, the averaged OG values (yellow curve), as computed from Eq. (26) and the adjusted OG obtained from the PSF (magenta curve), closely match those from the specific matrix, demonstrating the validity of this approach.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** (a) Modulated pyramid OG $r_{\text{mod}} = 3\lambda/D$. (b) Unmodulated pyramid OG. Average OG obtained for 20 residual phase screens. The equivalent Strehl ratio is 55% in the H band. The OG is obtained from the different matrices **OG**(0; δφ), **OG**(δφ; 0), and **OG**<sub>av</sub> and the focal plane camera **OG**<sub>GSC</sub>. The first 500 modes are shown for clarity. The shaded area represents the standard deviation.</figcaption></figure>

In the unmodulated case, shown in Fig. 4(b), the results are more nuanced. Although the average gains from Eq. (26) match the specific matrix-derived OG, discrepancies appear between **OG**<sub>GSC</sub> and the gains derived from the push-pull matrix around the residuals. This is attributed to the limitations of the convolutional model, which relies on linear assumptions and cannot fully account for the strong nonlinear behavior of the unmodulated pyramid. Nevertheless, **OG**<sub>GSC</sub> aligns with low-order modes of **OG**(0; δφ) and high-order modes of **OG**(δφ, 0). By adjusting **OG**<sub>GSC</sub> using Eq. (26), it is possible to achieve accurate estimations for low-order gains, though high-order estimations remain less precise.

### 4.1.2 Coupling terms

The second step toward constructing the hybrid reconstructor involves obtaining the full gain matrix, which includes the non-diagonal, i.e., coupling terms (cf. Sec. 2.3). The computation of the full gain matrix relies on simulations and a statistical approach. To achieve this, we consider a set $\mathcal{E}$ of residual phase screens representing various seeing conditions. Each set contains $N = 20$ phase screens, and for each one, a specific matrix is computed. The resulting gain matrices are derived using Eq. (24) from the averaged specific matrices over these phase realizations. In this manner, several gain matrices $T_i$, representative of different seeing conditions, are obtained.

An example of a full, dense gain matrix for a Strehl ratio of 42% in the H band is shown in Fig. 5(a). This matrix provides a detailed representation of the strength of the mode couplings under the given seeing conditions. To further illustrate this, Fig. 5(b) displays one line of the gain matrix, normalized by its diagonal term (the OG). The coupling strength varies between a few percent (~1% to 2%) and higher values (up to 10% of the OG), depending on the modes involved. Coupling is generally stronger for higher residual phases but remains proportional to the OG.

By averaging these gain matrices across the different phase screens in the set $\mathcal{E}$, it is possible to construct a representative gain matrix $\langle T \rangle_E$ for the given seeing conditions. This matrix accounts for the statistical coupling terms and provides a practical basis for constructing the hybrid reconstructor.

### 4.1.3 Hybrid reconstructor

Finally, for each residual amplitude, the average gain matrix $\langle T \rangle_E$ is obtained using the method described previously. This average matrix is subsequently updated with the OG corresponding to the current phase realization, derived from the GSC. As a result, the reconstructor used at each iteration closely approximates the specific matrix for the phase screen under consideration.
<figure><img src="image_7.png" /><figcaption>Fig. 5 (a) Dense gain matrix $T_1$. (b) Coupling of mode $KL_{30}$. Top: Full average gain matrix obtained for a phase set whose Strehl ratio is equal to 42% in H band. Bottom: Line 30 of gain matrices describing the coupling between the modes. The lines are normalized by the OG value. $T_1$ and $T_2$ are the average gain matrices obtained for a set of phases whose equivalent Strehl ratio is 85% and 42%, respectively, in the H band.</figcaption></figure>

This process involves normalizing the average gain matrix $\langle T \rangle_E$ by its diagonal terms, resulting in an intermediate matrix $T_{int}$ with ones on its diagonal and relative coupling strengths elsewhere. This intermediate matrix is then multiplied by the updated OG values extracted from the PSF image. The interaction matrix is subsequently updated using this dense gain matrix, as described in the following equations:

$$ 
\begin{aligned}
T_{\text{int}} &= \frac{\langle T \rangle_E}{\langle OG \rangle_E}, \\
T_{\text{updated}} &= \langle T \rangle_{\text{int}} \times OG_{\text{PSF}}, \\
R_{\text{updated}} &= (MT_{\text{updated}})^{\dagger},
\end{aligned}
\qquad (27) $$

where $\times$ denotes element-wise multiplication.

By integrating the statistical coupling terms from the averaged gain matrix with real-time OG estimations, this hybrid approach provides a practical method to dynamically adapt the reconstructor to the evolving turbulence. This ensures that the reconstructor remains aligned with the optimal performance of the specific matrix.

## 4.2 Frame-by-Frame Estimation of Phase Residual

This section focuses on the frame-by-frame estimation of residual phase screens with increasing amplitudes (the equivalent Strehl ratio in the H band is within the range SR $\in$ [30; 99]%). To evaluate this, we compare the reconstruction performance using full matrices ($M$, $M(\delta\phi; 0)$, $M(0; \delta\phi)$, $M_{av}$) and the hybrid matrices. The simulation parameters can be found in Table 1, we present the result for an unmodulated PWFS. The results are presented in Fig. 6(a), where different matrices are computed for each phase screen.

In Fig. 6(a), it is evident that the specific matrix (black curve) provides the most accurate results, as the WFE reaches the fitting error (cyan curve). Consistent with the results obtained in Sec. 2, in Fig. 3(b) and in Ref. 31, the matrix computed around the residual, $M(0; \delta\phi)$, offers an inaccurate phase reconstructor (blue dashed curve). To improve phase estimation, the average matrix $M_{av}$ (yellow curve) ensures a systematic reduction in WFE compared with the interaction matrix. The average matrix $M_{av}$ emerges as a viable alternative to the specific matrix, requiring less computational time and offering easier approximation through the gain formalism. As shown in Fig. 4, the average gain $OG_{av}$ is accurately estimated by adjusting the gains obtained from the PSF. These gains are closer to the real values derived from the specific matrix.

Another noteworthy reconstructor is the median matrix ($M_1$), computed around half the wavefront. This matrix represents the specific matrix of degree $n = 1$, and its WFE is plotted in magenta in Fig. 6(a). For small residuals, $M_1$ performs comparably to $M_{av}$. However, with
<figure><img src="image_8.png" /><figcaption>**Fig. 6** (a) WFE obtained with the full matrices. (b) WFE obtained with the hybrid matrix. WFE for residual phase screens, the equivalent Strehl ratio in the H band is within the range SR ∈ [30; 99]%. The shaded area represents the standard deviation.</figcaption></figure>

increasing amplitude, the median matrix outperforms $M_{av}$, reducing the WFE more effectively. The differences between $M_1$ and $M_{av}$ arise from the dependence of $M_{av}$ on the push-pull matrix around the phase, $M(0; \delta\phi)$. As the amplitude of the phase increases, $M(0; \delta\phi)$ becomes less representative of the linear relationship between intensity and wavefront, introducing errors into the estimation of $M_{av}$. Thus, in regimes of strong phase, $M_1$, corresponding to the push-pull matrix computed around half the wavefront, should be preferred over $M_{av}$.

Figure 6(b) depicts the results obtained using the hybrid approach. The green curve represents the WFE after estimating the phase screen with the hybrid reconstructor. In this case, the dense gain matrix is updated with the raw OG extracted from the GSC rather than the OG adjusted using Eq. (26). This choice was made because results without adjusting the gains were better. In regimes of strong turbulence or when using an unmodulated PWFS, the gain extracted from the GSC or its adjusted value either matches the specific gains of high-order modes or low-order modes, respectively. As a result, these gains are not fully representative of the true OG from the specific matrix due to the limitations of the convolutional approach [see Fig. 4(b)].

Although the hybrid approach does not achieve the same level of performance as the specific matrix, it allows for a 30% reduction in WFE compared with the classical approach. Moreover, its performance is comparable to that achieved using the median matrix $M_1$ or the average matrix $M_{av}$.

In addition, we performed simulations with a single gain matrix for the entire range of phase residuals. Updating the diagonal of this matrix produced results equivalent to those obtained with the average gain matrix at each amplitude. This demonstrates that the method robustly describes current seeing conditions using a single gain matrix while computing the current OG from the GSC.

The analysis of frame-by-frame estimation of phase residuals using hybrid methods inspires confidence in applying this approach to closed-loop systems with nonlinear WFS, such as the unmodulated PWFS. Updating the gain matrix enhances the robustness of the reconstructor, mitigating the risk of loop divergence.

## 4.3 Close-Loop Using the Hybrid Reconstructer

In this section, we evaluate the performance of the hybrid reconstructor, using the method described in the previous Sec. 4.1, in a closed-loop system with an unmodulated PWFS. Two scenarios are considered: (i) OG compensation only and (ii) full gain compensation using the hybrid approach. This analysis aims to understand the limitations of classical gain compensation and the capacity of our reconstructor to close the loop without modulation. In addition, it highlights the importance of the nondiagonal terms in the gain matrix, particularly for convergence under strong turbulence conditions. The performance of the unmodulated PWFS is compared with that of the classical modulated PWFS, which uses a standard interaction matrix. The AO loop parameters employed in these simulations are listed in Table 3.
<table><caption>**Table 3** Closed-loop parameters.</caption><tr><td>Frequency</td><td>1 kHz</td></tr><tr><td>Time of the loop</td><td>1 s</td></tr><tr><td>Frames delay</td><td>2</td></tr><tr><td>Integrator gain</td><td>0.4</td></tr><tr><td>Leak factor</td><td>0.99</td></tr><tr><td>Sensing wavelength</td><td>850 nm, I-band</td></tr><tr><td>Atmosphere</td><td>3 layers</td></tr><tr><td>Wind</td><td>v = 5 m.s<sup>-1</sup> ground layer</td></tr><tr><td>Noise</td><td>0</td></tr><tr><td>Science</td><td>H band</td></tr></table>

<table><caption>**Table 4** Statistical study of close-loop for three different $r_0$ made for 20 different atmosphere screens for each $r_0$. The average wind speed is $6.7 \text{ m.s}^{-1}$ considering the three atmospheric layers. The close-loop WFE and the standard deviation are specified in nm RMS. The average WFE is computed for the last 800 iterations.</caption><thead><tr><th rowspan="2">r<sub>0</sub> (cm)</th><th colspan="2">PWFS r<sub>mod</sub> = 0λ/D</th><th colspan="2">PWFS r<sub>mod</sub> = 3λ/D</th><th rowspan="2">Fitting</th></tr><tr><th>M</th><th>MOG</th><th>MT<sub>updated</sub></th><th>M</th></tr></thead><tbody><tr><td>10</td><td>144 ± 39</td><td>138 ± 38</td><td>94 ± 10</td><td>103 ± 11</td><td>86</td></tr><tr><td>7.9</td><td>390 ± 75</td><td>254 ± 66</td><td>145 ± 25</td><td>143 ± 21</td><td>110</td></tr><tr><td>6.9</td><td>602 ± 82</td><td>276 ± 50</td><td>210 ± 44</td><td>181 ± 36</td><td>124</td></tr></tbody></table>

The gain matrix used for the unmodulated PWFS simulations is built around residuals corresponding to a Strehl ratio of 30% in the H band. The matrix is updated using GSC images, as described in Eq. (27). Closed-loop simulations are conducted for three different $r_0$ values in V-band (@500 nm): $r_0 = \{6.9, 7.9, 10\}$ cm, corresponding to seeing conditions of $\{1.5''$, $1.3''$, $1.0''\}$ @ 500 nm, respectively. These harsh conditions are chosen to assess the robustness of the hybrid method.

The simulations are performed without noise, as the objective is to demonstrate the feasibility of closing the loop without modulation using the hybrid reconstructor. The star magnitude is therefore not a concern in these simulations. An in-depth analysis of noise impact has been conducted in Ref. 26, highlighting the benefits of the unmodulated PWFS for low-flux regimes due to its higher sensitivity.

For consistency, all simulations use a fixed loop gain of $g = 0.4$, allowing direct comparison between reconstructors under identical loop configurations. The results, summarized in Table 4, illustrate the residual WFE and standard deviations for 20 different atmospheric screens at each $r_0$ value. The average WFE is computed over the last 800 iterations of the loop.

For the median seeing condition 1", the unmodulated PWFS performs better than the modulated PWFS, achieving a 10-nm RMS reduction in WFE. This corresponds to a 3-point Strehl ratio improvement in the H band (from 85% to 88%). For seeing ", the hybrid method achieves performance comparable to the modulated PWFS. Under stronger seeing condition (1.5" $\le r_0 = 6.9$ cm), the modulated PWFS performs better than the unmodulated PWFS. This discrepancy is attributed to limitations in OG accuracy, which may inadequately reflect true specific gains, especially for strong turbulence. Indeed, they prove to be an approximation of the real specific gains, matching either the low-orders or the high orders or not matching at all
for strong turbulence (see Fig. 4). This mismatch may dominate here; however, it is important to highlight that with the proposed method it is possible to close the loop without modulation even under such harsh seeing conditions.

The hybrid reconstructor demonstrates the feasibility of closing the AO loop without modulation under challenging conditions. Although its performance in strong turbulence (1.5") does not surpass the modulated PWFS, it achieves a significant reduction in WFE under other conditions, highlighting the potential of the hybrid approach.

## 4.4 Discussion Around the Hybrid Reconstructer

This section summarizes the practical implementation of the hybrid reconstructor for on-sky operations. Furthermore, we explore the possibility of obtaining the specific matrix in practice.

### 4.4.1 Implementation of the hybrid reconstructor and OG tracking

For on-sky experiments, the dense gain matrix could be obtained using a simulation of the AO system to generate dense gain matrices under various seeing conditions. Moreover, to take into account the internal aberrations of the AO bench (optical elements defects), the computation of the different sets of specific matrices could be done on the AO bench directly if it is equipped with an internal calibration source. By applying AO residuals on the DM, all the additional error terms from the optics will be taken into account. Regarding the OG computation, using a GSC for frame-by-frame gain tracking can be computationally demanding. However, this study³² has shown that recalculating gains every few frames yields similar results to real-time calculations, thereby reducing computational requirements while maintaining system performance. With advanced real-time computers (RTCs), gain tracking can be conducted on millisecond time-scales, making this method feasible.

In Ref. 17, the authors demonstrate that the numerical complexity of their nonlinear iterative approach is comparable to that of linear methods, enabling a potential implementation at 1 kHz. Similarly, Ref. 14 explores the feasibility of running complex nonlinear reconstructors using deep learning and concludes that the numerical complexity remains within computational limits for operating an AO system at 1 kHz. Consequently, the proposed hybrid reconstructor—combing matrix-vector multiplication and convolutive products to compute the OG—can be implemented in a real system, as its complexity is lower than the aforementioned examples.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Closed-loop $r_0$ = 6.9 cm and average wind speed of $v$ = 6.7 m.s$^{-1}$, loop gain = 0.4. Comparison between the classical modulated PWFS $r_{\text{mod}}$ = $3\lambda/D$ blue-dashed curve and the different reconstructors using the nonmodulated PWFS. The classical interaction matrix in the red-dashed curve, the modal gain compensation in the green curve, and the full gain compensation in the yellow-dotted curve. The fitting error is plotted in cyan.</figcaption></figure>
### 4.4.2 Rainbow approach

It has been shown that obtaining the specific matrix requires access to push-pull matrices constructed around moving operating points from $\phi_r$ to $\phi_r + \delta\phi$. Two conditions are necessary: (i) performing push-pull calibrations faster than the input perturbation timescale. Because the specific matrix depends on the phase to be measured, it must be updated before the input changes. (ii) Accessing moving operating points: the calibration must be performed around the shifting operating points corresponding to the phase $\delta\phi$ with $t \in [0,1]$.

This paragraph focuses on an alternative approach to obtaining and computing the specific matrix by scaling the wavefront with the wavelength. Indeed, the computation of the specific matrix requires calculating push-pull matrices at different operating points. Using wavelength scaling, the moving operating points are defined as

$$t\delta\phi \text{ with } t \in [0,1]. \qquad (28)$$

By performing push-pull calibration of the WFS simultaneously at several wavelengths, the specific matrix can be computed. It should be noted that this approach is valid only for certain Fourier filtering WFSs, specifically those whose response is achromatic with respect to wavelength.

The unmodulated PWFS, operating around the zero-reference phase ($\phi_r = 0$), is an example of such an achromatic sensor. In this specific case, the specific matrix can be obtained from the push-pull matrices at various wavelengths using the following equations:

$$M(\delta\phi; 0) = \int_{0}^{1} dt M(0; t\delta\phi) = \int_{\lambda_0}^{\infty} \frac{\lambda_0}{\lambda^2} d\lambda M\left(0; \frac{\lambda_0}{\lambda} \delta\phi\right). \qquad (29)$$

where the parameter $t$ is linked to the ratio $\lambda_0/\lambda$, $\lambda_0$ is the sensing wavelength, and $\lambda$ represents the additional wavelengths chosen for push-pull calibration.

The use of a broadband detector, such as microwave kinetic inductance detectors (MKIDs) (e.g., Ref. 33), allows a single-phase residual to be applied to the DM. This setup enables the derivation of push-pull matrices across multiple wavelengths and would be made possible with the computation of the specific matrix.

The outlined approaches provide practical methods to construct the hybrid reconstructor or the specific matrix for real applications. On-sky, a combination of simulations and real-time gains tracking offers a feasible method for AO systems. By leveraging focal plane cameras (GSC) and dense gain matrices, it becomes possible to build an optimized reconstructor at each frame. In laboratory settings, push-pull matrices can be calibrated across a variety of phases and wavelengths. This allows stable closed-loop operation for nonlinear WFSs, such as the unmodulated PWFS, even under challenging seeing conditions.

## 5 Conclusion

This paper has demonstrated that the nonlinearity of Fourier filtering WFSs can be accurately described using a matrix formalism. The key lies in defining a specific matrix, which is constructed as the average of push-pull matrices around dynamic operating points. This matrix depends on the phase being measured. It accurately captures the loss of sensitivity resulting from the expansion of the electromagnetic field at the focal mask, making it possible to interpret nonlinearity as a self-modulation of the wavefront.

An algorithm has been presented for computing the specific matrix in numerical simulations, and its pseudo-inverse has been shown to provide accurate phase estimations over a wide range of input amplitudes. These results confirm that it is feasible to determine a matrix reconstructor even outside the WFS's typical dynamic range, which is particularly relevant for high-sensitivity WFSs, such as the unmodulated PWFS. The main limitation of this formalism is that it relies on the knowledge of the phase one tries to estimate.

For achromatic sensors such as the unmodulated PWFS, it is possible to obtain the specific matrix experimentally. This can be achieved by leveraging the scaling of phase with wavelength. In this approach, the push-pull interaction matrices are calculated at different wavelengths, corresponding to the dynamic operating points, and then used to compute the specific matrix. This method can even handle highly wrapped phases, thanks to its polychromatic basis.
However, practical implementation raises challenges due to the dependency of the specific matrix on the phase being measured, requiring real-time updates for each frame. To address these challenges, the OG formalism has been adapted. This redefinition of OG accounts for the self-induced over-modulation caused by the wavefront's amplitude. By extracting the accurate OG from the specific matrix and complementing it with cross-talk terms, the phase estimation accuracy is significantly improved.

Two practical approximations for the specific matrix have been proposed: the average matrix ($M_{av}$) and the median matrix ($M(0; \phi_r + \delta\phi/2)$). These approximations have been evaluated, demonstrating that the average matrix yields modal gains that closely match those of the specific matrix. The coupling terms, derived statistically for expected seeing conditions, enhance the performance of the reconstructor. By combining these gains and coupling terms, a hybrid reconstructor is created, blending statistical and frame-by-frame approaches.

Importantly, this hybrid method allows for stable closed-loop operation of an unmodulated PWFS under harsh seeing conditions, down to 1.3". Although compensating for nonlinearity might increase noise propagation in systems using faint guide stars, an optimal reconstructor that considers noise could be a valuable extension of this work. In the presented results, the noise was not included, allowing for simpler reconstructors to be implemented using data fusion with a focal plane camera to dynamically optimize the reconstructor during each loop iteration.

## 6 Appendix: Proof of the Specific Matrix Formula

For a Fourier filtering system with a tip/tilt modulation stage, the detector intensity may be obtained from the Fraunhofer diffraction laws

$$ I(\phi_r + \delta\phi)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'\mathbb{I}_P}|_{\vec{r}} e^{i(\phi_r|_{\vec{r}'}-\phi_r|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \hat{W}|_{\vec{r}'-\vec{r}} e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})}, \quad (30) $$

where $\mathbb{I}_P$ is the indicator function of the pupil, $m$ is the transparency function of the Fourier filtering mask, and $w$ is the modulation weighting function. $\vec{r}', \vec{r}$ and $\vec{R}$ are the 2D spatial variables. The reduced intensity is by definition

$$ \Delta I(\delta\phi; \phi_r) \equiv I(\phi) - I(\phi_r), \quad (31) $$

$$ =I(\delta\phi + \phi_r) - I(\phi_r), \quad (32) $$

which leads to

$$ \Delta I(\delta\phi; \phi_r)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'\mathbb{I}_P}|_{\vec{r}} e^{i(\phi_r|_{\vec{r}'}-\phi_r|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \hat{W}|_{\vec{r}'-\vec{r}} (e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})} - 1). \quad (33) $$

The last term is then modified as follows:

$$ e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})} - 1 = \frac{e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} 1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}}). \quad (34) $$

It allows us to split the integral of Eq. (33) into two parts

$$ \begin{aligned} \Delta I(\delta\phi; \phi_r)|_{\vec{R}} = & \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'\mathbb{I}_P}|_{\vec{r}} e^{i(\phi_r|_{\vec{r}'}-\phi_r|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \\ & \qquad \hat{W}|_{\vec{r}'-\vec{r}} \left( \frac{e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right) 1\delta\phi|_{\vec{r}'} - \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'\mathbb{I}_P}|_{\vec{r}} \\ & e^{i(\phi_r|_{\vec{r}'}-\phi_r|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}'} \tilde{m}|_{\vec{R}-\vec{r}} \hat{W}|_{\vec{r}'-\vec{r}} \left( \frac{e^{i(\delta\phi|_{\vec{r}'}-\delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right) 1\delta\phi|_{\vec{r}}, \end{aligned} \quad (35) $$

and then to invert $\vec{r}$ and $\vec{r}'$ in the first integral because they are dummy variables
$$
\begin{equation}
\begin{aligned}
\Delta I(\delta\phi; \phi_r)|_{\vec{R}} &= \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'} \mathbb{I}_P|_{\vec{r}} e^{i(\phi_r|_{\vec{r}} - \phi_{r'}|_{\vec{r}'})} \hat{m}|_{\vec{R}-\vec{r}} \hat{\tilde{m}}|_{\vec{R}-\vec{r}'} \\
&\quad \hat{w}|_{\vec{r}-\vec{r}'} \left( \frac{e^{i(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} - 1}{1(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} \right) 1\delta\phi|_{\vec{r}} - \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'} \mathbb{I}_P|_{\vec{r}} \\
&\quad e^{i(\phi_r|_{\vec{r}'} - \phi_{r}|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}} \hat{\tilde{m}}|_{\vec{R}-\vec{r}} \hat{w}|_{\vec{r}'-\vec{r}} \left( \frac{e^{i(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right) 1\delta\phi|_{\vec{r}}.
\end{aligned}
\tag{36}
\end{equation}
$$

It leads to

$$
\begin{equation}
\Delta I(\delta\phi; \phi_r)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'} (\mathbb{I}_P|_{\vec{r}} \delta\phi|_{\vec{r}}) 1 \left[ e^{i(\phi_r|_{\vec{r}} - \phi_{r'}|_{\vec{r}'})} \hat{m}|_{\vec{R}-\vec{r}} \hat{\tilde{m}}|_{\vec{R}-\vec{r}'} \right. \\
\left. \qquad \hat{w}|_{\vec{r}-\vec{r}'} \left( \frac{e^{i(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} - 1}{1(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} \right) - e^{i(\phi_r|_{\vec{r}'} - \phi_{r}|_{\vec{r}})} \hat{m}|_{\vec{R}-\vec{r}} \hat{\tilde{m}}|_{\vec{R}-\vec{r}} \right. \\
\left. \qquad \left. \hat{w}|_{\vec{r}'-\vec{r}} \left( \frac{e^{i(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right) \right]. \tag{37}
\end{equation}
$$

Moreover,

$$
e^{i(\phi_r|_{\vec{r}} - \phi_{r'}|_{\vec{r}'})} = \overline{e^{i(\phi_r|_{\vec{r}} - \phi_{r'}|_{\vec{r}})}}, \quad (38)
$$

$$
\hat{m}|_{\vec{R}-\vec{r}} \hat{\tilde{m}}|_{\vec{R}-\vec{r}'} = \overline{\hat{m}|_{\vec{R}-\vec{r}'} \hat{\tilde{m}}|_{\vec{R}-\vec{r}'}} \quad (39)
$$

$$
\left( \frac{e^{i(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} - 1}{1(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} \right) = \overline{\left( \frac{e^{i(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right)}. \quad (40)
$$

Finally, the tip/tilt weighting function w is a real function because it codes the time spent for
each tip/tilt. Consequently

$$
\hat{w}|_{\vec{r}-\vec{r}'} = \overline{\hat{w}|_{\vec{r}'-\vec{r}}}. \quad (41)
$$

As a conclusion

$$
\begin{equation}
\Delta I(\delta\phi; \phi_r)|_{\vec{R}} = \int_{\mathbb{R}^4} d^2\vec{r}' d^2\vec{r} \mathbb{I}_P|_{\vec{r}'} (\mathbb{I}_P|_{\vec{r}} \delta\phi|_{\vec{r}}) 1 \left[ e^{i(\phi_r|_{\vec{r}} - \phi_{r'}|_{\vec{r}'})} \hat{m}|_{\vec{R}-\vec{r}} \hat{\overline{m}}|_{\vec{R}-\vec{r}'} \right. \\
\left. \qquad \hat{w}|_{\vec{r}-\vec{r}'} \left( \frac{e^{i(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} - 1}{1(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} \right) - c.c. \right]. \tag{42}
\end{equation}
$$

The reduced intensity can therefore be written with the kernel formalism

$$
\Delta I(\delta\phi; \phi_r)|_{\vec{R}} = \int_{\mathbb{R}^2} d^2\vec{r} (\mathbb{I}_P|_{\vec{r}} \delta\phi|_{\vec{r}}) \mathbf{K}(\delta\phi; \phi_r)|_{\vec{R}, \vec{r}}, \quad (43)
$$

where the nonlinear kernel $\mathbf{K}(\delta\phi; \phi_r)$ equals

$$
\mathbf{K}(\delta\phi; \phi_r)|_{\vec{R}, \vec{r}} = 2 \mathrm{Im} \left[ \tilde{m}|_{\vec{R}-\vec{r}} \int_{\mathbb{R}^2} d^2\vec{r}' \mathbb{I}_P|_{\vec{r}'} \hat{m}|_{\vec{R}-\vec{r}'} \hat{w}|_{\vec{r}'-\vec{r}} \right. \\
\left. \qquad e^{i(\phi_r|_{\vec{r}'} - \phi_r|_{\vec{r}})} \frac{e^{i(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} - 1}{1(\delta\phi|_{\vec{r}'} - \delta\phi|_{\vec{r}})} \right]. \quad (44)
$$

It should then be noted that

$$
\frac{e^{i(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} - 1}{1(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})} = \int_{0}^{1} dt e^{it(\delta\phi|_{\vec{r}} - \delta\phi|_{\vec{r}'})}. \quad (45)
$$

Consequently, the nonlinear kernel may be written as an integral of a linear kernel with
a moving reference phase

$$
\mathbf{K}(\delta\phi; \phi_r) = \int_0^1 dt \mathbf{K}(0; \phi_r + t\delta\phi). \quad (46)
$$
This equation is exactly the same as Eq. (22), which gives the relation between the specific matrix and the push-pull interaction matrices. This is not surprising because kernels and interaction matrices are intrinsically linked. Indeed, any interaction matrix may be understood as the response of the kernel regarding a given phase basis

$$M_i|\bar{\mathbf{R}} = \int_{\mathbb{R}^2} d^2\bar{\mathbf{r}} (\mathbb{I}_P|\bar{\mathbf{r}}|\Phi_i|\bar{\mathbf{r}}) K|\bar{\mathbf{R}}, \quad (47)$$

In the specific matrix case, the previous equation becomes

$$M(\delta\phi; \phi_r)_i|\bar{\mathbf{R}} = \int_{\mathbb{R}^2} d^2\bar{\mathbf{r}} (\mathbb{I}_P|\bar{\mathbf{r}}|\phi_i|\bar{\mathbf{r}}) K(\delta\phi; \phi_r)|\bar{\mathbf{R}}, \quad (48)$$

which becomes with Eq. (46)

$$M(\delta\phi; \phi_r)_i|\bar{\mathbf{R}} = \int_0^1 dt \int_{\mathbb{R}^2} d^2\bar{\mathbf{r}} \mathbb{I}_P |\bar{\mathbf{r}}| \phi_i |\bar{\mathbf{r}}| K(0; \phi_r + t\delta\phi)|\bar{\mathbf{R}}, \quad (49)$$

$$= \int_0^1 dt M(0; \phi_r + t\delta\phi)_i|\bar{\mathbf{R}}. \quad (50)$$

This constitutes the mathematical proof of Eq. (22).

## Disclosures

The authors declare there are no financial interests, commercial affiliations, or other potential conflicts of interest that have influenced the objectivity of this research or the writing of this paper.

## Code and Data Availability

The simulations were done using the OOMAO Matlab toolbox.³⁴

## Acknowledgments

M.Cisse would like to thank Alexis Lau for her careful proofreading of the article. O. Fauvarque would like to thank L. de Raphélis-Soissan for his advice on the mathematical aspects and Romain Fétick and Cédric Taïssir Héritier for their very useful bibliographic suggestions. This work benefited from the support of the French National Research Agency (ANR) with Wave-front sensors for adaptive Optics on extremely Large telescope using Fourier filtering (Grant No. ANR-18-CE31-0018), APPLY (Grant No. ANR-19-CE31-0011), and LabEx Focal plane array for Universe sensing (Grant No. ANR-11-LABX-0013); the Programme Investissement Avenir F-CELT (Grant No. ANR-21-ESRE-0008); the Action Spécifique Haute Résolution Angulaire (ASHRA) of Centre National de la recherche scientifique/INSU co-funded by CNES, the Évaluation-Orientación de la Coopération Scientifique (ECOS) Comisión Nacional de Investigación Científica y Tecnológica (CONYCIT) France-Chile cooperation (Grant No. C20E02); the Opticon RadioNet Pilot-H2020 Framework Programme of the European Commission's (Grant No. 101004719); STIC AmSud (Grant No. 21-STIC-09); and the Région Sud and the french government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université -A*MIDEX (Grant No. AMX-22-RE-AB-151).

## References

1. R. Ragazzoni, "Pupil plane wavefront sensing with an oscillating prism," *J. Mod. Opt.* **43**(2), 289-293 (1996).

2. L. Schatz et al., "Three-sided pyramid wavefront sensor, part 1: simulations and analysis for astronomical adaptive optics," *J. Astron. Telesc. Instrum. Syst.* **7**(4), 049001 (2021).

3. O. Fauvarque et al., "General formalism for Fourier-based wavefront sensing," *Optica* **3**(12), 1440-1452 (2016).

4. B. Vohnsen, S. Castillo, and D. Rativa, "Wavefront sensing with an axicon," *Opt. Lett.* **36**(6), 846-848 (2011).

5. R. Clare et al., "Numerical evaluation of pyramid type sensors for extreme adaptive optics for the European extremely large telescope," in *Adapt. Opt. for Extremely Large Telesc.* **5** (2017).

6. M. N'Diaye et al., "Design optimization and lab demonstration of ZELDA: a Zernike sensor for near-coronagraph quasi-static measurements," *Proc. SPIE* **9148**, 91485H (2014).
7. V. Chambouleyron et al., "Variation on a Zernike wavefront sensor theme: optimal use of photons," *Astron. Astrophys.* **650**, L8 (2021).

8. M. Cisse et al., "The phase-shifted Zernike wave-front sensor," *Proc. SPIE* **12185**, 121850X (2022).

9. O. Guyon, "Limits of adaptive optics for high-contrast imaging," *Astrophys. J.* **629**(1), 592 (2005).

10. F. Rigaut and E. Gendron, "Laser guide star in adaptive optics-the tilt determination problem," *Astron. Astrophys.* **261**(2), 677-684 (1992).

11. V. Chambouleyron et al., "Optimizing Fourier-filtering WFS to reach sensitivity close to the fundamental limit," *Proc. SPIE* **12185**, 121852T (2022).

12. O. Fauvarque et al., "Variation around a pyramid theme: optical recombination and optimal use of photons," *Opt. Lett.* **40**(15), 3528-3531 (2015).

13. C. Weinberger et al., "Design and training of a deep neural network for estimating the optical gain in pyramid wavefront sensors," in *Adapt. Opt. and Appl.*, Optica Publishing Group, p. JF1B-6 (2022).

14. R. Landman and S. Haffert, "Nonlinear wavefront reconstruction with convolutional neural networks for Fourier-based wavefront sensors," *Opt. Express* **28**(11), 16644-16657 (2020).

15. J. Nousiainen et al., "Toward on-sky adaptive optics control using reinforcement learning-model-based policy optimization for adaptive optics," *Astron. Astrophys.* **664**, A71 (2022).

16. F. Archinuk et al., "Mitigating the nonlinearities in a pyramid wavefront sensor," *J. Astron. Telesc. Instrum. Syst.* **9**(4), 049005 (2023).

17. V. Hutterer and R. Ramlau, "Nonlinear wavefront reconstruction methods for pyramid sensors using Landweber and Landweber-Kaczmarz iterations," *Appl. Opt.* **57**(30), 8790-8804 (2018).

18. R. A. Frazin, "Efficient, nonlinear phase estimation with the nonmodulated pyramid wavefront sensor," *JOSA A* **35**(4), 594-607 (2018).

19. V. Chambouleyron et al., "Using the Gerchberg-Saxton algorithm to reconstruct nonmodulated pyramid wavefront sensor measurements," *Astron. Astrophys.* **681**, A48 (2024).

20. R. V. Shack, "Production and use of a lenticular Hartmann screen," in *Spring Meeting of Opt. Soc. of Am. Vol.* 656 (1971).

21. O. Fauvarque et al., "General formalism for Fourier-based wavefront sensing: application to the pyramid wavefront sensors," *J. Astron. Telesc. Instrum. Syst.* **3**(1), 019001 (2017).

22. O. Fauvarque et al., "Kernel formalism applied to Fourier-based wave-front sensing in presence of residual phases," *JOSA A* **36**(7), 1241-1251 (2019).

23. V. Korkiakoski, C. Vérinaud, and M. Le Louarn, "Improving the performance of a pyramid wavefront sensor with modal sensitivity compensation," *Appl. Opt.* **47**(1), 79-87 (2008).

24. V. Deo et al., "A telescope-ready approach for modal compensation of pyramid wavefront sensor optical gain," *Astron. Astrophys.* **629**, A107 (2019).

25. V. Chambouleyron et al., "Focal-plane-assisted pyramid wavefront sensor: enabling frame-by-frame optical gain tracking," *Astron. Astrophys.* **649**, A70 (2021).

26. G. Agapito et al., "Non-modulated pyramid wavefront sensor-use in sensing and correcting atmospheric turbulence," *Astron. Astrophys.* **677**, A168 (2023).

27. C. Vérinaud, "On the nature of the measurements provided by a pyramid wave-front sensor," *Opt. Commun.* **233**(1-3), 27-38 (2004).

28. C. M. Correia et al., "Performance limits of adaptive-optics/high-contrast imagers with pyramid wavefront sensors," *Mon. Not. R. Astron. Soc.* **495**(4), 4380-4391 (2020).

29. E. H. Moore, "On the reciprocal of the general algebraic matrix," *Bull. Am. Math. Soc.* **26**, 394-395 (1920).

30. R. Penrose, "A generalized inverse for matrices," *Proc. Cambridge Philos. Soc.* **51**(3), 406-413 (1955).

31. M. Cisse et al., "Description of the non-linear behaviour of Fourier filtering wavefront sensor using matrix formalism: the specific matrix," *Proc. SPIE* **13097**, 130970W (2024).

32. A. Striffling et al., "Toward the full control of NCPA with the pyramid wavefront sensor: mastering the optical gains," in *Adapt. Opt. for Extremely Large Telesc., 7th edition* (2023).

33. E. O'Connor, A. Shearer, and K. O'Brien, "Energy-sensitive detectors for astronomy: past, present and future," *New Astron. Rev.* **87**, 101526 (2019).

34. R. Conan and C. Correia, "Object-oriented Matlab adaptive optics toolbox," *Proc. SPIE* **9148**, 91486C (2014).

**Mahawa Cisse** is an optical engineer, who graduated from the Institut d'Optique Graduate School in Paris-Saclay in 2021. Recently, she obtained her PhD in 2024 at l'Université Aix-Marseille, working at ONERA and LAM in Marseille, France. Her PhD research focused on adaptive optics, wavefront sensing, and nonlinear wavefront sensors, such as the Fourier filtering wavefront sensor. She is a member of SPIE.

**Olivier Fauvarque** is a researcher in optical instrumentation for oceanography at Ifremer, France. At the same time, he continues his activities in wavefront analysis for adaptive optics, which he started during his thesis at the Laboratoire d'Astrophysique de Marseille (LAM).
**Vincent Chambouleyron** is a postdoc at UC, Santa-Cruz, USA, working in the adaptive optics groups. He did his PhD at LAM working on highly sensitive wavefront sensors.

**Nicolas Levraud** is a postdoc at LESIA, Meudon, France, working in the adaptive optics groups. He did his PhD at ONERA, LAM, and INAF Arcetri working on wavefront sensing for Extremely Large Telescopes.

**Charlotte Z. Bond** is an adaptive optics engineer at UKATC, Edinburgh, working on an adaptive optics system for the Extremely Large Telescope.

**Benoit Neichel** is a senior adaptive optics scientist currently employed at Laboratoire d'Astrophysique de Marseille. He obtained his PhD in 2008 between ONERA and Observatoire de Paris, then employed for 5 years at the Gemini South telescope to implement the first facility MCAO system and joined LAM in 2013.

**Jean-François Sauvage** is a senior adaptive optics scientist currently employed at ONERA and welcomed at LAM in the frame of the ONERA-LAM team. He obtained his PhD in 2007. He then worked on adaptive optics for defense for 5 years and then worked on the SPHERE extreme AO system and its installation on Paranal Observatory in 2014. He then joined LAM in 2015.

**Thierry Fusco** received his PhD (partial correction and anisoplanatism in AO) in 2000 from Nice University. He has been employed at ONERA since 2001 where he is now a research director. He is scientifically responsible for adaptive optics activities.