1. Introduction

Optical vortices (also known as branch points or phase singularities) occur when the phase of an optical field becomes undefined in a region of zero amplitude. They were comprehensively studied by Nye and Berry in the 70’s [1] and have continued to be a subject of interest in a number of different fields ever since, including optical tweezers [2] and optical communications [3]. There is a spiral phase change increasing from 0 to 2mπ around the undefined phase at the centre as demonstrated in Fig. 1, where m is an integer number corresponding to the charge of the vortex. Vortices naturally occur when a laser beam propagates through the atmosphere and they form in pairs [4]. These pairs consist of both a positive and a negative vortex, similar to an electric dipole, with the sign determined by the direction of the rotation of the spiral phase pattern around the centre of the vortex. The wavefront dislocation which joins the two oppositely signed vortices is along the line where there is 2mπ jump in the phase and is referred to as a branch cut. A paper by Fried and Vaughan in 1992 [5] showed that branch cuts are unavoidable in the phase function when a wave is distorted by the atmosphere and that they line up along the regions of low intensity.

 

Fig. 1. A 2π radian vortex in an atmospherically distorted phasefunction, the vertical axis is the phase in radians and the xy plane is the space where the phase is measured. This shows the spiral phase around the centre of the vortex as well as the wavefront dislocation created by the vortex.

Download Full Size | PDF

There is much interest in trying to correct for the distortions induced by the atmosphere on a propagating laser beam, with one the main techniques for accomplishing this being adaptive optics (AO). AO systems measure the wavefront and correct for its aberrations by using a controllable optical element and they generally consist of a wavefront sensor, a deformable element and a control system along with other static optical elements. It is now being used for a number of different applications, including imaging through the atmosphere [6] and line-of-sight optical communications [7]. Some work has highlighted a number of the difficulties that optical vortices can cause to adaptive optics (AO) systems for both the hardware [8–10] and control system [11].

While some proposals to overcome certain hardware issues have been put forward, such as which deformable mirror configuration is best [12], others remain. Two major difficulties in this field are the accurate and reliable detection of the vortices in the wavefront and the subsequent faithful reconstruction of the measured wavefront. Different reconstructors have been suggested [13,14] and some of these have been shown to work in laboratory conditions [15,16], while others work in theory but are impractical in reality [17, 18].

The issues of wavefront estimation and compensation are interlinked as correct reconstruction is impossible without proper sensing of the wavefront. In this paper we show that accurate detection of optical vortices is possible with a Shack-Hartmann wavefront sensor in conditions similar to those experienced in the atmosphere. The detection of the Shack-Hartmann spot positions, and thus the wavefront slopes, is improved by the use of a matched filter algorithm instead of a centroiding algorithm. The vortex detection is accomplished using a technique called the branch point potential method [19], which was first proposed by LeBigot and Wild [20–22] along with the slope discrepancy technique described by Tyler [23].

2. Shack-Hartmann wavefront sensing

2.1. Shack-Hartmann spot detection and centroiding using the Hotelling Observer

In this paper the Hotelling observer [24, 25] is used for the purpose of spot detection and centroiding. The Hotelling observer is sometimes referred to as a prewhitening matched filter [26]. In the process of prewhitening, the data, g, is divided by the data covariance matrix, K_g, with the aim of producing spatially stationary, uncorrelated noise. In this study an average point spread function (PSF) from a reference image of spots, produced from a flat wavefront, was used to obtain an estimate of the data covariance matrix. This PSF is then used to flatten the residuals (prewhitening) and to estimate the spot’s position via matched filtering. For the case of faint companion detection in astronomy, it has been shown that the spatial Hotelling observer outperforms the spatial matched filter both in terms of detection and signal to noise ratio [27]. The use of matched filters is also not new to Shack-Hartmann wavefront sensing [28, 29]. We will show that the Hotelling observer provides a better detection and centroiding method than a matched filter alone.

The observer/algorithm is charged with deciding upon the presence of a spot in an image. There are two possible outcomes to this decision: the spot present hypothesis, referred to as H ₁ or the spot absent hypothesis, denoted H ₀. The ideal observer for such a test is the likelihood ratio:

