Real-time Adaptive Optics with pyramid wavefront sensors: A theoretical analysis of the pyramid sensor model

We consider the mathematical background of the wavefront sensor type that is widely used in Adaptive Optics systems for astronomy, microscopy, and ophthalmology. The theoretical analysis of the pyramid sensor forward operators presented in this paper is aimed at a subsequent development of fast and stable algorithms for wavefront reconstruction from data of this sensor type. In our analysis we allow the sensor to be utilized in both the modulated and non-modulated fashion. We derive detailed mathematical models for the pyramid sensor and the physically simpler roof wavefront sensor as well as their various approximations. Additionally, we calculate adjoint operators which build preliminaries for the application of several iterative mathematical approaches for solving inverse problems such as gradient based algorithms, Landweber iteration or Kaczmarz methods.


Introduction
Ground-based telescope facilities suffer from degraded image quality caused by atmospheric turbulence. When light from a distant star passes the Earth's atmosphere, initially planar wavefronts get distorted due to turbulent air motions causing fluctuations of the index of refraction. Therefore, advanced Adaptive Optics (AO) systems [35,64] are incorporated in innovative telescope systems to mechanically correct in real-time for the distortions with deformable mirrors. The shape of the deformable mirrors is determined by measuring wavefronts coming from either bright astronomical stars or artificially produced laser beacons. The basic idea is to reflect the distorted wavefronts on a mirror that is shaped appropriately such that the corrected wavefronts allow for high image quality when observed by the science camera (see Figure 1). The according positioning of the mirror actuators implies the knowledge of the incoming wavefronts. Thus, in Adaptive Optics one is interested in the reconstruction of the unknown incoming wavefront Φ from available data in order to calculate the optimal shape of the deformable mirror. Unfortunately, there exists no optical device which is able to measure the wavefront directly. Instead, a wavefront sensor (WFS) measures the time-averaged characteristic of the captured light that is related to the incoming phase. The wavefront sensor splits the telescope aperture into many small, equally spaced subapertures and detects the intensity of the incoming light in each of the subapertures. The number of subapertures defines the spatial resolution of the AO system. The detector of the sensor provides the intensity data by integrating the light from the celestial object. From the detector's analog signal, digital sensor measurements are sampled with a time period T which is usually extremely short, i.e., 0.3 -2 milliseconds. This paper is focused on the pyramid wavefront sensor (PWFS), which has been invented in the 1990s [60]. It is gradually gaining more and more attention from the scientific community, especially in astronomical AO, due to its increased sensitivity, improved signal-to-noise ratio, robustness to spatial aliasing and adjustable spatial sampling compared to the other popular wavefront sensor choice -the Shack-Hartmann (SH) sensor. Several theoretical studies [8,21,22,48,61,62,70,80,82,81,85] including numerical simulations and laboratory investigations with optical test benches [4,33,52,59,78,84] have Figure 1: Principle of wavefront correction [2]. A deformable mirror reflects a perturbed wavefront and propagates the corrected, planar wavefront to the science camera yielding improved image quality.
confirmed the advantages of pyramid wavefront sensors while additionally promising surveys were operated on sky [24,23,25,27,32,56,57,58]. The current development of a new era of Extremely Large Telescopes (ELTs) with primary mirrors of 25 − 40 m in diameter brings new challenges to the field of Adaptive Optics. For ELTs, pyramid sensors show enhanced performance in real life settings, e.g., they provide the ability to sense differential piston modes induced by diffraction effects of realistic telescope spiders that support secondary mirrors and perform even under significant levels of non-common path aberrations [79,20]. Pyramid wavefront sensors are going to be included in many ELT instruments [5,12,79,20,31,30,45,51,53,54,80]. We would like to mention that, apart from astronomical applications, the pyramid wavefront sensor is also applied in adaptive loops in ophthalmology [10,13,15,43] and microscopy [42,44]. Therefore, wavefront reconstruction algorithms for pyramid wavefront sensors are in high demand. The main goal of this paper is to provide an extensive mathematical analysis of the PWFS operators in order to develop suitable wavefront reconstruction methods. So far, the existing algorithms are (with a small number of exceptions as, e.g., [11,14,28,29,37,46,47,48,83]) based on a linear assumption of the pyramid sensor model [36,41]. Nevertheless, the simplifications of the non-linear pyramid operator allow for acceptable wavefront reconstruction quality.
Basically, wavefront reconstruction from pyramid sensor data consists in solving two non-linear integral equations P Φ = [s x , s y ] with respect to the unknown wavefront Φ, where P = [P x , P y ] is a singular Volterra integral operator of the first kind.
In the following, we derive and analyze the mathematical model for both the non-and modulated pyramid sensor. The paper is organized as follows: In Section 2, we give a brief introduction into the physical background of pyramid and roof wavefront sensors. Afterwards, we derive the singular and non-linear forward PWFS operator P and roof sensor operators R. In order to simplify the problem of solving the WFS equation, we calculate linearizations R lin for roof wavefront sensors in Section 3. Furthermore, in Section 4, we evaluate the corresponding adjoint operators which are necessary for the application of linear iterative methods such as gradient based algorithms or Landweber iteration for wavefront reconstruction addressed in an upcoming second part of the paper.

