Using the Gerchberg–Saxton algorithm to reconstruct nonmodulated pyramid wavefront sensor measurements

,


Introduction
The pyramid wavefront sensor (PWFS; Ragazzoni 1996) is a Fourier-filtering wavefront sensor (Fauvarque 2017).These sensors are commonly used to measure aberrations in optical systems.Inspired by the Foucault knife test, the original PWFS consists of a four-sided glass pyramid located at an intermediate focal plane and a detector that captures images of the four pupils created by each beam passing through the different faces of the pyramid.This configuration efficiently converts phase into intensity, but it lacks the necessary dynamic range to accurately measure atmospheric turbulence aberrations, which can induce optical path differences on the order of several waves.To address this issue, the PWFS is often paired with a modulator that causes the electromagnetic (EM) field to circulate around the pyramid tip during the camera acquisition time.Modulation drastically improves the PWFS linearity, but it has three main drawbacks.First, the modulation leads to a strong decrease in PWFS sensitivity, especially for low-order modes.Second, modulation of the PWFS alters the nature of the signal (Vérinaud 2004;Guyon 2005), causing the response to resemble that of a slope sensor, which makes it difficult to detect phase discontinuities.This is particularly problematic for the next generation of giant segmented-mirror telescopes (GSMTs), where wavefront control will involve correcting not only for turbulence-induced aberrations, but also for those induced by the telescope itself (fragmentation, Schwartz et al. 2017;Bertrou-Cantou et al. 2022;Demers et al. 2022 and segmentation Chanan et al. 1998).Finally, adding the modulation mirror stage requires the use of more optics and leads to difficulties related to fast-steering components (stability issues, temperature constrains, speed limitations, failure risks, etc.).Therefore, extending the dynamic range of the nonmodulated PWFS (nPWFS) to remove the need for modulation could allow the use of the PWFS at its full potential while removing the requirement for moving parts, making it a less complex system.Many studies have been conducted to mitigate PWFS nonlinearities, and several options have been suggested.It was proposed to keep the matrix formalism and consider the PWFS as a linear varying-parameter system: the reconstruction matrix evolves according to the phase to be measured (Korkiakoski et al. 2008;Deo et al. 2019), but this technique usually requires some knowledge of the statistics of the measured phase (Chambouleyron et al. 2020).Gradient-descent methods have also been investigated (Frazin 2018;Hutterer et al. 2023).Finally, another approach that is more appealing today is to use a machine-learning approach to reconstruct the nPWFS signal (Landman & Haffert 2020;Nousiainen et al. 2022;Archinuk et al. 2023).
The goal of this paper is to present a novel approach to invert the nPWFS signal.This approach is based on the principle of reciprocity of light propagation and the Gerchberg-Saxton (GS) algorithm.In Sect.2, we present the principle of the reconstruction algorithm in detail.Section 3 assesses the basic performance of the reconstructor mainly in terms of dynamic range, but also in terms of noise propagation and broadband performance.Section 4 highlights the results of an experimental implementation of this new reconstructor on the Santa Cruz Extreme Adaptive optics Laboratory (SEAL) testbed (Jensen-Clem et al. 2021).Finally, Sect. 5 describes a possible way to push the dynamic range and convergence speed of this reconstructor even further by using sensor fusion.

Reciprocity of the light-propagation principle
The reconstructor presented in this paper relies on one of the most basic properties of light propagation: the principle of reciprocity of the light propagation.The idea is to construct a high-fidelity numerical model of the nPWFS and use the measurements to send the light backward in the numerical nPWFS.This technique can be applied to any Fourier-filtering wavefront sensor, but we focus on the nPWFS in this study.
The numerical model of the nPWFS was built by propagating the light assuming Fraunhofer approximation from one pupil plane to the next while it passes through the pyramid mask in an intermediate focal plane.In more detail, the EM-field in the nPWFS detector plane can be written as where Ψ p is the EM-field in the entrance pupil, m is the complex shape of the pyramid mask, ⋆ is the convolution product, and • denotes the Fourier-transform operator.Equation (1) allows us to simply simulate the light propagation from the entrance pupil plane to the WFS detector plane.The back-propagation of light from the detector plane to the entrance pupil plane is where • is the conjugate operator.This last equation can simply be understood in the nPWFS case: we write m as m = e i∆ , where ∆ is a 2D real function describing the phase corresponding to the pyramid shape.Passing the pyramid in the opposite direction means that the light propagates through the inverse phase mask, that is, m = e −i∆ .We can then write the entrance pupil and detector EM-fields in their complex form, where A p and A d are the amplitudes, and ϕ p and ϕ d are the phases of the electromagnetic fields.The goal of wavefront sensing is to recover the phase ϕ p in the entrance pupil.Light back-propagation cannot easily be performed from the detector measurements because we only have access to intensities in the detector plane The phase ϕ d is therefore missing in the measurements.Hence, we propose to use an iterative algorithm called the Gerchberg-Saxton (GS) algorithm to propagate the light back and forth in the numerical model of the nPWFS.