or its logarithm λ(g|r⃗_spot)) = ln[Λ(g|r⃗_spot)], where pr(g|H ₁,r⃗_spot) is the probability density function (PDF) of the data g under the H ₁ hypothesis with a spot located at the position vector r⃗_spot and pr(g|H ₀) is the spot absent PDF. The Hotelling observer takes on the forum:

where $w=K_g^{-1}{s}_{{\overrightarrow{r}}_{\mathrm{spot}}}$ , is the template to be matched to the data and $S_{{\overrightarrow{r}}_{\mathrm{spot}}=A.h\left({\overrightarrow{r}}_{\mathrm{spot}}\right)}$ is the signal of the spot, where A is the intensity of the spot (a scalar quantity) and h(r⃗) is an estimate of the PSF at the position vector r⃗. It can be shown [30] that when the data is normally distributed the Hotelling observer is equal to the log-likelihood ratio, λ(g|r⃗_spot)) = t_Hot(g|r⃗_spot), as both are linear functions of the data and optimal. The Hotelling observer therefore only requires knowledge of the modeled spot signal and the data covariance matrix. Therefore the goal is to find a template which maximises the correlation between the template and the data. If K_g|H_i for i ∊ {0,1} are the covariance matrices under the H ₀ and H ₁ hypotheses then the average covariance matrix is:

The Hotelling observer is now expressed as:

This expression shows how the Hotelling observer can be seen as a prewhitening matched filter, as the operation of multiplying the estimated spot signal by the inverse of the data covariance matrix is akin to the classical signal analysis process of prewhitening. Assuming the background level, b, of each spot image can be estimated along with its variance, σ ², and an estimate of the PSF is known it can be shown that K_g is diagonal with it elements given by [31]:

The cross-correlation of the template vector, w, and the data, g, can be computed either in the image plane or in the Fourier plane. In the Fourier plane the position of the spot is estimated by using a parabolic fit to interpolate the peak of the cross-correlation [28, 29]. In this study both the data and the template vector were zero padded to five times their original size to improve the resolution of the cross-correlation. In the following Hot_MF will refer to using the Hotelling observer in this manner and MF will refer to using a simple matched filter with this approach. In the image plane the peak of the cross-correlation is found, and hence the spot position, by shifting the position of the template spot until maximum is found. In practice this PSF fitting type approach is carried out by an unconstrained maximisation of the Hotelling observer where the value of the Hotelling observer is only dependent upon the position of the test spot in the template vector [25]. This image plane method will be referred to as Hot_ML.

To test the method, ten thousand noisy spot images (size 19 × 19 pixels), with random spot locations, were simulated using a Gaussian PSF profile with a mean full-with-half-max (FWHM) of 2.5 pixels varying normally with a standard deviation of 0.25 pixels from image to image. The noise in each image consisted of Poisson noise and normally distributed readout noise. For each image the three matched-filter type algorithms, as well as a classical centroiding algorithm, were tasked with locating a spot in the image. The PSF profile used in the template vectors was also Gaussian in shape and had the same FWHM as the mean of the data. The Hot_MF showed the lowest mean error on the estimation of the position of the spot (in pixels) followed by MF and Hot_ML, see Table 1. The centroiding algorithm faired worst; however, it should be noted that this centroiding algorithm was not optimised as in [29]. The centroiding algorithm used in the comparison is a simple, non-iterative centroiding algorithm. It thresholds the data and finds the centre of mass of a spot bigger than a certain size. The threshold is specific to each lenslet as it is intensity weighted for each subaperture. If two or more spots are detected above the threshold within one subaperture the brightest one is chosen as the focused Shack-Hartmann spot.

Table 1. Summary of testing the Hot_MF, Hot_ML, MF and centroiding algorithms on ten thousand simulated spot images. The Hot_MF showed the lowest mean error in spot position estimation. This algorithm was also the most efficient spot detector having the highest AUC of the four tested observers.

View Table

