J ul 2 01 0 Optical Quadratic Measure Eigenmodes

We report a mathematically rigorous technique which facilitates the optimization of various optical properties of electromagnetic fields. The technique exploits the linearity of electromagnetic fields along with the quadratic nature of their interaction with matter. In this manner we may decompose the respective fields into optical quadratic measure eigenmodes (QME). Key applications include the optimization of the size of a focused spot, the transmission through photonic devices, and the structured illumination of photonic and plasmonic structures. We verify the validity of the QME approach through a particular experimental realization where the size of a focused optical field is minimized using a superposition of Bessel beams.


Introduction
The decomposition of fields into eigenmodes is a well established technique to solve various problems within physical sciences. The most prominent example refers to Schrödinger's equation, within the field of quantum mechanics, where energy spectra of atoms are determined via the eigenvalue spectra and associate 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], and optical cavities [5]. 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 operators such as orbital and spin angular momentum [6].
In this paper, we report a novel method which we term "Quadratic Measure Eigenmodes (QME)" which represents a generalization of the powerful concept of eigenmode decomposition within the field of optics. Crucially, it is shown 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 QME 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.
From a theoretical perspective, the QME optimization method is mathematically rigorous and may be distinguished from the multiple techniques currently employed ranging from genetic algorithms [7] and random search methods [8] to direct search methods [9]. 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 in the problem. In contrast, our proposed QME 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 background of the QME theory and show its properties in a general context of optimizing the quadratic measures of interfering waves. In the second part, we apply the QME formalism to maximize the transmission through apertures and to minimize the focal spot size. For these applications, we describe, respectively, the electromagnetic field as a superposition of scalar Laguerre-Gaussian beams and vectorial Bessel beams. In the final part of the paper, we report a particular experimental implementation of the QME method using computer controlled spatial light modulators to squeeze the spot size of a superposition of Bessel beams. The paper concludes with a discussion of the particular results obtained and with general comments on the versatility of the QME method to a wide range of problems.

Fundamental concepts
The QME method is based upon two fundamental properties of the electromagnetic field and its interactions. Firstly, the approach relies on the linearity of the electromagnetic fields, i.e., the sum of two solutions of Maxwell's equations is itself a solution of them. As we consider free space propagation, this criteria is satisfied. The second property relates to the interaction of the electromagnetic field with its environment. All such interactions can be written in the form of quadratic expressions with respect to the electric and magnetic fields. Examples include the energy density, the energy flow, and Maxwell's stress tensor. This allows us to designate appropriate QME to various parameters (e.g. spot size) and subsequently ascertain the optimal eigenvalue which, in the case of a spot size operator, yields a sub-diffraction optical spot. In this section, we present the details of the theory underpinning our approach.

