## Abstract

We report a mathematically rigorous technique which facilitates the optimization of various optical properties of electromagnetic fields in free space and including scattering interactions. The technique exploits the linearity of electromagnetic fields along with the quadratic nature of the intensity to define specific Optical Eigenmodes (OEi) that are pertinent to the interaction considered. Key applications include the optimization of the size of a focused spot, the transmission through sub-wavelength apertures, and of the optical force acting on microparticles. We verify experimentally the OEi approach by minimising the size of a focused optical field using a superposition of Bessel beams.

© 2011 Optical Society of America

## 1. Introduction

The decomposition of fields into eigenmodes is a well established technique to solve various problems within physical sciences. The most prominent example is the Schrödinger’s equation within the field of quantum mechanics, where energy spectra of atoms are determined via the eigenvalue spectra and associated wavefunctions of the Hamiltonian operator. Indeed, electron orbits are eigenmodes of the energy, angular momentum, and spin operators [1] and as such they deliver fundamental insights into the physics of atoms. Within classical mechanics, modes of vibration of music instruments give, for example, their resonant frequencies while their spectrum is associated with the shape of the instrument [2]. In the optical domain, mode decomposition is used in order to describe light propagation within waveguides [3], photonic crystals [4], optical cavities [5], laser resonators [6], and the optical forces on Mie-sized particles [7]. In the case of waveguides and photonic crystals, for example, eigenmodes describe electromagnetic fields that are invariant in their intensity profile as they propagate along the fibre or crystal. Additionally, these modes are orthogonal and as such light coupled to one of these modes remains, in theory, in this mode forever. This optical mode decomposition can be expanded to include additional properties such as orbital and spin angular momentum [8]. All-together, the eigenmode expansion method is a well-established method for the representation of the propagation of optical fields.

In this paper, we report a novel method which we term “Optical Eigenmodes (OEi)” which represents a generalization of the powerful concept of eigenmode decomposition going beyond the propagation properties of light. Crucially, we show that eigenmode decomposition is applicable to the case of any quadratic measure which is defined as a function of the electromagnetic field. Prominent examples of optical quadratic measures include the energy density and the energy flux of electromagnetic fields. The OEi method makes it possible to describe an optical system and its response to incident electromagnetic fields as a simple mode coupling problem and to determine the optimal “excitation” for the given measure considered. Intuitively, a superposition of initial fields is optimized in a manner that the minimum/maximum measure is achieved. For instance, the transmission through a pinhole is optimized by maximizing the energy flux through the pinhole [9].

From a theoretical perspective, the OEi optimization method is mathematically rigorous and may be distinguished from the multiple techniques currently employed ranging from genetic algorithms [10] and random search methods [11] to direct search methods [12]. The major challenge encountered in any such approximate optimization and engineering of optical properties is the fact that electromagnetic waves interfere. As such the interference pattern not only makes the search for an optimum beam problematic but crucially renders the superposition found unreliable, as the different algorithms may converge on different local minima which are unstable with respect to the different initial parameters of the problem. In contrast, our proposed OEi method yields a unique solution to the problem and directly determines the optimum (maximal/minimal) measure possible.

In the first part of the paper, we introduce the OEi method and show its properties in a general context of optimizing the quadratic measures of interfering waves. In the second part, we apply the OEi formalism to minimize the focal spot size and discuss the appearance of superoscillating fields. For these applications, we describe, respectively, the electromagnetic field as a superposition of scalar Laguerre-Gaussian beams, vectorial Bessel beams or more general plane waves within the angular spectral decomposition representation of light. In the third part of the paper, we report a particular experimental implementation of the OEi method using computer controlled spatial light modulators to squeeze the spot size of a superposition of Bessel beams. In the last part of the paper, the method is applied in numerical 3D modelling to determine the OEi yielding the largest transmission through a sub-wavelength aperture and the largest optical force on a micrometer sized particle. The paper concludes with a discussion of the particular results obtained and with general comments on the versatility of the OEi method to a wide range of problems. A short annex compares the convergence properties of the OEi approach with standard phase front correction methods.

## 2. Method

The OEi method exploits both the linearity of Maxwell’s equations and the quadratic dependence of light-matter interactions on the electromagnetic field {**E**, **H**} where **E** and **H** denote the electric and magnetic field vectors, respectively. Table 1 provides a list of common examples of such interactions. These interactions may be written in a general quadratic matrix form

*A*) labels the light-matter interaction defined in Table 1. The vectors

