Alias-Free Basis for Modal Sensorless Adaptive Optics Using the Second Moment of Intensity
Oleg Soloviev
physics.app-ph


Oleg SolovievAlias-free basis for sensorless adaptive optics
















Alias-Free Basis for Modal Sensorless Adaptive Optics Using the Second Moment of Intensity

    Oleg Soloviev[DCSC, 3mE, Delft University of Technology,
		Mekelweg 2, 2628 CD Delft, the Netherands,

			o.a.soloviev@tudelft.nl] 

    December 30, 2023
========================================================================================================================================










In theory of optical aberrations, an aberrated wavefront is represented by its coefficients in some orthogonal basis, for instance by Zernike polynomials.
However, many wavefront measurement techniques implicitly approximate the gradient of the wavefront by the gradients of the basis functions.
For a finite number of approximation terms, the transition from a basis to its gradient might introduce an aliasing error.
To simplify the measurements, another set of  functions, an “optimal basis” with orthogonal gradients, is often introduced, for instance Lukosz-Braat polynomials.
The article first shows that such  bases do not necessarily eliminate the aliasing error and secondly considers the problem of finding an alias-free basis on example of second-moment based indirect wavefront sensing methods.
It demonstrates that for these methods any alias-free basis should be formed by functions simultaneously orthogonal in two dot-products and be composed of the eigenfunctions of the Laplace operator.

The fitness of such alias-free basis for optical applications is analysed by means of numerical simulations on typical aberrations occurring in microscopy and astronomy.







§ INTRODUCTION


In   optics, the quality of an optical system can be characterised by the rms of the  wavefront aberration.
Thus, in adaptive optics, the rms of the residual wavefront error is the measure of the correction:

    ϵ = ∫_A (ϕ() - φ_DM())^2  ,

where ϕ and φ_DM denote the incoming and correction wavefronts, and  is  the coordinate vector inside  aperture A.
On the other hand, wavefront sensors often provide information based on the gradient of the wavefronts, and control algorithms  use a quantity related (maybe implicitly) to the rms of the gradient of residual aberration:

    ϵ_∇ = ∫_A ∇ϕ() - ∇φ_DM()^2  .

While without limitation on φ_DM, minimisation of both errors  gives equivalent up to a constant phase term results, this is not the case when φ_DM is confined to some subspace, for instance, spanned by the response functions of a deformable mirror.


This article considers an abstract mathematical problem of finding such a basis that minimisation of both error given by <ref> over any subspace spanned by a finite subset of the basis functions provides equivalent results and discusses the advantages of using such bases in optical applications on example of indirect wavefront sensing method based on the second moment of intensity.
The method is briefly recapitulated in the following section.





§ INDIRECT WAVEFRONT SENSING




In feedback adaptive optics (AO) systems, the goal of wavefront correction is to minimise the rms of the residual aberration <cit.>,
which is measured by means of the wavefront sensing.
In applications with (quasi-)static aberrations, such as microscopy, indirect (also referred to as wavefront sensor-less) methods provide some advantage in comparison with direct wavefront sensing and at present can be considered as a common approach to aberration sensing and correction in AO microscopy <cit.>.
In these methods, the aberration is corrected by optimising some image quality metric as a function of the control signals applied to an adaptive element.
The chosen metric should  of course be related to the rms value of the residual aberration, so that optimisation of the metrics should minimise the wavefront error. 

In model-based indirect methods,  the correction wavefront is modelled as a linear  function of some parameters c = (c_1,…,c_N)^T, where N is the number of the degrees of freedom in the model;
the parameters  c are often (but not necessarily) linearly related to the control signals v of the adaptive element: v = Lc.
In this case one effectively controls the adaptive element via modes formed from the response functions by columns of matrix L.

If the image quality  metric J attains its optimum value J_o at some optimal point c_o = (c_1,o,…,c_N,o)^T,
it  can be considered as quadratic in the vicinity of c_o: 

    J(c_1,…,c_N) ≈ J_o ±∑_i,j =1^N m_ij (c_i-c_i,o)(c_j - c_j,o),   providedc-c_ois small,
 

where m_ij are the elements of some positive semidefinite matrix M≽ 0.
In this case, the metric can be optimised using one or more iterations consisting of P N +1 independent images each <cit.>, where P is some integer number, provided the coefficients m_ij of the quadratic form are known.

