© S. Oberti et al. 2022

1 Introduction

The term ‘super resolution’ (SR), as applied in this work, is a technique found to serve the needs of the computational imaging field, whereby the joint application of optical design and signal processing techniques is engineered to obtain higher-resolution (HR) data products from multiple low-resolution (LR) samples. Such SR techniques are relatively standard in vision science whenever image data is available with sub-pixel shifts between LR samples. The first work on image upscaling was published in 1984 (Tsai & Huang 1984) and the term ‘super resolution’ itself appeared around 1990 (Irani & Peleg 1990). Park et al. (2003) provides an excellent review, a sign of how, two decades on, SR imaging has been established as a prominent technique.

Such SR techniques are not new in the astronomical observation and instrumentation communities. One relevant example is the ‘drizzle’ method developed to enhance the resolution of under-sampled HST images from Fruchter & Hook (2002), where the high spatial frequency information of the HST images are recovered by combining sub-pixel dithered images. In this paper, we explore the application of SR to wavefront reconstruction for astronomical adaptive optics (AO), as proposed in the presentation from Oberti et al. (2017).

Reconstructing wavefronts instead of focal plane images is particularly interesting. Indeed, in the field of AO assisted astronomy instruments, the optical systems are typically designed to achieve a close to Nyquist-Shannon sampling of the diffraction limited image. In other words, the images are well sampled and the potential benefits from SR are moderate. On the other hand, the pupil image is typically not limited by diffraction in current AO systems. Specifically, the field stop is sufficiently opened to allow for the set of spatial frequencies of interest to be sensed by the WFS. Moreover, the bi-dimensional wavefront surface can be described by a theoretically unlimited set of spatial frequencies, hence, pushing the SR limit down to extremely small scales.

It turns out that the reconstruction problem is often ill-conditioned because of insufficient LR data or is based on an ill-posed system model. Consequently, some degree of regularisation is required to make it solvable, providing realistic and stable solutions. Even in the absence of multiple LR samples, the use of regularisation provides a means to estimate signal content beyond the native Nyquist-Shannon frequency. We call this strategy ‘statistical SR’, which has by now become well-established in AO. Indeed, the use of regularised reconstruction in the form of an minimum mean square error (MMSE) or maximum a posteriori (MAP) stochastic estimator (Fusco et al. 2001; Conan 2014; Correia & Teixeira 2014); additionally and more generally, methods that draw upon Tikhonov regularisation are now commonplace. Furthermore, SR can also be approached ‘geometrically’, which we aim to explore in this work. This approach is based on the idea that the samples from multiple wavefront samples can be combined to allow for the reconstruction of higher-resolution spatial frequencies than what could be achieved by each single sample alone. Multiple sampled wavefront measurements (mostly across space, however, we could similarly use video sequences across time or a mix of both) are thus combined from all available lower-resolution WFS measurements, provided that each contains some form of unique phase information. In practice, however, each SR facet can only be interpreted asymptotically: it is ‘geometric’ when we approach the limit of least-squares reconstruction and ‘statistical’ when fewer and fewer samples are used.

At this stage, the following observation is paramount: multiple frames are natively granted for SR in classically designed tomographic systems, where multiple WFSs are used, each looking at a single guide-star on a distinctive line of sight. The multitude of diverse WFSs back-projections at different ranges along those lines of sight sample the turbulence in altitude at offset locations (grid points), which is sufficient to perform SR at every altitude range. The exception is the pupil plane where all sampling grids perfectly overlap, making the relative phase information vanish. Depending on the geometric configuration, this may be the case as well at certain altitudes for which the relative grid offsets are integer multiples of one sampling step.

Building on the footsteps of Oberti et al. (2017), the SR concept has been implemented and investigated for several AO applications. For instance, the shape of a high-order DM has been estimated by reconstructing a posteriori a temporal sequence of LR measurements from a single SH WFS, with proper fractional sub-aperture shifts applied to each temporal sample (Woillez et al. 2019). Actually, such synthetic aperture wavefront sensing technique had already been applied in the field of eye aberrometry by Bara et al. (2013).

However, the main applications of SR relate to tomographic AO systems where the signal from multiple WFS are combined. A first prominent example include Ellerbroek & Rigaut (2001) and Wang et al. (2012) who tomographically reconstructed the 3D wavefront with an over-sampling factor of up to 2 from many identical WFS conjugated to the pupil-plane. However, the principles of SR (‘statistical SR’ in this case) were not drawn upon at that time.

In Oberti et al. (2017), numerical simulations were used to describe the application of SR to GALACSINFM (Oberti et al. 2018), the laser tomography AO (LTAO) mode in operation at ESO’s Very Large Telescope (VLT). In fact, SR is already built-in the tomographic reconstruction process as the four WFSs exhibit slight relative mis-registrations that are not compensated for by hardware; instead, these are calibrated and taken into account in the system model, the tomographic interaction matrix. Hence, SR is enabled through the MMSE reconstruction process. More recently, the ELT Multi Conjugate Adaptive Optics (MCAO) system formerly named MAORY now MOR-FEO, incorporated SR in its design. The concept has been analyzed via end-to-end simulations, with the first studies enabling SR in altitude layers without adding misregistrations in the pupil plane, such as in the works by Oberti et al. (2019) and Busoni et al. (2019). Another study taking advantage of the full ‘geometrical’ SR by introducing voluntary shifts and rotations between the sampling grids was carried out by Agapito et al. (2020). With such misregistrations (characterised by a diversity of sampling geometries), super-resolved wavefronts can be reconstructed in all layers, including the DM plane. As a result, it becomes possible to choose to either gain in aliasing reduction at the iso number of sub-apertures or to reduce the number of sub-apertures, while preserving a similar level of performance. In this case, SR opens the door to new possibilities for WFS design choices that allow for an increase in the sub-aperture size and fitting more pixels per sub-aperture, for example, when dealing with LGS spot elongation. Consistent results have been obtained for several other tomographic AO systems, as, for instance, the HARMONI Laser Tomography Adaptive Optics (LTAO) mode (Fusco et al. 2022) the Giant Magellan Telescope (GMT) GLAO and LTAO systems (Van Dam 2017, priv. comm.), the future VLT MCAO system MAVIS (Cranney et al. 2021; Rigaut et al. 2020), and KAPA, the Keck Observatory LTAO system (Wizinowich et al. 2020).

A meta-uniform sampling can be simply achieved with four sampling grids by introducing proper offsets of half a sub-aperture width, as depicted in Fig. 1. It appears clearly that with grids of finite size, the uniformity is not granted at the edge of the aperture. To enable SR at the ground but also over most conjugation altitudes, one obvious way is to break sampling and reconstruction symmetries (or regularities), and diversify the measurements (i.e. avoid redundancy) by adopting a non-uniform sampling. The most straightforward way is to offset and rotate the sampling grids with respect to one another. Other non-uniform sampling options can be envisioned, with some examples presented in Appendix B.

This method is meant to support the control of a N_act × N_act degree-of-freedom deformable mirror (DM) based on measurements of several WFS each of a number of sampling elements, N_sa × N_sa, much smaller than N_act, that is, typically N_sa ≈{1/2 − 2/3}N_act, depending on the number of WFS. SR greatly relaxes the requirements on WFS alignment. For our application domain, we break free from registering the DM to the WFS as accurately as possible, yet once we have completed this step, we need to know the registration precisely in order to forward model this information and estimate higher order data. In other words, SR allows us to trade accuracy for precision; thus, while this distinction may seem elusive, it is actually fundamental to the process. We note that the registration must indeed be calibrated precisely – but not more precisely than what classical AO systems require. Indeed, Agapito et al. (2020) suggests that the super-resolved cases are robust to model calibration errors: the performance is less sensitive to an error in the knowledge of the introduced misregistration amplitude than that of the classical MCAO system with co-aligned lenslet array grids in the pupil plane. Furthermore, for multi-WFS AO systems, asymmetric sampling may reduce the unseen or badly-seen modes (Neichel et al. 2009).

It does not go unnoticed that there are certain commonalities with other phase diverse techniques in imaging and AO. While classical phase-diversity relies on introducing a longitudinal phase encoding in the wavefront plane along the propagation path (Gerchberg 1974), the technique we explore here introduces it transversely with respect to the propagation. The SR technique operates on fused data from all available lower-resolution WFS measurements provided each contains some form of unique phase information – the basic prerequisite for achieving super resolution. The unique phase information can be used to separate the aliased high-frequency from the content in low-frequencies of interest, and the higher-resolution wavefront can be accurately reconstructed. Multiple samples can be acquired sequentially following the same principle, as long as the observed scene (the wavefront in our case) remains static.

In the following, the SR technique is first described from a theoretical point of view in the general and mono-dimensional case (Sect. 2). In particular, we highlight the provided benefits in terms of sensitivity to higher spatial frequencies and aliasing reduction. We show that the signal can be reconstructed without aliasing with uniform or non-uniform samples (under some conditions) at the Nyquist-Shannon average sampling rate. Then, the bi-dimensional case is investigated further through analytical models (Sect. 3.1). For SH WFSs, the sensitivity functions (Sect. 3.2) are studied in the diffraction limited and extended source cases. In Sect. 3.3, the application of SR to multi WFS AO systems is then analysed via numerical models by evaluating the reconstruction error and the noise propagation. Finally, we propose an extension of the SR concept to a single WFS, in the specific case of the pyramid WFS (Sect. 4).

Fig. 1 Illustration of four lenslet arrays, aligned to the left and shifted by 1/2 of a sub-aperture to the right, where the four colours used (red, green, blue, and black) are no longer superimposed.

2 Super-resolution principle

In this section, for the sake of clarity, we first recall the super-resolution principle in the case of a mono-dimensional signal sampled by multiple periodic grids. Then, the main results applicable to AO are highlighted, in the direct space and the Fourier domain. This chapter provides a general description in the sense that the described SR technique can be applied to any continuous signal. In the specific case of adaptive optics, it may be applied to any kind of WFS, measuring the phase, its derivative, Laplacian, or another linear transformation of the phase.

2.1 Super-resolution framework in the direct space