Pyramid and roof wavefront sensor
As it can be seen in Figure 2, the main component of the pyramid sensor is a four-sided glass pyramidal prism placed in the focal plane of the telescope pupil. The incoming light is focused by the telescope onto the prism apex. The four facets of the pyramid split the incoming light in four beams, propagated in slightly different directions. A relay lens, placed behind the prism, re-images the four beams, allowing adjustable sampling of the four different images I ij , i, j = {1, 2} of the aperture on the CCD camera. The two measurement sets s x , s y are obtained from the four intensity patterns as where I 0 is the average intensity per subaperture. For the sake of simplicity, we focus on the transmission mask modeling approach, which ignores the phase shifts introduced by the pyramidal prism and models the prism facets as transmitting only. The interference effects are neglected as well.
A dynamic circular modulation of the incoming beam allows to increase the linear range of the pyramid sensor [81] and is also used to adjust its sensitivity. The modulation can be accomplished in several ways: either by oscillating the pyramid itself [60], with a steering mirror [8,49], or by using a static diffusive optical element [49]. The circular modulation path of the focused beam on the pyramid apex is shown with a dashed circle in Figure 2. We consider the roof wavefront sensor as a stand-alone WFS and as a simplification of the pyramid sensor. For this type of sensor, two orthogonally placed two-sided roof prisms instead of the pyramidal one are used. The orthogonality of the roof leads to a decoupling of the dependence of the sensor data on the incoming phase in x-and y-direction. For roof wavefront sensors, a linear modulation path is additionally considered in the literature as an approximation to the circular one. The understanding of the physics behind the pyramid sensor has changed over time. In the beginning, the PWFS was introduced with dynamic modulation of the incoming beam and described within the geometric optics framework as a slope sensor similar to SH sensors, but having a higher sensitivity [19,60,62,82]. Later, the role of the beam modulation was questioned and the pyramid sensor without modulation was studied as well [47,50,61]. According to the Fourier optics based analytical model derived in [48], the non-modulated PWFS measures a non-linear combination of one and two dimensional Hilbert transforms of the sine and cosine of the incoming distorted wavefront. Then, it was recognized that the dynamic modulation of the beam allows to strengthen the linearity of the sensor and to increase its dynamic range. Taking modulation into account within the Fourier optics framework complicates the non-linear forward model even more. However, a linearization of both models, with and without modulation, is possible under certain simplifying assumptions [8,81]. Within the linearized model, it was shown that the modulated pyramid sensor measures both the slope and the Hilbert transform of the wavefront, depending on the frequency range and the amount of modulation [81].
The full continuous measurements s x (x, y) and s y (x, y) of the pyramid wavefront sensor are not available in practice. For the description of the discrete pyramid sensor we perform a division of the continuous two dimensional process into finitely many equispaced regions called subapertures. The data are then assumed to be averaged over every subaperture which corresponds to the finite sampling of the pyramid sensor. Hence, the data grid is predetermined by a subaperture size of d · d with d = D n , where D represents the telescope diameter, i.e., the primary mirror size, and n the number of subapertures in one direction. Additionally, we restrict the availability of measurements to the size of the region captured by the sensor. For several telescope systems the pupil is annular instead of circular since a secondary mirror shades the primary mirror, making the area of central obstruction hardly attainable for photons. Thus, the remaining light in the area of the central obstruction does not produce reliable measurements. Moreover, the incoming phases are defined on R 2 but for the control of the deformable mirror, we are only interested in the reconstructed wavefront shape on a restricted domain (bounded by the size of the telescope pupil).
In the following, we describe the annular telescope aperture mask by Ω = Ω y × Ω x ⊆ [−D/2, D/2] 2 as shown in Figure 3. Single lines of the annular aperture are represented by Ω x = [a x , b x ] and Ω y = [a y , b y ], with a x < b x , a y < b y being the borders of the pupil for fixed x and y correspondingly. The sensor provides measurement on the region of the CCD-detector D. Throughout the paper, we do not distinguish between pupil and CCD-detector and assume D = Ω. Further, the limitation onto the CCD detector (indicated as multiplication with a characteristic function of D) is not marked explicitly for the underlying operators for simplicity of notation. However, please keep in mind the compact support of the considered functions because of the restricted size of the aperture and the CCD-detector. For instance, we consider the norms in L 2 R 2 but since the operators map from and to functions with compact support on Ω it is equivalent to considering the norms in L 2 (Ω).

