1. Introduction

Curvature sensors are non-interferometric devices that provide the phase of an electromagnetic wave by means of measurements of intensity. They were proposed by Roddier at the end of the 80’s given their high use-potential in Adaptive Optics to correct stellar images affected by atmospheric turbulences [1–2]. The success levels reached together with the simplicity and robustness of the experimental devices has fostered the use of these sensors (and other variations for non-uniform intensity) in a number of different applications such as evaluation of the optical quality of telescopes and of optical components [3–7] and, more recently, in the detection of wave-front tilt for accurate pointing of laser beams [8] and phase imaging of thermal atomic beams [9].

Curvature sensing is based on the Transport of Intensity Equation (TIE) that associates the irradiance ce I(x, y; z) and phase ϕ(x, y; z) of a paraxial light beam propagated in a given z direction [10]. Curvature sensors are limited to the case in which the distribution of intensity on the plane for which phase estimation is desired is uniform, I ₀, and the TIE thus associates the local curvature of the phase with the axial derivative of its intensity:

Therefore, if the axial derivative of the intensity is known, the phase can be obtained by simply resolving the Poisson equation described in the above equation.

The experimental estimate of the axial derivative of intensity is acquired by subtracting the irradiance measured on two adjacent planes and dividing by the distance separating them; from a mathematical point of view, the approach is by means of finite differences. Here we will use the central difference approximation (two planes symmetrically placed from z=0) since it is the most common scheme of measurement, but all the mathematical treatment developed here can be applied to any other approximation, as forward or backward difference approximation. In any of the cases, the quality of the final result depends particularly on obtaining precise estimations. This is because if the planes are very close, although from a purely mathematical point of view the calculation of the axial derivative is more precise, the sensitivity is low [3]. On the other hand, if the distance between the planes is large, the signal is less affected by the noise but the calculation of the derivative is less precise. In short, the inevitable presence of noise in intensity measurements gives rise to the need for an analysis of the separation between planes that provides an ideal result. Despite the proven effectiveness of these sensors in a large number of applications, the problem of the proper and objective selection of the planes of intensity in function with the detection conditions has, in our view, yet to be resolved. Some authors have done short studies in this respect and have offered qualitative criteria for selecting the separation between planes as a function of the noise level present in the intensity. Specifically, references [10] and [11] propose separating the planes by an undetermined distance within a bounded interval determined by the characteristics of the noise and the detector used. It is then up to the researcher to select a distance from within this range. Moreover, the criteria used to position the planes is not normally made explicit in the literature and this means that there is no guarantee that the separation chosen is the one providing the best measurement [4–7,9]. In light of this situation it becomes essential to highlight the need to progress in the pursuit of objective criteria to serve as a guide in the proper selection of the planes that guarantee the most precise result.

This work tackles the problem of determining the separation between planes that allows for the best estimate of the axial derivative bearing in mind the random disturbance of the experimental data. The results depend on the statistics assumed for the noise. We will suppose that intensity measurements at each point of any plane are affected by uncorrelated random noise with zero mean. Different statistics for the noise must be studied separately.

This article is organized as follows: the second section sets the theoretical framework of the problem and the variance of the axial derivative of the intensity will be obtained as a function of the noise level of the measurements and the separation between planes. The third section proposes a way of choosing the distance between planes that optimizes sensor performance. And finally in the fourth section an illustration of the proceeding with two simple examples will be found.

In the future we will resolve the spacing between planes that optimizes the boundary signal. This issue is not addressed in this work because it requires a mathematical analysis different from that of the domain signal and we feel that is should be dealt with independently.

2. Theory

As was mentioned in the introduction, the phase of a paraxial beam propagating in the z direction can be obtained by resolving the Poisson equation