2.1.1 Combination of several sampling grids: general case

Let u := {u_n}_n∈z be a sequence of ℝ and f : ℝ → ℝ a real function. If f is discretely sampled at the points of the sequence u, then the result of the sampling is represented by the distribution:

$$\begin{array}{cc}{f}_u\left(x\right):={\displaystyle \sum _{n=-\infty }^{+\infty }f\left(x\right)\delta \left(x-u_n\right)}& \forall x\in ℝ.\end{array}$$(1)

It follows from Eq. (1) that if we own N sampling grids u¹, u², …, u^N, and that if w is any sequence whose set of points is the union of the points of u¹, u², …, u^N, then

$$\begin{array}{cc}{f}_w\left(x\right)={\displaystyle \sum _{\ell =1}^N{f}_{u^{\ell }}\left(x\right)}& \forall x\in ℝ.\end{array}$$(2)

In other words, the sum of N sampling of f at the points of the sequences u¹, u², …, u^N is a sampling of f at the union of the sequences, that is, f_w. For the AO scientist, each sampling sequence stands for the sampling grid corresponding to one WFS, composed of sampling points in the wavefront space, basically located at the center of each sub-aperture. In the following, without loss of generality and for any considered dimension, we use the term ‘sampling grid’ instead of ‘sampling sequence’. The quantity N stands for the number of WFS samples at play, namely, the number of WFSs or the number of temporal wavefront measurements that are combined. The sampling cell – which we call a ‘sub-aperture’ – is nominally the squared sub-aperture that averages the wavefront within its bounds.

2.1.2 Periodic non-uniform sampling

Following our brief introduction to the general framework, we go on to consider a case of particular interest for the AO community: the periodic non-uniform sampling (PNS) case. The PNS consists in combining several uniform sampling grids into a single non-uniform meta-sampling grid. The PNS is a scheme of deterministic sampling that was first introduced by Kohlenberg (1953).

Let us consider N uniform sampling grids, u¹,u², …, u^N, with the same sampling step, h. In AO, this corresponds to a set of N WFSs with h as the spatial distance between the centers of two neighbouring sub-apertures. The sampling step, h, is classically equal to the sub-aperture width. The first sampling grid u¹ is given by:

$$\begin{array}{cc}{u}_n^1=nh+u_0^1& \forall n\in \mathbb{Z},\end{array}$$

where $u_0^1\in ℝ$ is the reference origin, which can be for instance the center of the telescope pupil. For any ℓ ∈ {1, …, N}, the sampling grid, u^ℓ, is obtained by applying a translation of length ∆_ℓ to u¹. For the AO scientist, ∆ is equivalent to a shift of the WFS sampling grid number ℓ relative to an arbitrary reference WFS, whose ℓ = 1. So, we have ∆₁ = 0 and we suppose ∆_ℓ ∈ (0, h) if ℓ ≠ 1. Then, for every ℓ ∈ {1, …, N}, we have:

$$\begin{array}{cc}{u}_n^1=nh+u_0^1+{\text{Δ}}_{\ell }& \forall n\in \mathbb{Z}.\end{array}$$

Let w be any sequence whose points are the union of those of the sampling grids u¹, u², …, u^N. For example, w can be taken such that

$$\begin{array}{cccc}{u}_q^{\ell }=w_{qN+\ell -1}& \forall q\in \mathbb{Z}& \text{and}& \ell \in \left\{1,\dots ,N\right\}.\end{array}$$(3)

In this way, if ∆₁ < ∆₂ < ∙⋯ < ∆_N, the sequence {w_n}_n∈Z is sorted in increasing order (its points are labelled according to their spatial order) and if n − n′ = N, then w_n and w_n′ are two successive points of a same sampling grid, u^ℓ. It is interesting to note here that the operation described in Eq. (2) resembles a spatial reconstruction process, whereby the sampling points are re-arranged with respect to their spatial coordinates. In other words, the sampling sequences are not concatenated but, rather, organized according to the order of the sampling points. A simple toy example is easy to grasp: a sinusoid that aliases on each WFS, when samples from multiple WFSs are rearranged, is no longer seen to alias (depending on N, ∆_ℓ and the actual frequency of the sinusoid).

2.1.3 From fractional sample shifts to meta-uniform sampling

In the case where the shifts ∆_ℓ are multiple of $\frac{h}{N}$, that is:

$$\begin{array}{cc}{\text{Δ}}_{\ell }=\left(\ell -1\right)\frac{h}{N}& \forall \end{array}\ell \in \left\{1,2,\dots ,N\right\},$$(4)

then Eq. (3) can be expressed as:

$$w_{qN+\ell -1}=u_q^{\ell }=qh+u_0^1+\left(\ell -1\right)\frac{h}{N}=\left(qN+\ell -1\right)\frac{h}{N}+u_0^1,$$

for every q ∈ Z and ℓ ∈ {1, …, N}. Therefore,

$$\begin{array}{cc}{w}_n=n\frac{h}{N}+u_0^1& \forall n\in \mathbb{Z}.\end{array}$$

Then,

$${\displaystyle \sum _{\ell =0}^{N-1}{f}_{u^{\ell }}\left(x\right)=f_{\text{w}}\left(x\right)={\displaystyle \sum _{n=-\infty }^{+\infty }f\left(x\right)\delta \left(x-n\frac{h}{N}-u_0^1\right)=\left(f·{\text{III}}_{\frac{h}{N}}\right)\left(x\right),}}$$(5)

where

$${\text{III}}_{\frac{h}{N}}\left(x\right):={\displaystyle \sum _{n=-\infty }^{+\infty }\delta \left(x-\frac{nh}{N}-u_0^1\right)}$$

is the Dirac comb of period h/N centered at $u_0^1$.

Let us now consider a superposition of N sampling grids composed of sub-apertures of a size, h. In each sub-aperture, the function $f={\overline{s}}_h$ is the result of averaging the signal, s, by means of a convolution with the top hat function χ_h : ℝ → ℝ, which is the uniform density over the interval I_h := [−h/2, h/2],

According to Eq. (5), the signal sampled by the superposition of sampling grids, that we call ‘meta-uniform sampling’, is thus given by the distribution:

$$\begin{array}{cc}{s}_{N,h}\left(x\right)=\left({\overline{s}}_h·{\text{III}}_{\frac{h}{N}}\right)\left(x\right)=\left(s*{\chi }_h\right)·{\text{III}}_{\frac{h}{N}}\left(x\right)& \forall x\in ℝ.\end{array}$$(6)

We note that Eq. (6) suggests that the we can decouple the size of the sub-aperture, h, from the effective sampling step, h/N, that is, that we can generate a ‘meta-sampler’ with regular sampling points located every h/N and whose sub-apertures have a size that is different from the sampling step; this differs from classical AO designs, where the spatial sampling step and the sub-aperture width are the same.

2.2 Mono-dimensional representation in the Fourier domain

The Fourier transform of Eq. (6) is expressed as:

$${\widehat{s}}_{N,h}\left(x\right)=h\widehat{s}\left(k\right)\sin \text{c}\left(\pi hk\right)*{\text{III}}_{\frac{h}{N}}\left(k\right),$$(7)

with k as the spatial frequency in m⁻¹ when considering the wavefront sensing-related spatial sampling case. The Fourier transform of the function s is multiplied by the sinc function that depends only on the sub-aperture size, h. The sinc function can be understood as a smoothing function, standing for the individual sub-aperture transfer function. The sampling is achieved via the convolution by a Dirac comb ofperiod N/h. The latterconvolution leads to replicating the spectrum in the frequency domain due to the periodic spatial sampling. There is no loss of information due to the sinc function, except at the null points located at multiple integers of N/h in the Fourier space.

Let us also define the super-resolution factor F as the ratio between the maximum spatial frequency, k_max, that can be reconstructed without aliasing and the native Nyquist frequency v = 1/2h, that is:

$$F=\frac{k_{\max }}{v}.$$(8)

So, F stands for the gain in resolution achieved with N sampling grids, with respect to the case of a single grid with sampling step, h.

The meta-uniform sampling case (Eq. (7)) carries the following key features: Firstly, from the Nyquist-Shannon Theorem, we know that in the case of Eq. (7), the signal can be reconstructed up to a frequency of k_max = N/2h without aliasing. The super-resolution factor F is therefore equal to N. Thus, combining properly shifted sampling grids leads to increasing the spatial resolution by a factor, F, thus decreasing the aliasing. Secondly, the sub-aperture size is preserved, hence, the sensitivity function remains the same as the native one. Finally, the SR technique offers the capability to design a ‘meta-sensor’ whose sampling step is independent from the sub-aperture size. This independence allows us to decouple the desired spatial sampling period from the sensitivity function in multi-WFS AO systems, which is a new design paradigm.

Considering a super-resolution factor F of 2, it is interesting to note that there is no longer any overlap in the spectra – until the first null of the sinc function. Figure 2 highlights the benefit of such a feature in terms of aliasing mitigation.

The sensitivity is limited by the sinc function in the Fourier space (Eq. (7)) and by the aliasing of its replicated self. From Fig. 2, we can see that beating the v = 1/2h spatial resolution limit leads to aliasing reduction when introducing a SR factor of F = 2 (red curve). However, the sensitivity is greatly reduced with respect to a case with a twice larger native sampling frequency (black curve). With a SR factor of F > 2, it would be possible to go beyond a resolution of 2v = 1/h (not shown on this figure). Nevertheless, the sensitivity would again be limited by the sinc function that depends only on the sub-aperture size. For spatial frequencies higher than the first null of the sinc function, the sensitivity (red curve) is much lower than what would be achieved with a native sampling that is four times greater. This illustration indicates that noise propagation should be analysed with care in presence of SR. Indeed, on the one hand, the sensitivity is limited by the size of the sub-aperture. On the other hand, a larger sub-aperture size allows a higher signal-to-noise ratio to be reached. These two opposite effects ought to be taken into account to evaluate the noise propagation of a super-resolved system.