The structure of matrix M is influenced by the choice of the metric and 
the choice of the modes (that is by defining matrix L).

The particular case of the  diagonal  matrix is obviously the most advantageous, as it allows to adjust each of the modes independently, which increases the accuracy of the wavefront reconstruction and decreases the number of required iterations; it is achieved by careful selection of the metric function, or by selection of the modes that diagonalises the matrix m_ij:

    M = VSV^T,


where S is some positive semidefinite diagonal matrix and the rows of V define the new modes.
The shape and properties of V depend on the orthogonalisation procedure used: 
for instance, Gram-Schmidt orthogonalisation returns a  lower-triangular V; and in SVD  V is an orthogonal matrix formed by the eigenvectors of M, with S is formed by the eigenvalues of M.

The modal functions that diagonalise the metric J  are  usually referred to as “optimal” <cit.> and can be obtained analytically, numerically  or empirically <cit.>. 

A special class of metrics is based on measuring the change of the second moment (SM) of the image.
This method is applicable both to an image of a point source <cit.> and of an extended object <cit.>. This metric can be expressed via the squared gradient of the residual aberration; consequently it is quadratic also for larger aberrations and can be optimised in N+1 steps.




It has been shown that the Lukosz-Braat (LB) polynomials are the optimal modes for the second-moment methods <cit.>.
The LB polynomials are obtained from Zernike polynomials by Gram-Schmidt  orthogonalisation of their gradients <cit.> and thus minimise the transverse ray aberration (the second moment of a PSF). This is an example of an analytically obtained optimal basis.

For a low-order deformable mirror, the LB polynomials cannot be approximated well enough with the mirror response functions, and  the optimal modes are obtained by SVD orthogonalisation of their gradients <cit.>, or by SVD decomposition of an empirically obtained matrix M using <ref>. Modes obtained in such a way can be considered as numerical and empirical optimal bases respectively; this approach is more often used in practice. It was shown <cit.> that use of SVD modes is more advantageous over approximated LB modes for an adaptive element with a small number of actuators.

The subject of this paper is the derivation of an 
analytical 
optimal basis for the case when the number of degrees of freedom of an adaptive element is very large, such as for SLM or photo-addressed DM <cit.>. 
By “very large” or “practically infinite” we assume that the adaptive element is able to represent any analytical mode of interest accurately enough. 
We can mark the following differences with the case of a small number of degrees of freedom.

Firstly, because of the increased dimensionality, it becomes unpractical to find  the optimal modes via SVD  and they should be determined analytically — one needs to perform N^2 measurements to empirically obtain matrix M, and the results of this calibration procedure would be prone to errors, as a higher density of actuators usually results in decreased amplitude of individual response functions.


Secondly, unlike the small-dimensional case, where all optimal bases are related to each other via a (small-dimensional) isometric transformation <cit.>, in a large-dimensional case, the corresponding transformation would result in a computation of large series, prone to the numerical errors; and anyway this procedure is effectively equivalent to just defining  new modes point-wise.

At last, while in the case of a small number of actuators one can use all of them to obtain the best correction that can be achieved, in big-dimensional situations one needs to limit effective number of degrees of freedom to some smaller value  N' 
in order to avoid overexposure of the specimen;
in general, one wants to keep N' as small as possible. 
It is natural to think that  any incoming aberration ϕ  can be decomposed as a sum of two  functions, one from the space that an adaptive element (or only those its modes that are used) can correct and the residual which cannot be corrected, so one can  simply ignore it:

    ϕ = ∑_n=1^N' f_n +∑_n>N' f_n = φ_DM + φ_r,

where, to make the decomposition unique, the function φ_r is chosen in the least square-sense approximation, which is equivalent to the condition that φ_DM and φ_r are mutually orthogonal.
Then, following this reasoning, matrix M is obtained and diagonalised only for the first N' basis terms.
However, it is often forgotten that there are functions that cannot be corrected but which can still  be “seen” by the metric and, consequently, there are functions with a non-zero correctable part φ_DM that will not be seen by the method and remain under-corrected.
In this situation, the accuracy of the correction using only the first N' modes would not only depend on the statistical properties of the aberration, but also on the way the optimal basis was obtained. 
For instance, LB modes, obtained by sequential orthogonalisation of the gradients of the Zernike modes by Gram-Schmidt algorithm, introduce a systematic error  as is shown in Section <ref>.