Gerchberg-Saxton algorithm
The GS algorithm was first proposed by Gerchberg (1972) and is widely used to perform image sharpening from point spread function (PSF) images (Fienup 1982;Ragland et al. 2016).To perform this algorithm in our case, we assumed that we have access to a measurement of the entrance pupil amplitude A p (we show an easy practical way to obtain this quantity with the nPWFS below).In the nPWFS framework, we have two complex quantities (Eq.( 3)) for which the amplitudes are known and with a relation that links them (Eq.( 1)).The principle of the GS algorithm is to propagate the light back and forth in the numerical model of the nPWFS, injecting the knowledge of the amplitudes of the complex quantities we are trying to retrieve at each step.We propose to go through one iteration of the GS algorithm applied to the nPWFS based on the schematic given in Fig. 1.The amplitude for the entrance pupil and the detector plane used for this example are true measurements from the SEAL testbed (highlighted with yellow dots).One iteration of the GS algorithm can be split into four parts: 1. We compute the detector EM-field complex amplitude A p = I p .For the first iteration, the complex EM-field in the detector plane is built by using ϕ d = arg(A p ⋆ m), which corresponds to the phase in the detector plane when a flat wavefront is propagated through the nPWFS system.We therefore have a first estimate Ψ d that can be back-propagated in the system (through Eq. ( 2)). 2. A first estimation of Ψ p is then obtained.3.Because we already have access to A p , the amplitude found through back-propagation is discarded and replaced by the measurement of A p while keeping the estimated phase ϕ p .
The entrance pupil plane EM-field can then be propagated in the system (direct propagation, through Eq. ( 1)). 4. A new estimation of Ψ d is obtained.As previously done for the entrance pupil plane, we discard the amplitude and replace it by the measurement of A d given by the detector and keep the estimated phase ϕ d .We can then return to step 1 and iterate again.We call one iteration of the GS algorithm the numerical operation that consists of these four steps.The GS algorithm is therefore an iterative algorithm, which for one iteration performs four fast Fourier transforms (FFTs).It is important to note that this reconstructor assumes coherence of the EM-field.Therefore, this algorithm does not work for the modulated PWFS.It also raises the question of the impact of measurements with larger spectral bandpasses, which we discuss in the next section.We also emphasize that this GS principle can be applied to any FFWFS, but we focus on the nPWFS here.

Introduction
Detector plane Entrance pupil plane Amplitude Phase I p e i ‹ m " f p qe igrf p qs (2) " I ´Ip (5) " ´I ´Ip (6) August 2021 1 Introduction

Detector plane
Entrance pupil plane Amplitude Phase August 2021 1 Introduction

Detector plane
Entrance pupil plane Amplitude Phase I p e i ‹ m " f p qe igrf p qs (2) " I ´Ip (5) " ´I ´Ip (6) August 2021

Introduction
Detector plane Entrance pupil plane Amplitude Phase Discard Keep I p e i ‹ m " f p qe igrf p qs (2) " ´I ´Ip (6) I p e i ‹ m " f p qe igrf p qs (2) " ´I ´Ip (6) scope of this study, and therefore, we restrain ourselves to using a well-known simpler algorithm based on Ghiglia & Romero (1994) throughout this paper.

Performance of the GS algorithm reconstructor
This section gives a first overview of the basic performance of the GS algorithm.We focus on the performance of this reconstructor in terms of dynamic range, but also evaluate the number of iterations needed, the noise propagation, and the broadband impact on reconstruction.This study is not meant to be exhaustive, and parameters such as the minimum number of pixels needed on the detector with respect to the phase amplitude to be measured are not assessed, nor do we offer a detailed analysis of the impact of model errors.All the simulations simply match the SEAL testbed configuration.The main simulation parameters are the following ones: each pupil on the PWFS detector has 106 pixels across (matching the SEAL testbed configuration), which is a realistic case because the best low read-out-noise (RON) cameras can run at a few kHz with resolutions of about 250 px × 250 px.We use a Shannon sampling of 2 (4 pixels per λ/D) for the nPWFS model.For a more realistic simulation, the phase screens were first simulated and propagated with a resolution four times higher (4x106 px across) through a highresolution model of the nPWFS with a Shannon sampling of 2 as well.The signal was then binned to produce the nPWFS image with 106 px across.We worked with a modal basis of the first 500 Zernike modes (excluding piston), which is close to the control space that can be achieved with the deformable mirrors installed on the SEAL testbed.The linear reconstructor was created in a standard fashion by building an interaction matrix and computing its pseudo-inverse.No modes were filtered during the inversion because the measurements are oversampled.This led to a good conditioning number of the interaction matrix.The simulations tools we used were sourced internally, except for the turbulent phase screens, which were created using the HCIpy library (Por et al. 2018).