Fig. 2 Normalised sensitivity of a virtual averaging-and-sampler sensor used for illustration purposes only. The black curve corresponds to the sensitivity function of a sampling grid of step h/2. The dotted curve represents the sensitivity of a sampling grid of step h. The red curve is obtained by dealiasing the spectral sensitivity up to 1/h, that is, the first null of the sinc function with a super-resolution factor of 2; by combining two sampling grids of step h relatively shifted by h/2 with respect to each other.

2.3 Signal reconstruction: Number of sampling grids and optimal offsets

In this section, we wish to address signal reconstruction by analyzing the number of sampling grids required for a reconstruction without aliasing, as a function of their relative spatial offsets. From the Nyquist–Shannon sampling theorem, we know that the choice of N offsets according to Eq. (4) gives a meta-uniform sampling that provides a SR factor of F = N. We analyse here whether or not there are other ways to select the offsets in order to reach a higher SR factor, or at least to preserve it, without adding extra sampling grids.

In Appendix A, we analyze the reconstruction of band-limited signals by convolutions with sinc interpolation functions within the PNS framework. Under conditions that simplify the analysis (but that may be too restrictive), we show that with N sampling grids, F cannot be higher than N, independently of the offsets choice. Moreover, the meta-uniform sampling is the only solution with N sampling grids that reaches F = N, meaning that it allows us to reconstruct the signal without aliasing until the frequency Nv. With non-uniform configurations, that is, with samples that are not equally spaced, more than N sampling grids are necessary. Otherwise, any departure from the meta-uniform sampling scheme will generate some reconstruction error. We can show that this error term is bounded and converges to zero when tending towards the meta-uniform sampling. However, we see hereafter that other hypotheses and more elaborated reconstruction approaches allow us to show that some well-chosen non-uniform sampling schemes can reach the same performance as the meta-uniform sampling.

In practice, it may be tricky to place the sampling grids exactly at the positions required to generate a meta-uniform sampling. So, we considered how many sampling grids would be needed to achieve a given SR factor F, even if the sampling grids are disposed without particular care, possibly with random offsets. We show that with N = 2n₀ − 1 randomly positioned sampling grids, where n₀ = ⌈F⌉ is the ceiling function of F, we can always reach a SR factor at least equal to F.

To be specific and give some concrete examples, if we want to reconstruct a signal whose maximal frequency satisfies v < k_max ≤ 2v, we need 1 < F ≤ 2, and n₀ = 2. So, in the mono-dimensional case, we may use two sampling grids providing a meta-uniform sampling, or three sampling grids that are randomly positioned, excluding the trivial case where some offsets are multiple of the native sampling step. In the bi-dimensional case, we can show that four sampling grids are needed with a meta-uniform sampling scheme, while nine randomly positioned sampling grids are sufficient.

As previously mentioned, with a more elaborated approach based on Lagrange polynomials, it is possible to reconstruct non-uniformly sampled signals without aliasing, even outside the PNS framework. Indeed, the Paley–Wiener–Levinson theorem (Marvasti 2001; Wiener 1930; Paley & Wiener 1934; Levinson 1940), which generalizes the Whittaker–Shannon–Kotelnikov sampling theorem (Marvasti 2001) from uniform to non-uniform samples, states that a band-limited signal can be reconstructed from its samples if the two following sufficient conditions are respected.

Firstly, according to Beutler (1966), the average sampling rate k_e, that is the inverse of the average sampling step h_e, satisfies the generalized Nyquist-Shannon condition:

$$k_{\text{e}}=\frac{1}{h_{\text{e}}}\ge 2k_{\max }.$$(9)

In the PNS framework, h_e is always equal to h/N. To reach a SR factor F = N, at least N sampling grids are necessary.

Secondly, the Kadec condition (Kadec 1964) states that the non-uniform sampling should be constructed by relocating the points of a uniform sampling of step 1/2k_max no further than 1/8k_max from their original location.

We have seen that under some conditions, both uniform and non-uniform sampling enable a theoretically perfect reconstruction of the sampled band-limited signal. In the PNS framework, the N sampling grids’ offsets shall be close to but not necessarily equal to a multiple of h/N, in order to achieve optimal SR (F = N). We show in Sect. 3.3 that numerical simulations confirm the results from the first part of this paper. The meta-uniform sampling is indeed optimal, in the sense that it provides the best reconstruction performance. Nevertheless, it is important to note that the performance is improved as soon as any non-zero shift is introduced. Moreover, there is a shift range around the optimal shift that provides stable performance, which indicates its robustness. This shift range is consistent with the Kadec condition (Kadec 1964) for a perfect reconstruction with non-uniform sampling.

3 Super-resolution wavefront reconstruction with Shack-Hartmann WFSs

3.1 Bi-dimensional model

We go on to transpose the previous models to a SH-WFS, whose measurements s(x) are well approximated by the geometrical-optics linear model (Rigaut et al. 1998):

$$s\left(x\right)=G\phi \left(x\right)+\eta \left(x\right),$$(10)

where G is a phase-to-slopes linear operator mapping aperture-plane guide-star wavefronts φ(x) onto WFS measurements over a bi-dimensional space indexed by x = [x, y]; η(x) represents white noise due to photon statistics, detector read noise and background photons. Both φ and η are zero-mean functions of Gaussian probability distributions and known covariance matrices ∑_ϕ and ∑_η respectively. Noise is assumed both temporally and spatially uncorrelated.

It can be shown for the nominal case when the sub-aperture size d coincides with the sampling step h (Correia et al. 2017) that:

$$G=\text{III}\left(\frac{x}{h}\right)\times \left[\text{Π}\left(\frac{x-1/2}{d}\right)\otimes \nabla \right],$$(11)

where ⊗ is a two-dimensional (2D) convolution product, × is a point-wise multiplication, ∏(·) is the top-hat separable function, ∇ is the gradient operator, and III(x) is a comb function (a bi-dimensional sum of Dirac delta functions) that represents the sampling process. In the following, we consider the classical case of d = h.

Now the following argument stands: if we properly combine four sets of measurements provided by convolutional operators:

$$G=\left[\begin{array}{c}{G}_1\\ G_2\\ G_3\\ G_4\end{array}\right]$$(12)

with a relative offset between them of h/2, corresponding to the meta-uniform sampling in 2D, that is:

$$\begin{array}{lll}{G}_1\hfill & =\hfill & \text{III}\left(\frac{x}{h}\right)\times \left[\text{Π}\left(\frac{x-\left[1/2,1/2\right]}{h}\right)\otimes \nabla \right],\hfill \\ G_2\hfill & =\hfill & \text{III}\left(\frac{x-\left[0,1/2\right]}{h}\right)\times \left[\text{Π}\left(\frac{x-\left[1/2,1/2\right]}{h}\right)\otimes \nabla \right],\hfill \\ G_3\hfill & =\hfill & \text{III}\left(\frac{x-\left[1/2,0\right]}{h}\right)\times \left[\text{Π}\left(\frac{x-\left[1/2,1/2\right]}{h}\right)\otimes \nabla \right],\hfill \\ G_4\hfill & =\hfill & \text{III}\left(\frac{x-\left[1/2,1/2\right]}{h}\right)\times \left[\text{Π}\left(\frac{x-\left[1/2,1/2\right]}{h}\right)\otimes \nabla \right],\hfill \end{array}$$(13)

it leads to a meta-SH system with half the sampling step and the same averaging cell size or sub-aperture width as in the native case:

$$\overline{G}=\text{III}\left(\frac{x}{h/2}\right)\times \left[\text{Π}\left(\frac{x-1/2}{h}\right)\otimes \nabla \right],$$(14)

therefore lifting the sub-aperture imposed spatial cut-off frequency, which becomes, in the case of a SR factor of F = 2:

$$k_{\max }=2v=\frac{1}{h}.$$(15)

Based on the fact that the SH-WFS measurement model is a set of convolution integrals, the treatment in the spatial-frequency domain is straightforward. Let the Fourier-domain representation of Eq. (10):

$$\tilde{s}\left(\kappa \right)=\tilde{G}\tilde{\phi }\left(\kappa \right)+\tilde{\eta }\left(\kappa \right),$$(16)

with κ = (k_x, k_y) ∈ ℝ² as the frequency vector and symbol $\underset{·}{~}$ used for Fourier-transformed variables. The time-dependence could be added in with recourse to the frozen-flow hypothesis. Using common transform pairs for the individual operations (Oppenheim & Willsky 1997), the Fourier representation of the measurements in Eq. (10) given in Eq. (16) can be expanded to:

$$\begin{array}{lll}\tilde{s}\left(\kappa \right)\hfill & =\hfill & 2i\pi h{\displaystyle \sum _m(\kappa +m/\left(h/N\right)\tilde{\phi }\left(\kappa +m/\left(h/N\right)\right)}\hfill \\ \hfill & \hfill & \times \widetilde{\text{Π}}\left(\kappa h+m\right)\times {\text{e}}^{i\pi \left(\kappa h+m\right)}+\tilde{\eta }\left(\kappa \right).\hfill \end{array}$$(17)

From Eq. (17), it is apparent that the wavefront sensing operator $\tilde{G}$ is not purely a spatial filter except for functions within the SH-WFS pass-band (Ellerbroek 2005).

3.2 Trade-off between sensitivity and noise propagation

In practical terms, any interested AO practitioner may wonder (quite rightly) whether SR signal reconstruction above the nominal Nyquist–Shannon cut-off is associated with a somewhat greater noise propagation, which could limit the prospective SR advantages. We first address this key point in the spatial-frequency domain and later using the Karhunen-Loeve (KL) modal basis. We show in either case that SR grants comparable noise propagation without sacrificing effective sensitivity.

3.2.1 Sensitivity functions revisited: Fourier S/N

We start by revisiting commonly-accepted expressions for sensitivity functions, often called the Fourier signal-to-noise ratio (S/N; Vérinaud 2004). The sensitivity defined as a S/N (the ratio between the produced signal and the noise affecting the measurement) is a function of the effective number of photons and spot size on each sub-aperture, as follows

$$S/N={\text{Ξ}}_{\text{SH}}^2/{\sigma }_{\text{ph}}^2.$$(18)

The SH-WFS sensitivity can be readily computed from Eq. (17) for functions in the pass-band assuming the compact form from Correia & Teixeira (2014):

$${\text{Ξ}}_{\text{SH}}^2={\left(2i\pi h{k}_x\right)}^2{\left(\frac{\sin \left(\pi k_x\right)}{\pi k_{x}h}\right)}^2={\text{Ξ}}_{\text{grad}}^2{\text{Ξ}}_{\text{aver}}^2,$$(19)

where we assume the effective SH signal to be the angle of arrival on each sub-aperture instead of the phase difference at opposite sub-aperture edges (there is a normalisation factor, h, between both). The 2D case is now straightforward to obtain. Let the following photon and detector noise expressions (Roddier 1999; Hardy 1998):

$${\sigma }_{\alpha ,\text{ph}}^2=\frac{1}{8}\frac{1}{n_{\text{ph}}}{\theta }_{\text{spot}}^2,$$(20)

$${\sigma }_{\alpha ,\text{det}}^2=\frac{1}{12}\frac{{\sigma }_{\text{e}}^2}{n_{\text{ph}}^2}\left(\frac{N_{\text{s}}^2}{N_{\text{D}}}\right){\theta }_{\text{spot}}^2,$$(21)

with n_ph as the number of incident photons per frame and per sub-aperture, θ_spot the spot size, σ_e the rms number of detector read-out photo-electrons, and N_s and N_t as the window used for processing and sub-aperture diffraction sizes respectively. In the remainder of this paper, we limit the analysis to photon noise only, since it is not far-fetched to consider practically zero-noise detectors today.

In diffraction-limited cases with point-sources, then ${\theta }_{\text{spot}}=\frac{\lambda }{h}$. It it straightforward to check that if we chose to split the aperture into n_SA-times as many sub-apertures (linearly across the pupil), the photon-noise is increased by a factor of $n_{\text{SA}}^4$ that results from the product of $n_{\text{SA}}^2$ fewer photons per sub-aperture and a diffraction spot size that is twice as large, resulting in another $n_{\text{SA}}^2$ factor for the bi-dimensional case. Considering all the measurements over the whole pupil, an average over $n_{\text{SA}}^2$ times more measurement points is performed, resulting in a final $n_{\text{SA}}^2$ sensitivity decrease factor. This factor is at the origin of the so-called ‘full-aperture gain’, which leads to the use of a single full aperture spot measurement for tip and tilt modes. A comparative plot shown in the top panel of Fig. 3 illustrates this fact. This would lead to the conclusion that using a super-resolved meta-SH-WFS setting leads to sensitivity gains over a twice-resolved SH and even to a four-times-resolved SH over the original sampling band $\left[0:\frac{1}{2h}\right]$ and some portions of the $\left[\frac{1}{2h}:\frac{1}{h}\right]$. Therefore, SR appears very well suited for AO applications with diffraction limited WFS, for instance, those sensing with a corrected near-infrared natural guide star (NGS).

We do, however, expect to deploy multi-WFS mostly on laser-tomographic AO cases. In this case, splitting the aperture into fewer sub-apertures leads to small (or none at all) sensitivity gains because the spot is no longer a function of the diffraction but, instead, it is of a fixed size. In the bottom panel of Fig. 3, we can inspect the relative sensitivities: they are roughly the same until half the control radius of the nominal case of $\left[0:\frac{1}{4h}\right]$; by this time, the larger sub-aperture size averaging function (the sinus cardinal) starts bending down the sensitivity curve. In the SR case, we can only be as sensitive as the red curve, yet with the aliasing cut-off frequency being a function of the over-sampling used. The periodic nature of the sensitivity function with an initial zero at κ = 1/h m⁻¹ tells us that a SR factor of up to two is of practical relevance but not much beyond. Indeed, the high noise amplification at the frequencies with close-to-zero sensitivity is to be traded-off for increased sensitivity at frequencies located beyond the native cut-off frequency.

If we opt to evaluate the noise propagation instead, we can see from Fig. 4 that SR would not lead to higher figures at low spatial-frequencies and that it would extend the effective maximum reconstructed frequency avoiding the zeros of the sine² function for which the noise propagation is infinite. The sensitivity evaluated on a single sub-aperture tells us little about how it is distributed over modes spanning the aperture. For tip-tilt, it is easy to compute since the global estimate is a straight average of each sub-aperture’s measurements, resulting in an effective $1/n_{\text{SA}}^2$ gain (diffraction-limited sensing) or no gain (seeing or LGS spot limited sensing). In the next section, we present a more general case of modal reconstruction for multi-WFS AO systems.

Fig. 3 Fourier S/N functions. Top: Fourier S/N for a diffraction-limited wavefront sensing case (Eq. (6)). Bottom: seeing-limited or laser spot-sensing case. The number of incident photons is considered the same for all configurations. The bottom panel assumes a spot size with one unit.

Fig. 4 Noise propagation coefficient as a function of spatial frequency. This plots are computed as the inverse of the S/N functions from Fig. 3. The fluxes are scaled by the collecting area, namely, the square of the sub-aperture width.

3.2.2 Noise propagation: modal analysis

In order to properly treat the general case, the noise propagation is considered on a per-mode basis, following Rigaut & Gendron (1992):

$$\text{n.p.}=\text{diag}\left\{R{R}^T\right\},$$(22)

where R stands for the reconstructor matrix. We use this formulation instead of the original one provided by Eq. (14) in Rigaut & Gendron (1992) due to the use of regularised reconstruction. If the least-squares solution were adopted, we would find Rigaut’s (D^TD)⁻¹, with D the system modal interaction matrix concatenating the sensor’s response to each mode.

Figure 4, which is essentially the inverse of Fig. 3 as can readily be seen, depicts the noise propagated as a function of spatial frequency. When accounting for the sub-aperture size in the computation of the noise propagation, we can confirm that SR configurations do not propagate more noise than classical configurations, along lower spatial frequencies.

3.3 Numerical application to multi-WFS-AO systems

In this section, we numerically simulate the SR concept applied to a multi-SH-WFS AO system. First of all, we can easily confirm that SR can be promoted by building a multi-SH-WFS model and measuring the Interaction Matrix (IM) between the phase and the meta-SH-WFS, combining several LR WFS, within the configuration described in Table 1.

We go on to compute a set of modal IMs, between the 1364 KL modes as input and the four SH-WFS slopes as output, this for all configurations (40 × 40, 20 × 20, 10 × 10, co-aligned, or shifted). The rank of the IM indicates the number of degrees of freedom that the super-resolved system is sensitive to. The singular value decomposition (SVD) of the IM highlights the sensitivity of the system to the singular modes. In other words, the singular value associated to each singular mode indicates how well such a mode is seen by the system. Figure 6 highlights the benefit provided by the shifted configurations. In the classical cases (solid lines), the sensitivity drops to 0 for a KL index equal to the number of slopes measured by each individual SH-WFS. On the other hand, with SR (dashed lines), the sensitivity is enhanced and the sensitivity limit is pushed towards that of the co-aligned cases with double native resolution. Concretely, the super-resolved 10 × 10 (resp. 20 × 20) configuration provides some additional sensitivity up to the singular modes located at the sensitivity limit for the 20 × 20 (resp. 40 × 40) co-aligned configuration.

We then perform an open-loop reconstruction of the phase by inverting the previously computed IMs in a MMSE fashion. We note that each KL mode pastthe LR SH-WFS cut-off frequency v = 1/(2h) aliases onto its corresponding low-frequency spectrum, yet the inversion of the meta SH-WFS IM plays the role of a numerical de-scrambler as if the samples were spatially reordered and, therefore, de-aliased. Figure 7 indicates that the SR configurations significantly reduce the reconstruction error with respect to the co-aligned cases. The cut-off frequency is twice higher than the respective classical case with as many reconstructed modes as number of sensing element per WFS $\frac{\pi }{4}{n}_{\text{SA}}^2$. Moreover, the aliasing features vanish close to the cut-off frequency for the super-resolved cases. Recent end-to-end closed-loop simulations performed in the context of ELT tomographic AO systems (MORFEO and HARMONI LTAO) confirm the benefits of the SR technique in terms of aliasing mitigation at an iso number of sub-apertures, as well as the potential reduction of the number of sub-apertures, while preserving a similar level of performance as the classically dimensionned AO system (Fusco et al. 2022). Summing the modal residual variances shown in Fig. 7 provides an estimation of the total reconstruction wavefront error (WFE) for the shifted cases.

In Fig. 8, we evaluate the residual WFE as a function of shift amplitude in percentage of a subaperture width. A first observation is the performance gain that can be obtained when reconstructing four times more modes than the number of degrees of freedom of a single WFS (colored squares). It can be explained by the ‘statistical SR’ properties. In this case, the MMSE process extrapolates the measurements statistically towards larger spatial frequencies, by assuming the power spectrum described in the Kolmogorov covariance priors. Next, we focus on the ‘geometric SR’ properties. We can observe a periodicity of the achieved WFE with a periodic shift that is a multiple integer of the sub-aperture width, as expected. However, with any other fractional shift amplitude, the WFE is improved with again an optimal value of one-half (1/2) of a sub-aperture, which corresponds to a meta-uniform sampling generated by the super-resolved system. At this optimal offset of half of a sub-aperture width, it is interesting to note that the 10 × 10 (resp. 20 × 20) shifted configuration achieves a WFE that tends towards that from a 20 × 20 (resp. 40 × 40) co-aligned system with a number of reconstructed modes equal to the native number of degrees of freedom $\left(\frac{\pi }{4}{n}_{\text{SA}}^2\right)$.

As already pointed out, the ‘geometric SR’ provides an improvement of the WFE as soon as there is a small offset introduced. So, even non-optimal shifts induce a dramatic improvement of the reconstruction quality. Furthermore, the stability of the WFE around the optimal shift values is remarkable. In summary, Fig. 8 confirms on one hand that the meta-uniform sampling is optimal, as then demonstrated in Appendix A. On the other hand, it appears unnecessary to accurately tune the system with such optimal offsets since the performance appears weakly sensitive to the offset value. Actually, this observation is consistent with the Kadec condition (Kadec 1964), which corresponds to a range of shifts of ±12.5% of a sub-aperture width around the 1/2 sub-aperture offset, within which the samples can be perfectly reconstructed. In our simulation, the MMSE reconstructor was not optimized for each shift value, hence, there is a slight performance degradation when going away from the meta-uniform sampling case. The Kadec range can be interpreted as a tolerance to uncontrolled misalignments, as long as these ‘misregistrations’ are calibrated or at least known. Such feature opens the door to a diversity of geometrical transformation breaking the symmetries of the wavefront sampling topology. Furthermore, Agapito et al. (2020) has suggested that SR is more robust to model calibration errors than ‘classical’ AO system of the same dimension.

In the case of a multi-WFS tomographic AO system, super resolution is ‘built-in’ when it comes to the turbulent layers conjugated in altitude. Indeed, when projecting the lenslet array grids onto the altitude layers, the off-axis WFS grids are shifted proportionally to the conjugation altitude. It is straightforward to assume that there are some altitude layers where the resulting shift is a multiple integer of the sampling step, hence, no super resolution is granted in this case. An optimized geometrical arrangement, possibly including a rotation of the lenslet array grids with respect to each other, would allow for this effected to be mitigated.

Finally, we evaluate how SR behaves in terms of noise propagation. Figure 9 (essentially a numerical confirmation of Fig. 4) depicts the noise propagated as a function of the KL mode number for a combination of four WFS that are: (1) co-aligned and (2) optimally shifted by half a sub-aperture width of SH-WFS. The noise propagation follows the expected (n + 1)⁻¹ power law with the shifted system (on which SR is applied) not exceeding the co-aligned case. We note that for each of the tested cases, the noise propagation within the original cut-off frequency range is almost indistinguishable from that of the co-aligned system; whereas when it goes past the original Nyquist-Shannon cut-off frequency, it rolls off, asymptotically saturating to the (n + 1)⁰ regime.

Fig. 5 SR configurations illustrated via a projection in the pupil plane of one sub-aperture from each WFS (1 color for each): (a) 2D surface is improperly paved with shifts along the sampling grid’s axes; (b) 2D surface is properly paved with the addition of shifts at 45° with respect to the sampling grid’s axes, as in Eq. (13).

Fig. 6 Singular value decomposition of a multi-WFS mono layer modal IM. The singular values are plotted as a function of the singular modes index. The solid lines corresponds to 4 co-aligned SH-WFS while the dashed lines represent the shifted configurations shown in Fig. 5b.

Fig. 7 Open-loop reconstruction of an incident phase corresponding to a Kolmogorov turbulence with r0 = 12.6 cm. The dashed lines represent the co-aligned cases. The solid lines stand for the super-resolved cases achieved by combining four SH WFS that are properly shifted by half (1/2) of a sub-aperture. The same number of modes is reconstructed for the co-aligned cases and their respective shifted cases, that is, up to $4.\frac{\pi }{4}{n}_{\text{SA}}^2$ modes, which is the spatial resolution that would be provided by a single WFS with twice as many sub-apertures per diameter.

Fig. 8 Residual wavefront error from open-loop reconstruction as a function of shift amplitude for the shifted and super-resolved configurations. In each case, up to four times the native number of modes $\left(4.\frac{\pi }{4}{n}_{\text{SA}}^2\right)$ are reconstructed. The dashed black lines represent the Kadec condition boundaries within which the non-uniform samples can be perfectly reconstructed. The colored squares stand for the reconstruction error from the co-aligned case when $\frac{\pi }{4}{n}_{\text{SA}}^2$ modes are reconstructed.

4 SR WFR with pyramid WFSs

We argue that the pyramid wavefront sensor (P-WFS) also benefits from an effect that is comparable to SR in the Shack-Hartmann WFS (SH-WFS). Unlike the latter, the SR principle can be directly applied to a single P-WFS when we allow a more general formalism than Eq. (5), whereby the samples are obtained from transformations of the initial function. In the case of the P-WFS, each quadrant applies a different transformation, each of which is required for proper signal reconstruction. However, two diagonally opposite quadrants measure the same information for the spatial frequencies larger than the modulation radius (Guyon 2005). Thus, the mono-dimensional SR formalism presented in this paper (see Sect. 2) can, in principle, be applied to the corresponding modes for each pair of quadrants whose pixel grids are shifted in one direction by a fraction of a pixel.

Here, the multiple samples are nominally four diffraction patterns produced at the corresponding re-imaged pupil planes, however, with the sampling grid pixels shifted with respect to one another in such a way that their sampling becomes interlaced by a fractional pixel. By adjusting the P-WFS apex angles, one such configuration can be readily obtained. The same process can, in principle, be carried out simply by rotating the pyramid optic with respect to the detector. This relaxes the specification and manufacturing constraints, yet in order to obtain the optimal SR-enabling configuration, we may end up re-introducing anew such design constraints that we wanted to relax in the first place. Clearly, more investigation on the optimality and robustness are required and we invite the community to contribute to this endeavour.

To exploit the SR enhancement however, the use of the four diffraction intensity patterns (i.e. intensity maps) instead of the customary slope maps appears to be required. It is known that the ‘x’ and ‘y’ slopes-maps contain the complete information on perfectly aligned systems, that is, exactly registered quadrants and a perfect four-sided prism with homogeneous and stationary pupil illumination (Deo et al. 2018). However, additional slope maps are required when these conditions fail. In our opinion, the ‘extension’ proposed in Deo et al. (2018) is unnecessary since it represents a full-rank linear combination of diffraction intensity patterns. This intermediate step adds nothing to the use of intensity maps directly, as it neither adds nor removes useful information to that originally contained in the pixels’ intensities directly. Furthermore, it can only be detrimental to the pipelined pixel processing by creating additional idle time before matching pixels become available during detector read-out.

Mathematically, out of the eight possible linear combinations (we can take each intensity map either positive or negative and there are four of these), one full-rank slopes-map is:

$$\left[\begin{array}{c}{s}_x\\ s_y\\ s_{xy}\\ s_{\text{t}}\end{array}\right]=\left[\begin{array}{cccc}1& -1& 1& -1\\ 1& 1& -1& -1\\ 1& -1& -1& 1\\ 1& 1& 1& -1\end{array}\right]\left[\begin{array}{c}{i}_1\\ i_2\\ i_3\\ i_4\end{array}\right]/F_{\text{t}},$$(23)

where F_t is a scalar value that represents the total flux captured on all the four re-imaged pupils. We have made a slight change to the last entry, s_t, which in Deo et al. (2018) is the total flux s_t = [1,1,1,1] , [i₁, i₂, i₃, i₄]^T but leads essentially to the same results for the rank of either transformation matrix is 4. As a matter of fact, we do go beyond the recommendations of Deo et al. (2018) in that we do not restrict ourselves to developing a misalignment mitigation strategy. Under the SR framework, we claim that by expressly misaligning the PWFS system, we can enhance sensitivity and therefore gain in the delivered optical performance.

The use of intensity-maps has obvious implications for realtime processing (i.e. twice as large reconstruction matrix) that needs to be carefully counter-balanced against the performance gains. To illustrate this behaviour, Fig. 10 shows the SVD decomposition from a 16 × 16 lenslet based P-WFS system with a 33 × 33 regular DM. We clearly see that the introduction of a relative offset (translation) in the pixels sampling grid leads to sensitivity gains provided that the intensity-maps are processed. The opposite is obtained when the slopes-maps are processed instead. The mathematical demonstration and interpretation will be provided in a forthcoming paper.

Fig. 9 Noise propagation coefficients, computed as diag{RRT}, shown in logarithmic scale and scaled by the flux collecting surface of the sub-aperture. The shifted cases are plotted with solid lines. For the sake of visibility, dashed lines are overlayed at the level of the spatial frequencies where the aliasing starts being prominent for the co-aligned cases. Along the low order modes, the noise propagation is quasi identical for all cases and fit the expected law in (n + 1)−1. Then, when exceeding the original Nyquist–Shannon cut-off frequency, the noise propagation coefficient converges towards a constant asymptotic limit.

Fig. 10 SVD decomposition on a 16 × 16 P-WFS system with a 33 × 33 DM. Top: intensity map processing. Bottom: slopes-maps processing. We clearly see that using the intensity-maps on an offset pixel sampling system leads to sensitivity gains, whereas the same is not true if the slopes-maps were to be processed.

5 Conclusion and outlook

Super-resolution is a computational imaging technique that is aimed at upscaling a signal’s spatial resolution from multiple lower resolution samples. Wavefront reconstruction – a process at the very heart of astronomical adaptive optics systems – provides a particularly fertile setting for the application of SR. Unlike in imaging systems where the ultimate spatial resolution at the focal plane hits the fundamental diffraction limit, operating in the wavefront domain in or near the pupil plane offers the potential to significantly increase the sensing resolution, within the technological ability to sample the wavefront (spatially and temporally) subject to the inescapable concerns around adequate S/N.

To make SR viable, we have identified several ways to introduce unique phase signatures in the transverse plane to the beam propagation, the simpler one being a straight offset. As it stands, SR promises to lift instrument design assumptions, for instance, that n × n measurement samples are needed to drive a n × n actuators DM, that the registration between the DM and the WFS sampling grid needs to be finely controlled, and that the guide-star asterism be regular (along with a few other aspects), leading to a major overhaul in how the AO community may approach the designs of future instruments. By decoupling the WFS sampling step from the averaging cell size (sub-aperture width), we are led to a paradigm shift where the same AO control performance can be achieved with fewer WFS sampling points. This, in turn, leads to fewer detector pixels, smaller systems and ultimately a lighter computational load. Another strategy is to optimize the performance while keeping the same number of sensing degrees of freedom, by reconstructing modes beyond the cut-off frequency and mitigating the aliasing. Several tomographic AO systems in operation (GALACSI NFM) or in the design phase (MORFEO, HARMONI LTAO, MAVIS, KAPA, GMT’s GLAO and LTAO, etc.) could benefit from the implementation of SR in one way or the other.

With N regular sampling grids of step h, the optimal SR sampling scheme is obtained by introducing fractional offsets that are multiple of $\frac{h}{\sqrt{N}}$ in the bi-dimensional case. In this way, we can generate a ‘meta sampler’ corresponding to the so-called ‘meta-uniform’ sampling scheme. However, the non-uniform sampling case is appealing because it relaxes the absolute alignment constraints. Concretely, the location of the grids is flexible. The optimal reconstruction performance can even be reached within the Kadec condition, that is, the departure from the meta-uniform sampling by less than a quarter of the Nyquist–Shannon period. Furthermore, we are yet to establish optimal, or at least practical, SR sampling geometries for future tomographic AO systems with a non-square number of WFSs. With six WFS (MORFEO, HARMONI LTAO) or eight WFS (MAVIS), the solution is not obvious. Preliminary geometry optimization results can be found in Cranney et al. (2021). Some other ways to break the system symmetries are proposed in Appendix B.

We have shown that SR does not entail larger noise propagation coefficients and that is applicable to multi-Shack–Hartmann WFS systems as well as to single pyramid WFS. In either case, the relative alignment accuracy is traded-off for the precision with which the registration is known to inform the model-based SR reconstruction. We do not exclude – in fact, we very much welcome – data-driven, machine-learning algorithms that can further optimize this trade-off. The very fast developments of these techniques in the field of SR applied to imaging may be inspiring. A review of deep learning approaches for imaging SR is, for example, provided in Wang et al. (2021).

Acknowledgements

This work benefited from the support of the WOLF & APPLY project ANR-18-CE31-0018 & ANR-19-CE31-0011 of the French National Research Agency (ANR), the project Fortalecimiento del Sistema de Investigación e Innovación de la Universidad de Valparaíso UVA20993 and the project ANID/REDES 190071. It has also been prepared as part of theactivities of ORP H2020 Framework Programme of the European Commission’s (Grant number 101004719). Authors are acknowledging the support by the Action Spécifique Haute Resolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES, the support of LABEX FOCUS ANR-11-LABX-0013 and the support of the programme d’investissement avenir ANR-21-ESRE-0008. We would like to thank Christophe Verinaud and Guido Agapito for reading parts of the manuscript and providing pertinent comments. Finally, we wish to acknowledge the fruitful discussions with Pierre-Yves Madec, Michel Tallon, Miska Le Louarn and Jason Spyromilio at the inception of this idea back in 2016. Thank you for your critical review and wise feedback.

Appendix A Signal reconstruction and optimal sampling grid positions

For the sake of completeness, we provide in this appendix mathematical arguments that help us estimate the number of sampling grids required for a reconstruction without aliasing, as a function of their relative spatial offsets. This study can be done with simple mathematical tools, if we consider a specific context that we describe below.

We propose to reconstruct a band-limited function f: ℝ → ℝ sampled by a set of N grids with an interpolation formula of the form:

$$f\left(x\right)={\displaystyle \sum _{\ell =1}^N{\displaystyle \sum _{n=-\infty }^{+\infty }f\left(u_n^{\ell }\right)s_{\ell }\left(x-u_n^{\ell }\right),}}$$(A.1)

where for each ℓ ∈ {1, …, N}, the sequence $u^{\ell }={\left\{u_n^{\ell }\right\}}_{n\in \mathbb{Z}}$ gives the points where the function, f, is sampled. So, we have to find a suitable set of reconstruction functions s₁, …, s_N. The existence of such functions depends on the spectral properties of fand on characteristics of the set of sampling grids, such as their number, their step and their positions. One restriction lies in the fact that formula (A.1) is a convolution. It is a classical hypothesis, but it reduces the set of the possible reconstruction functions, s₁, …, s_N, and therefore the admissible positions of the grids. However, we can see that formula (A.1) leads to sinc-type reconstruction functions, which as such can be implemented at low numerical cost.

In the context of the PNS framework, we then suppose that all the grids have the same step h. However, we are free to set them at any position. Thus, they take the following form:

$$\begin{array}{ccc}{u}_n^{\ell }=nh+u_0^{\ell }& \forall n\in \mathbb{Z}& \ell \in \left\{1,\dots ,N\right\},\end{array}$$(A.2)

where $u_0^{\ell }$ indicates their initial position, which is supposed to be different for each of them.

In the following, we establish the conditions, relating the maximal frequency of f to the characteristics of the sampling sequences, which ensure the existence of the functions s₁, …, s_N allowing us to reconstruct f with the formula (A.1). For a fixed maximal frequency, there are many possible options for the positions of the sampling grids, provided that we own a sufficient quantity of them. However, we are interested in finding the positions that need the smallest number of sampling grids as possible. As we go on to see, there is no surprise in this regard. With a set of uniform sampling grids, the meta-uniform sampling (Eq. (4)) is optimal: it allows us to reach a given maximal frequency with the lowest number of sampling grids.

A.1 Reconstruction conditions

To find the reconstruction functions s₁, …, s_N, we solve Eq. (A.1) in the Fourier space. In this context, we assume the hypothesis that all involved functions admit a Fourier transform. So, we are looking for conditions that ensure the existence of some functions s₁, …, s_N such that:

$$\begin{array}{cccc}\widehat{f}=\widehat{g}& \text{where}& g\left(x\right):={\displaystyle \sum _{\ell =1}^N{\displaystyle \sum _{n=-\infty }^{+\infty }f\left(u_n^{\ell }\right)s_{\ell }\left(x-u_n^{\ell }\right)}}& \forall x\in ℝ.\end{array}$$

Forthe sampling grids that we consider, the Fouriertransformof g is given by:

$$\widehat{g}\left(k\right)=C_0\left(k\right)\widehat{f}\left(k\right)+{\displaystyle \sum _{n\ne 0,n=-\infty }^{+\infty }{C}_n\left(k\right)\widehat{f}\left(k-v_{h}n\right),}$$(A.3)

where

$$\begin{array}{ccc}{C}_n\left(k\right):=\frac{1}{h}{\displaystyle \sum _{\ell =1}^N{\widehat{s}}_{\ell }\left(k\right)e^{-2i\pi v_{h}n{u}_0^{\ell }}}& \text{and}& v_h:=\frac{1}{h}.\end{array}$$(A.4)

This computation is not direct, but classical, and it can be found, for instance, in Kohlenberg (1953).

In what follows, we suppose that f is band-limited, hence $\widehat{f}\left(k\right)$ vanishes outside some symmetric interval K. That is, there exists a frequency, k_max > 0, such that $\widehat{f}\left(k\right)=0$ for all k ∉ K, where K := (−k_max, k_max). In most practical situations f does not necessarily have a bounded support. But, $\widehat{f}\left(k\right)$ may decrease rapidly as |k| goes to infinity, as does, for instance, the power spectrum density of turbulent wavefronts following a Kolmogorov power law.

The translated copies $\widehat{f}\left(k-v_{h}n\right)$ of $\widehat{f}$, that appear in Eq. (A.3), may have some non zero values in the set K overlapping $\widehat{f}$ in this set, and hence producing aliasing. This happens if for some n ≠ 0 a set K_n := (−k_max + v_hn, k_max + v_hn) satisfies K ∩ K_n ≠ Ø. It is the case if and only if |n| < n₀, where

$$n_0:=\lceil \frac{2k_{\max }}{v_h}\rceil =\lceil \frac{k_{\max }}{v}\rceil ,$$(A.5)

and ⌈·⌉ is the ceiling function and v = 1/2h is the native Nyquist frequency. This can be reformulated in the following way: if for some integer n₀ the frequency k_max satisfies (n₀ − 1)v < k_max ≤ n₀v, then K ∩ K_n ≠ ∅, if and only if |n| < n₀.

It follows, that in K, but also in the bigger set:

$$I:=\left[-n_{0}v,n_{0}v\right]\supset K,$$

the formula (A.3) is expressed as:

$$\begin{array}{cc}\widehat{g}\left(k\right)=C_0\left(k\right)\widehat{f}\left(k\right)+{\displaystyle \sum _{n\ne 0,n=-n_0+1}^{n_0-1}{C}_n\left(k\right)\widehat{f}\left(k-v_{h}n\right)}& \forall k\in I,\end{array}$$(A.6)

since $\widehat{f}\left(k-v_{h}n\right)=0$ for every k ∈ I if |n| ≥ n₀, see Figure A.1.

Fig. A.1 Fourier transform $\widehat{f}$ (blue) and its translated copies (red). Here, 2v < kmax ≤ 3v and therefore n0 = 3. The translated copies $\widehat{f}\left(k-v_{h}n\right)$ overlap f for n ∈ {−2,−1,1,2} (dashed curves). On the interval I = [−3v, 3v], there is no contribution of $\widehat{f}\left(k-v_{h}n\right)$ for any |n| ≥ n0.

To ensure that $\widehat{g}=\widehat{f}$ in I, and therefore in K, we can impose for every k ∈ I that C₀(k) = 1 and C_n(k) = 0 for all

$$\begin{array}{l}n\in \left\{-n_0+1,\dots ,n_0-1\right\}\backslash \left\{0\right\},\text{\hspace{0.17em}\hspace{0.17em}that\hspace{0.17em}is,}\hfill \\ \begin{array}{cc}{\displaystyle \sum _{\ell =1}^N{\widehat{s}}_{\ell }\left(k\right)=h}& \forall k\in I,\end{array}\hfill \end{array}$$(A.7)

$$\begin{array}{ccc}{\displaystyle \sum _{\ell =1}^N{\widehat{s}}_{\ell }\left(k\right)e^{2i\pi v_{h}n{u}_0^{\ell }}=0}& \forall n\in \left\{-n_0+1,\dots ,n_0-1\right\}\backslash \left\{0\right\}& \forall k\in I.\end{array}$$(A.8)

Now, for k ∉ I, we have $\widehat{f}\left(k\right)=0$ and therefore we want $\widehat{g}\left(k\right)=0$. To cancel the contributions of all the translated copies of $\widehat{f}$ outside of K, we can simply set:

$$\begin{array}{cc}{\widehat{s}}_1\left(k\right)=\cdots ={\widehat{s}}_N\left(k\right)=0& \forall k\notin I.\end{array}$$(A.9)

Remark: Eq. (A.7) and Eq. (A.8) are sufficient conditions to obtain $\widehat{g}=\widehat{f}$ in the set I and to reconstruct f using formula (A.1). However, it may be possible to propose weaker but more elaborated conditions. For the sake of simplicity, we can restrict the study to Eqs. (A.7) and (A.8).

We can sum up our results and introduce more suitable notations for the forthcoming study of the system given by Eq. (A.7) and Eq. (A.8). We consider a band-limited function f, whose Fourier transform vanishes outside of (−k_max, k_max), and N sampling sequences of step h and initial positions $u_0^1,\dots u_0^N$ (see Eq. (a.2)). If we introduce the quantities:

$$\begin{array}{cc}{\beta }_{\ell }:=e^{2i\pi v_h{u}_0^{\ell }}=e^{2i\pi u_0^{\ell }/h}& \forall \ell \in \left\{1,\dots ,N\right\},\end{array}$$

then the system defined by Eq. (A.7) and Eq. (A.8) is of the form AX = Y, where X is a vector of size N, A is the (2n₀ − 1) × N matrix defined by:

$$\begin{array}{ccc}A:=\left(\begin{array}{cccc}{\beta }_1^{-n_0+1}& {\beta }_2^{-n_0+1}& \cdots & {\beta }_N^{-n_0+1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 1\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_1^{n_0-1}& {\beta }_2^{n_0-1}& \cdots & {\beta }_N^{n_0-1}\end{array}\right)& \text{and}& Y:=\left(\begin{array}{c}0\\ ⋮\\ h\\ ⋮\\ 0\end{array}\right).\end{array}$$(A.10)

If X = (x₁, …, x_N) is a solution of the system A.10, and ŝ₁ …, ŝ_N are defined for any ℓ ∈ {1, …, N} by:

$$\begin{array}{cccc}{\widehat{s}}_{\ell }\left(k\right)=x_{\ell }& \forall k\in I& \text{and}& {\widehat{s}}_{\ell }\left(k\right)=0\forall k\notin I,\end{array}$$(A.11)

then $\widehat{g}=\widehat{f}$. It follows that Eq. (A.1) holds, with the reconstruction functions:

$$\begin{array}{cc}{s}_{\ell }\left(x\right)=\frac{x_{\ell }{n}_0}{h}\text{sinc}\left(\frac{\pi n_{0}x}{h}\right)& \forall \ell \in \left\{1,\dots ,N\right\}.\end{array}$$(A.12)

We can offer a simple example and confirm the Nyquist–Shannon Theorem anew at the same time.

Example: By the Nyquist-Shannon Theorem, we know that, if k_max ≤ 1/2h = v, a sampling sequence of step, h, allows us to perfectly reconstruct f. We can re-prove that result. Indeed, if k_max ≤ 1/2h, then formula (A.5) gives n₀ = 1 and the system (A.10) to solve is expressed as:

$$x_1+\cdots +x_N=h.$$

If N ≥ 2, we can produce an infinite number of solutions of this equation, and we have as many possible choices for ŝ₁ …, ŝ_N.

However, a single sampling grid (N = 1), is enough to ensure the existence of a solution, which is x₁ = h. So, we let

$$\begin{array}{cccc}{\widehat{s}}_1\left(k\right)=h& \forall k\in I& \text{and}& {\widehat{s}}_1\left(k\right)=0\forall k\notin I,\end{array}$$

where I = [−1/2h, 1/2h]. By applying Eq. (A.12) and adopting Eq. (A.1), we obtain s₁(x) = sinc(πx/h) and

$$f\left(x\right)={\displaystyle \sum _{n=-\infty }^{+\infty }f\left(nh+u_0^1\right)\text{sinc}\left(\frac{\pi }{h}\left(x-nh-u_0^1\right)\right),}$$

as expected.

In the following, we study the conditions on N and $u_0^1,\dots u_0^N$ that ensure the existence of a solution to the system (A.10). As explained, if such a system has a solution, then we can reconstruct f using formula (A.1). If the system does not have solution, then there is no functions ŝ₁ …, ŝ₁, that satisfies Eq. (A.7) and Eq. (A.8). This does not necessarily mean that f cannot be reconstructed with formula (A.1), since Eq. (A.7) and Eq. (A.8) are sufficient, but not necessary, conditions to obtain $\widehat{g}=\widehat{f}$.

A.2 Optimality of the uniform sampling

The system (A.10) has 2n₀ − 1 equations with N unknowns. The number of equations is imposed by the highest frequency kmax of the signal f and by the step h of the sampling grids, through the formula (A.5). The number of unknowns is equal to the number N of sampling grids that we own. The coefficients of the system depend on h and can be tuned by our choice of the initial positions $u_0^1,\dots u_0^N$ of the sampling grids. Here, we are interested in the initial positions which ensure the existence of a solution to Eq. (A.10), with the smallest number of sampling grids N as possible.

A.2.1 Minimal number of sampling grids

If N = 2n₀ − 1, then h is a square matrix with a Vandermonde determinant. This determinant vanishes if and only if β_ℓ = β_l′ for some ℓ ≠ ℓ′. This occurs if and only if $u_0^{\ell }=u_0^{{\ell }^{\prime }}+qh$ for some q ∈ ℤ, that is, if the two grids ℓ and ℓ′ have the same points (see Eq. A.2). Therefore, if the sampling grids are all different, then the matrix A is invertible. In such a case, the system (A.10) has a solution that is unique for each possible position of $u_0^1,\dots {\mathrm{,u}}_0^N$. We can then reconstruct our signal choosing the functions ŝ₁ …, ŝ_N according to Eq. (A.11). We conclude that with N = 2n0 - 1 grids, we can reconstruct a signal such that (n₀ − 1)v < k_mαx ≤ n₀v, even if the initial positions of the grids are chosen randomly.

However, for aptly chosen positions of $u_0^1,\dots ,u_0^N$, we can reduce the number of independent equations of the system and obtain a solution even if N < 2n₀ − 1 . Indeed, we may suppose that we have N = n₀ sampling grids of the form:

$$\begin{array}{cc}{u}_n^{\ell }=nh+u_0^1+\left(\ell -1\right)\frac{h}{N}& \forall \ell \in \left\{1,\dots ,N\right\}.\end{array}$$(A.13)

Then, for every ℓ ∈ {1, …, N} and j ∈ {−n₀ + 1, …, −1}, we have

$${\beta }_{\ell }^{j+n_0}=e^{2i\pi \left(j+N\right)u_0^{\ell }/h}={\beta }_{\ell }^j{e}^{2i\pi N{u}_0^1/h}{e}^{2i\pi \left(I-1\right)}={\beta }_1^N{\beta }_{\ell }^j.$$

This implies that to multiply the n₀ − 1 first lines of AX = Y by ${\beta }_1^N$ gives a system where each of the n₀ − 1 first equations is equal to one of the n₀ − 1 last equations. Therefore, if X = (x₁, …, x_N) is a solution for:

$$\left(\begin{array}{cccc}1& 1& \cdots & 1\\ {\beta }_1& {\beta }_2& \cdots & {\beta }_N\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_1^{n_0-1}& {\beta }_2^{n_0-1}& \cdots & {\beta }_N^{n_0-1}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{c}h\\ 0\\ ⋮\\ 0\end{array}\right),$$(A.14)

then X is also a solution of AX = Y. The square system (A.14) has a Vandermonde determinant. Since the points of the grids (A.13) are all different, the system (A.14) is invertible. We can verify that its unique solution is given by x₁ = x₂ = ⋯ = x_N = h/N. So, we re-prove that with N = n₀ grids chosen according to Eq. (A.13), we can reconstruct a signal such that (n₀ − 1)v < k_max ≤ n₀v.

Now, we may wonder if there is some smarter way to select the positions of the sampling grids, which would allow us to further reduce their number. The answer is negative. Indeed, if N ≤ n₀ − 1, we can extract, from Eq. (A. 10), the subsystem:

$$\left(\begin{array}{cccc}{\beta }_1& {\beta }_2& \cdots & {\beta }_N\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_1^N& {\beta }_2^N& \cdots & {\beta }_N^N\end{array}\right)\left(\begin{array}{c}{x}_1\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{c}0\\ ⋮\\ 0\end{array}\right),$$(A.15)

which must be satisned by any solution of Eq. (A.10). as Eq. (A.15) is an invertible homogeneous system, its unique solution is X = (0, 0, …, 0), which is not a solution of Eq. (A.10). Thus, Eq. (A.10) has no solution if N ≤ n₀ − 1.

A.2.2 Unicity

We have shown that N = n₀ is the minimal number of sampling grids for the existence of a solution to the system (A.10). In particular, if we choose our grids according to Eq. (A.13), then (A.10) has a solution. At this point, we can verify that for any other choice of the position of n₀ grids, the system (A.10) has no solution.

To that aim, we look for a necessary condition for the rank of A to be equal to the rank of A|Y, where A|Y is the matrix A with an extra row equal to Y. First, we have rank(A) = N, since, on the one hand, the number of rows N = n₀ of A is smaller than its number of lines 2n₀ − 1 ; and, on the other hand, A contains a N × N Vandermonde invertible matrix. Thus, it is enough to look for a necessary condition for rank(A|Y) = N.

As N = n₀, A|Y is a (2N − 1) × (N + 1) matrix, from which we can extract a collection of N − 1 matrices B₁, …, B_N−1 of size (N + 1) × (N + 1) and defined by:

$$\begin{array}{cc}{B}_i:=\left(\begin{array}{ccccc}{\beta }_1^{-i}& {\beta }_2^{-i}& \cdots & {\beta }_N^{-i}& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 1& \cdots & 1& h\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\beta }_1^{N-i}& {\beta }_2^{N-i}& \cdots & {\beta }_N^{N-i}& 0\end{array}\right)& \forall i\in \left\{1,\dots ,N-1\right\}.\end{array}$$

We can have rank(A|Y) = N, only if the det(B_i) = 0 for all i ∈ {1, …, N − 1}. A Laplace expansion along the last row gives:

$$\det \left(B_i\right)={\left(-1\right)}^{N+i+2}h\det \left(C_i\right){\displaystyle \prod _{\ell =1}^N{\beta }_{\ell }^{-i}}$$

where

$$\begin{array}{cc}{C}_i:=\left(\begin{array}{cccc}1& 1& \cdots & 1\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_1^{i-1}& {\beta }_2^{i-1}& \cdots & {\beta }_N^{i-1}\\ {\beta }_1^{i+1}& {\beta }_2^{i+1}& \cdots & {\beta }_N^{i+1}\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_1^N& {\beta }_2^N& \cdots & {\beta }_N^N\end{array}\right)& \forall i\in \left\{1,\dots N-1\right\}.\end{array}$$

Let P(X) = a_NX^N + ⋯ + a_iXⁱ + ⋯ + a₀ be the polynomial defined by

$$P\left(X\right):=\left|\begin{array}{ccccc}1& 1& \cdots & 1& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\beta }_1^i& {\beta }_2^i& \cdots & {\beta }_N^i& X^i\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\beta }_1^N& {\beta }_2^N& \cdots & {\beta }_N^N& X^N\end{array}\right|.$$