The rest of the article is organised as follows. 
First, we show on an example the problem of aliasing if suboptimal modes are used. 
Second, we show that second-moment metric infers a new dot product in the space of modes and show that diagonalisation of M is equivalent to orthogonalisation of the modal basis. 
Then  we derive the condition for the aliasing-free basis. 
Finally, we present and discuss the results of numerical simulations using optimal modes on examples of typical aberration in microscopy and in astronomy.




§ ALIASING ERROR IN THE SECOND-MOMENT METHODS


In this section we demonstrate how use of the first terms of a basis, although optimal in the sense of diagonalisation of matrix M, can result in a  systematic aliasing error. 
Although the aliasing or coupling error is familiar to all using adaptive optics with direct or indirect wavefront sensing, we dwell a little bit on it to demonstrate that the error is caused completely by the choice of basis and can be avoided by using an “even more optimal” basis.
We name the bases correspondingly as suboptimal and optimal (or aliasing-free) bases.

In the examples, we use indirect wavefront sensing based on  any image metric related to the computation of the second moment of the intensity in the focal plane. 
The concise description of such a metric for images of a point source is given in the next section.
The same reasoning, however, can be applied and similar examples can be constructed for other metrics.




 §.§ Second-moment based wavefront sensor-less method

For an image of a point source obtained with control signals v,  let metric J(v) be defined as a second moment of its normalised intensity  I(x,y) = I(x,y|v):

    J(v) = ∑_x,y I(x,y) (x^2 + y^2)/∑_x,y I(x,y) ,

where the origin (x,y) = 0 is given by the intersection of the optical axis with the imaging plane.

Under some natural limitations, the sums in the equation above can be considered as discrete approximations of two-dimensional integrals over ℝ^2, and by denoting as φ = φ(v) and ϕ the phase  generated by the adaptive element and  that of the initial aberration, in the case of uniform illumination, metric J can be shown <cit.> to be linear dependent on the averaged square gradient of the phase of light field in the aperture:

    J(v) ≈ J_d + 
    
    1/4 π^ 2 ∫_^2P^2()∇ (ϕ- φ(v))^2 /∫_^2 P^2() 
    ,

where = (ξ, η) is the coordinate vector in the pupil plane, P() is the pupil function, and J_d is the second moment of the diffraction-limited PSF.

Assuming linear dependence of the correction phase on the control signals

    φ (v) = ∑_n=1^N f_n() v_n,

where f_n(),  n =1, …, N are the response functions of the adaptive element,  expanding the modulus square in the right hand side of  <ref> and completing the full square, the metric J can now be written in the form of <ref> as 

    J(v) ≈ J_d + m_o + (v -v_o)^T M (v -v_o),

where m_o ≥ 0 is some constant dependent on φ, such that the minimum value J_o of <ref> is given by  J_o = J_d + m_o, and the elements m_i,j of the matrix M are given by the averaged inner product of the gradients of the response functions f_i and f_j:

    m_i,j = 1/4 π^ 2 ∫_^2P^2()  (∇ f_i, ∇ f_j)  /∫_^2 P^2() .

If all  elements m_i,j are known (either from the calibration procedure or by diagonalisation of M in a theoretical, numerical, or empirical way), there remain N+1 unknowns — m_o, v_o — in <ref>; they can be extracted using K ≥ N+1 measurements  J_k = J(v_k),  k=1,…, K of the metric in K test points v_k.
Thus, in the case of diagonal M, the test points c_k,  k=0,…, N can be chosen as 

    c_k = (0,…, 1,…, 0),

where 1 is present only on the k-th position.
The calculation of the optimal point c_o from measurements J_k is then straightforward.




 §.§ Example of aliasing in suboptimal basis





