Wave characterisation and aberration correction using hybrid direct search

There are a number of applications of shaped beams that require high accuracy reproductions of a precise wavefunction of light, such as: beam modes [1], knots [2], perfect vortices [3], and optimised trapping fields [4], vector polarized beams [5, 6] and a number of others [7]. Accurate retrieval of phase and amplitude is essential for appropriate and precise construction of structured light. To achieve this, the point spread function of the light must be determined [8–14]. One scheme to achieve this at the target plane is to consider separate regions of uniform phase and amplitude on a mode converter, such as a spatial light modulator (SLM), which enables an approximation of the phase and amplitude of a wavefunction to be found [11, 15–17]. These methods leverage the fact that correlations between simple closed regions ('patches') of the projected wave function are easy to measure. However, the strength of a correlation depends on the integrated amplitude of each patch which can depend on amplitude and phase of the wavefunction. Other schemes such as IMPACT iterate towards optimisation of the point spread function at scanned points within a sample [18]. Whatever method that is used to sample the projected light must thus be sensitive to potentially greatly varying light levels. The main problems we want to address are related to the fact that singular points in the field are a non-local phenomenon [19–21]. The schemes [11, 15–17] we have previously mentioned, perform poorly in simple regions where the correlation strength is reduced by low light levels, such as around vortexes. Assuming that the transformation (co-ordinate change or inverse) back to the mode converter apparatus can be achieved, these problems can be addressed by simultaneously finding a large part of the wavefunction. Here we have chosen to use Fourier analysis of a continuous time sequence of connected, but irregularly sized correlation regions to find the a large portion of the wavefunction of light from a single set of measurements.

In order to determine a wavefunction it must first be sampled in either space or time. Modal decompositions are one way the wavefunction can be characterised. This can take the form of a fundamental beam mode basis, such as Laguerre–Gaussian or Hermite–Gaussian modes [22–26], or it can be carried out in terms of a decomposition into momentum or position subspaces as we have previously discussed in the introduction. In an optical system these components are coupled through diffraction and can be measured as spatial correlations. One of the fundamental problems when using spatial correlations to perform tomography of a field is that the measured outputs are always inseparable product states. If the wavefunction is approximated in each region as a single complex number then a correlation measurement can be performed to find it [11]. These patches can be made irregular in size and shape to allow correlations where the component fields have similar amplitudes. To investigate the correlations between only two patches, a pattern can be added to the rest of the plane to allow spatial filtering of the targeted modes. There is an alternative strategy—simultaneously use the entire modulator and filter correlations in the temporal frequency domain over space.

Figure 1 depicts regions of a wavefunction emitted by a beam shaping device labelled with (0,1,3,5,7) to represent number of times the phase is swept by 2π. If the phase of region 0 (figure 1(a)) was kept constant in time and 1, 3, 5, and 7 would correspond to the frequency, then in the optical Fourier plane the complex amplitudes of the integer temporal Fourier components $\left(n = \left\{1,3,...,N\right\}\right)$ would be given by:

$$\begin{equation} f_n = E_0^*E_n, \end{equation} \tag{ 1 }$$

for complex amplitudes E₀ and E_n of the Fourier transformed wavefunction. Odd couplings in the Fourier spectrum can only occur with the zero frequency component patch. The even frequency difference components can be found as sums of conjugate products of the odd components of the spectrum multiplied by $E_0^*E_0$. There is a caveat; this scheme, by default, only works under the condition that there is sufficient light in the reference (zero frequency) region. If $E_0^*E_0 = 0$ then the odd frequency components disappear at these points in space of the Fourier spectrum. Theoretically, it remains possible to calculate by setting one of the remaining components to have $E_n^* = E_n$. Given this, we still need to determine the precise wavefunctions comprising the product states. There is an extremely large number of combinations of phase with amplitude that reproduce the correlations measured in an experiment. Without more information many schemes that could find the spatial properties of amplitude and phase will likely fail to converge to a sensible result. To avoid this, we would typically assume uniform properties, such as the average amplitude and phase across each patch. In principle, this seems to be a good compromise, and it works well with tiled rectangular patches that are precisely in the Fourier plane. This does not work so well for irregularly tiled patches or for devices in some other plane because both the width and peak position of the interference can shift. Central to our problem is the ability to consider the wavefunction in the target and device planes. To do that we will employ our own inverse transform. The inverse transform problem space greatly reduces in size and will more likely converge to the wavefunction at the device plane if we use the fact that each patch occupies a definite space with a definite relative position to the zero frequency component. If there is too little information or if the inverse transform is not well approximating the system the optimisation may not be possible.