Then, a Laplace expansion along the last row gives a_i = (−1)^i+N+1 det(C_i) for all i ∈ {1, …, N − 1}. Therefore, we have det(B_i) = 0 for all i ∈ {1, …, N − 1} if and only if

$$P\left(X\right)=a_N{X}^N+a_0.$$

The coefficients a_N and a₀ are obtained by Laplace expansion, and if we denote V_N(β₁, …,β_N) the Vandermonde determinant of size N, then:

$$P\left(X\right)=V_N\left({\beta }_1,\dots ,{\beta }_N\right)\left(X^N+{\left(-1\right)}^{N+2}{\displaystyle \prod _{\ell =1}^N{\beta }_{\ell }}\right).$$

On the other hand, we have P(β_ℓ) = 0 for all ℓ ∈ {1, …, N}, which gives:

$${\beta }_1^N={\beta }_2^N=\cdots ={\beta }_N^N={\left(-1\right)}^{N+1}{\displaystyle \prod _{\ell =1}^N{\beta }_{\ell }.}$$

It follows that (β_ℓ/β₁)^N = 1 for all ℓ ∈ {1, …, N}. In other words, each is a N^th root of the unity. Since all the β_ℓ are supposed different, we have

$$\begin{array}{cc}{\beta }_{\ell }={\beta }_1{e}^{2i\pi \left(\ell -1\right)/N}=e^{2i\pi \left(u_0^1/h+\left(\ell -1\right)/N\right)}& \forall \ell \in \left\{1,\dots ,N\right\}\end{array}.$$