**a**and

**a**

^{†}are comprised of the superposition coefficients

*a*and their complex conjugates, respectively. The elements are constructed by combining the respective fields {

_{j}**E**

*,*

_{j}**H**

*} and {*

_{j}**E**

*,*

_{k}**H**

*} for*

_{k}*j*,

*k*= 1 . . .

*N*. More precisely, we have:

**M**

^{(}

^{A}^{)}defines a spectrum of real eigenvalues ${\lambda}_{k}^{(A)}$ and associated eigenvectors ${\mathbf{v}}_{k}^{(A)}$. Each of these eigen-vectors corresponds to a superposition of fields {

**E**

*,*

_{j}**H**

*} termed here Optical Eigenmode (OEi). Crucially, we may now extremize the light-matter interaction considered; that is, the extremal eigenvalue ${\lambda}_{\text{ext}}^{(A)}$ and the associated eigenvectors ${\mathbf{v}}_{\text{ext}}^{(A)}$ deliver the superposition of fields $\{{\mathbf{E}}_{\text{ext}},{\mathbf{H}}_{\text{ext}}\}=\left\{{\sum}_{j=1}^{N}{v}_{\text{ext},j}^{(A)}{\mathbf{E}}_{j},{\sum}_{j=1}^{N}{v}_{\text{ext},j}^{(A)}{\mathbf{H}}_{j}\right\}$ which extremizes the interaction*

_{j}*m*

^{(A)}while keeping the total field contributions constant

**a**

^{†}

**a**= 1.

In this paper, we apply the OEi concept to minimize the size of a laser spot within a surface ROI. One way to define the spot size of a laser beam is by measuring (whilst keeping the total intensity constant) the second order moment *w* of its intensity distribution [13]. Crucially, *w* can be expressed in terms of *m*^{(0)} and *m*^{(2)} (see Table 1) or the respective matrix representations Eq. (1) as follows:

**M**

^{(0)}and

**M**

^{(2)}are termed the

*intensity operator (IO)*and

*spot size operator (SSO)*, respectively. According to Eq. (3), the minimum spot size is obtained by the OEi associated with the smallest eigenvalue

*λ*

^{(2)}of the SSO provided that the IO is simultaneously diagonalized and normalized to 1. Direct evaluation shows that this is precisely achieved by the combined OEi

**M**

^{(2)}and in the intensity normalised eigenbase ${v}_{k,j}^{(0)}/\sqrt{{\lambda}_{k}^{(0)}}$ of

**M**

^{(0)}.

For our proof-of-principle studies described in the remainder of the paper, we also applied the scalar version of the OEi method where a set of scalar fields *E _{i}* is considered in order to determine the IO and SSO as

#### 2.1. Smallest focal spot using Laguerre Gaussian beams

Using a superposition of LG beams we can minimize the size of a focal spot using the representation of the SSO (6) in cylindrical coordinates. It is important to note at this point that we only retain the intensity OEis whose eigenvalues are within a chosen fraction of total intensity. This is equivalent to considering only beams that have a significant intensity contribution in the ROI. Intuitively, the optimization procedure may perform so well that a spot of size zero is finally obtained if no intensity threshold is applied. Figure 1 shows the smallest spot superposition where we observe the appearance of sidebands just outside the ROI. These sidebands are a secondary effect of squeezing the light below its diffraction limit. It is these sidebands that decrease the efficiency of the squeezed spot with respect to the maximal possible intensity in the ROI as calculated via the IO. Using the ratio between these two intensities we can define the intensity Strehl ratio [14] for the SSO (see Fig. 2b). We remark that both, the spot size and the Strehl ratio, show resonances as a function of the ROI size. This can be explained by considering the number of intensity eigenmodes used for the spot size operator. Indeed, as the ROI size decreases, so does the number of significant intensity eigenmodes. Each time one of these modes disappears (step in Fig. 2), we have a sudden increase in the minimum spot size achievable accompanied with an enhanced Strehl ratio as we drop the most intensity inefficient mode. Overall, the Strehl ratios determined in our studies predominantly exceeded values of 1% even when spots were tightly squeezed. Therefore, the observed decrease of intensity is not to severe in terms of potential applications of squeezed beams such as optical manipulation and imaging.