Zoom In Zoom Out Reset image size

To retrieve a wavefunction a convergent scheme needs to be chosen. In our case, each patch must independently contribute to the correlation measurements. We have almost no information about one of the planes, in this case, one within the sample [11, 15–17]. The computational methods for either wavefunction or amplitude retrieval require that information about the device and measurement regions are known [27–31]. We want to construct a general algorithm to bridge the gap between the first and second case so that we can resolve fine details of the field, such as the amplitude and phase near any singularity present.

The correlation measurements give us only the product of amplitude states as expressed in the previous section. To isolate the absolute value of the amplitudes, we first divide through by the magnitude of the zero frequency component. This yields the amplitude in whichever plane we chose to optimise. However, we only know the phase of each frequency component as a difference from the zero frequency components' phase. The task for the algorithm is to find which value of phase needs to be added to both the zero and non-zero frequency components correspond to the appropriate region of the device. Represented as equations, the auxiliary product of each patch, n, is:

$$\begin{align} f_n(x,y)& = E_0^{*}(x,y)E_n(x,y)\nonumber\\[7pt] f_n(x,y)& = |E_0(x,y)||E_n(x,y)|\exp{\mathrm{i}(\arg{E_n(x,y)}-\arg{E_0(x,y)})} \end{align} \tag{ 2 }$$

$$\begin{align} \psi^{(a)}_n(k_x,k_y;\phi)& = \hat{\mathcal{O}}\left(\frac{f_n(x,y)}{|E_0(x,y)|}\right)\nonumber\\ & \quad + f_n(x_a,y_a)\left(\exp{\mathrm{i}\phi}-1\right)\exp{\mathrm{i} \mathbf{k}\cdot{\mathbf{r}_a}} \end{align} \tag{ 3 }$$

where E_n are complex field amplitudes that are present at the detector, $\hat{\mathcal{O}}$ is the Fourier transform operator (${x,y}\rightarrow{k_x,k_y}$) that returns the coordinates back to the device, normalised to a point spacing of one, the second term of equation (3) can be added to a pre-computed Fourier transform to check the effect of a change of phase on the device, φ, at point $\mathbf{r}_a = (x_a,y_a)$. $\psi^{(a)}_n$ is the unnormalised wavefunction at the device plane as function of changing the phase difference. To measure the effect of changing φ we define a cost (or value) function for the patch:

$$\begin{equation} c^{(a)}_n(\phi) = \frac{\int\int \left|A_n(k_x,k_y) \psi^{(a)}_n(k_x,k_y,\phi)\right|^2\mathrm{d}k_x\mathrm{d}k_y}{\int\int \left|\psi^{(a)}_n(k_x,k_y,\phi)\right|^2\mathrm{d}k_x\mathrm{d}k_y}, \end{equation} \tag{ 4 }$$

that should be maximized, where $\psi(k_x,k_y,\phi)$ is the Fourier transform of the current patch of the wavefunction for a change of phase at a point and $A_n(k_x,k_y)$ is a mask function that is 1 inside and 0 outside the patch, n. Given that only a single point is modified the cost function must result in a sinusoidal variation about a mean value as a function of φ, e.g. of the form: $c^{(a)}_n(\phi) = a\cos{(\phi+b)}+c$. It can thus be accurately fitted in a Fourier series to guarantee uniqueness a minimum of five equal steps of φ are needed to determine the maximum and minimum phase in the domain $[-\pi,\pi-2\pi/N]$:

$$\begin{align} z& = \frac{1}{N}\sum_k c^{(a)}_n(\phi_k) \exp{-\mathrm{i} \phi_k}, \end{align} \tag{ 5 }$$

$$\begin{align} \bar{c}^{(a)}_n(\phi) & = 2|z|\cos(\phi+\arg{z})+\frac{1}{N}\sum_k c^{(a)}_n(\phi_k), \end{align} \tag{ 6 }$$

where $\bar{c}_n$ is an accurate fit of c_n that resolves the cost function for all phase differences. The optimisation is one where the cost function should increase as each patch approaches the value of the wavefunction. The most straight-forward way of doing this is to find the product of all c_n : $C^{(a)}(\phi) = \prod_n c^{(a)}_n(\phi)$. The maximum of this function can be found numerically.