The signal-to-noise ratios (SNR) reported in Table 1 were defined as the ratio of the mean to the standard deviation of the measured data [6]. In practice the SNRs were calculated by comparing the peak value of the cross-correlation (or the data) to the standard deviation of the values of the cross-correlation located in a ring around the edge of the image, i.e.

 

Fig. 2. Plots of the LROC curves for: Hot_MF, Hot_ML and MF.

Download Full Size | PDF

where the overline denotes average.

For the three matched-filter type observers it was possible to calculate a localisation receiver operating characteristic curve (LROC) [30]. The LROC curve is a plot of the true positive fraction of correct decisions made by an observer, in this case the number of correct spot detections coupled with localising a spot to within a given tolerance ±ε = 4 × 10⁻⁴ pixels, versus the false positive fraction of decisions made, i.e. the false alarm rate. The area under the LROC curve, denoted AUC [30], is used as the figure of merit for the performance of the algorithms, see Fig. 2.

We see from Table 1 and Fig. 2 that the cross-correlating Hotelling observer outperforms all the other tested algorithms in terms of mean centroiding error and in terms of the number correct spot detections, i.e. Hot_MF has the highest area under the LROC curve, AUC = 0.82. Furthermore the Hot_MF has a much higher SNR compared to the straight matched filter, MF, and centroiding algorithm. Hence this increased SNR coupled with a high performance in spot detection and localisation shows that the Hot_MF is capable of detecting spots with very low intensities which is critical when analysing spot images in close proximity to a vortex. It should be noted that computing the Hotelling observer in the image plane could be improved by interpolation of the PSF so that the resolution of the cross-correlation could be improved, as was done with the Fourier plane method.

2.2. Vortex detection

The slopes (or phase gradients) calculated from the Shack-Hartmann wavefront sensor are related to the optical wavefront by the following equation

where A is a matrix which relates the phase gradients also referred to as the geometry matrix, s, to the wavefront phase, ϕ.

The least mean square error reconstruction, ϕ _lmse, as described by Fried [11], can be written in the following form

This however does not take account of all the phase in the reconstruction process: some of the phase is, in essence, ignored and treated as noise. This is what Fried [11] called the hidden phase and this is the part of the phase which is discontinuous when optical vortices are present. It can be defined as

where ∇ϕ is the gradient of the phase of the field, ∇ϕ_lmse is the gradient of the least-squares phase and ∇ϕ_hid is the gradient of the hidden phase.

The hidden phase is the part of the phase which is the curl of the vector potential (the least squares reconstruction only looks at the scalar potential) and thus can be defined as

where the H(r) above is the Hertz potential, which is defined as:

where h(r) is the Hertz function. It will be this Hertz potential that we will take advantage of when using the branch point potential method to detect vortices in a distorted phase function.

Although the curl component of the phase is hidden to the least squares reconstruction method, it can be calculated from the slope discrepancy Eq. (16) [23] of the reconstruction. This is done by calculating slopes from a wavefront reconstructed by the least squares method and subtracting them from the slopes originally measured by the wavefront sensor. If a system is completely noise free, then the slope discrepancy is exactly equal to the curl component of the phase. However, experimentally the slope discrepancy will be noisy and have fitting errors which were filtered out by the least squares reconstruction.

The most commonly used method to examine the curl component of the phase to find the position and sign of optical vortices is the contour sum method of gradient circulation. This involves summing the phase gradients calculated from the wavefront sensor around in a closed loop and can be described mathematically by the following equation

where Δ^x() are the horizontal phase differences and Δ^y() are the vertical phase differences.

If there are no vortices present, and thus the phase is continuous, then this sum will be equal to zero (ignoring noise). It will also equal zero if there are two vortices of opposite sign enclosed by the loop of phase gradients, as they will cancel each other out. However if the sum around the loop is equal to, or a multiple of, ±2π then there is a phase singularity enclosed in the loop. By convention the (i,j)th pixel is the top left of the four pixels and the vortex is placed there, a more exact position of the branch point could be established by further numerical methods but we know of no standard technique. The overall robustness of the contour sum method has been questioned under heavy noise conditions [18, 21,32] making it problematic to use this with the slope discrepancy technique and therefore we have decided to use a different method [21, 22]. The issue of noise with respect to the contour sum method is raised because it only does a gradient circulation around a closed loop of four phase gradients, whereas the branch point potential method described below employs all of the available slope data by using a least squares reconstruction method to calculate the potential function.