On a final note, we remark that squeezing light below its diffraction limit may be associated with the effect of super-oscillations [15]. This refers specifically to the ability to have a local *k*-vector (gradient of the phase) larger than the spectral bandwidth of the original field. To visualize this effect, in the case of OEi spot size optimized beams, we have calculated the spectral density of the radial wave-vector for the smallest planar spot [16]. As shown in Fig. 3, this spectral density clearly identifies a spectral bandwidth (white background in Fig. 3). Regions of the beam which exhibit locally larger wave-vectors than the ones supported by this spectral band width correspond to super-oscillating regions. The local wave vector is defined as *∂ _{r}* arg(

*u*(

*r*)) where arg(

*u*) defines the phase of the analytical signal

*u*. In this particular case, we observe that super-oscillations occur in the dark region of the beam. Additionally, when the ROI is large compared to the Gaussian beam waist

*w*

_{0}, there are no super-oscillating regions. These only appear when the beam starts to be squeezed.

#### 2.2. Smallest focal spot using Bessel beams

The paraxial approximation employed above in the case of LG beams can be used to describe sub-diffracting beams but breaks down when beams are tightly focused. As a consequence we must consider full vectorial solutions of Maxwell’s equations. Here, we have chosen Bessel beams as a basis-set and determined the superposition of Bessel beams which minimized the spot size in a planar finite ROI. Note that the problem of the finite intensity of Bessel beams [17] is easily circumvented here due to the finite ROI size considered. The monochromatic electric vector field of the vectorial Bessel beam may explicitly be expressed as [18]

*k*=

_{t}*k*

_{0}sin(

*θ*) and

*k*=

_{z}*k*

_{0}cos(

*θ*) are the transversal and longitudinal wave vectors with

*θ*the characteristic cone angle of the Bessel beam.

**e**

*,*

_{x}**e**

*and*

_{y}**e**

*are the unit vectors in the Cartesian coordinate system. The parameter*

_{z}*ℓ*corresponds to the azimuthal topological charge of the beam while

*α*and

*β*are associated with the polarization state of the beam. The magnetic field

**H**was deduced according to Maxwell’s equations. Figure 4 shows a comparison between the Airy disk, the Bessel beam and the OEi optimized spot considering a numerical aperture of NA= 0.1. As in the case of the LG beams, squeezing the focal spot is accompanied by side bands and a loss in efficiency shown by the Strehl ratio (see Fig. 5).

## 3. Experimental OEi

#### 3.1. Experimental implementation of the OEi concept

To perform an experimental OEi optimization we have used the following setup: A HeNe laser beam is expanded and subsequently amplitude modulated by a spatial light modulator (SLM) display operating in conjunction with a pair of crossed polarizers. Analogously to a liquid crystal display on a computer or laptop monitor, the liquid crystal SLM display rotates the polarization of the incident light by an angle depending upon the voltage applied to the display pixels. The amplitude modulated beam is then imaged onto a second SLM display through a pair of lenses. This second SLM display along with a subsequent Fourier lens and aperture served to modulate the phase of the laser beam in the standard first order configuration [20]. The field modulations of interest were encoded as RGB images where the blue channel represented the amplitude and the green channel the phase modulation. The SLM controller extracted these information and applied the two channels to the respective panel. We have performed calibration measurements to ensure that both the amplitude and phase modulation exhibited a linear dependence on the applied 8-bit color value between 0 and 255. A CCD camera allowed us to record images of laser fields in the Fourier plane of lens 5.

The experiment consist in determining the OEi in the CCD camera plane whilst shaping and superimposing the test fields *E _{j}* in the SLM planes. In the following, we indicate the plane of interest by a

*z*-coordinate along the optical axis where

*z*=

*z*

_{1}and

*z*=

*z*

_{2}refer to the SLM and CCD camera plane, respectively. According to this convention we shape a set of test fields

*E*(

_{j}*z*

_{1}) =

*A*(

_{j}*z*

_{1})

*e*

^{iϕj (z1)}both in amplitude

*A*and phase

_{j}*ϕ*in the SLM plane, and the associated intensities

_{j}*I*(

_{j}*z*

_{2}) ∝ |

*E*(

_{j}*z*

_{2})|

^{2}are detected in the CCD camera plane. The amplitudes

*A*(

_{j}*z*

_{2}) were determined from these intensities by simply taking the square root. We used the three-step phase retrieval algorithm described in Ref. [21] to retrieve the phase modulations

*ϕ*(

_{j}*z*

_{2}). The determination of the phase and amplitude of the beam in the CCD plane allows us to numerically vary the ROI without redisplaying the test fields. Using these fields, the IO and SSO are determined according to Eqs. (5) and (6), respectively.