For the test case shown in figure 1 the iterative algorithm worked on a sequence of points sorted in descending order of amplitude. The largest amplitude components have the biggest contribution to the wavefunction. The algorithm analyses an increasing fraction of the total number of elements until every one is optimised on the last iteration. The results from this test case are shown in figure 2. The first step in the reconstruction is to obtain the amplitude of the zero frequency patch, as seen in figure 2(a). An initial phase, figure 2(b) is chosen. Figure 2(c) depicts the change in the aggregate cost function as the algorithm is iterated, showing an overall increase in fitness. After the algorithm has finished, the phase of the zero frequency patch is obtained, shown in figure 2(d), which gives an excellent degree of correspondence to the ground truth depicted in figure 2(e). The results show an excellent reconstruction of the original wavefunction (figures 1(a) and (b)) in figures 2(f) and (g). Importantly, it captures the presence of the singularity in the wavefunction which would be difficult to capture for a patch with assumed uniform amplitude and phase.

Zoom In Zoom Out Reset image size

In summary the proceedure is:

• (a)  
  Perform intensity measurements (e.g. trace in figure 1(e)) on the temporally modulated regions (e.g. with frequencies depicted in figure 1(a)) on a pixelated spatial light modulator, taking a selection of pixels or points of detection. Create a trial wavefunction at the beam shaping device.
• (b)  
  Calculate the device region cost function for the current point on the detector (e.g. equation (4)). Find the optimal phase (stationary point) from the global cost function.
• (c)  
  Calculate the new trial wavefunction using a inverse transform/basis. Goto (b) until the cost function stabilises.

With the initial low noise conditions used in the simulation the method had a high probability of success. In experiments there are limits to our ability to accurately detect signals due to the presence of various types of noise which will in turn affect our ability to recover the wavefunction. The experimental apparatus used to examine limits of the method is shown in figure 3. The experiment consists of an infrared laser (IPG Photonics, YLR-10-1064-LP) with a roughly Gaussian beam profile and polarisation optics designed to control the intensity of the beam. Light enters the experimental section consisting of a low reflectivity pelicle beamsplitter so that the incident laser power can be monitored. The light reflects of a SLM (Medowlark optics, Spatial Light Modulator—Large 512 × 512, dvi) through a lens onto a camera (Teledyne LUMENERA, Lw115) to measure the Fourier plane of the SLM.

Zoom In Zoom Out Reset image size

Our first experiment was designed to find a wavefunction reconstruction similar to the model that could be obtained in a lab setting. To perform this comparison, we first found the phase function of the SLM to produce a look up table for the SLM using the methodology described in the appendix. We created a similar initial condition in our experiments to that used in the example shown in the iterative scheme section. We split up the spatial regions of the spatial light modulator such that the zero order patch is large but not completely dominating the projected light. The most important part of the process in the optimisation procedure is to find the reference phase for the zero frequency patch. As this experiment uses similar parameters to the computational example a reasonable measure of success would be to obtain a similar amplitude and phase of the modelled zero frequency component. Figure 4(a) shows the phase-only diffraction pattern send to the SLM. Figure 4(b) shows the phase only component projected to the position targeted by the pattern in (a). As was the case with the simulation example, a charge one phase singularity can be seen which would make it difficult to reconstruct if uniform phase and amplitude was assumed for each discrete element.

Zoom In Zoom Out Reset image size

Figure 5 shows results from this experiment. We obtained a time series long enough to distribute the random frequency contributions of air currents and light saturation of the camera. Using an aggregate of the even coupled frequency terms we found the average intensity of the zero frequency component. The algorithm converged in a similar manner as to how it did in the simulated model, the results can be seen in figures 5(b)–(d). Figures 5(e) and (f) show the wavefunction reconstructed by the algorithm yielding excellent agreement with the incident beam including the phase singularity. The other wavefront curvature on the wavefunction shown in figure 5(f) can be attributed to the aberration of the light accrued as it propagates to the camera.

Zoom In Zoom Out Reset image size

Using the recovered wavefunction we applied the conjugate phase to the SLM to remove the singularity from the original wavefunction and perform aberration correction, producing a symmetric and tight focus. The method not only corrects the point spread function of the incident beam (figures 6(a) and (d)), it also removes the vortex (figures 6(b) and (e)) without any significant additional artefacts (figures 6(c) and (f)).

Zoom In Zoom Out Reset image size

An alternative arrangement of device regions are shown in figures 7(a) and (b). Using the same criteria in the data processing as before we obtained the amplitude and phases shown in figures 7(c)–(f). To increase the speed of the analysis and to improve convergence, low intensity pixels at the noise level of the camera were not included in the calculation. As a result, the agreement between the two patterns is excellent. In both instances the algorithm produced wavefunctions with different global phases. To allow a simpler direct comparison between the two recovered wavefunctions the phase of a pixel near the maximum of intensity was zeroed in both figures 7(e) and (f). This data shows that the method converges to very similar results for different patch configurations. Coupled with the demonstration of the aberration correction containing a phase singularity in figure 6, there is support to the claim that the the reconstruction of patch regions using a global cost function is sufficient to reconstruct the wavefunction projected on to the camera.