Convergence speed and first comparison with a linear reconstructor
The first test we present was meant to provide an idea of the convergence speed of the GS reconstruction.It is also be a way to draw a first comparison between the linear model and our reconstructor.To do this, we studied the reconstruction of four different turbulent phase screens following a Von-Karman spectrum and corresponding to four different configurations of D/r 0 , where D is the telescope diameter, and r 0 is the Fried parameter.To produce a fair comparison between the linear reconstructor and the GS reconstructor, which gives the phase with a much higher resolution (106 px across), these screens were projected on the 500 Zernike modes basis before propagation (no aliasing).The reconstruction error was estimated at each GS iteration (the estimated phase was systematically unwrapped) and is expressed as the ratio of the rms error in the entrance pupil plane (error = rms rec /rms input ).
The results are plotted in Fig. 2. The GS reconstruction accuracy increases with the number of iterations before it stabilizes.The shorter the phase, the fewer iterations are needed to reach the best reconstruction.Hence, for D/r 0 = 1, only 10 iterations are required, whereas more than 1000 iterations are required for D/r 0 = 32 (corresponding to a value of r 0 = 25 cm for a D = 8 m telescope).This figure also shows that in any case and after a few iterations, the GS reconstructor performs better than the linear reconstructor.In the small-phase regime (in our context, this would correspond to regimes smaller than D/r 0 ≈ 0.5), both perform equally well because the nPWFS would be working in its linear range.To illustrate reconstruction products, the input phase and reconstructed phases (linear and GS after 2000 iterations) for D/r 0 = 32 are shown in Fig. 3.It is clear that even if the GS reconstruction still underestimates the phase (due to nPWFS saturation; this is discussed in more detail in Sect.4.3), the shape of the retrieved phase is much more similar than the shape produced by the linear reconstructor.As a side note, the phase of the shape given by the linear reconstructor is highly underestimated, and the pattern seems quite different from the A48, page 3 of 11 Detector plane Entrance pupil plane Amplitude Phase Discard Keep input phase with higher spatial frequencies, suggesting that it would be hard to start a closed loop with this reconstruction.As the nPWFS is sensitive to phase discontinuities, the linear reconstruction shape might partially be explained by the impact of phase wrapping on the measurements.A detailed analysis of this phase-wrapping effect on the linear reconstruction is beyond the scope of this paper, but remains an important point to understand nPWFS behavior.
To illustrate how the GS-reconstructed phases evolve with the iteration, the reconstructed phases for iterations 1, 200, and 2000 are shown in Fig. 4 for the strongest turbulent case D/r 0 = 32.The top row presents the wrapped phase given by the GS algorithm, and the bottom row shows the corresponding unwrapped phases.The algorithm appears to quickly converge toward the good overall wrapped shape and then improves the estimation by scaling the phase closer from the input amplitude (therefore increasing the phase wrapping).For Fig. 2, we used seeing-limited phase screens as examples for a first insight into the behavior of the GS reconstructor.In a true AO system, the nPWFS would typically work in closed loops around residual phases, so that typical AO residual phases could also have been used for this previous analysis.We instead decided to use a full power-law turbulence phase screen because (i) the closed-loop bootstrapping is performed on full turbulence in any   case, and (ii) using AO residual phase screens would have been required to choose a specific system configuration.In order to refine the comparison between linear and GS reconstructors, it was decided to build linearity curves for Zernike modes.The results are presented in the next section.Finally, from this first analysis, the GS reconstructor apparently outperforms the linear reconstructor, but it requires tens of iterations to be effective.Each propagation back and forth requires four 2D FFTs (Fig. 1).In the current implementation, numpy.fft.fft2 was used to compute the Fourier transforms for a problem size of [424,424] (106 pixels across the pupil, with a sampling of twice Shannon), averaging at 7.1 ± 0.4 ms for each FFT on the SEAL control computer.Although the 2D FFTs algorithm scales as O(N 2 log N) overall, problem sizes that are a factor of 2 larger or that are divisible by large powers of 2 or 3 can run significantly faster.Padding to 448 = 2 6 × 7 improves the numpy.fft.fft2runtime to 3.2 ± 0.07 ms.With the downsampled problem size of [212,212] (which could easily be reached by decreasing the sampling and reducing the number of subapertures) and a more efficient FFT algorithm called FFTW (Frigo 1999), the runtime decreases to 640 ± 27 µs, which is an order of magnitude compared with the initial runtime.This represents an efficiency of about 5500 mflops, which is below the stated achievable performance for similar CPUs of ∼12 000 mflops.It is likely that an improved CPU would be able to achieve this performance, enabling us to run each FFT at the full problem size in less than 1 ms.Furthermore, it is possible that a lowerlevel implementation of the GS algorithm would be able to make use of the repeated FFT per iteration, for example, by optimizing memory access or by varying the FFT problem size per iteration.Based on prior benchmarking conducted for the CUDA-based cuFFT implementation (Kunkel et al. 2017), GPU computation could reduce this further by up to an order of magnitude.This could allow us to run a few iterations of the algorithm within a millisecond, and would open the path for closed-loop operation.

