Optical Eigenmodes; Exploiting the quadratic nature of the energy flux and of scattering interactions

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.


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]. Alltogether, 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.

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 where we considered a superposition of fields {E, H} = ∑ N j=1 a j E j , ∑ N j=1 a j H j and where (A) labels the light-matter interaction defined in Table 1 Given the Hermitian form of (2), we remark that the light-matter interaction M (A) defines a spectrum of real eigenvalues λ k . Each of these eigenvectors corresponds to a superposition of fields {E j , H j } termed here Optical Eigenmode (OEi). Crucially, we may now extremize the light-matter interaction considered; that is, the extremal eigenvalue λ ext, j H j which extremizes the interaction 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 (1) as follows: where M (0) and M (2) are termed the intensity operator (IO) and spot size operator (SSO), respectively. According to (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 where v (2) min,k is the eigenvector associated with the smallest eigenvalue of M (2) and in the intensity normalised eigenbase v Table 1. Time averaged quadratic measures m of common light-matter interactions. The integration either over a volume V or a surface S which in general corresponds to the Range of interest = ROI of the measure. In the optical momentum case, it corresponds to a closed surface surrounding the scattering object with F · u representing the optical force in the direction defined by the unit vector u. For surface integrals, n is the normal unit vector to the surface considered. The definition of the electromagnetic energy density is 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 M and M respectively. These scalar expressions are equivalent to the respective vector versions listed in Table 1 and determined through equation (2). The scalar version of the optimimum OEi (5) explicitly reads

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.

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] 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 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).

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 (z 1 ) both in amplitude A j and phase φ j in the SLM plane, and the associated intensities 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 threestep 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 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.

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.

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 equation (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. 9 a-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·u z ) 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.

Discussion and Conclusion
We have experimentally and theoretically demonstrated an approach based on optical eigenmodes 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.