∀(x, y) within the domain. Thus, the quality of the final result depends strongly on the precise experimental estimation of the right hand side of the above equation. We will assume that the intensity i(x, y; z) measured at each point (x, y) of any plane z is affected by noise that is supposed to be random with zero mean value, negligible covariance between the different measurements and with variance σ ². Therefore, at each point (x, y) of a given plane z, the intensity measurements can be described by:

where n is the noise of the measurement.

The estimate of the axial derivative at each point (x, y) of the plane z=0, that we express as Î′(x, y;0), can be obtained easily by applying the central finite differences approximation:

In order to obtain an accurate estimate of the intensity derivative at a given point (x, y), we will search for the position of the plane z (z=z(x, y)) so as to verify the following two conditions:

where the upper bar denotes mean value.

In doing so we will keep Taylor’s formula about z=0 in mind [12]. Thus, by assuming that I(x, y; z) is a continuous, single valued function of z having continuous first derivative and having second derivative, then:

where the third addend on the right side of the Eq. (6) represents the remainder term of the Taylor’s series, in which ξ ₊(x, y) is a value within the interval [0, z] and ξ _-(x, y) within [-z,0].

Substituting Eq. (6) in Eq. (4) and taking Eq. (3) into account:

Under the assumption that there is no noise correlation between the points of the two planes, it can be deduced from the above equation that the variance of the estimated axial derivative is given by:

where the first addend is due to the noise of the intensity measurements and the second one depends on the specific nature of the wave front being tested. Mindful of the fact that the curvature sensors measure local wave-front curvatures that are nearly flat, it can be assumed that the second derivative of the intensity at any point (x, y) on any plane within the interval [0, z] is a bounded magnitude, the bounds of which depend on the specific application in which the sensor is used and, therefore, the above equation demonstrates that for every point (x, y) a plane separation, z _min (x, y), could be found allowing for minimization of the curvature measurement variance at this specific point. Nevertheless, it is clearly impossible to put Eq. (8) into practice due to two fundamental problems: the fact that the remainder term of the Taylor’s series for each point of the space is unknown and the fact that the separation between planes providing more exact measurements depends on each point (x, y).

3. Selecting the ideal measurement plane

The basic aim of this section is to propose a selection method for the distance between planes that will provide the best results. First of all we do a local analysis of the propagation of a wave front allowing us to overcome the two problems referred to above; then we will propose a method for choosing the separation between planes and finally, in the next section, we will illustrate the method with a numerical simulation and the application for the retrieval of phases affected by Kolmogorov turbulences.

3.1 Local phase analysis

On the one hand, we must consider that we are dealing with paraxial approximation that treats Huygens wavelets as if they were emitted in a narrow cone with a parabolic approximation for the wave-front shape of each wavelet (geometrical optics approximation). It can therefore be assumed that in a neighborhood of a point (x ₀, y ₀), Ω(x ₀, y ₀), of the plane z=0, an approximation to the phase, ϕ(x, y), can be calculated by a second order power series:

where a(x ₀, y ₀)=ϕ(x ₀, y ₀) represents a local piston; $b(x_0,y_0)={\left(\frac{\partial \varphi }{\partial x}\right)}_{(x_0,y_0)}$ and $c(x_0,y_0)={\left(\frac{\partial \varphi }{\partial y}\right)}_{(x_0,y_0)}$ local tilts assumed to be negligible [11]; and $d(x_0,y_0)={\left(\frac{1}{2}\right)\left(\frac{{\partial }^2\varphi }{\partial x^2}\right)}_{(x_0,y_0)}$ , $e(x_0,y_0)={\left(\frac{1}{2}\right)\left(\frac{{\partial }^2\varphi }{\partial y^2}\right)}_{(x_0,y_0)}$ and $f(x_0,y_0)={\left(\frac{{\partial }^2\varphi }{\partial x\partial y}\right)}_{(x_0,y_0)}$ are related to the local curvature and the Hessian of the phase. Equation (9) means that the phase on a plane can be considered as a covering of paraxial spherocylindrical wave fronts, i.e., for every point (x ₀, y ₀) of the plane z=0 there exists a neighborhood, Ω(x ₀, y ₀), where the local phase can be accurately described by a quadratic approach.