This method, called the branch point potential method, inspects the curl of the phase by determining the Hertz potential of the field. It is done by first rotating the phase gradients calculated from the wavefront sensor by 90 degrees, which is accomplished by a simple matrix multiplication.

where (R _π/2) is a 90 degree rotation, s are the slope values measured by the wavefront senor with s_x and s_y the slopes x-values and y-values respectively.

Then the pseudo-inverse of the geometry matrix is calculated:

where A ^† is the adjoint of the geometry matrix and the m ² value is inserted to prevent the matrix from becoming singular. The pseudo-inverse could also be calculated using a singular value decomposition.

The potential function, V is then calculated by multiplying the rotated slopes by the pseudo-inverse of the geometry matrix.

The potential function, V, is then reformed into a square array of the same size as the Shack-Hartmann lenslet array and the position and sign of the vortices is found by locating the peaks or valleys in the potential function. The positive vortices are designated by the peaks in the potential plot and the negative vortices are designated by the valleys. This potential function can be used to locate the optical vortices because it is the Hertz potential Eq. (11) and therefore it is able to see the hidden phase of the field.

We have found that this method can be improved upon by using it in conjunction with the slope discrepancy technique. As mentioned previously the slope discrepancy is the slope field obtained by subtracting the Shack-Hartmann slopes from the slopes ascertained from a least squares reconstruction of phase, as shown

where ∇ϕ _slpdis are the phase gradients of the slope discrepancy, ∇ϕ_lmse are the phase gradients of the least-squares phase and s are the phase gradients calculated from the wavefront sensor.

The slope discrepancy phase gradients are then used as the inputs into the branch point potential detection method. As our results will show, when these are used instead of the Shack-Hartmann slopes we get improved vortex detection. This happens because tilts and other large aberrations that are present in the wavefront have a strong contribution to the slope data and can mask the relatively small circulation of the slopes when a vortex is present. By removing the continuous part of the phase from the slope data the vortex circulation becomes more apparent. Since these large aberrations are essentially noise to our vortex detector we get increased levels of detection. It should also be noted that the areas where we most want to evaluate the phase gradients (areas near the centre of the vortex) are the most difficult to measure due to low light levels. This also means that they are the slopes which are effected most by the large, turbulence-induced aberrations.

3. Laboratory experiments

3.1. Experimental method

An experimental setup was built in the laboratory to conduct practical tests of the theory outlined above. A schematic of the setup is shown in Fig. 3. The system is configured with two parallel beams which, after coming through a spatial filter from the laser, are split by a beam-splitter (BS1). One of these beams is a reference beam which is sent through a 4f system to the second set of beamsplitters (BS3,BS4), through the optical trombone and back into the wavefront sensing arm. The purpose of the optical trombone is to allow for different propagation distances along what is otherwise the same optical path. When Shack-Hartmann wavefront measurements are being carried out this beam is blocked off somewhere in the path before it rejoins the Spatial Light Modulator (SLM) beam.

The second beam going from BS1 goes through another beamsplitter (BS2) and is directed onto the SLM. After being reflected from the SLM the beam comes back through BS2 and into a 4f system, which has an iris placed at the focus to select only one diffraction order of the SLM. It then travels to BS4, where it rejoins the reference beam, and through the optical trombone. It then proceeds back through into the wavefront sensing arm. In the wavefront sensor itself a lens (L8) is used to expand the beam slightly to ensure it passes through the active area of the Shack-Hartmann lenslet array. Due to the relatively short focal length of the lenslet array a relay lens (L9) is needed to focus the beamlets on to the camera.