Consider the following  hypothetical example of wavefront correction with a deformable mirror able to reproduce any Zernike polynomial, and suppose we want to limit ourselves only to the correction of the first 6 polynomials[usually in microscopy, tip, tilt, and defocus (and of course piston) are discarded from the modal space as they correspond to the translation of the image in 3D space; here we kept them for the sake of simplicity—the example can be easily modified to include any other first N' Zernike modes as the modal space and a higher order mode as the aberration] using the SM-based method, that is, using metric J given by <ref>.
Suppose also that we can calculate without error the second moment of PSF, i.e. there is no noise, and the domain of (x,y) (the camera size) is infinite.

Lukosz-Braat (LB) polynomials diagonalise the matrix given by <ref> for a circular aperture for any N, so m_i,j = δ_ij,  i,j = 1,…,∞, where δ _ij is the Kronecker symbol; consequently they are often used as a theoretical basis for SM-based methods. 

The first 6 LB polynomials span
the modal space F   given by

    F= ⟨ 
    1, x, y, x^2, xy, y^2
    
    ⟩,

where ⟨·⟩ stands for the linear span.
This is  the same space as spanned by the first 6 Zernike polynomials, 
so any aberration that belongs to this space would be corrected by the mirror and SM method with zero error.





In our example, let's take the incoming wavefront ϕ outside of space F, for instance a LB polynomial of a higher degree  B_4^0 

    ϕ = B_4^0 =  6 r^4 - 8  r^2  +2.
 
As LB basis also diagonalises the metric J for any number of the first  terms, the value of metric J will not change when applying any test aberration from the first 6 modes.
Thus, the second moment methods will not see the aberration and leave it without change.

On the other hand, in the sense of <ref>, ϕ contains a non-zero part φ_DM = -2 r^2 +1 that should be corrected by the mirror. 
Thus, B_4^0 is not seen by the SM method and will be under-corrected. 

Symmetrically,
the spherical aberration Z_4^0 from the space {φ_r} of  uncorrectable terms of <ref> can be expressed as

    ϕ' = 6 r^4 - 6 r^2 +1 = B_4^0 +(2r^2 -1),

will be seen by the SM method as the defocus term Z_2^0 (see <ref>):

    ϕ≈φ' = 2 r^2-1.
 and will be over-corrected by applying  the corresponding correction.








This example can be considered as an illustration of aliasing, a manifestation of a “high-order” term Z_4^0 not belonging to the modal space F, as a non-zero “low-order”  term Z_2^0 ∈ F, or equivalently, the presence of a non-zero correctable term in a aberration “not-seen” by the metric. 

In other words, the aliasing here was caused by the choice of space F of the first N' mirror responses, that was not orthogonal  to a high-order function from the basis that diagonalises the metric (B_4^0).

From the dot-product matrix  of the LB polynomials (see second line of Fig. <ref>), it is clear that any subspace F formed by a finite number of LB modes (not necessarily the first ones) will have an infinite number of non-orthogonal to it LB modes.

Moreover, the aliasing will remain for any other orthogonal basis in F (e.g. a basis obtained by SVD of the first N' modes with tip, tilt, and defocus excluded), as choice of basis does not change the space itself. Thus we can conclude that the nature of Zernike (or equivalently LB) polynomials does not allow the use of the decomposition of <ref> without introducing an aliasing error. 

From our example, with the space F  formed by the first N' modes of the LB basis, it is clear that the necessary condition for a basis diagonalising SM metric J to be  aliasing-free   is that the functions from  its tail should be orthogonal to all the functions from its start. 
This can be generalised further to conclude that bases that diagonalise the metric but of which the dot-product matrix is not block-diagonal are not aliasing-free and are thus sub-optimal. 

Therefore for use in indirect wavefront sensing methods, it would be more appropriate to call a basis “optimal” if it a) is  formed by a set of orthonormal functions and b) diagonalises the metric.
In this redefined optimal basis, decreasing the number of decomposition terms will only introduce a corresponding truncation error caused by the non-zero norm of φ_r; in a sub-optimal basis, an  additional aliasing error will also be present. 

The natural question is whether a defined in such a way optimal, aliasing-free, basis does exist for a given metric. 
The next section provides an example of a positive answer.





 §.§ Example of an aliasing-free basis


It's relatively easy to find an example of an optimal basis for the second-moment metric and a square aperture. This is provided by the trigonometrical basis, defined on a unit square (x,y) ∈ [0,1] × [0,1] as

    f_k l () = ^2 π (k x + l y),    k, l ∈ℤ,

or its real-valued equivalent {sin (k x + l y), cos( k x + l y)}.
The basis is orthogonal  on the unit square, and it is easy to check that it diagonalises the matrix given by <ref> as its gradients are also orthogonal on the unit square, so the truncated Fourier series of the decomposition of ϕ in this basis 

can be used for alias-free indirect wavefront sensing.