On the other hand, we may consider that the intensity propagation of a paraxial beam along the axis can be described for each point (x, y) of a plane z in the following way:

where ∇²ϕ(x_c, y_c ) and H[ϕ(x_c, y_c )] represent the Laplacian and the Hessian of the phase on the plane z=0 evaluated at first kind critical points (x_c =x_c (x, y; z), y_c =y_c (x, y; z)) [13].

By taking this into account, we can deduce that at each point (x, y) within the neighborhood Ω(x ₀, y ₀) of a given point (x ₀, y ₀), the evolution of the intensity along the axis is given by

and it does not depend on the critical points given that ∇²ϕ_loc(x, y)=2(d+e) and H[ϕ _loc (x, y)]=4de-f ² ∀(x, y)∊Ω(x ₀, y ₀) (we have dropped the dependence on (x ₀, y ₀) of the expansion coefficients d, e and f).

If we now substitute Eq. (11) and the TIE (Eq. (1)) in Taylor’s series of Eq. (6), the remainder terms of both Taylor series can be easily evaluated:

and therefore the variance in the measurements of the derivative within the neighborhood Ω(x ₀, y ₀) depends on the local geometric characteristics of the phase only through its Laplacian and Hessian:

Given that both the Laplacian and the Hessian are invariants and in order to simplify the problem, we can rotate the coordinate axis in such a way that the equation of the conic defining Eq. (9) can be rewritten in its canonical form:

Given that ∇²ϕ_loc(x, y)=2(d+e)=2(D+E) and H(ϕ _loc (x, y))=4de-f ²=4DE Eq. (13) can be rewritten:

The minimization of the above equation with respect to variable z provides the value of the separation between planes so that the estimation of I′_loc(x, y; 0) is more exact; it should be pointed out that it does not depend on the coordinates of the point (x, y) within the region Ω(x ₀, y ₀) where the approximation to the phase given by Eq. (9) is accurate.

In order to study the behavior of s _loc(z) we will carry out below a detailed analysis for different values of the principal curvatures D and E, which will allow us to establish a method for choosing the measurement planes.

The following notation will be used:

and, due to the symmetries of Eq. (15), only the following cases will be analyzed:

Besides, we will only consider cases where 2(|D|+|E|)≤U. The purpose of taking an upper bound U for the sum of the absolute value of the principal curvatures and not for the absolute value of the Laplacian is to eliminate, for example, cases of hyperbolic wave fronts with very high principal curvature values that could produce zero or very small values for the Laplacian. This implicitly means putting a higher limit on the absolute value of the Hessian at each point which guarantees that the wave front being tested is nominally flat. Bearing this in mind, and taking into account the validity conditions of the Transport of intensity Equation it can be even assumed that if the maximum absolute value of the Laplacian for a point within the measurement domain of any phase is approximately equal to U, then at any other point of the domain it can be ensured that 2(|D|+|E|)≤U.

In Fig. 1 we plot the variation of s _loc(z) for the four main kinds of local approximations to the phase: Figure 1a shows the behavior for different cylindrical wave fronts, 1b for elliptical wave fronts, 1c for spherical wave fronts and 1d for hyperbolic wave fronts. Given that the qualitative behavior of the standard deviation of the different cases does not depend on the values taken for the variance in the intensity measurements, σ ², and for I ₀, we have considered σ=0.01 and unit intensity. Moreover, we have normalized the bound U to unity, so that 2(|D|+|E|)≤1. In any of these four cases it can be observed that:

We must point out that in all cases z _min moves away from z=0 as σ increases.

 

Fig.1. Standard deviation of estimated axial derivatives within Ω versus measurement planes position for (a) cylindrical wave fronts, (b) elliptical wave fronts, (c) spherical wave fronts and (d) hyperbolic wave fronts. In all cases σ=0.01.

