Focal-plane-assisted pyramid wavefront sensor: Enabling frame-by-frame optical gain tracking

Aims. With its high sensitivity, the pyramid wavefront sensor (PyWFS) is becoming an advantageous sensor for astronomical adaptive optics (AO) systems. However, this sensor exhibits signiﬁcant non-linear behaviours leading to challenging AO control issues. Methods. In order to mitigate these e ﬀ ects, we propose to use in addition to the classical pyramid sensor a focal plane image combined with a convolutive description of the sensor to fast track the PyWFS non-linearities, the so-called optical gains (OG). Results. We show that this additional focal plane imaging path only requires a small fraction of the total ﬂux while representing a robust solution to estimating the PyWFS OG. Finally, we demonstrate the gain that our method brings with speciﬁc examples of bootstrapping and handling non-common path aberrations.


Introduction
The pyramid wavefront sensor (PyWFS), which was proposed for the first time in 1996 by Ragazzoni (1996), is an optical device used to perform wavefront sensing. Inspired by the Foucault knife test, the PyWFS is a pupil plane wavefront sensor performing optical Fourier filtering with a glass pyramid with four sides that is located at the focal plane. The purpose of this glass pyramid is to split the electromagnetic (EM) field into four beams producing four different filtered images of the entrance pupil. This filtering operation allows the conversion of phase information at the entrance pupil into amplitude at the pupil plane, where a quadratic sensor is used to record the signal (Vérinaud 2004, Guyon 2005. Recently, the PyWFS has gained the interest of the astronomical community because it offers a higher sensitivity than the classical Shack-Hartmann wave-front sensor (WFS) that is commonly used in adaptive optics (AO) systems (Esposito & Riccardi 2001). However, the PyWFS exhibits non-linearities that prevents a simple relation between the incoming phase and the measurements, leading to control issues in the AO loop. Previous studies (Korkiakoski et al. 2008, Deo et al. 2019a) have demonstrated that one of the most striking effect of this undesirable behaviour is a time-averaged frequency-dependent loss of sensitivity when the PyWFS is working in presence of non-zero phase. This detrimental effect can be mitigated by providing an estimation of the so-called optical gains (OG), which are a set of scalar values encoding the loss of sensitivity with respect to each component of the modal basis. The goal of this paper is to present a novel way of measuring the OG. In the first section we introduce the concept of the linear parameter-varying system (LPVS) to describe the PyWFS, which opens the possibility of estimating the OG frame by frame instead of considering a time-averaged quantity. In the second section, we present a practical implementation of the method, enabling frame-by-frame OG tracking. Finally, we illustrate this OG tracking strategy in the context of closed-loop bootstrapping and handling non-common path aberrations (NCPAs).

PyWFS non-linear behaviour and optical gains
In the following, we call s the output of the PyWFS. This output can be defined in different ways. The two main definitions are called 'full frame' or 'slope maps'. In the first case, s is obtained by recording the full image of the WFS camera, for which a reference image corresponding to a reference phase has been removed (Fauvarque et al. 2016). In the second case, the WFS image is processed to reduce the useful information to two pupil maps usually called 'slope maps' (Ragazzoni 1996). The work presented here remains valid for full-frame or slope-map computation, and we decided to use the full-frame definition throughout.
When described with a linear model, the PyWFS outputs are linked with the incoming phase φ through an interaction matrix called M. This interaction matrix can be built with a calibration process that consists of sending a set of phase maps (usually a modal or a zonal basis) to the WFS with the deformable mirror (DM) and then record the derivative δs(φ i ) of the PyWFS response for each component of the basis (Fig.1). This operation is most commonly performed with the so-called push-pull method, consisting of successively sending each mode with a positive and then negative amplitude a to compute the slope of The interaction matrix (also called Jacobian matrix) is then the collection of the slopes recorded for all modes, M = (δs(φ 1 ), . . . , δs(φ i ), . . . , δs(φ N )).
( 2) In this linear framework, we can then link the measured phase with the output of the PyWFS by the relation This matrix computation formalism has interesting properties that are required in the AO control loop. However, the PyWFS exhibits substantial non-linearities that make the equation above only partially true. Mathematically, the deviation from linearity is expressed with the following inequality: where φ is a non-null given phase. When working around φ, the slope of the linear response of the sensor is therefore modified, During AO observation, the sensor is working around a nonnull phase φ corresponding to the residual phase of the system. As a consequence of Eq. (4), the response of the system is modified. Previous studies suggested updating the response slopes to mitigate this effect by relying on two main concepts. The first concept is the stationarity of the residual phases (Rigaut et al. 1998). For a given system and fixed parameters (seeing, noise, etc.), we can compute an averaged response slope for each mode. It has been proven (Fauvarque et al. 2019) that under this stationarity hypothesis, the averaged response slope depends on the behaviour of the statistical residual phases through their structure function (D φ ), δs φ(φ i ) = δs D φ (φ i ).The second concept is the diagonal approximation (Korkiakoski et al. 2008). This approximation implies considering no cross-talk between the modes, which means that the response slopes are only modified by a scalar value for each mode. This value is known as the OG. We then have δs D φ (φ i ) = t i D φ .δs(φ i ), where t i D φ is the OG associated with the mode i for a given residual phase perturbation statistics characterised by the structure function D φ . In this approximation, the shape of the response is left unchanged.
Finally, the interaction matrix is updated by multiplying by a diagonal matrix T D φ called the OG matrix, whose diagonal com- We used the scalar product presented in Chambouleyron et al. (2020) to calculate the diagonal components of this matrix, .
Several approaches to practically compute this matrix can be found in the literature. They can be split into two categories: those that are invasive for the science path, consisting of sending some probe modes to the DM to return to the OG (Esposito et al. 2015, Deo et al. 2019a, and those that rely on the knowledge of the statistics of the residual phases through the telemetry data to estimate the OG (Chambouleyron et al. 2020). In all the proposed methods, the OG can be seen as an evaluation of a timeaveraged loss of sensitivity of the sensor. Being able to accurately retrieve OG allows compensating for the sensitivity loss.

LPVS approach
As described by Eq. (4), the PyWFS outputs are affected by the incoming phase. The time-averaged definition of the interaction matrix M D φ , is limited to a statistical behaviour of the PyWFS, even though it has good properties. We propose a framework that addresses the non-linearities in real time, with an interaction matrix that is updated at every frame. To do so, we first assumed that the diagonal hypothesis holds. Then, and inspired by the automatic field domain, the PyWFS is now considered as an LPVS (Rugh & Shamma 2000): Its linear behaviour encoded by the interaction matrix is modified at each frame according to the incoming phase. Under this framework, the new expression of the PyWFS output can be written as where T φ is the OG matrix for the given measured phase φ.
Assuming the diagonal approximation holds, we can extract T φ from the interaction matrix computed around φ, For a given system, repeating this operation on a set of different phases will eventually lead to the time-averaged definition of the OG matrix, To illustrate the difference between the time-averaged response and a single realisation, we performed the simulation presented in Fig. 2. These simulations were made with parameters consistent with an 8m telescope and for two seeing conditions. All results showed in this paper rely on end-toend simulations performed with the OOMAO Matlab toolbox (Conan & Correia 2014). The exact conditions and parameters are summarised in Table 1. In the simulation, we can compute the exact PyWFS response by freezing the entrance phase and performing a calibration process around this working point. We therefore computed T φ for 1000 residual phase realisations, and show the OG variability for two seeing conditions in Fig. 2. This represents an optimistic context where the Fried parameter r0 is fixed through the complete simulation. By estimating the T φ with a time-averaging strategy, the errors on the OG corresponding to a given residual phase can reach more than some dozen percent (OG exhibiting a maximum deviation from the averaged value are highlighted. This result illustrates the potential gain Variability of closed-loop OG. For given system parameters we compute T φ for 1 000 phase realisations in two seeing configurations: r0 = 18 cm and r0 = 12 cm. The variability of the frame-by-frame OG is shown in the histogram in the right panel and by the highlighted extreme OG curves for each r0 case. of performing a frame-by-frame estimation of the OG instead of a time-averaged one. In the next section, a practical means for performing this frame-by-frame gain-scheduling operation is presented.

Principle
Obtaining an estimate of the OG values (the diagonal of T φ ) requires obtaining additional information describing the working point of the PyWFS at each moment, independently of the PyWFS measurements themselves. To this end, a specific sensor called gain-scheduling camera (GSC) is implemented. Empirically, it is well known that the PyWFS sensitivity depends on the structure of the EM field when it reaches the pyramid mask. For instance, the more this field is spread over the pyramid mask, the less sensitive the PyWFS. In addition, because sensitivity and dynamic range are opposing properties, a well-known technique used to increase the PyWFS dynamic range consists of modulating the EM field around the pyramid apex. In order to keep track of the sensor regime, we therefore suggest probing this EM field by acquiring a focal plane image synchronously with the Pyramid WFS data. This can be achieved by placing a beam splitter before the pyramid mask and recording the signal with a focal plane camera that has the same field of view as the pyramid (Fig. 3).
In this configuration, the focal plane camera, hereafter called the GSC, records the intensity of the modulated EM field seen by the pyramid. By using the same exposure time and frame rate as the WFS camera, the signal observed is then an instantaneous AO-corrected point-spread function (PSF) convolved with  the circle of modulation. This is illustrated in Fig. 4, where the modulation circle is shown on the left, and the replicas of this modulation circle by the focal plane speckles are shown on the right. By denoting Ω φ the GSC signal, we can therefore write where ω is the modulation weighting function. This function can be thought of as a map of the incoherent positions reached by the EM field on the pyramid during one integration time of the WFS camera. This function is thus a circle for the circularly modulated PyWFS (Fig. 5 right). Ω φ has to be understood as the effective modulation weighting function: The phase to be measured produces its own modulation, leading to PyWFS loss of sensitivity, and the GSC is therefore a way to monitor this additional modulation. The next step is now to link this focal plane information with the PyWFS optical gains and merge the GSC and PyWFS signal in one final set of WFS outputs. In a previous work (Chambouleyron et al. 2020), we demonstrated that the convolutive model of the PyWFS developed by Fauvarque et al. (2019) can be used to predict the averaged OG if the statistical behaviour of the residual phases (through the knowledge of their structure function) is known. In Eq. (11) we recall the expression of the PyWFS output in this convolutive framework, where IR is the impulse response of the sensor and the star denotes the convolutive product. In the framework of the infinite pupil approximation, the impulse response around a flat wavefront can be expressed through two quantities, the mask complex function m and the modulation function ω (Fig. 5), We propose here to combine this model with the signal delivered by the GSC in order to compute the impulse response IR φ of the PyWFS around each individual phase realisation. To do this, we replaced ω by the GSC data as described in Eq. (13), This new way to compute the impulse response can be considered as using the impulse response given for an infinite pupil system (Eq. 12) for which we replaced the modulation weighting function by the energy distribution at the focal plane, including both the modulation and the residual phase. Now that we are able to compute IR φ at each frame, we can estimate the OG matrix T φ through the following computation of its diagonal components as described in Chambouleyron et al. (2020), where IR calib is the impulse response computed for the calibration state, most commonly for φ = 0 (Fig. 4, left).

Accuracy of the estimation
It is now possible to test the accuracy of our estimator by comparing T φ and T φ . To do this, we computed the true T φ through end-to-end simulations by proceeding through the ideal way described in section above: An interaction matrix was computed around each given residual phase, from which the OG matrix was derived (Eq. 8). This provides the ground truth to which the gains estimated with the GSC are compared. First results are shown in Fig. 6 for different seeing and modulation conditions. As illustrated in Fig. 6, the real and estimated OG agree well, demonstrating the accuracy of the proposed method.
For the parameters used in our simulations, the estimation remains accurate regardless of whether we are in open loop or closed loop. The ripples seen in the ground-truth OG curves are smoothed in the convolutive framework. The convolutive product given in Eq. (11) tends to smooth the output of the PyWFS even when the impulse response is computed around a nonzero phase. Figure 6 also shows a slight deviation for low-order modes for a low-modulation regime and a strong entrance phase (open loop here).

Robustness to noise
The GSC has shown to be a reliable way to perform a fast OG tracking, but it requires using a fraction of the photons available in the sensing path. This inevitably competes with the gain of sensitivity provided by the PyWFS. The goal of this section is then to demonstrate that our GSC approach is only weakly affected by photon noise and therefore requires only a small number of photons while performing an accurate frameby-frame OG estimation. To this end, we propose to inject noise in the data delivered by the GSC and to probe the effect on the OG estimation.
We ran simulations with the same parameters as described above. The sensing path works around the central wavelength λ c = 550 nm with the given bandwidth ∆λ = 90 nm and an ideal transmission of 100%. The exposure time of the GSC is 2 milliseconds (frame rate of the loop), and 10% of the photons are used by the GSC camera. The GSC pixel size corresponds to Shannon sampling of the diffraction-limited PSF. In this given configuration, the data recorded by the GSC for a given closedloop residual phase (r0 = 14 cm, r mod = 3 λ/D) are presented in Fig. 7 (top) for (a.) a noise-free system, (b.) a guide-star magnitude equal to 8 (c.) a guide-star magnitude equal to 10, and (d.) a guide-star magnitude equal to 12. For these three noise configurations (mag = 8, 10, and 12), we estimated the OG for 500 realisations of the noise. The results are given Fig. 7 (bottom  part). The introduction of noise leads to an increased OG estimation error, which logically scales with the signal-to-noise ratio (S/N) according to √ n ph . However, the GSC approach also still performs a satisfactory OG estimation even for low-magnitude guide stars. For even fainter guide stars, the noise effect might be mitigated by integrating the GSC data over several frames. A trade-off between noise propagation and OG error would then be required. These results are crucial because they demonstrate that the GSC can be used with only a small fraction of WFS photons, leading to a limited repercussion on the S/N on the PyWFS. We therefore have a way to estimate the OG, and to some extent increase the linearity of the sensor while having a reduced effect on its sensitivity.

GSC spatial sampling
Another aspect is the sampling of the GSC detector with respect to the modulated PSF. If an under-sampling could be In the chosen configuration, the exposure time is 2 ms, and we collect 10% of the photons in the sensing path. a: infinite number of photons. b: Guide-star magnitude = 8 (n ph = 55 000 on the GSC). c: Guide-star magnitude = 10 (n ph = 9000 on the GSC). d: Guide-star magnitude = 12 (n ph = 1400 on the GSC). Bottom: OG estimate for the noise-free system compared with the three noisy configurations. considered, it would reduce the number of pixels required by the GSC, and consequently reduce the practical implementation complexity. To test this, we ran our algorithm for various samplings of the GSC in order to see the effect on the OG estimation. The results for a given closed-loop residual phase (r0 = 14 cm, r mod = 3) are given in Fig. 8. The sampling of the PSF can go below the Shannon sampling (2 pixels per λ/D) without significant effect on the estimate. This result depends on the modulation radius r mod used, and we note that the OG estimate is not affected as long as the pixel size d px satisfies the Shannon criterion for the modulation radius, When this criterion is not respected, the undersampled modulation circle is seen as a disc (Fig. 8), which affects the OG estimate for low-order modes. As a concrete example, a PyWFS for the Extremely Large Telescope (ELT) working at λ = 800 nm with a field of view of 2 arcsecs and with a sampling of Shannon/4 on the GSC would require a GSC camera with no more than 250 × 250 px. This limited size allows for the use of low-readout noise cameras such as GSC, and remaining in a photon-noise limited regime.
To conclude this section, we have shown that it is possible to perform OG fast-tracking by using an image of the modulated EM field at the focal plane. Our method uses a so-called GSC providing non-biased information on the working point of the PyWFS, and the subsequent OG estimate using a convolutive model. We demonstrated that the GSC can work with a limited number of photons and pixels, which makes the practical implementation fully feasible. The next section is dedicated to quantifying the performance benefits of OG fast tracking with the GSC.

Application to specific AO control issues: Bootstrapping and NCPA handling
As shown in the previous sections, the GSC allows tracking the PyWFS OG frame by frame and compensating for these nonlinearities. We illustrate here two possible situations in which the GSC can significantly improve the performance: bootstrapping and NCPA handling.

Bootstrapping
During the AO loop bootstrap, the PyWFS faces large amplitude wavefronts (due to uncorrected turbulence), leading to significant non-linearities that may prevent the loop from closing. Therefore this step is critical because it corresponds to the moment at which the OGs are the most important. Monitoring them frame by frame in order to update the reconstructor helps closing the AO loop. Because of the timescales involved in the AO loop bootstrap, this problem cannot be tackled by other OG handling techniques that were previously studied in the literature. The best solutions already proposed endures necessarily delays of a few frames (Deo et al. 2019b). Here, we can estimate the OG corresponding to the current measurement frame: This is an unprecedented feature. We show different images delivered by the GSC during the bootstrap operation in Fig. 9. The corresponding estimated OGs are also plotted, compared with the end-to-end computation giving the true OG values. While the loop is closing, the OG varies from low values to higher values, indicating that the residual phases reaching the PyWFS decrease: The loop is closing, and the DM is starting to correct the atmospheric aberrations. Our technique performs a precise OG follow-up during all the steps of the process, at the frame rate of the loop. We can use our frame by frame OG estimation to update the reconstructor while the loop is closing. The reconstructor is the pseudo-inverse of the interaction matrix. We can therefore relate it to the OG matrix and the calibration interaction matrix through the following formula: By doing so, we show that it is possible to close the loop faster. A simulation example is presented by comparing a loop bootstrap with and without OG compensation by the GSC camera (Fig. 10). This example, with a limited benefit in practice, shows how a fast OG tracking combined with the corresponding update of the reconstructor can be applied to mitigate all types of short-timescale residual variations, such as seeing bursts.

NCPA handling
Handling NCPA is emerging as one of the main issues due to PyWFS OG, as was demonstrated for instance on the Large Binocular Telescope (Esposito et al. 2015). How to handle this issue while having an accurate OG estimation was discussed in a previous paper (Chambouleyron et al. 2020). We briefly recall the main problem. The NCPA reference measurements are recorded around a diffraction-limited PSF and need to be rescaled by the OG while working on sky: s(φ NCPA ) ← s φ (φ NCPA ). To compute s φ (φ NCPA ), we need to have estimate T φ , We show here the results of a simulation in which we used the GSC to handle NCPA in the AO loop. We retained the same simulations parameters as before (caption of Fig. 2). The PyWFS modulation radius was r mod = 3 λ/D and r0 = 14 cm. The interaction matrix was computed around a flat wavefront. We injected 200 nm rms of NCPA into our system, distributed with a f −2 power law on the first 25 KL modes (except for tip-tilt and focus). In this configuration and for a flat wavefront in the science path (H band), the PSF in the wavefront sensing path (V band) is given in Fig. 11a and the signal Ω φ NCPA seen by the GSC is shown in Fig. 11b.
We then proceeded in the following way: We closed the loop on the turbulence, and after 5 s of closed-loop operation, the NCPA was added to the system. These NCPAs were then handled with different configurations, and the results were compared with the NCPA-free case. Figure 12 illustrates the results. The main conclusions from Fig. 12 are listed below. 1. When the NCPA is not compensated for (orange plot), the loop converges toward a flat wavefront in the sensing path. This induces a high loss of Strehl ratio (SR) in the science path, corresponding to the NCPA. 2. When a reference map s(φ NCPA ) in the PyWFS measurement is used without updating it by the OG, this leads to a divergence of the loop (so-called NCPA catastrophe, yellow plot). This can be explained by the fact that because of the OG, the PyWFS introduces too much NCPA, creating an even stronger aberrated wavefront. This aberrated wavefront increases the OG in the next frame, which continues to increase the aberration, and so on. This quickly causes the loop to diverge. 3. When the reference map is compensated for by the timeaveraged OG computed in the first 5 s of the loop by a Fig. 11. a: PSF on the pyramid apex when a flat wavefront is set in the science path. b: GSC signal when there are no residual phases and for a flat wavefront in the science path. c: GSC signal during closed-loop around NCPA. long-exposure image of the GSC (purple plot), no NCPA catastrophe appears, and the final performance reaches an averaged SR of 82%. 4. When the reference map is compensated for by the OG computed at each frame, using the GSC camera (green plot), the final performance reaches an averaged SR of 86%. This solution is better than the previous one because we monitor the OG at each frame, and we also take the effect of the NCPA themselves on the OG into account. To illustrate this, we show the GSC image for a given closed-loop residual when the NCPA is compensated for in Fig. 11c. This study is a clear demonstration that our strategy can solve the AO control issue due to PyWFS OG. It also shows that even if the OGs are compensated for on a frame-by-frame basis, the ultimate performance (without NCPA) cannot be reached. This limitation is mainly due to the LPVS approach, which is characterized by a linear description of the whole sensing problem. Improving the performance further would probably mean starting to consider other non-linear (second-or third-order description) solutions, which goes beyond the computation framework of a simple matrix.

Conclusion
The PyWFS is a complex optical device exhibiting strong nonlinearities. One way to deal with this behaviour while keeping a matrix computation formalism is to consider the PyWFS as a LPVS. To probe the sensing regime of this system at each measurement, a gain scheduling loop needs to be implemented that gives information on the sensor regime at every moment. With this perspective, the OG compensation can be deployed on a frame-by-frame basis. We provided here an innovative solution to this end: the GSC combined with a convolutive model. As such, the PyWFS data synchronously merged on a frameby-frame basis with GSC data can be thought of as a single WFS combining images from different light-propagation planes. It therefore provides an efficient way to compensate for nonlinearities at each AO loop frame without any delay, and it significantly improves the final performance of the AO loop in terms of sensitivity and dynamic range as well as robustness. It also allows unambiguously disentangling the effect of OG from the full AO loop gain, which is a fundamental advantage for NCPA compensation. The GSC solution has now to be implemented on the AO facility bench LOOPS at the LAM for an experimental demonstration (Janin-Potiron et al. 2019).