During the course of our experiments we verified the linearity of our optical system by performing a comparison between what we term the “experimental superposition (Exp-S)” and the “numerical superposition (Num-S)”. The Exp-S refers to the case where the set of OEi optimized superposition coefficients *a _{i}* is used to encode the optimized superimposed field onto the SLM. The CCD camera then detected the intensity

*I*

_{Exp-S}(

*z*

_{2}) corresponding to this encoded optimized field. The Num-S utilizes the fields

*E*(

_{j}*z*

_{2}), which were individually measured to assemble the OEi operators, in order to

*numerically*determine the intensity distribution as ${I}_{\text{Num}-\text{s}}({z}_{2})\propto {\left|{\sum}_{i=1}^{N}{a}_{j}{E}_{j}({z}_{2})\right|}^{2}$. Crucially, linearity is verified if

*I*

_{Exp-S}(

*z*

_{2}) =

*I*

_{Num-S}(

*z*

_{2}). This is indeed observed in our experiments as demonstrated in the following subsection which features a comparison of experimental and numerical intensity distributions.

#### 3.2. Results and discussion

In our experiments, we used *N* = 11 non overlapping amplitude ring masks with a constant phase modulation as fields of interest *E _{i}*(

*z*

_{1}). After propagation through the Fourier lens 5 (see Fig. 6) the resulting fields

*E*(

_{i}*z*

_{2}) form a set of Bessel beams. Figure 7(a) shows the largest ring modulation encoded onto the SLM with the resulting Bessel beam shown in Fig. (b). As this particular Bessel beam comes along with the highest NA compared to the Bessel beams created with smaller ring modulations, the beam shown in Fig. (b) exhibits the smallest central spot of all beams realized in our experiments. The spot size of the Bessel beam featuring the smallest core is denoted as

*w*

_{B}and used as reference for the measurements presented below. For comparison Figure 7(c) depicts a circular aperture which is encoded onto the SLM in order to observe the Airy disk (see Fig. (d)). The spot size of the Airy disk is approximately 1.5 times larger than the core of the reference Bessel beam as expected [22].