Download Full Size | PDF

Furthermore, of all of the combinations of D and E values such that 2(|D|+|E|) has the aximum value, that which features the worst deviation for any z value is the one corresponding to a cylindrical wave front followed, in decreasing order, by the elliptical, spherical and hyperbolic (see Fig. 2). This is because the cylindrical wave fronts are the ones with the least linear behavior of intensity as is shown in Fig. 3 for different principal curvature values leading to the maximum value of the Laplacian.

 

Fig. 2. Error curves corresponding to four paraxial spherocylindrical phases with 2(|D|+|E|)=1. z _opt denotes the position of the plane which minimizes the standard deviation of the estimated intensity derivative for the cylindrical wave front.

Download Full Size | PDF

 

Fig. 3. Evolution of the intensity along the z axis. The curves plotted show the behavior for combinations of D and E for which 2(|D|+|E|)=1.

Download Full Size | PDF

We will now see how to use these local results to choose a proper separation between planes that optimizes the operation of the sensor.

3.2 Selection of the best position

In general terms, in the design of the set up of any type of wave-front sensor, one should have a certain degree of a priori knowledge of the phase to be measured in order to properly choose its principal components. For example, in the case of Shack Hartman wave-front sensors, prior information of the phase is needed in order to design arrays of micro-lenses with suitable sampling and numerical aperture. In the case of curvature sensors, once the measurement noise of the available detector is determined, the next step in the design is the selection of the separation between measurement planes.

In Fig. 2, z _opt denotes the bearing of the planes that minimizes the variance of the axial derivative for a cylindrical wave front with the maximum curvature (U=1) of all of those analyzed in Fig. 1. Let us now turn our attention to a phase to be measured that can be locally represented by a covering of any kind of spherocylindrical wave fronts with local curvatures varying between -1 and 1.

If we were to position the measurement planes at z<z _opt then, as can be deduced from Fig. 1, the estimated derivatives in any part of the domain will always show larger errors than would be the case if the planes were positioned at z _opt. On the other hand, if the planes were to be situated at z>z _opt then the estimated derivatives for the neighborhoods that can be approached by cylindrical wave fronts with strong curvatures would worsen. Therefore, if there were a significant density of neighborhoods within the measurement domain in which the phase could be approached by means of maximum curvature cylindrical wave fronts (U=1), the overall estimation of the phase would worsen.

Thus, in practice, if we only know that the phase to be measured features maximum curvatures of a certain value V, minimizing Eq. (15) for a cylindrical wave front of curvature V, i.e., D=V/2, E=0, we would obtain the optimum positioning of the measurement planes compatible with the available information.

Since usually one has or can deduce more information from the statistics of the phases to be measured, a more refined positioning can be performed. Logically, the more a priori information that one could have the more refined the final results would be.

For instance, to improve the above result, let us suppose that the a priori knowledge of the phases to be measured allows us:

If no further information is available about the phase to be measured, we propose situating the measurement planes at the distance, z _opt, that minimizes Eq. (15) for a cylindrical wave front with D=V/2, while V=L _max+ΔL _max.

In both situations, s _(D=V/2,E=0) (z _opt) represents the least upper error bound of the estimated derivative at any point of the domain and the evaluation of the standard deviation at z _opt for a plane wave front, s _(D=E=0) (z _opt), provides the greatest of the lower bounds, i.e.:

In short, depending on the a priori knowledge about the phase, a value V can be estimated, so that the optimum separation is determined by calculating the position z _opt for a cylindrical wave front with D=V 2, i.e., minimizing:

or equivalently, minimizing:

where h=Vz.

In Fig. 4 we show the behavior of s _{D=V/2, E=0}(h)/V for different values of σ together with the behavior of s _D=0,E=0(h)/V for plane wave fronts:

As σ increases z _opt=h _opt/V moves away from z=0. For any given value of σ, s _(D=V/2,E=0)(h _opt)/V provides the upper bound of error, while s _D=0,E=0(h _opt)/V corresponds to the lower one. Moreover, we can also reach the conclusion that if the intensity measurement error is in the range of 70% or higher, the error made in the estimate of the local curvatures is greater than or equal to 100% regardless of the position of the planes.

 

Fig. 4. Upper bound (solid line) and lower bound (dashed line) of error in the estimation of the axial derivative of the intensity for several values of the noise level σ.

Download Full Size | PDF

Finally, Fig. 5 shows for I ₀=1 the dependence of the optimum normalized position, h _opt, versus σ. This graph fits accurately to a second order polynomial in ln(σ) in the form:

with α=0.6278, β=0.1438 and γ=0.00946. From this formula the optimum separation of the planes can be easily calculated by simply taking the corresponding value of h _opt and dividing by the value of V obtained as explained above throughout the a priori knowledge of the phase to be measured.

 

Fig. 5. Optimum positioning versus noise level plotted within the operational range of the curvature sensor. The points in blue are evaluated by finding the minimum value of Eq. (18) for intensity normalized to unity. In red the linear curve fit is represented.

Download Full Size | PDF

4. Examples

4.1 Spherical aberration and defocus

This first illustrative example represents the retrieval of the local curvatures of a simple rotationally symmetric phase defined by ϕ=(r ⁴-4r ²)/16 where r∊[0,1].

 

Fig. 6. Standard deviation of the estimated derivative (blue line) and upper error bound (red line) and lower error bounds (green line) predicted by the theory.

Download Full Size | PDF

4.2 Atmospheric distortion

Let us assume that we have a priori knowledge of the variances and covariances of the coefficients of phase expansion as an infinite sum of orthogonal polynomials as is the case with phases affected by Kolmogorov turbulences. In this case, if the five first orders are considered, the temporal mean values are independent of the point within the pupil and they are given by [14]:

By taking λ=0.63 µm, Fried parameter r ₀=30 cm, telescope radius R=4 m and I ₀=1, then $V=\sqrt{〈{\left({\nabla }^2\varphi \right)}^2〉}=4.58·{10}^{-6}{\text{m}}^{-1}$ . By following the procedure explained above that assumes that there is a significant density of neighborhoods within the pupil in which the phase can be locally approached by means of cylindrical wave fronts (with D=V/2 and E=0), we evaluated the optimum positioning for different values of σ and the results are shown in the first column of Table 1.

Nevertheless the previous result can even be improved since there is more information available about phase when dealing with atmospheric turbulences. Thus, it can be thought that the isotropic behavior of the statistics of the phase supposes that the most probable covering of the phase is by means of spherical wave fronts, then the optimum positioning could be obtained by minimizing Eq. (15) not for a cylindrical wave front but for a spherical wave front, i.e., D=E=V/4. The second column of Table 1 shows the optimum positioning for the same parameters used in the previous case.

5. Conclusion

By considering statistics of intensity measurement noise random, zero mean value with standard deviation σ and by means of a local geometric analysis of the phase on the measurement plane, we have proposed a method to estimate the position optimizing the operation of the curvature sensors in the sense that said position gives rise to the best estimate of the axial derivative of the intensity when calculated as a central difference approximation. The steps presented allow for the simple calculation not only of an optimum measurement position but also for an upper and a lower bound to the error of each local estimated derivative. It must be stressed the fact that the evaluation of an optimum positioning depends on the a priori knowledge of the phases to be measured. We provide a formula to locate the measurement planes when only a minimum knowledge of the phase is available. In any situation the fact is that the most pronounced local curvature of a wave front with high statistical weight is what determines the separation between planes. In any case, as could be expected for small values of σ and/or presence of strong curvatures of the phase the distance between the planes becomes smaller. Finally, the operational limit of the sensor is established around noise measurements errors greater than 70% regardless of the local curvature values.