Linearity curves and plots of the dynamic range
To achieve a more quantitative comparison in terms of the dynamic range between the linear and the GS reconstructors, linearity curves for some Zernike modes were analyzed.Zernike modes within a full range of amplitudes were sent through the PWFS and were then reconstructed for four different cases: (i) nPWFS with a linear reconstructor (ii) modulated PWFS with a modulation radius of 3λ/D as a comparison point in terms of dynamic range, (iii) nPWFS with a GS reconstructor without the unwrapping step, and (iv) nPWFS with a GS reconstructor with phase unwrapping.The results for four different Zernike modes ranging from low to high spatial frequency (Z 6 , Z 19 , Z 150 , and Z 490 ) are plotted in Fig. 5.In this case, the GS reconstructor was used with an arbitrary number of 25 iterations (the impact of the number of iterations is assessed below).These linearity curves show the well-known behavior of the modulated PWFS with respect to the nPWFS: the increase in dynamic range for loworder modes located within the modulation radius is significant, and the dynamic range for high-order modes composed of spatial frequencies outside the modulation radius is comparable.The GS algorithm without the unwrapping algorithm shows an extended linear range for all the modes, but the response curves drop steeply after 1 rad rms of aberrations.This corresponds to the amplitude for which the phase starts to wrap (PtV greater than 2π), and therefore, the reconstructed phase starts to have a different shape from the unwrapped phase and the projection on Zernike modes is modified as more high-order modes are introduced.However, adding the unwrapping step to the reconstruction allows us to improve the linearity curves after 1 rad rms of input phase.It is clearly demonstrated here that the GS reconstructor with the unwrapping step allows us to significantly increase the nPWFS dynamic range for all modes.
For the plots in Fig. 5, it was arbitrarily chosen to run 25 iterations for the GS algorithm.The impact of the number of iterations on the linearity curve of the Zernike mode Z 19 is assessed in Fig. 6, where the linearity curve for this mode is shown using 1, 25, and 200 iterations when the GS algorithm is combined with the unwrapping step.For one iteration, the linearity curve does not follow the y = x slope around a null phase.This shows that even in a small-phase regime, only one iteration of the algorithm is not enough to accurately recover the phase.We note that two 2 iterations are enough in the small-phase regime (in the ideal case of these simulations, where there is no model error).As expected, the linearity curves then improve with the number of iterations, especially for the strong-aberration regime.
The linearity curves provided in Figs.reconstruction.Another approach to show the improved dynamics offered by the GS reconstructor is to generate dynamic range plots as proposed in Lin et al. (2022).In this method, we randomly introduced sampled aberrations based on Von Karman atmospheric power spectral density (PSD) into the PWFS system and vary the total wavefront error from 0 to 3 rad rms.For each wavefront error value, we injected and reconstructed 100 phase screens.This process was repeated for the nPWFS with the linear reconstructor, the GS algorithm with unwrapping (25 and 200 iterations), and the modulated PWFS (r mod = 3 λ/D).The mean values of the 100 reconstructions for each input amplitude for all configurations are plotted in Fig. 7.The figure reveals that the GS algorithm outperforms the linear reconstructor for the nPWFS.When we consider a satisfying reconstruction threshold as the point where reconstruction error reaches 10% of the input rms, then the GS algorithm extends the linearity range by a factor of approximately 3.Moreover, it achieves a dynamic performance comparable to the 3 λ/D modulated PWFS for input wavefront errors of up to 1.5 rad rms for 25 GS iterations and up to 2.3 rad rms for 200 GS iterations.For this test, the modulated PWFS shows the best linearity as expected from the analysis of the linearity curves presented Fig. 5: the modulated PWFS has the best dynamic range for low-order modes (which are affected by the modulation), and the input turbulent phases used to produce the curves in Fig. 7 have a higher amplitude for low-order modes.
We demonstrated that the GS algorithm combined with an unwrapping step can drastically improve the nPWFS dynamic range.Hence, it is possible to use this technique for reconstructing larger amplitude phases more accurately.Still, it is important to analyze the noise propagation through this reconstructor with respect to the linear framework because one of the motivations to use the nPWFS is the better sensitivity with respect to noise.