Another arm is used for interferometry; this arm starts with a flip mirror (FM1). When interferometric images are needed this flip mirror is dropped into the path of the beam before the wavefront sensing arm. It redirects the beam into lenses L10 and L11 which demagnify the beam size onto the CCD used for interferometry. When measurements involving the Shack-Hartmann sensor are needed FM1 is flipped out of the beam path and it continues into the wavefront sensing arm.

 

Fig. 3. Schematic of experimental setup used in the laboratory.

Download Full Size | PDF

To detect optical vortices in atmospheric like conditions we use the SLM to generate intensity scintillation and phase variations in a field typical of those encountered after the propagation of a laser beam through the atmosphere. Turbulence degraded optical fields were first generated by numerical simulations and the final field after the numerical propagation was encoded and applied to the SLM. In the simulations an initially plane wave is passed through a number of phase screens randomly generated using Kolmogorov statistics, equally spaced over the propagation path. The plane wave was propagated between phase screens by multiplying the Fourier transform of the optical field by the Fresnel transfer function [33, 34].

To ensure that an optical vortex is present in the final optical field, which is applied to the SLM, the initial plane wave is seeded with an optical vortex of known sign. This vortex will stay in the field over the total distance of the waves simulated propagation through the atmosphere. In general, a phase map with canonical vortices at locations (X_k,Y_k) with charge m_k can be generated using the following product wave function, where K is the index of the vortex. In this experiment we seeded the initial plane wave with a single vortex, (X ₁,Y ₁), at the centre of the field:

with Q being the resulting phase and X and Y are the coordinates of an array on a cartesian grid and X_k and Y_k defining the positions of the zero crossings. These simulations and all of the other analysis code is being implemented using MATLAB.

The desired amplitude and phase variations in the field are encoded into a binary phase mask of 0 to π radians and applied to the SLM [35, 36]. The SLM used is a ferroelectric liquid crystal, binary phase Spatial Light Modulator from Boulder Nonlinear Systems. It has 512 × 512 pixels with an active area of 7.7mm × 7.7mm and toggles between the true and inverse images imprinted once every cycle to prevent hardware damage. The phase and amplitude information is carried in the first diffraction order of the phase mask. In our system Fig. 3 the iris, between L4 and L5, spatially filters the +1 order for further propagation through the system.

The position of the optical vortices produced by the simulations are known and therefore we can test the efficacy of the branch point potential method against the interferometry technique of viewing vortices. The interferometric fringes of a vortex beam combined with a reference beam will have a fork (or a branch) present. It will be seen as one fringe splitting into two separate fringes, by pinpointing where this happens on the interferometry image and comparing it to the location given by the branch point potential method we can be certain of the veracity of the branch point potential method.

3.2. Experimental results

In these set of results we will demonstrate how the branch point potential methods capacity to detect vortices is improved by using the slope discrepancy technique alongside it. The differences in detection rates between two particular levels of turbulence will be highlighted.We will also show interferometric images of optical vortices in turbulent fields and their accompanying detection in the potential plots.

One of the initial results we noticed in the course of the experiments was that the use of the slope discrepancy greatly increased the chances of a correct detection by the branch point potential method in the presence of turbulence. We must first define what a correct detection of a vortex entails. A correct detection is when the vortex is seen in the correct position in the field and with the correct sign along with providing no other false detections in the field. Since we have seeded the turbulence simulations with an optical vortex and then used this as the pattern imprinted on the SLM, we know both the position and sign of the vortex. By examining the simulations we can also check if any other vortices are present and therefore rule any other detections as false.

As we can see from Fig. 4 the use of the slope discrepancy technique improves the detection rate approximately threefold for the scintillation indices of both σ ² _I = 0.154 and σ ² _I = 0.652, where σ ² _I is the scintillation index given by

 

Fig. 4. Bar chart showing the ability of the branch point potential method to detect vortices using both the raw Shack-Hartmann slopes and the slopes discrepancy technique.

Download Full Size | PDF

and I is the irradiance of the beam at the receiver.