§ DERIVATION OF AN ALIASING-FREE BASIS FOR ANY APERTURE


Now we're going to find an optimal basis for the SM metric for  a circular aperture, for which
both  the modes  and their gradients should form two orthogonal sets on a unit disk.
Moreover, we find a sufficient and necessary condition for a basis to be optimal

for second-moment based methods 
on any aperture.
For this, consider the following mathematical formalisation of the problem.



Consider the Hilbert space L = L_2[P], where P is the phase domain, that is, the pupil of the system, with the “standard” L_2 inner product defined as

    (f,g) = ∫_P f() g() ,

and let  F ⊂ L be a subspace, spanned by some N' functions or modes f_n(),  n=1,…, N' 


    F = ⟨ f_1, …, f_N'⟩.

Let ϕ∈ L be some function[for the sake of brevity, from now on the dependence on  will be omitted if it is not causing any problems, and tensor notation will be used, so a variable with an upper letter index c^n denotes a column vector c^n = [c^1, …, c^N']^T, and a variable with a lower index denotes a row vector; a repeating upper and low index denotes summation]
that is approximated by an  element φ of F:

    ϕ() ≈  c^0 + φ() = c^0 + ∑_n=1^N' c^n f_n () ≡ c^0 + c^n f_n,

where the piston term c^o is usually separated as it has no influence on the PSF.  



With ϕ representing the incoming wavefront, and f_1,…,f_N' representing the decomposition modes, the 
modal wavefront reconstruction 
may now be represented as an optimisation problem in the space L_2[P]  

    min_c^0, c^nϕ - c^0 - c^n f_n _2^2.



Without loss of generality, we can consider the set  f_n to be the first N' functions of a full and orthonormal basis f_n, n=1,…, ∞,

    (f_n, f_n') = δ_n,n',
 
where  δ is Kronecker's symbol.
Then the best approximation φ  of ϕ is given by truncation of its  decomposition by the basis functions

    ϕ = ∑_n=1^∞ c^n f_n  ⇒φ = ∑_n=1^N' c^n f_n,

where the decomposition coefficients c^n are obtained by the inner product of φ with the basis functions:

    c^n = (f^n,ϕ),

and are known as the (generalised) Fourier coefficients. 
From the orthonormality of the basis, we get the Euclidean norm of the approximation error ε = ϕ - φ as

    ε^2= ϕ^2 -φ^2 = ∑_n=N'+1^∞ (c_n)^2 = ϕ^2 -  ∑_n=1^N' (c_n)^2  .


From <ref> it follows that the second-moment-based methods approximate the gradient of the incoming wavefront with the gradients of the basis function, solving thus another optimisation problem

    min_ c^n∇ϕ  - c^n ∇ f_n _2^2.

Here, the piston term c_0 is cancelled by the gradient operator ∇.

Disregarding the piston term, the optimal (alias-free) basis should provide the same solutions for the problems <ref> and <ref>.




Let us remove the piston term from both problems by redefining L as a factor space by the gradient kernel: L=L_2[P]/{ϕ:ϕ=const}.
In the factor space L, the optimisation problem <ref> is equivalent to the decomposition of  L  into the direct sum of F and its orthogonal complementation O= F^⊥:

    L = F ⊕ O,
 
where  φ∈ F and ε∈ O are the orthogonal projections of ϕ on the modal space and its orthogonal complementation, and

    ϕ^2 = φ^2 + ε^2.





Let us also introduce another inner product in L, defined as the “standard” inner product[the inner product (∇ f, ∇ g) is defined here for the vector-functions from L× L in a natural way: (∇ f, ∇ g) =  (∂ f/∂ x, ∂ g/∂ x) + (∂ f/∂ y, ∂ g/∂ y)]
of the gradients:

    (f,g)_∇def=  (∇ f, ∇ g) =
    ∫_P ∇ f()·∇ g() ,

where a · b denotes the “standard” inner product in ℝ^2. One can easily verify that in the factor space L, relation (<ref>) satisfies all properties of an inner product,  which we will call the “gradient-dot product”, and therefore we can define a corresponding “gradient-norm” as

    f_∇ = √((f,f)_∇)≡√(∫_P ∇ f^2 ).


With the gradient-norm, the minimisation problem of <ref> can be written as 

    min_c^nϕ - c^n f_n_∇,

and solving this problem is equivalent to finding the orthogonal complementation of F  with respect to the gradient-dot product: 

    L = F ⊕_∇ O_∇.
 



Obviously, as <ref> and <ref> use different norms to minimise the residue, their solutions do not necessarily  match. 
It is easy to prove (see Fig. <ref>) that they would be the same if and only if the two orthogonal complementations of F are the same, that is 
    O=O_∇ .





Indeed, if for any ϕ its projections on F in both inner products match, i.e. φ  = φ ', then consequently its approximation errors match too, i.e. ε = ε'. Then for any vector x∈ O_∇, considering its approximation errors, one has  x= ε ' = ε∈ O, and thus O_∇⊆ O. In the same way,  O ⊆ O_∇ and thus O_∇ = O.

Vice versa, suppose that O_∇ = O. Then from Eqs.(<ref>) and (<ref>) and from the definition of the direct sum, we get that φ = φ '.

This explains why optimisation of the Zernike basis by changing to LB or SVD modes may introduce an aliasing error.
By performing this optimisation, we re-orthogonalise the basis modes in F with respect to the dot-gradient product and thus also implicitly introduce a new inner product, which changes the orthogonal complementation from O to O_∇.
For instance, no first N terms of the Zernike modes considered in the example above can form an optimal basis for the second-moment-based moments. 
Noll <cit.>  has calculated the coefficients of decomposition of the x- and y-components of the gradient of Zernike polynomials through themselves, and

from his calculation it follows that

    (Z_n^0,Z_n+2^0)_∇≠ 0  ∀ n,  n  is even.

Hence, if n is the maximal order of spherical aberration from F, then Z_n+2^0 ∉ O_∇. But as Z_n+2^0 ∈ O, this means O ≠ O_∇.

Let us call the basis  f_n, n=1,…,∞ in L “optimal” if its elements are both orthogonal and gradient-orthogonal:

    (f_i, f_j) = (f_i, f_j)_∇ = 0  if i≠ j.

For the modal space F formed by any N' functions of the basis F= ⟨ f_i_1,…,f_i_N'⟩, condition (<ref>) is obviously satisfied, and thus limiting the number of decomposition modes of the optimal basis will not introduce an aliasing error.
In this basis, the solution to the optimisation problem (<ref>) is given in a way similar to <ref>:

    c'_n = (f_n,ϕ)_∇/(f_n,f_n)_∇.




Now we are going to find a necessary condition for a basis to be optimal  

for an arbitrary  aperture.
Suppose F has matching orthogonal complements O and O'. Then for any vector ϕ, its orthogonal projection on F should match in both inner products

    _F (ϕ) = _F,∇ (ϕ).

From this it follows that the L_2-scalar product and the gradient-dot product of any test function ϕ and any modal vector f_n should be proportional to each other, or more accurately, if k_n is the gradient-norm of a basis vector f_n, then

    (ϕ, f_n) = 1/k_n^2 (ϕ, f_n)_∇.

Let us take a Dirac δ  as a test function, ϕ() = δ( -_0) for some _0 ∈ P. 
By definition of the gradient-dot product  <ref> and by Stokes theorem, we get

    f_n(_0) = ( δ( -_0), f_n) = 1/k_n^2 ( δ( -_0), f_n)_∇
     
    =  1/k_n^2 ( ∇δ( -_0), ∇ f_n) = - 1/k_n^2 (  δ( -_0), ∇^2 f_n) 
     
    = -1/k_n^2Δ f_n (_0)   ∀_0 ∈ P,

that is to say, f_n should be the eigenfunction of the Laplace operator, and is given by the solution of the Helmholtz equation in the aperture P:

    Δ f_n + k_n^2 f_n = 0

for some k_n ∈.
This equations comprises the necessary condition for {f_n} to form an optimal basis. 


As the eigenfunctions of the Laplace operator satisfying to the homogeneous boundary condition form an orthogonal basis, we have just proved that any optimal modal space for the gradient operator is given by the solutions to <ref> with some homogeneous boundary conditions.
The trigonometric basis for the square aperture of <ref> is a particular case of such a  solution
to <ref>, with periodic boundary conditions.
For a circular aperture, this is a classical problem of mathematical physics, with the analytical solution  bounded at =0 phase given in the polar coordinates = (r,θ) by

    f_m,n(r, θ) = J_m(k_m,n r) 
    	 0pt1cossin (m θ),

where J_m(r) are the Bessel functions of the first kind, and k_m,n some constants depending on the boundary conditions.


It's interesting to note that this set of functions, referred as the Bessel circular functions or the membrane vibration modes, has already popped up recently in the literature<cit.> as being more suitable than  Zernike polynomials for some applications, but, to our knowledge, was not yet considered as optimal for wavefront sensor-less methods.







§ NUMERICAL SIMULATION


















In this section we compare the accuracy of aberration correction using sub-optimal and optimal bases  with a numerical simulation.


It is not evident whether a basis satisfying the optimality criteria would suit a particular application, because optimality as defined in this paper means only that the result of two optimisation problems would be the same and no aliasing will be introduced by using the gradient norm. 
In other words, the basis might be optimal for a given wavefront sensing method, but not for representing a typical aberration of an application.



The magnitude of the aliasing error is also defined by the statistical properties of the aberration and  by the speed of convergence of the generalised Fourier series. 
For instance, the eigenfunctions defined with the Dirichlet boundary conditions seem a priori to be less suitable for  aberrations with non-constant values on the boundary.


As test cases, we have used a) wavefront aberrations produced by mouse oocyte cell as an example of specimen induced aberration in microscopy <cit.> and b) wavefront aberrations produced by Kolmogorov's turbulence, covering two important applications of AO[see Ref.  for the fitness analysis for ophthalmic applications].
Although the iterative nature of the indirect wavefront sensing methods makes them less suitable for the  case of the atmospheric turbulence, we have included this example for comparison purposes and to get an impression of the extent to which the improvement given by the optimal modes depends on the wavefront statistics.


The oocyte phase screens were obtained by the procedure described in Ref.  with the same parameters as its authors used for the results demonstrated in their Figure 5. 
The Kolmogorov phase screens were calculated by the method described in Ref.  with 20 phase screens for the ratios D/r_0 = 25/ i  with integer i in the range from 1 to 25. 
The typical phase screens for both test cases are shown in the left columns of 

<ref>.
The screens are shown wrapped for the presentation purposes only; unwrapped phase screens were used for approximations.
500 random phase screens for each of the cases were generated; several aberration-free phase screens were generated 

by the method of Ref. ; they were discarded, resulting in a total of 459 phase screens for the oocyte-case.







For each of the phase screens,
we compared the numerical results of its approximation 
by the gradients of the first 10 orders (= 65 terms) of the 
a) Zernike polynomials (denoted in the plots as Z);
b) LB polynomials  (LB);
c) SVD modes obtained from the first 4 orders of the Zernike modes with respect to the gradient-dot product (SML);
d) SVD modes obtained from  all 10 orders of the Zernike modes with respect to the gradient-dot product (SM), and eigenmodes of the Laplacian satisfying  
either to
d) Neumann boundary conditions (OMN) or
e) Dirichlet boundary conditions (OMD)[the results for the OMD case are not present in this paper, but are available in the supplementary material <cit.>]. 