Theoretical background
Electromagnetic waves To demonstrate our method, we consider monochromatic solutions of the free space Maxwell's equations: where E and H are the spatial part of the electric and magnetic vector fields and where ε 0 and µ 0 denote the vacuum permittivity and permeability. The time dependent carrier wave is given by exp(iωt). These monochromatic solutions of Maxwell's equations can be written in an integral form linking the electromagnetic fields on the surface A with the fields at any position r, where √ 2F = ( √ ε 0 E, √ µ 0 H) is a shorthand for the two electromagnetic fields having six scalar components F i . The integration kernel P i j corresponds to a propagation operator giving rise to different vector diffraction integrals such as Huygens, Kirchhoff [10], and Stratton-Chu [11].
Quadratic measures Crucially all "linear" and measurable properties of the electromagnetic field can be expressed as quadratic forms of the local vector fields and are therefore termed quadratic measures. For instance, the time averaged energy density of the field is proportional to F * · F = 1/2(ε 0 E * · E + µ 0 H * · H) while the energy flux to 1/2(E * × H + E × H * ). The asterisk * stands for the complex conjugate. Integrating the first quantity over a volume determines the total electromagnetic energy in this volume, and integrating the normal energy flux across a surface yields the intensity of the light field incident on this surface. All the quadratic measures m κ can be represented in a compact way by considering the integral where the kernel κ i j = κ † ji is Hermitian with † the adjoint operator including boundary effects for finite volumes. Table 1 enumerates some operators associated to common quadratic measures. The integrand part of most of these quadratic measures corresponds to the conserving densities, which together with the associated currents are Lorentz invariant [6]. The volume, over which the integral is taken, does not need to be the whole space and can be a region of space, a surface, a curve, or simply multiple points. To account for this general integration volume, we broadly term it the region of interest (ROI) in the following. Quadratic measures eigenmodes Finally, using the general definition (3) of the quadratic measure it is possible to define a Hilbert sub-space, over the solutions of Maxwell's equations, with the energy operator defining the inner product. Furthermore, any general quadratic measure defined by (3) can be represented in this Hilbert space by means of its spectrum of eigenvalues and eigenfunctions defined by: Depending upon the kernel κ i j or operator κ, the eigenvalues λ form a continuous or discrete real valued spectrum which can be ordered. This gives direct access to the solution of Maxwell's equations with the largest or smallest measure. The eigenfunctions are orthogonal to each other ensuring simultaneous linearity in both field and measure.
In the following, we study the case of different quadratic measure operators and their spectral decomposition into the QME. The convention for operator labeling we adopt is to use the shorthand QME followed by a colon and a shorthand of the operator name. In the following, we discuss some examples of different quadratic measure operators. QME: Intensity operator (QME:IO) The quadratic measure corresponding to the QME:IO measures the electromagnetic energy flow across a surface ROI: where κ is chosen such that it corresponds to the projection of the Poynting vector on the normal n to the surface. The eigenvector decomposition of this operator can be used to maximize the optical throughput through a pinhole or to minimize the intensity in dark spots. Considering a closed surface ROI surrounding an absorbing particle, the QME approach gives access to the field that either maximizes or minimizes the absorption of this particle. QME: Spot size operator (QME:SSO) One way to define the spot size of a laser beam is by measuring the second order intensity moment (SOIM) w [12]. w can be expressed as with m (0) as defined in Eq. (4) and the QME:SSO defined by where r is the position vector and r 0 the centre of the beam. Accordingly, the square root of the QME:SSO eigenvalues multiplied by two gives direct access to the respective beam size provided that the intensity is normalized to one within the ROI, i.e., m (0) = 1; QME: Optical force operator (QME:OFO) The optical force acting on a scattering object can be calculated by considering the momentum flux, given by Maxwell's stress tensor, across a surface ROI, surrounding the scattering particle. There are three quadratic measures associated with the optical force acting on the particle, one for each orthogonal direction: where corresponds to the measure of the force. The eigenvector decomposition of this operator can be used to maximize the optical scattering and optical trapping force on microparticles. This approach can be extended through the use of the angular momentum operator, defined locally by ir × ∇.

Applications
In this section, we put the QME concept into practice and provide a couple of examples which demonstrate the striking applicability of the QME formalism to answer different questions within the field of optics. One such question is determining the largest intensity that may be transmitted through a given optical structure such as metallic apertures [13]. Another question is what is the smallest optical spot one may achieve. Before we discuss these questions, we show how the operator formalism presented in section 2 is rendered into a practical formalism which we have also applied in the experiments described in section 4.

Practical implementation of the QME concept
For practical purposes, in particular in terms of the experimental realization of the QME concept, the optimization procedure is spatially separated; that is we consider 1) an initial plane at the propagation distance z = z 1 where we can superpose a set of i = 1 . . . N fields E i (x, y, z) and H i (x, y, z) (N > 1) shaped both in amplitude and phase and 2) a target plane at the propagation distance z = z 2 where we have the ROI within which the optimization is actually carried out. Due to linearity a superposition of fields in the initial plane is rendered into a superposition in the target plane both featuring the same set of superposition coefficients: .
Based on this, the QME:IO defined in (4) can be rewritten as The matrix M (0) is a N × N matrix with the elements given by the overlap integrals This matrix is equivalent to the QME:IO on the Hilbert subspace defined by the fields max will deliver the initial plane (z = z 1 ) or target plane (z = z 2 ) superposition which maximizes the intensity within the ROI. We remark, that the case of a superposition of scalar fields u i (x, y, z), the matrix operator (10) takes on the simpler form an approach used in section 3.2 to maximize the transmission through a pinhole. Similar to the QME:IO, the QME:SSO defined in Eq. (6) can be rewritten as where M (2) and b must be represented in the intensity normalised base in order to fulfill the requirement of unity intensity within the ROI (m (0) ! = 1, see paragraph "QME: Spot size operator (QME:SSO)" in section 2.1). M (2) is a N × N matrix with the elements given by After transforming back to the original base we denote the eigenvalues of M (2) as λ (2) k and the eigenvectors as v (2) k . Then, the eigenvector v (2) min associated with the smallest eigenvalue λ corresponds to the smallest spot achievable within the ROI. The respective superposition is for z = z 1 and z = z 2 . We employed this vectorial definition to minimize the size of a focused spot in section 3.4 using a superposition of vector Bessel beams. The spot size minimization performed both in section 3.3 using a superposition of LG beams and in the experimental section 4 is based on the simpler scalar version ot the QME:SSO matrix defined as In this subsection, we consider a superposition of LG beams propagating in the z-direction and modulating the carrier wave u c = exp(−i(k 0 z − ωt)). In cylindrical coordinates, LG beams are defined by [14]:

Maximizing transmission through apertures using Laguerre Gaussian beams
with z r = k 0 w 2 0 /2, q(z) = z + iz r , and w 2 (z) = 1 + z 2 /z 2 r where w 0 , k 0 , ω are the Gaussian beam waist, vacuum wave vector, and optical frequency, respectively. This beam profile is a solution of the paraxial equation for integer values of P and L corresponding respectively to the radial and azimuthal index of the beam. Within the paraxial approximation, the intensity of the beam transmitted through a planar ROI defined by a disk centered on r = 0 is proportional to R 0 2πr|u L P (r, φ , z)| 2 dr where R is the radius of the disk. The coefficient C L P = P!/(P + |L|)! is the normalization factor such that the total intensity of the beam, for an infinite ROI, is unity.
The transmission is maximized using the practical implementation of the QME concept described above in section 3.1; that is the QME:IO was assembled according to Eq. (12) using the respective representation in cylindrical coordinates. We only considered the radial family of LG beams (L = 0) and performed the optimization in the plane at z = 0. Figure 1 shows the final superposition u 0 max,P u 0 P (x, y, 0) in the case of a ROI with a radius equal to the waist of the Gaussian envelope.
In Fig. 2, we observe that the maximal transmission achievable via the QME intensity optimized beams depends on the number of LG beams considered in the superposition and on the size of the ROI. Indeed, the QME for a smaller ROI needs a larger number of LG modes to achieve 100% transmission.

Smallest focal spot using Laguerre Gaussian beams
Using a superposition of LG beams we have also minimized the size of a focal spot using the representation of the QME:SSO (17) in cylindrical coordinates. It is important to note at this point that we only retain the intensity eigenmodes 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 be performing so well that a spot of size zero is finally obtained if no intensity threshold is applied. Figure 3 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 QME:IO. Using the ratio between these two intensities we can define the intensity Strehl ratio [15] for the QME:SSO (see Fig. 4b).
We remark that both, the spot size and the Strehl ratio, show resonances as a function of the ROI size. This can be explained 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. 4), 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 for optical manipulation and imaging.
Analogously to the approach to calculate the smallest spot size in a planar ROI, we can determine the QME:SSO for a volume. In this case, the sidebands appear in the outside of the ROI in both the lateral and longitudinal directions. The different volumes considered in Fig. 5 suggest that there is an intrinsic link between squeezing light in the lateral direction and in the longitudinal direction.
On a final note, we remark that squeezing light below its diffraction limit may be associated with the effect of super-oscillations [16]. 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 QME spot size optimized beams, we have calculated the spectral density of the radial wave-vector for the smallest planar spot [17]. As shown in Fig. 6, this spectral density clearly identifies a spectral bandwidth (white background in Fig. 6). 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. We observe, that super-oscillations occur only 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.

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 base-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 [18] 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 [19] E = E 0 exp (iLφ + ik t z)) (αe x + β e y )J L (k t r) where 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 y and e z are the unit vectors in the Cartesian coordinate system. The parameter L 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 and the QME:IO and QME:SSO operators are assembled according to the expressions (10) and (15), respectively. Figure 7 shows a comparison between Airy disk, Bessel beam, and QME 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. 8). According to the theoretical foundations of the QME concept, the successful experimental implementation within the field of optics requires a linear optical system along with the ability to shape laser fields in both amplitude and phase. We have achieved this by using the experimental setup shown in Fig. 9. 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.