Recalling that ${\beta }_{\ell }=e^{2i\pi u_0^{\ell }/h}$, we obtain

$$\begin{array}{ccc}{u}_0^{\ell }=u_0^1+\left(\ell -1\right)\frac{h}{N}& \mathrm{mod}h& \forall \ell \in \left\{1,\dots ,N\right\}.\end{array}$$

We conclude that, when N = n₀, the sampling grids must satisfy Eq. (A.13) to ensure the existence of a solution of the system (A.10).

A.2.3 Conclusion

To sum up our findings, with N = 2n₀ − 1 = 2 ⌈k_max/v⌉ − 1 sampling grids of step h, we can always reconstruct a signal, f, with a maximal frequency of kmax. The sampling grid positions can be chosen at random, provided that the grids are all different. However, the same performance can be reached with only N = n₀ sampling grids, if their union is a uniform sampling grid of step h/N, that is, if they satisfy Eq. (A.13) modulo a translation o a size h. On the other hand, if we have fewer than n₀ sampling grids, or n₀ sampling grids that do not satisfy Eq. (A.13), then there is no set of reconstruction functions s₁, …, s_N which satisfy Eq. (A.7) and Eq. (A.8). Since these conditions might be too strong, we cannot draw a firm conclusion that the signal cannot be reconstructed, since it is possible that more elaborate approaches and reconstruction functions could offer a more promising result. Our study suggests nevertheless that the meta-uniform sampling Eq. (A.13) is optimal, because it allows us to reconstruct the signal without aliasing with the smallest number of sampling grids as possible. It is therefore the safest way to make use of a limited number of sampling grids.