Noise propagation through reconstruction
To investigate the way in which noise spreads across different reconstructors, the four configurations outlined in the preceding section were used: an nPWFS with a linear reconstructor, a modulated PWFS, an nPWFS with a GS reconstructor without phase unwrapping, and an nPWFS with a GS reconstructor combined with phase unwrapping.For this study, only photon noise " ✏ " ´✏ " IM ´1.sp q pokes " M modes2pokes .modes GS reconstructor after 12 iterations Linear reconstructor after 12 iterations propagation was analyzed because it is a fundamental noise that cannot be mitigated by technological improvement.The effect of noise propagation was simply evaluated by introducing noise into the measurements and reconstructing the signals considering the different reconstructors.Once again, the idea is not to run an extensive simulation study on noise propagation through the GS algorithm in order to assess a large parameter space, but rather to highlight the basic sensitivity performance of this reconstructor with respect to the nPWFS and the modulated PWFS.To keep the analysis relatively simple, we studied two configurations: noise impact on the reconstruction of a flat wavefront (often referred to as reference intensities), and in the case of a turbulent phase screen with an amplitude of 0.75 radians rms (this is beyond nPWFS linearity range, but low enough to ensure that phase is not wrapped).To assess the noise propagation relative to these configurations, the following procedure was used in the simulation: for a given number of photons available for the measurement, the mean variance of each mode after reconstruction was estimated by averaging the variance of the reconstructed modes over 200 noise realizations.This operation was repeated for different numbers of photons, setting the number of iterations in the GS algorithm to 25.This number of iterations was once again arbitrarily chosen because we found that the number of iterations does not change the impact of the noise propagation.
The phase reconstruction errors due to photon noise propagation in the case of a flat wavefront and a 0.75 radians rms turbulent screen with respect to the incident flux are plotted in Fig. 8.The incident flux is expressed as the total number of photons available for measurements assuming a prefect transmission and a quantum efficiency of one for the detector (as a reminder, the reconstruction was performed over 500 Zernike modes).
For the cases of the linear reconstructor and a flat wavefront, Fig. 8 (top) shows a well-known behavior: the phase estimation error due to noise propagation through the reconstructor is lower for the nPWFS than for the modulated one.In the case of the GS algorithm, noise propagation behaves differently whether the unwrap step is performed or not.
In the context of noisy measurements for a flat wavefront, when unwrapping is omitted, the GS algorithm still propagates more error than the nPWFS and modulated PWFS for the highflux regime, and it seems to perform slightly better than the modulated PWFS for the low-flux regime.For even lower flux Input amplitude (rad rms) regimes (not shown here), GS would even appear to perform better than the linear nPWFS reconstructor, but this behavior is misleading, in the sense that the GS algorithm always provides wrapped measurements, which prevents the amplitude of the reconstructed phase from diverging for extremely noisy measurements.Another way to explain this fact is to say that the GS has an intrinsic regularization that biases the results in this analysis.When the wrapping step is added, the GS reconstructor builds more error in the presence of noise compared to all the other configurations we studied.
Figure 8 provides the same analysis of the noise propagation for the case of a turbulent phase screen (0.75 nm rms).The same trend is visible, but there is an offset for the linear nPWFS reconstructor that corresponds to the reconstruction due to nonlinearity error (which is absent for GS cases or the modulated PWFS).
To better understand the noise propagation behavior, it is useful to examine how the error propagates along the modes for one flux regime.An example is given Fig. 9 for the case of 400 photons per mode.The unwrapping step clearly drives the noise to drastically propagate on the low-order spatial frequencies.
To conclude this brief study of noise propagation, the GS reconstructor combined with an unwrapping algorithm performs poorly in the presence of noise because it propagates more A48, page 6 of 11 noise than the modulated PWFS itself.Therefore, the sensitivity advantage of the nPWFS is lost while using this reconstructor.However, this implementation is currently in its most basic form.Noise propagation could be mitigated in the reconstructor through two aspects: on the one hand, by using a GS algorithm that is more robust against noise (Levin & Bendory 2019), and on the other hand, by using an unwrapping algorithm with a better behavior with respect to noise (Estrada et al. 2011).As mentioned in the introduction, using a nPWFS brings more advantages than only the gain in sensitivity with respect to noise, hence this current implementation of the reconstructor is still interesting for uses at high signal-to-noise ratio (S/N; application examples are given in the conclusion).It is also important to note that this drawback for the GS algorithm was not raised in previous studies that proposed its implementation to reconstruct the curvature of wavefront sensor signals (Guyon 2010).