Each of the basis functions and aberrations was discreetized on a 128 × 128 grid, with a circular aperture of a 125 pixels diameter (see <ref>). 
The Zernike polynomials were calculated using their analytical formulae and normalised with respect to their ℓ_2 norm on the aperture. 
Their gradients were approximated by first-order finite differences inside the aperture; the values on the aperture boundary were zeroed. 
The Lukosz-Braat polynomials were calculated numerically through Cholesky decomposition of the gradient-dot product matrix M'_Z of the Zernike modes (because  the gradients are not defined on the boundary, this provided more accurate results for the chosen resolution than calculation  using analytical formulae). 
The SVD modes were obtained by SVD decomposition of M'_Z. 
The optimal modes were obtained in Mathematica software with the  function with the corresponding boundary condition.  
All modes were, if necessary, first renormalised with respect to the aperture and then reordered by increasing their gradient-norm. 

The dot product and the gradient dot product matrices M and M' were calculated numerically for the discreetized bases for illustration purposes (see right column of <ref>); the matrices allow to see the key difference between seemingly similar basis function shapes. 
Please note how SVD modes diagonalise both matrices only for the orders from which they were calculated (the top left minor formed by the red lines in the matrix plots). 

For each of the phase screens and each of the bases, we have calculated
[a)]
	
  * its rms error by approximating it by the first N' Zernike terms, 
	
  * the rms error of its approximation by the gradients of the first N' terms of the basis,
	
  * the Strehl ratio, and
	
  * the second moment of the PSF of the remaining aberration (normalised to that of the diffraction-limited PSF).


	The Strehl ratio was defined as the maximum intensity of the central pixel, normalised to the diffraction-limited maximum. 
	This definition results in a specific dip of the plots for the oocyte-produced small aberrations, like in the top rows of <ref>: for instance, correction of tilt only would move the brightest pixel further away  from the image centre. 
 