Pyramid and roof wavefront sensor forward operators
The pyramid WFS sensor model is non-linear and extensive [48,69]. The measurements (s n x , s n y ) of the non-modulated sensor are connected to the wavefront Φ via a combination of 1d and 2d Hilbert transforms of sin Φ, cos Φ and their multiplications. Circular modulation complicates the model even more which makes the task of inverting the operator P , i.e., the full Fourier optics based model of the sensor rather difficult.
In the following P {n,c} = [P {n,c} Atmospheric turbulence is usually introduced by either using the von Karman or Kolmogorov turbulence model which describe the energy distribution of turbulent motions [64,65]. Following the Kolmogorov model, considerations in [16] initiate to expect smoother functions for representing atmospheric turbulence than just L 2 with a high probability. This leads to the assumption of Φ ∈ H 11/6 R 2 as already suggested in [18].
We start with the derivation of the analytic pyramid wavefront sensor model without taking interference effects between the four pupils on the CCD-detector into account. This pyramid sensor model is known as the transmission mask model. For the modulated sensor, we define a modulation parameter with α = rλ/D for a positive integer r representing the modulation radius and λ the sensing wavelength.

Pyramid forward model without interference (transmission mask model)
For the pyramid wavefront sensor, we only consider the non-modulated and the sensor with circular modulation since they make sense from the physical point of view.

Definition 1. We introduce the operators P {n,c}
x in x-direction given by and P {n,c} y in y-direction given by Ωy Ωx Ωy cos(x cos t) dt with the modulation parameter α λ defined in (2).
Note that the kernels of the involved integral operators are strongly singular. Hence, they are defined in the p.v. (principal value) meaning, i.e., the integrals above are meant in the sense of for any wavefront Φ. In the following, p.v.
Ωy Ωx Ωy is always meant as abbreviation of p.v.
Ωy in the context of the pyramid sensor operator.

Theorem 1.
The relation between the pyramid wavefront sensor data and the incoming wavefront is given by where P {n,c} denote the pyramid sensor operators defined in (3) -(4).
Proof. The proof is given in the Appendix.

Roof sensor forward models (transmission mask model)
The roof WFS constitutes a part of the pyramid WFS. In the roof sensor, the pyramidal prism is replaced by two orthogonally placed two-sided roof prisms, resulting in a decoupling of x-and y-direction.  [48,68,81]. Due to the physical setup of the roof WFS, linear modulation induces characteristics which are of interest especially for roof wavefront sensors. Hence, we additionally investigate the linear modulated roof wavefront sensor model. In case of circular modulation the amplitude of modulation is assumed to be α λ , already introduced in (2). For linear modulation, 2α denotes the angle of ray displacement along the desired direction, which is equivalent to the assumed circular modulation. We substitute the zero-order Bessel function by a sinc-term in order to describe a linear modulation instead of a circular one. For further investigations, we again do not take interference between the four beams into account and consider the roof WFS operator R {n,c,l} on the one hand as a standalone wavefront sensor and on the other hand as an approximation to the pyramid WFS P {n,c} .

Definition 2. We introduce the operators R {n,c,l}
x and R {n,c,l}

Theorem 2.
Using the operators defined in (6) -(7) the measurements of the roof WFS in the transmission mask model are written as Proof. See [8,70,81].
The operators P {n,c}  Proof. From the proof of Theorem 1 it follows that the pyramid sensor operators and (further) the roof sensor operators are non-linear and well-defined for any wavefront Φ ∈ H 11/6 R 2 . It remains to show R {n,c,l} Φ ∈ L 2 R 2 . The proof uses the boundedness (1) |sin (Φ)| ≤ |Φ| and the Hölder continuity (2) with α = 5/6 and Hölder constant C > 0 in one direction of any function Φ ∈ H 11/6 R 2 (cf Sobolev embedding theorem, e.g., [1,Theorem 5.4]). We start with showing that the integrand of (6) in fact is As the proper integral exists, the Cauchy principal value exists as well. It follows that (8) and further sup Hence, the roof sensor operators maps incoming wavefronts in H 11/6 R 2 to measurements in L 2 R 2 with compact support on the telescope pupil/detector.
From the considerations in the above proof it follows the following statement: Proposition 2. The non-linear operators P {n,c} : H 11/6 R 2 → L 2 R 2 , representing pyramid wavefront sensors, are well-defined operators between the above given spaces.
Proof. As already shown in the proof of Theorem 1, the pyramid sensor operators are non-linear, welldefined operators for any wavefront Φ ∈ H 11/6 R 2 . In order to verify that P {n,c} Φ ∈ L 2 we split the corresponding operators into two parts: with the roof sensor operators R {n,c} x defined in (6) and the second term With Proposition 1 and First, we focus on the non-modulated sensor and consider the operator S n x . The proof uses the L p -boundedness of the classical Hilbert transform for 1 < p < ∞ as found in, e.g., [9]. For our purposes, we define the Hilbert transforms H x in x-direction and H y in y-direction by Theorem 3 (Theorem 8.1.12, [9]). For Φ ∈ L p R 2 , 1 < p < ∞, the Hilbert transform defined in (10) exists almost everywhere, belongs to L p R 2 and satisfies with constants c p , d p > 0. (11) is often referred to as the Marcel Riesz inequality for Hilbert transforms.
The 2d Hilbert transform H xy is given as an operator H xy : Using trigonometric formulas, we rewrite S n x into Ωy Ωx Ωy Ωy Ωx Ωy and obtain for 2 ≤ p < ∞ by using the generalized Hölder inequality. The multiplication with the characteristic functions is omitted above due to simplicity of notation. It follows that In order to prove (P c Φ) ∈ L 2 we use relation (29) of the proof of Theorem 1 between non-modulated and modulated pyramid data, i.e., Let T denote one full time period and t ∈ [−T /2, T /2]. For deriving the time-dependent non-modulated pyramid sensor data, we introduce an operator M mod t given by for the periodic tilt Φ mod inducing modulation. As in (28), this tilt is represented by Due to the structure of Φ mod and its compact support on the aperture, it holds Φ mod (·, ·, t) ∈ H 11/6 R 2 Using s n x (·, ·, t) = P n M mod t Φ , (12), and the generalized Minkowski's integral inequality (1) (cf, e.g., [34, Theorem 202], [76]), we obtain which shows P c Φ ∈ L p R 2 , 1 ≤ p < ∞ for the modulated pyramid sensor operators. Note that we omitted the factor − 1 2 from (5). Merely the light which is captured on the telescope pupil Ω influences the pyramid sensor response. We consider only the light falling on the aperture (X Ω · Φ) but use the notation Φ ∈ H 11/6 R 2 for both considered variants of wavefronts with and without compact support on the telescope pupil. Alternatively, one can also define the pyramid sensor operators as P {n,c} respectively. However, since we apply, e.g., Plancherel's theorem we define the operators on the whole R 2 but keep in mind that one can always restrict to the region of the telescope aperture and CCD-detector. More precisely, the pyramid sensor operators are applied to functions with compact support on the pupil Ω and map to functions with compact support on the detector Ω.

Linearization of the roof sensor operators
Linear approximations of WFS operators around the zero phase are sufficient in closed loop AO in which the wavefront sensor measures already corrected and very small incoming wavefronts. The linearization of the operators can be obtained by different ways, e.g., by replacing Linearizations R {n,c,l},lin based on these approximations were already considered in [8,48,73,81]. We concentrate on linear approximations for the roof wavefront sensor operators by means of the Fréchet derivative. We start by calculating the Gâteaux derivatives. Then, we show that the Gâteaux derivatives coincide with the Fréchet derivatives and finally, we evaluate the corresponding linearizations.  (14) and R {n,c,l} y (Φ) ∈ L H 11/6 , L 2 respectively.
Proof. We utilize the representation (cf Taylor's theorem with Lagrange form of the remainder) for a θ = θ (Φ, ψ) ∈ (0, 1). For simplicity of notation, we omit the multiplication with the characteristic function of the aperture in front of every integral and the multiplication with the modulation kernels k {n,c,l} as these functions are independent of the phase Φ anyway.
The Gâteaux derivatives R {n,c,l} Obviously, the Gâteaux derivatives are linear in ψ.
By application of the Cauchy-Schwarz inequality, k {n,c,l} ≤ 1, as well as the Hölder continuity (1) with α = 5/6 and Hölder constant C > 0 of any function in H 11/6 R 2 , the statement results from Theorem 5. The Gâteaux derivatives (14) coincide with the Fréchet derivatives.
Proof. For the proof we use the following assertion stated in, e.g., [3,86]. Hence, it suffices to show that for any Φ 1 , Under the Lipschitz continuity (2) of the cosine function with the Lipschitz constant L > 0 we obtain

Simplified linearized operators L {n,c,l}
x Variations of the finite Hilbert transform operator allow to simplify the linear approximations of the roof sensor model. Let us, thus, consider the following operators. x − x dx (17) and L {n,c,l} y accordingly.
As derived in [26,76], L {n,c,l} are bounded operators on L p R 2 for 1 < p < ∞ due to the structure of the functions k {n,c,l} introducing modulation. According to the pyramid and roof sensor model, we consider the operators as L {n,c,l} : H 11/6 R 2 → L 2 R 2 . In case of no modulation the operator L n x coincides with the finite Hilbert transform -a singular Cauchy integral operator. Using the above defined operators, the linearized roof sensor measurements read as Equation (18) provides two possibilities for wavefront reconstruction. As it is shown in [68,69], when an AO system enters the closed loop, the first term in the forward model gains in importance and an assumption of neglecting the second term in the reconstruction procedure is justifiable. Therefore, one could either use the full linearized roof sensor model or ignore the second term L {n,c,l} x 1 as it is done in several existing algorithms for pyramid wavefront sensor [41], e.g., the Preprocessed Cumulative Reconstructor with Domain decomposition (P-CuReD) [72,73], the Pyramid Fourier Transform Reconstructor [71,74], the Finite Hilbert Transform Reconstructor [70], or the Singular Value Type Reconstructor [38].

Adjoint operators
Several iterative algorithms for solving inverse problems (e.g., the Landweber iteration, steepest descent method or conjugate gradient method for the normal equation) include the application of adjoint operators. In order to make these methods suitable for wavefront reconstruction, we derive the adjoints of the underlying operators. First, we will evaluate the Fourier transforms of the one-term assumptions L {n,c,l} defined in (17) and afterwards use Plancherel's theorem to calculate the corresponding adjoints.
The underlying operators L {n,c,l} are defined from H 11/6 R 2 into L 2 R 2 . In order to calculate the corresponding adjoint operators from L 2 R 2 into H 11/6 R 2 , we introduce the embedding operator and derive the adjoints according to [63] as, e.g., with c n (ξ) = i sgn (ξ) (21) for the non-modulated sensor, for the circularly modulated sensor, and for the linearly modulated sensor. The Fourier transforms of L {n,c,l} y are represented accordingly.
Proof. Similar considerations are contained in [73,81]. Since we examine the 1d Hilbert transform L 2 (R) → L 2 (R) we fix y ∈ Ω x and investigate the operators L {n,c,l} x : H 11/6 (R) ⊆ L 2 (R) → L 2 (R) defined according to (17) By the convolution theorem, the 1d convolution in (24) is a multiplication in the Fourier domain, i.e., As already used in [73,81], the Fourier transforms of the kernel functions are calculated as

Adjoint operators L {n,c,l}
x * in L 2 (R 2 ) As previously mentioned, it is sufficient to derive the adjoints as operators from L 2 into itself and use the embedding operator (19) in order to obtain adjoint operators from L 2 into H 11/6 [63]. where we consider the inner product with respect to L 2 R 2 or L 2 (Ω) respectively.

Conclusion
In this paper, we have considered the mathematical background for the pyramid wavefront sensor which is widely used in Adaptive Optics systems in areas such as astronomy, microscopy, and ophthalmology. The theoretical analysis of the forward operators of pyramid and roof wavefront sensors presented in this paper was aimed at the subsequent development of fast and stable algorithms for wavefront reconstruction from pyramid sensor data. The interference effects as well as the phase shift introduced by the pyramidal prism were neglected for simplicity. The analysis allows any kind of modulation (no, circular, linear) to be applied to the sensors. Due to the closed loop operation assumption, we can linearize the initially non-linear forward operators. Additionally, the closed loop setting allows to simplify the linearized operators further, ending up with measurements described by one term comparable to a finite Hilbert transform in case of the non-modulated sensor. Moreover, for the considered operators we derived the corresponding adjoint operators. These are used in the upcoming second part of the paper [39] which will be devoted to the application of iterative algorithms to the problem of wavefront reconstruction from pyramid wavefront sensor data. An extension of the analysis (linearization and calculation of adjoints) to the full pyramid sensor operator as well as the adaption of the aperture mask to segmented pupils on ELTs [17,40,55,66,67] is dedicated to future work.
Proof of Theorem 1. First, we consider the non-modulated case.
Let Φ ∈ H 11/6 R 2 denote the phase screen in radians coming into the telescope. Using wave optics based models and assuming constant amplitude over the full telescope pupil Ω = Ω y × Ω x , the complex amplitude (wave) ψ aper corresponding to this incoming phase reads as (x, y)) .
Note that ψ aper ∈ L p R 2 for all 1 ≤ p ≤ ∞ due to the compact support of ψ aper , the continuity of Φ, and further the continuity of ψ aper on Ω y × Ω x . In order to assume the continuity of ψ aper on R 2 we slightly modify the telescope aperture mask X Ω and approximate it by X Ω ∈ C ∞ 0 R 2 fulfilling X Ω = X Ω on Ω. The idea of the extended mask is to smoothen the sharp edges of the telescope pupil in order to guarantee ψ aper ∈ H 11/6 R 2 using (x, y)) .
The above assertion is fulfilled for an approximation of the aperture mask denoted by X Ω ∈ C ∞ 0 R 2 .
One possible representation of the modification X Ω can be constructed utilizing the following Lemma for the sets Ω = Ω y × Ω and Ω = (a y − , b y + ) × (a x − , b x + ) for a small > 0, i.e., Ω ⊂ Ω .
Outside of Ω we extend X Ω with zeros and obtain X Ω ∈ C ∞ 0 R 2 .
We investigate the construction of the new smooth aperture mask X Ω in more detail. As in [77], we consider the smooth function Additionally, we introduce f (z) := 1 2 f z for > 0. Using this function, S. L. Sobolev established a method [75] which is utilized for the smoothing of the characteristic function describing the telescope pupil. With the coordinate transformation for the function X Ω ∈ L p R 2 , 1 ≤ p ≤ ∞ we build the average function X Ω ∈ C ∞ 0 R 2 by Altogether, ψ aper is an element of H 11/6 R 2 using that from Φ ∈ H 11/6 it follows e −iΦ ∈ H 11/6 as stated in [6,7]. The adaption of the aperture mask is necessary to guarantee a well-defined mathematical derivation of the pyramid wavefront sensor model. Please note that and therefore X Ω is an arbitrarily good approximation of the aperture mask X Ω . The following physical argument supports the usage of a smoothed aperture mask Ω instead of the one with the sharp edges. Both masks have a compact support. Therefore, on the Fourier domain they are both represented with infinite spectra. Since < d (with d denoting the telescope subaperture size) is small, the difference between the two masks, when looking in the Fourier domain, appears only in the very high frequency components. Because the pyramidal prism is a physical device of finite size, it anyway cuts off high frequency components of the input. Additionally, the sensor brings spatial discretization in the model due to subaperture averaging. As a result, the spectra of the resulting sensor data contain frequencies only up to a given cut-off frequency defined via the subaperture size d as f cut = 1/(2d).
Therefore, in practice there is no difference between the smoothed and the sharp aperture mask. In the following and throughout the paper we write X Ω but keep in mind that we always mean X Ω for a small > 0 for the pyramid wavefront sensor model to be well-defined.
As a next step, we consider the point spread function (PSF) of the glass pyramid. The PSF is the inverse Fourier transform of the optical transfer function (OTF) of the pyramidal prism Within the transmission mask approach [48], the OTF only takes splitting of the light into account and ignores the phase shifts introduced by the pyramid facets. It is represented as a sum of 2d Heaviside functions The 2d Heaviside function is the product of two 1d Heaviside functions and the 1d Heaviside function can be represented as Therefore H 2d (ξ, η) = 1 4 [1 + sgn(ξ) + sgn(η) + sgn(ξ) · sgn(η)] , which for g m (ξ) := sgn ((−1) m ξ) results in Please note that in the above notation, g m is always meant as function of the first variable ξ and g n as function of the second variable η for 2d considerations. Furthermore, F x will denote the 1d Fourier transform in the first variable and F y the 1d Fourier transform in the second variable. With (27) and due to the linearity of the Fourier transform, the PSF of the pyramid is represented as a sum of four PSFs The OTF is a sum of products of functions depending either on ξ or η. Therefore, the inverse 2d Fourier transform reduces to products of 1d inverse Fourier transforms. For P SF mn pyr , we obtain P SF mn pyr (x, y) = The Fourier transforms of the involved constant and signum functions do only exist in a distributional sense. For test functions ϕ, we introduce the delta distribution δ defined as δ, ϕ = ϕ(0). This application is well-defined for continuous functions, e.g., ϕ ∈ H 1/2+ (R), and, on account of this, δ ∈ H −1/2− (R) for > 0. The distribution (p.v. 1 x ) is defined via the Cauchy principal value by for the 1d Hilbert transform defined according to (10). As the Hilbert transform H x : L 2 (R) → L 2 (R) is a well-defined operator (see, e.g., [9,Theorem 8.1.12]), the evaluation of (p.v. 1 x ), ψ is well-defined for ψ ∈ H 1/2+ (R) ⊂ L 2 (R), and thus (p.v. 1 x ) ∈ H −1/2− (R). Specifically, and in the same way Using the notations δ x in case the delta distribution is only applied in x-direction and δ y accordingly, v x := (p.v. 1 x ) and v y := (p.v. 1 y ) as well as v xy := (p.v. 1 xy ), the 2d P SF mn pyr ∈ H −1− R 2 is given by According to the standard description of optical systems, the complex amplitude ψ det coming to the detector plane is the inverse 2d Fourier transform of the complex amplitude after the pyramid ψ aper · OT F pyr which results (by application of the convolution theorem) in a convolution of the incoming complex amplitude ψ aper with the point spread function of the glass pyramid as described, e.g., in [48]. This step can now mathematically be written in the sense of distributions as a shifted PSF distribution applied to the complex amplitude of the incoming phase ψ det (x, y) = 1 2π P SF pyr ((x, y) − (·, ·)) , ψ aper .
Then, by linearity, the four independent beams ψ mn det , m, n ∈ {0, 1}, falling onto the detector are given by ψ mn det (x, y) = 1 2π P SF mn pyr ((x, y) − (·, ·)) , ψ aper = 1 2π The four complex amplitudes are explicitly formulated as , ψ aper (·, y) Now, the intensities on the detector are computed as If we abbreviate ψ aper by ψ and omit the arguments for simplicity of notation the four intensities are evaluated as det (x, y) ψ 10 det (x, y) Taking the sums according to (1), we obtain the non-modulated (indicated by the superscript n) pyramid sensor data s n x in x-direction as which simplifies to This can further be formulated as With Euler's and trigonometric formulas we obtain Taking the sums according to (1), the non-modulated pyramid sensor data s n y in y-direction are written as which is equivalent to ψ (x, y) v y (y − ·) , ψ (x, ·) − v y (y − ·) , ψ (x, ·) ψ (x, y) ψ (·, y) and results in I 0 · s n y (x, y) = i 4π   X Ω (x, y) · exp (−iΦ (x, y)) p.v.

R
X Ω (x, y ) · exp (−iΦ (x, y )) 1 y − y dy Using trigonometric formulas, we separate the time-dependent parts where we used the substitution t = 2πt/T and the definition of the zero-order Bessel function J 0 (x) = 1 π π 0 cos(x sin t) dt cos(x sin t) dt.
All steps of the proof can by performed for the data s {c} y analogously.