Broadband impact on the performance
It is crucial to take into account the potential influence of larger spectral bandpass measurements on the GS algorithm when reconstructing the EM field, as this technique assumes the coherence of the field.The nPWFS is achromatic in phase (Fauvarque 2017), meaning that it measures the phase of the incoming wavefront independently of the wavelength of the light (dispersion effects are neglected here).However, the amplitude of the wavefront does depend on the wavelength because aberrations present in the system and the turbulence usually introduce a fixed optical-path difference (OPD) in all wavelengths.Hence, the measured signal will scale proportionally with the wavelength.This leads to two advantageous properties: for a flat wavefront, all the wavelengths will give the same measurements.In the right conditions, the closed loop can therefore converge towards a null phase.Second, by choosing the central wavelength for the reconstruction, it allows scaling of the phase to compensate for larger and smaller wavelengths, which returns a monochromatic measurement in the case of the linear range (providing a flat spectrum).Overall, it is known that the bandwidth has a limited impact on nPWFS measurements in general, and therefore, it seems reasonable to expect the same for the GS reconstructor, even though it assumes monochromatic measurements.
A thorough investigation of the broadband impact on reconstruction would require an extensive study assessing various factors such as the amplitude of the measured phase, the sampling used for the measurements, and the bandwidth.Because the scope of this paper is limited, a brief example is given to illustrate the impact of a broadband measurement on the GS algorithm using the same configuration as in the previous section, but providing the algorithm measurements recorded by a polychromatic light with a flat spectrum ranging from 550 nm to 800 nm (bandwidth about 35%, sampled with 50 points in our simulation).To run the GS algorithm, the nPWFS model was simulated as working at a monochromatic wavelength, set as the broadband central wavelength (675 nm).Linearity curves for Zernike mode 19 in the case of different numbers of GS iterations are again displayed, but adding the one corresponding to the reconstruction of a broadband measurement (Fig. 10).Figure 10 demonstrates the limited impact of the broadband measurement on the reconstruction, even though the chosen bandwidth is large.
To confirm the GS algorithm robustness to broadband measurement, the reconstruction error as a function of iteration number for a turbulent screen in the case D/r 0 = 4 is plotted in Fig. 11.The reconstruction error stabilizes to a reconstruction error that is slightly larger than in the monochromatic case.
Overall, this short study points out that broadband measurements should not jeopardize the GS reconstruction scheme.We also mention that a polychromatic implementation of the GS algorithm could be imagined (Fienup 1999), but it will come at the expense of the computational time.
As the GS reconstruction scheme proposed in this paper is highly model dependent, it is important to show that our technique is robust enough to model errors, so that it can be implemented on a real experimental setup.This is shown in the next section, in which we provide an experimental demonstration of our reconstructor.

Experimental setup
The goal of this section is to present a laboratory demonstration of the GS algorithm on SEAL, an extreme adaptive optics testbed composed of several deformable mirrors (DM), wavefront sensors, and coronagraphic branches.A schematic layout of the SEAL testbed for the nPWFS subsystem alone is presented in Fig. 12.The SEAL components relevant for the experiment described here are the source at λ = 635 nm, spatial light modulator (SLM) with 1100 pixels across the pupil diameter (van Kooten et al. 2022), an IRISAO segmented deformable mirror (DM) has six segments across the pupil diameter, a low-order ALPAO DM has nine actuators across the pupil diameter and a A48, page 7 of 11 high-order BMC DM has 24 actuators across the pupil diameter.The focal-plane camera is sampled at 2.3 Shannon and the nPWFS is composed of a double-rooftop mask (Lozi et al. 2019) with 106 pixels across the pupil diameter on the camera (same sampling as the simulations reported in this paper).
To build a reliable model of the SEAL nPWFS and run the GS algorithm, two quantities are required: an image of the pupil, and the shape of the pyramid mask.To obtain an image of the pupil, large tip-tilt offsets were applied on the ALPAO DM in order to move the PSF away from the pyramid tip (∼20λ/D) and place it successively in each nPWFS quadrant.By doing so, four images of the pupil were obtained through each side of the pyramid mask.The SEAL pupil was then computed by recentering and averaging these four images, with an estimated accuracy below one pixel.An image of the SEAL pupil measured through this method is shown in the top left panel of Fig. 1.The segments gaps created by the IRISAO DM are clearly visible.In order to compute the pyramid shape, the four pupil images were simply registered in order to produce the corresponding tip-tilt for each face of the pyramid mask.The phase-to-DMs (ALPAO and BMC) registration was made the following way: the central actuator of the DM was pushed and reconstructed, and then again pulled and reconstructed.The difference of the images corresponds to the phase of this actuator.Then a waffle pattern was sent to the DM and was reconstructed in order to register the actuator positions.Finally, the phase of the central actuator was fitted with a Gaussian function and was duplicated at the other actuatos positions.This calibration process requires only four images (two for the central actuator, and two for the waffle).For the ALPAO and BMC DM, the reconstructed phase is largely oversampled (106 pixels across versus 23 actuators for the BMC).The nPWFS and associated GS algorithm are routinely used on the SEAL testbed with both DMs as the standard way to flatten the wavefront in order to correct for quasi-static aberrations, achieving a wavefront error of about 17 nm rms (as measured by the nPWFS).
The nPWFS image after the closed loop is shown in Fig. 13 (left).A high-frequency pattern caused by the BMC DM is clearly visible (a well-known effect, called the print-through or quilting effect).The corresponding residual phase was reconstructed from the reference image (best flat after the closed loop on static aberrations) and propagated in the model.The obtained image is displayed in Fig. 13 (right), showing that the model exhibits high fidelity with the measurements.The small differences that can be spotted include the print-through effect, which seems slightly underestimated in the simulated image most likely because the field of view of the simulated nPWFS is smaller than the real one, which acts like a spatial low-pass filter), and a faint ring pattern can be seen between the two bottom pupil images of the true image (most likely caused by a dust particle on the glass pyramid).It is worth mentioning that all the results presented in the following subsections were obtained for measurements made at high S/N.