The results for 5 representative aberrations of various strengths for each set are shown in Figs. <ref> and <ref>.


Figure <ref> presents the box-plot for the relative rms error for the ensembles of the phase screens; it is shown only for the first third of the modes for better visibility.






From the results shown in Fig. <ref>, it is clear that the  optimal modes obtained from solution of Helmholtz equation with Neumann boundary conditions (denoted as OMN) provide on average better results for smaller dimension of the modal space (first 4 orders of Zernike, denoted as SML), and
similar approximation quality to that of the numerically calculated SVD modes obtained from the first 10 orders of Zernike polynomials (SM).
This result can be explained by the fact that for large enough N, the part of the aberration not covered by the span of the basis function is small, and correspondingly the aliasing error introduced by it is small as well.



















§ DISCUSSION AND CONCLUSIONS



While the results of <ref> describe the best analytical basis for SM-based methods, this basis can be difficult to implement in practice because of alignment errors, necessity of exact calibration of the adaptive element and so on, so SVD modes are used commonly.
Numerical results of <ref> show that
singular modes obtained from the large number of Zernike polynomials look similar to the optimal modes (<ref>), which can be seen as the limiting case of the SVD-ing with the number of modes going to infinity. 
This is also confirmed by the similarity of the results (see <ref>, cases SM and OMN). 

 