Appendix B Non-uniform sampling schemes

We have shown that both uniform and non-uniform sampling schemes allow for the reconstruction of a band-limited signal without aliasing, under some conditions presented in this paper. In addition to offsets and rotations, the AO design offers several other sampling schemes that facilitate SR. A non-exhaustive list of recommendations is as follows:

• differently rotate and offset (possibly stretch) the WFS sampling arrays by non-integer multiples of each other (Figure B.1);

• use asymmetric guide star asterisms projected on-sky (Figure B.2);

• use WFSs with different numbers and sizes of sub-apertures (Figure B.3);

• allowing the sub-aperture sizes in each lenslet array to vary.

In either case, the unique phase information collected by each WFS can be used to separate the aliased high-frequency from the low-frequency content of interest and the higher-resolution wavefront can be accurately reconstructed.

Fig. B.1 Six differentially rotated WFS configurations shown at three different arbitrary conjugation altitudes creating growing separations (highly exaggerated for illustration only).

Fig. B.2 Six regular WFS configurations with a spiral shaped asterism shown at three different conjugation altitudes creating growing separations.

Fig. B.3 Six WFS configurations with two kinds of sub-aperture, one being the double in size of the other, shown at three different conjugation altitudes creating growing separations.

References

All Tables

All Figures

	Fig. 1 Illustration of four lenslet arrays, aligned to the left and shifted by 1/2 of a sub-aperture to the right, where the four colours used (red, green, blue, and black) are no longer superimposed.