Linearity curves
As a first demonstration of the GS algorithm performance on the SEAL testbed, some of the linearity curves obtained in the simulations in the previous section are reproduced.This study was made in the following way: Zernike modes were sent with the SLM to gain a good knowledge value of the input phase.Following the same procedure as before, the linear reconstructor was calibrated with the Zernike modal basis, hence the reconstruction directly gives the value of each of the modes.For the GS algorithm, the phase was reconstructed for each pixel and was then projected on the Zernike basis.
The linearity curves obtained for the Zernike mode 19 and 150 are shown in Fig. 14 for the 200 iterations we used for the GS algorithm.The conclusions from the simulations are confirmed in this experiment: the GS demonstrates a higher dynamic range than the linear reconstructor.However, for higher-order modes, the GS algorithm without phase unwrapping seems to under-perform compared to the simulation (Fig. 5).This could be explained by model errors that affect the performance slightly.

Example of phase reconstruction and close loop
As a demonstration of the GS capabilities in reconstructing strong phase aberrations on the bench, a D/r 0 = 32 phase screen was generated on the SLM (bottom left panel of Fig. 15).The corresponding recorded nPWFS image is shown in the top left panel of the same figure, and the reconstructed unwrapped phase after 200 iterations is displayed in the bottom right corner.Figure 15 also shows a simulated image of the nPWFS signal, produced by simply propagating the reconstructed phase in the nPWFS model.Although the simulated and the real image are almost identical, the phase is largely underestimated (by more than a factor of 2; this effect was also observed in the simulations presented in the previous section).This is explained by the fact that the nPWFS saturates: for an increasing input phase, the signal reaches a point at which it almost does no longer change.A simple example of this saturation is the case of a tip-tilt aberration: when the PSF is moved several tens of λ/D in one quadrant, only one pupil image of the nPWFS is illuminated.Displacing the PSF even farther away will almost have no impact on the A48, page 8 of 11  measurements because the illuminated pupil image already concentrates the entire flux.The saturation is a strong effect that will limit the reconstruction range for any type of linear but also nonlinear reconstructor because it implies that two different phases can lead to the same measurements.Hence, the saturation seems to be an intrinsic limitation in the nPWFS measurement, and inverting large-amplitude aberrations more accurately would require additional knowledge of the phase to be measured.A potential solution to even further improve the dynamic range of GS reconstruction is sketched in the next section.
Reconstructed phases at different iteration steps are displayed in Fig. 16.As shown before, the phase estimation improves with iterations, but the estimation appears to stop improving after just aa few hundred iterations instead of a few thousand iterations in the simulations.Hence, the reconstructed phases for 200 and 2000 iterations are highly similar.Once again, this could be explained by the differences between the nPWFS on the SEAL testbed and its model used for reconstruction.
It is also possible to compare the GS reconstruction with the linear reconstruction.To do this, a push-pull interaction matrix was measured on the bench by sending the first 500 Zernike modes on the SLM.The command matrix was then computed by taking the pseudo-inverse of the interaction matrix, and it was used to reconstruct the signal.The comparison between the linearly reconstructed phase and the GS reconstructed unwrapped phase projected on the first 500 Zernike modes is shown in Fig. 17.Our GS reconstruction clearly also shows a significant improvement compared to the linear reconstruction on the bench.To conclude the testbed demonstration, we present results of a closed-loop run using the BMC DM (controlling all the modes) on the static turbulent screen presented in Fig. 15.The linear reconstructor uses a push-pull zonal interaction matrix, and the GS reconstructor uses 25 iterations with phase unwrapping.For the closed-loop AO control, a simple integrator controller with a loop gain of 0.5 was used.The PSFs obtained after 12 closedloop steps are presented in the PSF corresponding to the best flat after closing the loop with the nPWFS on bench aberrations (with the BMC DM), and the uncorrected PSF corresponding to the atmospheric phase displayed on the SLM.The best flat PSF exhibits a faint and vertical light pattern that crosses its core.This comes from diffraction effects due to the SLM.The closed-loop test clearly confirms the extended dynamic range of our reconstructor by showing that the closed loop with the GS reconstructor outperforms the closed loop in the linear framework (which finally diverges after a few dozen steps).It also shows that despite saturation effects, the AO closed-loop scheme helps to correct the high-amplitude phase with the nPWFS combined with the GS algorithm.As the wavefront error drastically improved after bootstrap, we could imagine a closed-loop scheme for which the number of GS iterations used for the reconstruction is reduced in the residual-phase regime.
These experimental tests demonstrate the GS reconstructor performance on an optical bench.They show that it provides an improvement over the linear reconstructor.Hence, the f-GS could be a powerful tool for drastically increasing the convergence speed of the algorithm and performance in the case of high amplitude aberrations.These advantages come with a more complicated practical implementation, requiring splitting the flux (an operation that could introduce noncommon-path aberrations) between two synchronized cameras.The details of this implementation will be considered in future studies, with the main motivation of building a demonstration of the f-GS on the SEAL testbed.We note that f-GS may not be needed in the case of a classic AO closed-loop scheme because the nPWFS would work around residual phases and therefore far from the saturation regime after bootstrap is completed.Nevertheless, this f-GS setup could open the nPWFS to a wider range of applications that would require performing wavefront sensing on large phases.