Zoom In Zoom Out Reset image size

We have designed a scheme to take correlation measurements from a small number of patches that represent a wavefunction. We reconstruct the amplitude and phase of the wavefunction traversing an optical system to the point of measurement with the aid of a global cost function. The method works best when good coverage of detection in the target region is possible. When low-signal high noise data is removed from the reference data for reconstruction the detail of the recovered wavefunction reduces. All data and codes are available online [32].

This method shows promise in reconstructing the wavefunction of light travelling through an optical system. There are some good and useful features of the method, the first of which is that compared to the method [11], a smaller number of time points are needed to characterise the wavefunction. Second, the method developed here can also find singularities in the wavefunction without introducing extra aberrations from the assumption of constant phase and amplitude in each patch.

The pattern displayed on the beam shaping device has a maximum frequency, $f_\mathrm{max}$, and thus the theoretical minimum number of sample points is $N_\mathrm{min} = 2f_\mathrm{max}+1$. Running the experiment at this number of time points did not result in an accurate reconstruction of the wavefunction. Review of the captured data revealed two features of the data that could produce inconsistent results. The first is that despite the apparent linearity of the camera, the minimum signal was hardware limited to a pixel brightness level of 3 or 4 out of the 255 available. In a Fourier transform each discontinuity results in ringing in the transformed data set. Second, we saw random movement of the pattern caused by air currents within the lab. This too will add noise to the Fourier transform. A combination of two improvements were used to overcome these limitations. First, the quality was improved by picking a capture time about ten times the maximum frequency determine from the number of patches, which reduces statistical noise by about a factor of ≃ 3. The second was to remove contribution of noise dominated pixels when performing the optimisation. This speeds up processing at the expense of reducing the number of resolvable features at the device. In our particular implementation we removed all points in the Fourier transform with complex amplitudes smaller than the dark level of the camera.

In our application of the method we found that increasing the number of scanning regions would typically make optimisation more difficult in experiments. One cause of this is that the optimisation cost function becomes weaker with increasing number of elements. There is no limit to the extent of cost functions to choose and further study to find more appropriate cost functions for this problem should improve convergence. Secondly, the noise in all of the amplitude products will increase as the light present in each patch reduces. Thus, a careful balance of size and light level needs to be struck to maintain the advantage of optimisation over differently shaped correlation regions.

There are a number of areas where further developments are possible for wavefunction retrieval using this algorithm. Like other extant methods, our method does not account for complex scattering interactions that change as a function of time as it explicitly assumes a simple transfer function (Fourier transform). An adaptation of the method using filtering methods to measure particular sub-spaces of scattering may allow for good characterisation of these interactions so that it can be applied to imaging through scattering media and determination of other transfer functions.

It may be possible to produce a variation of this technique for use in the far-field zone to correct focal spots in microscopes using light scattered by small particles or nano-antennas by using a sparse reconstruction algorithm [33, 34] which will work for low-order aberrations. We would also like to see if the cost function we have chosen has uses for other systems and for similar problems, such as wavefunction design. Potentially, pixel crosstalk and fringe effects may influence the measurement when larger number of correlation regions are used [35, 36] and may need to be investigated. The scheme outlined here could be implemented on GPU to improve computation time using a number of light modulation packages, for example, Bowman et al [37] and Ploschner et al [38].

The method used in our work reduces the number of time steps needed to acquire good high detail of the wavefunction because it assumes the fact the field at all points on the device are related. The reduction to the number of measurements is offset by an increase in the amount of processing. Taking computation into account means that the combined time needed to perform a correlation measurement and process it on a desktop computer (ca. 2014) is similar to the assumption of patches with uniform complex amplitude. This method will gain practical performance benefits as computing power continues to improve.

The data that support the findings of this study are available upon reasonable request from the authors.

This research was supported under Australian Research Council's Discovery Projects funding scheme (Project Number DP180101002).

The authors declare that there are no conflicts of interest related to this article.

The scheme used to find the look up table on the high-efficiency SLM used in the experiment utilizes the fact that for a regularly alternating phase pattern of two values the intensity of the reflected spot and one of the first/second order diffracted have the relationships [39, 40]:

$$\begin{align}I_0& = \cos^2{\frac{\phi}{2}}+I_\mathrm{off}, \end{align} \tag{ A.1 }$$

$$\begin{align}I_1& = \left(\frac{2}{\pi}\right)^2\sin^2{\frac{\phi}{2}}, \end{align} \tag{ A.2 }$$

$$\begin{align}I_2& = \left(\frac{2}{3\pi}\right)^2\sin^2{\frac{\phi}{2}}, \end{align} \tag{ A.3 }$$

where the phase difference between the two phase levels is φ, I₀ is the reflected intensity, and $I_\mathrm{off}$ is component of light that is reflected off the SLM. Assuming that $I_\mathrm{off}\ll I_0$ for two unknown phase levels x₁ and x₂, the phase difference can be found from the ratio of intensities (equations (A.1) and (A.2)):

$$\begin{align} \frac{I_1}{I_0}& = \left(\frac{2}{\pi}\right)^2\tan^2{\frac{\phi}{2}}, \end{align} \tag{ A.4 }$$

$$\begin{align}\phi& = 2~\arctan\left({\frac{\pi}{2}\sqrt{\frac{I_1}{I_0}}}\right). \end{align} \tag{ A.5 }$$

Digital cameras are known to have non-linearities in the infra-red region. To get an accurate quantification of the camera response an intensity map to pixel signal level may be required. Further, if the zero level is known, then when $I_1 = I_0$, we obtain a target phase of $\phi = 2~\arctan{\frac{\pi}{2}}$.

This does not completely solve the look up problem as we also need to obtain a reference to represent the component of light that is reflected off the SLM. Not taking account for this offset results in a bias on the target phase difference which takes the form of some unknown factor. $I_\mathrm{off}$ should be determined to get an accurate scan. However, a scan over the region can be approximated by a single iteration of the scheme derived from equations (A.1) to (A.5). A second iteration with this offset can then be used to find a more accurate version of the look up table, potentially saving a significant amount of time.

If the mapping to linear detector response is known then it is straight forward to take a scan of the diffracted light as a function of pixel level. The relationship of intensity between the difference in pixel level and phase is

$$\begin{equation} \frac{\mathrm{d}I}{\mathrm{d}\phi}\simeq\frac{1}{k}\frac{\mathrm{d}I}{\mathrm{d}x}. \end{equation} \tag{ A.6 }$$

Importantly, we note that these expressions are first order for small perturbations. Thus the intensity fluctuation is proportional to both the pixel level and phase. By taking the cumulative differences in intensity to the peak value we obtain an explicit relationship between the actual phase level (π) and cumulative but otherwise unnormalised phase shift as a function of pixel level. From this we can invert the expressions to find a look up table between the pixel level and phase. Figure A1 shows how a calibration can be found from the equations and data presented. The calibration uses two facts: (a) If

$$\begin{equation} I_0^{\mathrm{max}} = \left(\frac{\pi}{2}\right)^2 I^{\mathrm{max}}_1, \end{equation} \tag{ A.7 }$$

then the phase shift between two pixel levels, x₁ and x₂, is π (represented by the black arrows), and (b) for a small variation in phase, $x_2 = x_1+\Delta$ the difference in intensity is $\Delta I_1\ll I\,_1^{\mathrm{max}}$, the phase difference caused by Δ is close to proportional to difference of intensity (represented by the green arrows in the inset of figure A1).

Zoom In Zoom Out Reset image size

Using the data in figure A1 and a starting level, x₀, we project from the diagonal, where phase is equal, to the curve defined by the maxima of I₁ (or minima of I₀) twice to obtain a scale for 2π phase difference. The precise phase function curve is then found from the differences in intensity when the pixel level changes by a small constant amount. To find the look up table we add the intensities as a function of x,

$$\begin{equation} S(x) = \int_{x_0}^x \frac{\mathrm{d} I_1}{\mathrm{d}x^{^{\prime}}}\mathrm{d} x^{^{\prime}}, \end{equation} \tag{ A.8 }$$

so the phase function becomes:

$$\begin{equation} \psi(x) = \pi\frac{S(x)-S(x_0)}{S(x_b)-S(x_a)}, \end{equation} \tag{ A.9 }$$

where x₀ is the starting pixel level, x_a is the phase difference defined earlier, φ = π, and x_b is φ = 2π. If ψ(x) is monotonic then the look up is Φ(x) = ψ⁻¹(x). The inversion can be performed numerically though not analytically in general. From the data obtained we find a normalised phase function as shown in figure A2.

Zoom In Zoom Out Reset image size