In the text

Fig. 2 Normalised sensitivity of a virtual averaging-and-sampler sensor used for illustration purposes only. The black curve corresponds to the sensitivity function of a sampling grid of step h/2. The dotted curve represents the sensitivity of a sampling grid of step h. The red curve is obtained by dealiasing the spectral sensitivity up to 1/h, that is, the first null of the sinc function with a super-resolution factor of 2; by combining two sampling grids of step h relatively shifted by h/2 with respect to each other.

In the text

	Fig. 3 Fourier S/N functions. Top: Fourier S/N for a diffraction-limited wavefront sensing case (Eq. (6)). Bottom: seeing-limited or laser spot-sensing case. The number of incident photons is considered the same for all configurations. The bottom panel assumes a spot size with one unit.
In the text

	Fig. 4 Noise propagation coefficient as a function of spatial frequency. This plots are computed as the inverse of the S/N functions from Fig. 3. The fluxes are scaled by the collecting area, namely, the square of the sub-aperture width.
In the text

	Fig. 5 SR configurations illustrated via a projection in the pupil plane of one sub-aperture from each WFS (1 color for each): (a) 2D surface is improperly paved with shifts along the sampling grid’s axes; (b) 2D surface is properly paved with the addition of shifts at 45° with respect to the sampling grid’s axes, as in Eq. (13).
In the text

	Fig. 6 Singular value decomposition of a multi-WFS mono layer modal IM. The singular values are plotted as a function of the singular modes index. The solid lines corresponds to 4 co-aligned SH-WFS while the dashed lines represent the shifted configurations shown in Fig. 5b.
In the text

Fig. 7 Open-loop reconstruction of an incident phase corresponding to a Kolmogorov turbulence with r0 = 12.6 cm. The dashed lines represent the co-aligned cases. The solid lines stand for the super-resolved cases achieved by combining four SH WFS that are properly shifted by half (1/2) of a sub-aperture. The same number of modes is reconstructed for the co-aligned cases and their respective shifted cases, that is, up to $4.\frac{\pi }{4}{n}_{\text{SA}}^2$ modes, which is the spatial resolution that would be provided by a single WFS with twice as many sub-apertures per diameter.

In the text

Fig. 8 Residual wavefront error from open-loop reconstruction as a function of shift amplitude for the shifted and super-resolved configurations. In each case, up to four times the native number of modes $\left(4.\frac{\pi }{4}{n}_{\text{SA}}^2\right)$ are reconstructed. The dashed black lines represent the Kadec condition boundaries within which the non-uniform samples can be perfectly reconstructed. The colored squares stand for the reconstruction error from the co-aligned case when $\frac{\pi }{4}{n}_{\text{SA}}^2$ modes are reconstructed.

In the text

Fig. 9 Noise propagation coefficients, computed as diag{RRT}, shown in logarithmic scale and scaled by the flux collecting surface of the sub-aperture. The shifted cases are plotted with solid lines. For the sake of visibility, dashed lines are overlayed at the level of the spatial frequencies where the aliasing starts being prominent for the co-aligned cases. Along the low order modes, the noise propagation is quasi identical for all cases and fit the expected law in (n + 1)−1. Then, when exceeding the original Nyquist–Shannon cut-off frequency, the noise propagation coefficient converges towards a constant asymptotic limit.

In the text

	Fig. 10 SVD decomposition on a 16 × 16 P-WFS system with a 33 × 33 DM. Top: intensity map processing. Bottom: slopes-maps processing. We clearly see that using the intensity-maps on an offset pixel sampling system leads to sensitivity gains, whereas the same is not true if the slopes-maps were to be processed.
In the text

	Fig. A.1 Fourier transform $\widehat{f}$ (blue) and its translated copies (red). Here, 2v < kmax ≤ 3v and therefore n0 = 3. The translated copies $\widehat{f}\left(k-v_{h}n\right)$ overlap f for n ∈ {−2,−1,1,2} (dashed curves). On the interval I = [−3v, 3v], there is no contribution of $\widehat{f}\left(k-v_{h}n\right)$ for any |n| ≥ n0.
In the text

	Fig. B.1 Six differentially rotated WFS configurations shown at three different arbitrary conjugation altitudes creating growing separations (highly exaggerated for illustration only).
In the text

	Fig. B.2 Six regular WFS configurations with a spiral shaped asterism shown at three different conjugation altitudes creating growing separations.
In the text

	Fig. B.3 Six WFS configurations with two kinds of sub-aperture, one being the double in size of the other, shown at three different conjugation altitudes creating growing separations.
In the text