Generally a scintillation index of σ ² _I > 1 is described as strong turbulence [37] and therefore the scintillation encountered here can only be described as being in either the weak or moderate turbulence regimes. This means that few vortices will be created naturally by the numerical simulation of propagation through the atmosphere. Those that have been created are very close together (ie. within one Shack-Hartmann lenslet) and do not effect the Shack-Hartmann slopes. The reason that the simulations are initially seeded with one vortex is to ensure that there is a detectable vortex in the final optical field, post propagation. At higher turbulence levels (ie. in the strong regime) one would expect to see more naturally occurring vortices and with the pairs spaced further apart, this has not been investigated in this paper.

For σ ² _I = 0.154 the detection rate goes from 0.23 for the raw Shack-Hartmann slopes to 0.63 for the slope discrepancy technique. For σ ² _I = 0.652 the detection rate goes from 0.07 for the raw Shack-Hartmann slopes to 0.21 for the slope discrepancy technique. This is a result for a comparison of 30 fields, using a numerical simulation of 15 phase screens employing Kolmogorov statistics with a λ of 633 nm and a 2.8 km propagation path, for both σ ² _I = 0.154 and σ ² _I = 0.652. A failure is defined as either a non-detection or as a detection along with false positives, in the case of σ ² _I = 0.154 most of the failures are due to false detection. The reason for the drop in detection rates when going into a higher turbulence regime can be explained by the existence of more phase aberrations and a higher degree of scintillation.

Using the interferometer arm of the system we can obtain images of how the interferometric fringes for a turbulent vortex beam appear. Here we show two different fringe images of optical vortices placed at the same location but with different turbulence levels, Fig. 5(a) and Fig. 6(a). We also have the maps that are output from the vortex detection code (using the slope discrepancy technique) for both of these, Fig. 5(b) and Fig. 6(b). It is worth pointing out that the lenslet array does not image all of the SLM, it only images the central portion of it. They show an obvious valley indicating the presence of a negative vortex at that location. Once a valley (peak) goes below (above) a certain threshold value then it is deemed to be a vortex, however since this threshold is set by heuristic methods it is not ideal. An obvious improvement to the branch point potential method would entail finding a more rigorous method to determine vortex detection.

 

Fig. 5. The interferometric fringes of a turbulent vortex beam for σ ² _I = 0.154 (a) and the detected vortex shown as a contour plot (b). The vortex is detected as it is the only peak/valley above/below the threshold used to separate a vortex from noise.

Download Full Size | PDF

 

Fig. 6. The interferometric fringes of a turbulent vortex beam for σ ² _I = 0.652 (a) and the detected vortex shown as a contour plot (b). The vortex is detected as it is the only peak/valley above/below the threshold used to separate a vortex from noise.

Download Full Size | PDF

4. Conclusion

In this paper we have shown that it is possible to experimentally locate optical vortices in a turbulent optical field using a Shack-Hartmann wavefront sensor. The branch point potential method is used to detect the vortices and it’s effectiveness is improved approximately threefold by using the slope discrepancy technique along with it, as seen in Fig. 4. From the same figure it is also obvious that the detection rate falls with increased turbulence, this is presumably due to the higher scintillation and larger phase aberrations across the optical field at these turbulence levels.

We also use the Hotelling observer, a prewhitening matched filter, to give a better evaluation of Shack-Hartmann spot positions. This improves the wavefront slopes and therefore aids in the ability to detect vortices. It is also shown, Table 1, to be more capable at spot detection and centroiding than a straight matched filter or a common centroiding algorithm. It also retains the ability to do good spot detection at low light levels, which is important when dealing with spots close to the vortex centre.

We feel it has been necessary to improve the Shack-Hartmann spot positioning to improve vortex detection in turbulent optical fields. However it is equally obvious that work should continue to improve the ability of the branch point potential method to detect optical vortices. This is especially the case when higher turbulence levels are considered, where more vortex pairs are created. These improvements should include a more rigorous and enhanced approach to vortex determination from the potential maps instead of the thresholding method currently used. Other methods to reduce the levels of noise in the potential maps are also being investigated, particularly in relation to edge effects.