The results of the performed OEi spot size minimization are shown in Fig. 8 for different sizes of the ROI. To begin with, the comparison of the Num-S intensity distribution *I*_{Num-S}(*z*_{2}) (top row) and the Exp-S intensity distributions *I*_{Exp-S}(*z*_{2}) (bottom row) clearly reveals good agreement and thus verifies the linearity of our optical system as elucidated above. For completeness, the central row shows the Exp-S superposition in RGB format as encoded onto the SLM. The color code features a blue channel representing the amplitude modulation from 0 (black) to 1 (blue) and a green channel corresponding to phase modulations from 0 (black) to 2*π* (green). Next, we conclude from the measured relative spot size *w*/*w*_{B} that the spot size decreases if the ROI size is reduced. The reduced spot size is achieved at the expense of the spot intensity which is redistributed to a ring outside of the ROI similar to the theoretical results presented in section 2.2 and Fig. 4. Referring to the Exp-S data, for *R* = 7 pixel the spot size is reduced to 72 % of the size of the reference Bessel beam’s core and even further to 50 % for *R* = 4 pixel. The latter result is somewhat vague, though, due to the low spot intensity which may be truncated by the sensitivity threshold of the CCD detector and thus may appear smaller. However, our experimental results overall clearly verify the OEi concept applied to spot size minimization. Moreover, the results strongly suggest that the OEi optimization may indeed squeeze spots to the subdiffractive regime since the optimal superposition of Bessel beams not only beats the Airy disk but also the reference Bessel beam diameter.

## 4. Applicability of the OEi method to scattering interactions

In this section, we demonstrate, on the basis of a numerical study, how the OEi method can equally be applied in some scattering interaction processes. Indeed, the OEi method presented above is applicable to free space propagation and can be directly extended to linear scattering processes where the optical interaction is expressed as a quadratic form of the field. As in the case of the smallest spot operator, the OEi of these light matter interactions can be used to determine the electromagnetic field profile delivering the largest or the smallest interaction strength. In this section, we show two numerical examples illustrating the cases presented in Table 1. The numerical modelling is performed using a finite element method (Comsol) and solving the fully vectorial monochromatic Maxwell’s equations in 3D. The structures considered here are embedded in a larger computational domain surrounded by perfectly matched layers. Figure 9 shows the electric field amplitude |**E**| for the different cases considered. To implement the OEi method, we determine the matrix operator with the help of Eq. (2). Here, we use angular spectral decomposition [19] of the incident light field corresponding to a numerical aperture of NA=0.8.

In the first example (Fig. 9a–c), we use the intensity operator associated to the measure *m*^{(0)} to determine the largest transmission through a sub-wavelength aperture (diameter=200nm) in a thin layer of silver (thickness=200nm). The incident light field considered is linearly polarised and the transmission is determined across the output surface of the aperture.

The electric amplitude |**E**| of most efficient transmission OEi is shown in Fig. 9a illustrating scattering of the aperture. The OEi on its own in Fig. 9b. The transmission enhancement factor, with respect to the tightest Bessel beam achievable for a numerical aperture of NA=0.8 (Fig. 9c), is 2.1 and 1.55 with respect to the Airy diffraction limited disk with the same numerical aperture.

In the second example, we use the optical momentum operator associated with the quadratic measure *m*^{(F·uz)} as defined in Table 1. The momentum OEi with the largest positive eigenvalue corresponds to the field profile (Fig. 9 e–f) giving the largest optical force on the microparticle. Figure 9d shows the field amplitude |**E**| of this OEi on its own and scattering of the microparticle (Fig. 9e). The optical force enhancement factor, with respect to the plane wave illumination (Fig. 9f), is of 49.3 and of 1.33 with respect to the Airy diffraction limited disk with the same numerical aperture.

## 5. Discussion and Conclusion

We have experimentally and theoretically demonstrated an approach based on optical eigen-modes that enables the minimization of the free space spot size of a beam. Using full vectorial simulations in 3D, we have shown how the OEi approach can be used to optimise light-matter scattering interactions in the case of transmission through sub-wavelength apertures and optical forces on micro-particles. The generic nature of our approach means that it can be applied to other cases where the measure has a quadratic form and the propagation is linear. In the present paper we have verified the rigor of the method by demonstrating the experimental spot size operator and intensity operator optimization using Laguerre-Gaussian and Bessel light modes using a dual SLM to implement the technique. Future work will aim to extend this method to optimise the size and contrast of optical dark vortices, the linear Raman scattering or the fluorescence of different samples, the optical trapping force, and the angular momentum transfer in optical manipulation.

## Annex: Dimensionality study

In this appendix, we verify the convergence properties of our optical eigenmode method as a function of the dimensionality of the probing base used. Here, we consider a plane wave basis corresponding to the angular spectral decomposition of any incident field. The propagation is modelled using direct Fourier transform optics between the SLM plane and the target plane. Additionally the phase front of the incident beam is randomly changed to simulate the effect of high aberrations within the optical train. The dimensionality for the optical eigenmode method corresponds to the number of plane waves, *N*, taken into account, whilst retaining a fixed numerical aperture. The convergence behaviour is compared to standard phase front correction methods which can also be used to achieve focalised spots in the case of highly aberrated light fields. More precisely, these methods are based on the variable partitioning of the SLM to create *N* beamlets whose phases are individually changed such that all beamlets constructively interfere in the focal target point [23,24]. This approach delivers a final correction phase mask that is able to correct for aberrations in the SLM incident light field additionally to the pre-correction of the propagation aberration between the SLM and the target focal plane [24]. This correction mask delivers an enhanced focal intensity *η* shown in black in Fig. 10a for both standard methods [23, 24]. For comparison purposes, we define a phase front correction mask using the phase part of the intensity operator eigenmode (red in Fig. 10a) which shows that, for all three methods, the enhancement factor scales linearly with the dimension *N*.

Further, Fig. 10b shows that the standard approaches cannot beat the diffraction limit, given by the Airy disk. In stark contrast, the smallest spot size eigenmode delivers sub-diffraction over almost the whole dimensional range considered with increasing efficiency (blue curve in Fig. 10a) for increasing dimensionality.

## Acknowledgments

We thank the UK EPSRC for funding through the Nanoscope Basic Technology grant. Mark Dennis is acknowledged for the introduction to super-oscillations. KD is a Royal Society-Wolfson Merit Award Holder.

## References and links

**1. **C. Cohen-Tannoudji, *Quantum Mechanics* (Wiley, New York, 1977).

**2. **M. Kac, “Can one hear the shape of a drum?” Am. Math. Mon. **73**, 1–23 (1966). [CrossRef]

**3. **A. Sudbo, “Film mode matching: a versatile numerical method for vector mode field calculations in dielectric waveguides,” Pure Appl. Opt. **2**, 211 (1993). [CrossRef]

**4. **P. Bienstman and R. Baets, “Optical modelling of photonic crystals and VCSELs using eigenmode expansion and perfectly matched layers,” Opt. Quantum Electron. **33**, 327–341 (2001). [CrossRef]

**5. **J. Reithmaier, M. Röhner, H. Zull, F. Schäfer, A. Forchel, P. Knipp, and T. Reinecke, “Size dependence of confined optical modes in photonic quantum dots,” Phys. Rev. Lett. **78**, 378–381 (1997). [CrossRef]

**6. **H. Kogelnik and T. Li, “Laser beams and resonators,” Appl. Opt. **5**, 1550–1566 (1966). [CrossRef] [PubMed]

**7. **J. Barton, D. Alexander, and S. Schaub, “Theoretical determination of net radiation force and torque for a spherical particle illuminated by a focused laser beam,” J. Appl. Phys. **66**, 4594–4603 (1989). [CrossRef]

**8. **M. Mazilu, “Spin and angular momentum operators and their conservation,” J. Opt. A **11**, 094005 (2009). [CrossRef]

**9. **F. García-Vidal, E. Moreno, J. Porto, and L. Martín-Moreno, “Transmission of light through a single rectangular hole,” Phys. Rev. Lett. **95**, 103901 (2005). [CrossRef] [PubMed]

**10. **A. Assion, T. Baumert, M. Bergt, T. Brixner, B. Kiefer, V. Seyfried, M. Strehle, and G. Gerber, “Control of chemical reactions by feedback-optimized phase-shaped femtosecond laser pulses,” Science **282**, 919 (1998). [CrossRef] [PubMed]

**11. **M. R. Dennis, R. P. King, B. Jack, K. O. Holleran, and M. J. Padgett, “Isolated optical vortex knots,” Nat. Phys. **6**, 118 (2010). [CrossRef]

**12. **L. C. Thomson, G. Whyte, M. Mazilu, and J. Courtial, “Simulated holographic three-dimensional intensity shaping of evanescent-wave fields,” J. Opt. Soc. Am. B **25**, 849–853 (2008). [CrossRef]

**13. **R. Paschotta, *Encyclopedia of Laser Physics and Technology* (Wiley-VCH, 2008).

**14. **T. Sales and G. Morris, “Fundamental limits of optical superresolution,” Opt. Lett. **22**, 582–584 (1997). [CrossRef] [PubMed]

**15. **M. Berry, “Faster than Fourier,” in “*Quantum Coherence and Reality; in celebration of the 60th Birthday of Yakir Aharonov*,”, J. S. Anandan and J. L. Safko, eds. (Singapore, 1994), World Scientific, pp. 55–65. [PubMed]

**16. **M. R. Dennis, A. C. Hamilton, and J. Courtial, “Superoscillation in speckle patterns,” Opt. Lett. **33**, 2976–2978 (12). [PubMed]

**17. **J. Durnin, “Exact solutions for nondiffracting beams. I. the scalar theory,” J. Opt. Soc. Am. A **4**, 651–654 (1987). [CrossRef]

**18. **K. Volke-Sepulveda, V. Garcés-Chávez, S. Chavez-Cerda, J. Arlt, and K. Dholakia, “Orbital angular momentum of a high-order Bessel light beam,” J. Opt. B-Quantum S. O. **4**, S82–S89 (2002). [CrossRef]

**19. **B. Richards and E. Wolf, “Electromagnetic Diffraction in Optical Systems. II. Structure of the Image Field in an Aplanatic Systems,” Proc. R. Soc. London, Ser. A. **253**, 357–379 (1959).

**20. **R. Di Leonardo, F. Ianni, and G. Ruocco, “Computer generation of optimal holograms for optical trap arrays,” Opt. Express **15**, 1913 (2007).

**21. **D. Malacara, *Optical Shop Testing* (Wiley-Interscience, 1992), 2nd ed.

**22. **H. Wang, L. Shi, B. Lukyanchuk, C. Sheppard, and C. T. Chong, “Creation of a needle of longitudinally polarized light in vacuum using binary optics,” Nat. Photonics **2**, 501–505 (2008). [CrossRef]

**23. **I. Vellekoop and A. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. **281**, 3071–3080 (2008). [CrossRef]

**24. **T. Čižmár, M. Mazilu, and K. Dholakia, “In situ wavefront correction and its application to micromanipulation,” Nat. Photonics **4**, 388–394 (2010). [CrossRef]