Experimental implementation of the QME concept
To conform this experimental section to the conventions introduced in section 3.1 we remark that we shaped a set of test fields both in amplitude and phase in the initial plane at z = z 1 which coincided with the two SLM panels and subsequently minimized the size of a focal spot in the target plane at z = z 2 which coincided with the CCD camera chip. For our proofof-principle experiments we ignored polarization effects and considered a set of scalar fields 1 . . . N) where the CCD camera detected the associated intensity I i (x, y, z 2 ) ∝ |E i (x, y, z 2 )| 2 . Both the QME:IO and the QME:SSO were assembled from these fields according to the scalar expressions (12) and (17), respectively.
The amplitudes A i (x, y, z 2 ) were determined by simply recording an associated intensity image I i (x, y, z 2 ) and subsequently taking the square root. We used the three-step phase retrieval algorithm described in Ref.
[21] to retrieve the phase modulations φ i (x, y, z 2 ). This algorithm is based on interference with a reference field E R (x, y, z) = A R (x, y, z)e iφ R (x,y,z) . The reference field's intensity was distributed uniformly over the ROIs considered in our experiments by adding a square phase φ R ∝ (x 2 + y 2 ) to the reference field using the SLM phase panel. Moreover, a constant phase term φ k = 2πk/3 was added for k = 0, 1, 2, and the superimposed fields E i (x, y, z 1 ) + E R (x, y, z 1 )e iφ k were then encoded onto the SLM. The associated intensity distributions explicitly read where the spatial coordinates (x, y, z) were omitted for brevity. These three intensity distributions represent a 3-dimensional equation system with three unknowns I bg,i , γ i and ∆φ i . The latter is explicitly obtained as where atan2{ζ , ξ } is the two argument arctangent function corresponding to the argument of the complex number ξ + iζ . Note that the reference phase φ R cancels out when calculating the operator elements due to multiplication of a complex conjugate field ∝ e −iφ R (x,y,z) with a complex field ∝ e iφ R (x,y,z) . Therefore Eq. (21) yields the adequate phase modulation required to assemble the QME operators. 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 QME 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 (x, y, z 2 ) corresponding to this encoded optimized field. The Num-S utilizes the fields E i (x, y, z 2 ), which were individually measured to assemble the QME operators, in order to numerically determine the intensity distribution as I Num-S (x, y, Crucially, linearity is verified if I Exp-S (x, y, z 2 ) = I Num-S (x, y, 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. In our experiments, we used N = 11 non overlapping amplitude ring masks with a constant phase modulation as fields of interest E i (x, y, z 1 ). After propagation through the Fourier lens 5 (see Fig. 9) the resulting fields E i (x, y, z 2 ) form a set of Bessel beams. Figure 10(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 SOIM 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 10(c) depicts a circular aperture which is encoded onto the SLM in order to observe the Airy disk (see Fig. (d)). The SOIM of the Airy disk is approximately 1.5 times larger than the core of the reference Bessel beam as expected [22].

Results and discussion
The results of the performed QME spot size minimization are shown in Fig. 11 for different sizes of the ROI. To begin with, the comparison of the Num-S intensity distribution I Num-S (x, y, z 2 ) (top row) and the Exp-S intensity distributions I Exp-S (x, y, 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 SOIM 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 3.4 and Fig. 7. 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 QME concept applied to spot size minimization. Moreover, the results strongly suggest that the QME 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.
We have performed an additional comparison of the experimental results to the theoretically predicted ones: The relative SOIM w/w B was evaluated for the Num-S for ROI radii ranging from R = 1 pixel to R = 50 pixel in steps of one pixel (see Fig. 12). In comparison to the corresponding graph shown in Fig. 8 we observe agreement in terms of both the general decrease of w/w B with decreasing R and the peaks occurring periodically along the R-axis. Surprisingly, for R ≥ 10, we do not observe values of w/w B ≈ 1; that is the QME spot minimization does not deliver, as one might naively expect, the reference Bessel beam featuring the smallest SOIM w B of all Bessel beams considered. This is due to the ring structure of Bessel beams which significantly adds to the SOIM and therefore to the spot size in the case of large ROIs. As a consequence, the QME optimization aims to superimpose the set of Bessel beams in a manner that the rings are destructively interfered within the ROI, only retaining the central core.
Given that the ring structure is essential for the reduced core size of Bessel beams compared to the diffraction limited Airy disk we observe values of w which are larger than w B . Indeed, the values of w are fairly close to the size of the Airy disk which was determined as w ≈ 1.5w B . This strongly suggests that for large ROIs the QME spot minimization yields a superposition of Bessel beams which closely resembles the diffraction limited Airy disk. Finally, we remark that for R < 7 the relative SOIM w/w B becomes very small and has to be taken with care due to the possible truncation of the spot mediated by the CCD detector's sensitivity threshold as mentioned above.

Discussion and Conclusion
We have theoretically and experimentally demonstrated a novel approach based on quadratic measures eigenmodes that enables the optimization of different optical measures. The theory that we employ is rigorous and based on considering the light-matter interaction as a quadratic measure originating from the fields described by Maxwell's equations. Excitingly, we can define many quadratic measure operators to which our approach is applicable (see Table 1). The method is thus very powerful and the generic nature of our approach means that it may be applied for example to optimize the size and contrast of optical dark vortices, the Raman scattering or fluorescence of any samples, the optical dipole force, and the angular/linear momentum transfer in optical manipulation. In the present paper we have verified the rigor of the method by demonstrating experimental spot size operator and intensity operator optimization using Laguerre-Gaussian and Bessel light modes using a dual SLM approach to implement the technique. We envisage the QME approach as providing a powerful and versatile theoretical and practical toolbox. Our generic approach is applicable to all linear physical phenomena where generalized fields interfere to give rise to quadratic measures.