The optimal basis for a given aperture is not unique (for instance, both OMN and OMD are optimal), 
and thus the question of a proper optimal basis for a given application remains open. 
Similarity of SVD and optimal bases suggests that an approximation of the optimal basis can be used. 
For instance, for a piecewise smooth aberration like that of the oocyte-induced aberrations, SVD orthogonalisation of Gaussian radial base functions might give better results.
Moreover, as one can see by the first row of <ref>, for a sparse aberration present only at one edge of the aperture and a small number of correction modes, the Lukosz-Braat functions give better results.
This can be explained by the specific shape of the aberration, which has 70 per cent of its power uniformly contributed by the first 30 Zernike polynomials.

An important conclusion on the practical use of SM-based wavefront sensorless methods can be drawn from   <ref>.
To minimise the aliasing effect, the response functions of  the adaptive element should be a good approximation of some optimal basis, i.e. approximately satisfy <ref>.
Since the response functions of a membrane deformable mirror, a bi-morph mirror, and a free-edge mirror with push-pull actuators  satisfy either Poisson or biharmonic equations with a piece-wise constant right-hand part <cit.>, none of them, nor any basis obtained by orthogonalisation of their gradients, can form a true optimal basis, and thus the use of low-order deformable mirrors for correction with a second-moment based method will suffer from aliasing effects. 
Again, some approximation to the optimal basis can still be obtained with these types of mirrors with a large number of actuators, e.g. a photo-controlled deformable membrane mirror <cit.>.
For a low-order adaptive element, the layout of the actuators should be optimised in order to get response functions closer to the optimal modes (and not to Zernike modes, for instance).









To conclude, in this paper
we have derived the necessary condition for an analytically defined basis to be optimal, that is, aliasing-free, for wavefront sensorless methods based on the second moment of intensity, namely that each function of the basis should satisfy the Helmholtz equation, and we have provided examples of such a basis for square and circular apertures. 
The analytical optimal basis can be considered as a limiting case of a numerically computed optimal basis via SVD of the gradient-dot product when the number of basis functions is sufficiently large. 
The proposed analytical modes can be useful when controlling an adaptive element with a large number of actuators, as probe modes (starting point) for building an empirical optimal basis, and when designing low-order adaptive elements for SM-based methods.





§ ACKNOWLEDGEMENTS


The author is grateful to G. Vdovin, M. Verhaegen, and D. Wilding (TU Delft) for their attention to the work and valuable comments and to D. Soloviev for proofreading the manuscript. 




plain
 