Conclusion
We introduced a new way to invert the nPWFS measurements that relies on the numerical model of the sensor and on the use of the GS algorithm.We demonstrated that the GS-based reconstructor along with an unwrapping algorithm can drastically improve the nPWFS dynamic range of the reconstruction compared to the linear framework and can extend its linearity range by a factor of approximately 3.Moreover, it can achieve a dynamic performance comparable to the 3 λ/D modulated PWFS up to 2 rad rms.This technique was successfully demonstrated on the SEAL testbed at UCSC, on which it was used to close the loop on high D/r 0 turbulent phase screens.This reconstructor has two drawbacks, however: it is highly computationally complex, which prevents it from being used for real-time control purposes (in its current implementation), and the noise propagation is worse compared to the linear reconstructor.However, noise propagation for only the most basic implementation of GS and phase unwrapping was studied in this paper.There might be a better-suited GS algorithm combined with a noiserobust phase-unwrapping algorithm that could improve the noise propagation.The GS reconstructor already has several areas that can be useful in an AO system even though they were only tested at high S/N at slow speeds.For example, calibration purposes (as done routinely on the SEAL testbed), segment or fragment phasing, or highly sampled phase reconstruction.It could also support a second-stage nPWFS running in real time with a linear reconstructor through soft real-time reconstruction of the phase in order to compute optical gains, and it might also be used to reconstruct telemetry data with higher fidelity.
The most promising path toward a real-time implementation of this technique is to use focal-plane-assisted GS, which we showed in the simulation in Sect. 5.This technique uses the fact that the GS algorithm is well suited for sensor fusion and seems to significantly increase the convergence speed of the algorithm while allowing us to measure high-amplitude phases with better accuracy.The next steps for this approach is to push our understanding of the f-GS and start to implement it on the SEAL testbed, using a PSF image as an additional sensor.
Finally, the GS algorithm was applied in this paper for nPWFS alone, but we argue that it can be used for all Fourierfiltering WFS and maybe beyond.In a more general way, this reconstructor belongs to a wider class of reconstructors that use a numerical model of the sensor and iterative algorithms to increase the dynamic range.
A48, page 1 of 11 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe to Open model.Subscribe to A&A to support open access publication.
Principle of the GS algorithm applied to nonmodulated PWFS.The yellow dots highlight the fact that these are true data obtained on the SEAL testbed.

ChambouleyronFig. 2 .
Fig. 2. Convergence speed of the GS algorithm.The GS algorithm was applied for input turbulent screens for four different seeing conditions.The dashed lines represent the corresponding reconstruction error in the linear framework and the x-axis is given in log-scale.

Fig. 5 .
Fig. 5. Linearity curves for different Zernike modes, ranging from low to high order.The GS algorithm is run with 25 iterations in this case.

Fig. 6 .
Fig. 6.Linearity curves for the Zernike mode 19 for different numbers of iterations while performing the GS algorithm.The number of iterations is written in the top right part of the curves.

Fig. 7 .
Fig. 7. Dynamic range plots produced by sending 100 randomly sampled phase screens following a Von Karman PSD for each input amplitude and reconstructing them.

Fig. 8 .Fig. 9 .
Fig. 8. Reconstruction errors due to photon noise.Top: case of an input flat wavefront.Bottom: case of an input turbulent wavefront of 0.75 radians rms.

Fig. 12 .
Fig. 12. Simplified layout of the SEAL testbed for the nPWFS subsystem.

Fig. 13 .
Fig. 13.Flat wavefront on the SEAL nPWFS.Left: nPWFS image on SEAL after the closed loop.The nPWFS estimated wavefront error is 17 nm rms.Right: corresponding simulated image.The same scaling is used for both images.

Fig. 14 .
Fig. 14.Top: linearity curve for Z 19 obtained on SEAL.Bottom: linearity curve for Z 150 obtained on SEAL.The GS algorithm is used with 200 iterations.The yellow dots highlight the fact that these curves were obtained on the SEAL testbed.

NFig. 18 .
Fig. 18.Closing the loop with the BMC DM on a static atmospheric screen.Top left: PSF corresponding to best flat on the SEAL testbed.Top right: uncorrected PSF for an input phase D/r 0 = 32 on the SLM.Bottom left: PSF after 12 closed-loop steps using the GS reconstructor.Bottom right: PSF after 12 closed-loop steps using the linear reconstructor.

Fig. 21 .
Fig. 21.Reconstructed phases for the f-GS compared to the standard algorithm for one iteration and 4000 iterations.All phases are displayed with the same scaling.