Optimising superoscillatory spots for far-field super-resolution imaging.

Optical superoscillatory imaging, allowing unlabelled far-field super-resolution, has in recent years become reality. Instruments have been built and their super-resolution imaging capabilities demonstrated. The question is no longer whether this can be done, but how well: what resolution is practically achievable? Numerous works have optimised various particular features of superoscillatory spots, but in order to probe the limits of superoscillatory imaging we need to simultaneously optimise all the important spot features: those that define the resolution of the system. We simultaneously optimise spot size and its intensity relative to the sidebands for various fields of view, giving a set of best compromises for use in different imaging scenarios. Our technique uses the circular prolate spheroidal wave functions as a basis set on the field of view, and the optimal combination of these, representing the optimal spot, is found using a multi-objective genetic algorithm. We then introduce a less computationally demanding approach suitable for real-time use in the laboratory which, crucially, allows independent control of spot size and field of view. Imaging simulations demonstrate the resolution achievable with these spots. We show a three-order-of-magnitude improvement in the efficiency of focusing to achieve the same resolution as previously reported results, or a 26 % increase in resolution for the same efficiency of focusing.


Introduction
In recent years, it has been comprehensively shown that the resolution of optical imaging can be pushed beyond Ernst Abbe's diffraction limit. In biological super-resolution imaging, the widely used techniques all require fluorescent labelling of the sample (for a review see [1]). Numerous techniques for achieving super-resolution using evanescent waves have also been suggested and demonstrated [2][3][4][5][6][7], but these all require a probe within tens of nanometres of the sample, meaning they are not generally applicable to thick samples such as biological cells. The Pendry-Vesselago superlens [8,9] theoretically allows unlabelled imaging from the far-field, but still faces many challenges in practical realisation. Superoscillatory imaging [10,11] promises non-invasive unlabelled super-resolution imaging from the far-field. Superoscillatory imaging relies on producing an optical diffraction pattern with localised spots that are spatially much smaller than the conventional diffraction limit allows. The fact that this is possible was noted by G. Toraldo di Francia in 1952 [12], in applying concepts from super-directive antennas in the microwave community to optical instruments. Closely linked to this phenomenon is the mathematical concept that a function can locally oscillate much faster than its fastest Fourier component [13,14]. The relevance of this mathematical concept to several physical systems was noted by Berry [15], inspired by the work of Aharonov et al. on weak measurements in quantum systems [16]. Berry and Popescu [17] then proposed a method for creating optical superoscillations using a complex diffraction grating; at the same time Huang et al. [18] provided the first demonstration of optical superoscillations by diffraction from a quasiperiodic array of nanoholes. The cost of these small superoscillatory spots is that energy is directed away from the superoscillatory features [19]. In optics, this manifests itself as high energy sidebands to any superoscillatory spot.
Developing this further, Rogers et al. [10] built an optical superoscillatory microscope and demonstrated its super-resolving capabilities. A superoscillatory lens of concentric rings produces a diffraction pattern with the required superoscillatory spot. The width and number of rings is optimised for spot size using binary particle swarm optimisation; appropriate constraints ensure an experimentally usable field of view (FoV) and sideband intensity. We use the term FoV throughout to refer to the region in which the superoscillatory function is controlled, as this has a clear meaning in the optical applications we propose. The term is equivalent to the region sometimes called the 'superoscillation region' or the 'region of silence' in the superoscillation literature. There is no theoretical limit to the size of spot or FoV that can be produced [20], but it is clear that there are limits in resolution imposed by the proportion of energy in the spot relative to the sidebands and the size of the FoV. Therefore, techniques for optimising superoscillatory spots are crucial to developing practical instruments and probing the limits of superoscillatory optical imaging.
Superoscillatory spots have been produced in a number of ways, for example by a quasi-periodic array of nanoholes in a metal screen [18], binary ring masks [21], using micro/nanofibre arrays [22] and in random speckle patterns [23]. Theoretical techniques for designing superoscillatory functions with particular properties have also been demonstrated, for example by manipulating zeros [24], by prescribing values of the superoscillatory function and minimising total energy [25,26] or energy yield [27], and by using variational techniques to find minimum energy solutions that satisfy specified constraints on the amplitude and/or derivative of the desired superoscillatory function [19,28].
Methods combining theoretical design of superoscillatory functions with techniques for producing superoscillatory spots in optical instruments have also been demonstrated [29][30][31]. Zacharias et al. [32] achieve superoscillatory focusing in their beams, created by the superposition of optical Bessel beams and generated with a spatial light modulator (SLM).
By adapting techniques from superdirective antenna theory, Wong and Eleftheriades [33] produce superoscillatory functions which are realised using a SLM in their super-microscope. This is an extension of the 1D technique in [34], which uses an expansion in Tschebyscheff polynomials and determination of zeros to find the narrowest peak for a given grass level and FoV. The authors also extend this technique to 3D in [35].
Canales et al. [36] optimise the design of an annular binary pupil filter for specified restrictions on spot size, Strehl ratio and/or intensity ratio of the central spot to the first sideband. This is done by expanding the point spread function using the first-order Bessel functions as a basis, and then solving equations for the restrictions to obtain the mask parameters. This extends the methodology in [37,38] by removing the second order approximation for the point spread function and by combining amplitude and phase in the mask design.
Optical eigenmodes are used to directly determine the smallest spot of maximal intensity within a given region of interest in [39,40]. This method is used to find a range of spots tailored to different measurement tasks by considering a number of FoV in [41]. For each FoV, the relative spot intensity (to that of the nearest side-band) is found. These spots are created using a SLM and focusing lens, and imaging is demonstrated with a confocal system.
Complete control of the shape of a superoscillatory feature is achieved in [42], extending work done in [43] on diffraction-free beams. The authors use a superposition of high order Bessel beams that share the same transverse wave number, so that their superposition is propagation-invariant. A Taylor expansion of the desired superoscillatory shape determines the coefficients of the Bessel beams to achieve the desired superoscillatory feature.
These methods introduce varying levels of control over different aspects of the constructed superoscillations. However, none of these provide simultaneous optimisation of spot size, relative spot intensity and the FoV, which are crucial to realising the potential of superoscillatory imaging. We present a method which introduces a significant extra degree of control over the superoscillations, by simultaneously optimising the size of a superoscillatory spot and its relative intensity to the sidebands for any given FoV. There is no theoretical limit on the size of a superoscillatory spot or its isolation from the sidebands [20], so it is the intensity of the sidebands relative to the spot and the distance from the spot to the sidebands that limit the practically achievable imaging resolution. We find the optimal trade-off between spot size and relative intensity of the spot for a number of FoV. Including the relative intensity of the spot to the sidebands in our optimisation means we can achieve a three-order-of-magnitude improvement in focusing efficiency for the same imaging resolution, or a 26 % improvement in resolution for the same focusing efficiency.
Our numerical method uses the circular prolate spheroidal wave functions (CPSWFs) as a basis set on the FoV [44]. Huang and Zheludev [20] note the equivalent property of the prolate spheroidal wave functions (PSWFs) [45] in 1D, but the practical applications of 1D focusing are limited. 2D focusing, on the other hand, is of great practical importance in many areas of optics. The CPSWFs allow full control of the intensity in the FoV for radially symmetric 2D systems: any possible intensity profile in the FoV can be constructed as a linear combination of CPSWFs. As in the 1D case, there is no theoretical limit on the spot size or the FoV, as shown in 2D in Section 2. A genetic algorithm and a semi-analytic optimisation method (verified by comparison to the genetic algorithm) are used in Sections 2.1 and 2.2 to obtain a set of optimal compromises between the spot size and relative intensity for a number of different FoV. This introduces new capabilities in practical spot design. We show in Section 3, through imaging simulations, that these spots offer significantly improved resolution and focusing efficiency compared to those in the literature.

Constructing superoscillatory spots
Before demonstrating how the CPSWFs can be used to construct any radially symmetric shape of superoscillatory spot in 2D, some definitions and specification of notation are required. Assuming radial symmetry, the 2D intensity profile of light in the focal plane of a conventional lens can be described by a function, I(r). Here r ∈ [0, ∞) is the radial distance from the centre of the focal plane. An ideal I(r) for super-resolution imaging has a central peak (when reflected in the axis of r = 0) with vanishingly small width, and is fully isolated from its sidebands. Assuming also monochromatic light of wavelength λ, and that r is normalised to that wavelength, then a spot in I(r) of width less than approximately 0.5 is smaller than the conventional diffraction limit allows.
Considering a central peak at r = 0, the full width at half maximum (FWHM) provides a measure of the width of the peak and is defined to be the width of the spot at half the value at its peak (at r = 0). This is shown in Fig. 1(a). Also shown is the FoV, defined to be a number of wavelengths, D, such that appropriate restrictions can be placed on I(r) for r ∈ [0, D] (usually specifying the region of low intensity around the peak, for example a region where the intensity does not exceed a small fraction of the central intensity). This region is sometimes referred to elsewhere as the 'superoscillation region' or 'region of silence'. Note that in our later numerical results the intensity may increase rapidly at the very edge of the FoV, effectively shrinking the FoV by a few percent. A reflection in the axis of r = 0 (as shown in Fig. 1) is always assumed when referring to the FWHM (representing a planar slice of the 2D intensity profile).
An ideal intensity profile, I(r), thus has a central peak with vanishingly small FWHM for any size of FoV. Before considering limitations inherent in imaging systems, we construct such an ideal I(r) mathematically by extending the technique in [20] to 2 dimensions. As such, it is necessary to also consider the electric field strength E(r) such that |E(r)| 2 = I(r). As in the 1D case, the function E(r) must be band-limited so that the superoscillatory spot can be created by the diffraction/focusing of propagating light. Now the (band-limited) CPSWFs form an orthogonal basis on the interval [0, D] and an orthonormal system on the interval [0, ∞) for all functions of finite energy (i.e. of L 2 ([0, D]) and L 2 ([0, ∞)), respectively) [44]. They also form an orthonormal basis on [0, ∞) of all such functions that are band-limited under the Hankel transform. Thus for r ∈ [0, D] any sufficiently smooth, real function with finite energy can be constructed from a series of (band-limited) CPSWFs, with the resulting series extending to r ∈ [0, ∞) and clearly also being band-limited. So a band-limited function approaching the delta function can be constructed for any choice of D; that is the FWHM can be made arbitrarily small for any choice of the FoV. (Note that [44] give the above result for the interval [0, 1] rather than [0, D]. However, re-scaling of the system to an axis where r is normalised to D is trivial.) To construct such a series of CPSWFs numerically, it is first necessary to construct the individual CPSWFs. This is done using the numerical procedure in [44], and the alternative, more usable form of the CPSWFs, φ k (r) (with integration weighting r). These are hereafter referred to as the basis functions, and the squares of the first four |φ k (r)| 2 , k = 1, 2, . . . , 4 for r ∈ [0, 1] are plotted in Fig. 1(c) (normalised to have equal energy) as well as their power spectra in Fig. 1 The functions are reflected in the axis of r = 0 to show a planar slice of the 2D function. When applied to light, the |φ k (r)| 2 represent the intensity profile of each basis function. Note that the numerical procedure of [44] returns the basis functions in order of decreasing eigenvalue which coincides with a decreasing proportion of energy in the FoV and a decreasing FWHM (seen most easily in the inset to Fig. 1(c)). The bandlimited nature of the CPSWFs is seen in Fig. 1 Since the CPSWFs form a basis set on L 2 ([0, D]), the shape of the spot inside r ∈ [0, D] can be controlled without consideration of the region where r > D, and since this is valid for any D then the FoV can be made arbitrarily large. An obvious target function to construct using this band-limited basis set is a sinc function on r ∈ [0, 1] with vanishingly small FWHM. For example, the function f t (r) = sinc(arπ) has vanishingly small FWHM at r = 0 as a → ∞. An approximation to this function f N (r) which approaches f t (r) as N → ∞ can be constructed within the FoV: (1) Figure 1(b) shows the target function, f t (r) = sinc(arπ) for a = 8, which has a FWHM of 0.1λ. Any value can be specified for the FoV; in this figure λ = 1 and the FoV is 1λ. Also shown are f N (r) for N = 3, 6, 9, 12 and 15 (where N is the number of basis functions used in the approximation). N = 15 gives a very good approximation to f t (r) and the corresponding intensity profile has a FWHM of 0.1λ. This illustrates how an arbitrarily small band-limited spot can be constructed with an arbitrarily large FoV in 2D. If a spot with E(r) = f t (r) could be produced and used successfully in imaging there would be no limit to the achievable resolution. However, imaging systems are limited by the intensity outside the FoV, commonly called the sidebands. A high sideband intensity is problematic for two reasons: first, forming a faithful image of a sample requires suppressing the contribution of the image from the sidebands [10]. This is achievable for moderate sideband intensities, but becomes impossible if they become too intense, with the precise intensity determined partly by the FoV. The second reason is that, for a spot with limited total energy (as in experiment), as the sidebands grow in intensity, the central spot intensity must necessarily reduce. For very large intensity sidebands the intensity of the central spot will therefore drop below the detectable limit. Considerations beyond the FWHM and FoV are therefore necessary to produce spots that can be  Intensity /arb. used in practice in imaging. A useful metric for this is the intensity ratio, IR, defined to be the ratio of the maximum intensity outside the FoV and the intensity of the spot at r = 0: This measure is deliberately chosen in preference to the commonly used energy measures such as total energy or energy yield [19,27,47], as the solution with optimal energy yield may not be optimal for optical imaging. The quality of the image obtained is strongly dependent on the distribution of intensity in the sidebands, with intensity nearer the FoV being more problematic. A lower-energy-yield solution, but where the sideband energy is distributed further from the FoV, may enable finer resolution than a higher-energy-yield solution.
Consider the spot with electric field profile E 15 (r) = f 15 (r) so that I 15 (r) = | f 15 (r)| 2 . This has FWHM 0.1λ and an IR of 4.5 × 10 41 . Considering instead the less ambitious electric field profile of E 6 (r) = f 6 (r) with an (intensity) FWHM of 0.27λ, the IR is approximately 7.3 × 10 5 . Changing the target function to an electric field profile with a = 3 in f t (r) gives a FWHM of 0.30λ and is well approximated by 6 basis functions, giving a FWHM of 0.29λ and IR approximately 2700. All these IRs are too high to be practical in imaging systems, so the procedure above, while demonstrating an interesting theoretical point (on which we will later build), does not produce practical spots for imaging.
This problem is not specific to this particular target function. In general, high-order basis functions are necessary to approximate the precise shape of a given target function well, and these contribute significantly higher intensity sidebands. The basis functions φ k have an exponentially decreasing proportion of energy inside the FoV with increasing k [44], thus if higher-order basis functions are used, the high intensity sidebands increase rapidly. Since the detailed shape of the target function is largely irrelevant (a narrow central peak is all that is required), a more appropriate optimisation problem, minimizing the width while ignoring the precise shape, is considered below.

Optimal spots for imaging
In the previous section, superoscillatory functions of theoretically ideal electric field profiles are constructed. To produce results applicable to imaging, we now optimise the function, f (r), taking the intensity of the sidebands into account. As above, we seek a series of basis-function coefficients that will give an optimal electric field profile. However, rather than using the standard technique of specifying a target function and then solving for the coefficients of the basis functions, we take a more general approach. Considering radially symmetric 2D systems with a specified FoV, D, we define a multi-objective optimisation problem as follows: The problem space is all possible electric field profiles in the imaging plane; here represented by the set of all linear combinations of the basis functions φ k (r). As shown above, there is no limit on the first objective, i.e. the FWHM can be made vanishingly small. However, it is well known (see, for example, [19]) that there is a trade-off between the first objective and the other two. The addition of the second objective provides a solution set that is usable for imaging, so that the minimum FWHM can be found for a range of acceptable IRs. As noted earlier, this measure (as opposed to energy measures), allows solutions which may be better suited to optical imaging yet have lower energy yield. Apparent optimal solutions with small FWHM but negligible energy are not a concern despite the point measure used, since the requirement that the functions are band-limited causes such solutions to have high IRs. The third objective is necessary as the intensity profile within the FoV is no longer controlled by approximation to a target function. The IR within the FoV, IR in , is defined to be the ratio of the maximum intensity beyond the central peak and within the FoV to the intensity of the central peak. Some regions of notable intensity within the FoV are acceptable in imaging, and adding this as a third objective enables the best compromises to be found. The result of the optimisation is a Pareto front, which is the set of points which form the best compromises between all these objectives. Thus the set of optimal points for a given particular set-up can be read from the (3D) Pareto front.
To allow a numerical exploration of the infinite problem space, the series of basis functions is truncated to include the first N = 5 basis functions. This truncation is justified because, as mentioned earlier, the proportion of energy inside the FoV decreases exponentially over the sequence of basis functions. We find that basis functions of order 6 and higher have negligible effect within the FoV for the range of FWHM considered here, when a reasonable IR is enforced. The initial population of the genetic algorithm is a set of randomly generated sequences of coefficients of length N, with each sequence representing the function obtained by adding the basis functions using that sequence of coefficients.
Three fitness functions corresponding to the three objectives compute the FWHM, IR and IR in for each member of the population. We use the multi-objective genetic algorithm gamult, available in the MATLAB Global Optimization toolbox, which runs until the relative change in the movement of the Pareto front is less than a set tolerance. For 5 basis functions, a population of 5000 and a tolerance of 1 × 10 −4 the code typically runs for 100-110 generations and takes approximately 3.5 minutes on a standard desktop PC. There is negligible difference in the result when using different starting populations.
The 2D Pareto fronts, i.e. the set of best FWHM-IR compromises for three different FoV, are shown in Fig. 2 (FoV = 1.6λ in Fig. 2(a), FoV = 1λ in Fig. 2(b) and FoV = 0.5λ in Fig. 2(c)). For easier visualisation, these points have been selected from the 3D Pareto front such that they have an IR in of no more than 0.25. Considering a FWHM of 0.30λ and a FOV of 1λ, this can now be achieved with an IR of approximately 10, compared to the IR of 2700 obtained with the earlier sinc approximation. A FWHM of 0.41λ can be created with an IR of 0.2 (note this point has an IR in of 0.002). This is a notable improvement on the conventional diffraction limit with an extremely low IR, making it very easy to use for imaging with particular benefits during investigations of new imaging systems.
It is easy to calculate the total energy of the optimal solutions by integrating the intensity over all space, allowing the energy cost of the superoscillations to be quantified. In doing so, we find that the total energy depends polynomially on FWHM and exponentially on the FoV, in line with what is known in the literature [19]. So, although the Pareto fronts are not guaranteed to show the global optimal solutions due to the nature of genetic algorithms, they agree with previously observed trade-offs in superoscillatory functions and give significant improvements in achievable imaging resolution, while also being easily produced (see section 3).

Two-function optimisation
An optimisation method that can be used in real time (runs in around 20 ms on a standard desktop PC) is more practical for adjusting for different FoV while a superoscillatory imaging system is in use, since the optimisation has to be run for each FoV. Such a method, termed the two-function method, is presented here, where the optimal ratio of the coefficients of only two carefully selected basis functions is found. By considering linear combinations of only two basis functions, we reduce the problem space significantly, and make searching it much faster. IR in is restricted to be less than 0.25, but in most cases is significantly smaller. In contrast to the genetic algorithm, which considers the full problem space, many possible solutions are disregarded in this method. The results are, however, comparable to those of the genetic algorithm, and the optimisation problem is reduced to a simple search of a decreasing function.
This method uses the observation that the basis functions themselves are good solutions, which is unsurprising since the PSWFs are known to be energy-yield optimal under certain conditions in 1D [27]. The careful selection of two basis functions relies on the ordering of the basis functions by decreasing FWHM (see Fig. 1(c)). We specify a target FWHM and FoV, and the algorithm selects the two basis functions φ l and φ g with FWHM just less than, and just greater than the target FWHM. Since only the proportion of φ l relative to φ g is needed, the coefficient of φ l is set to 1. The optimal coefficient of φ g is then a search for one parameter over which the corresponding FWHM is strictly decreasing, allowing a function with the target FWHM (within a set tolerance) to be found very quickly. Note that using exactly two basis functions enables a continuous range of FWHM (rather than the discrete set obtained if using single basis functions) while keeping the optimisation simple: adding further basis functions would add unnecessary complexity without significantly improving the result.
The resulting Pareto fronts obtained using repeated runs of this method for different FoV are shown together with the Pareto fronts from the genetic algorithm in Figs. 2(a)-2(c). This shows sufficiently close agreement for this sub-second optimisation method to be a practical alternative to the more general approach of using a genetic algorithm, despite the reduced problem space. It also supports the observation that the basis functions are themselves near optimal solutions. (The apparent minor discrepancy for a FoV of 0.5λ and FWHM of approximately 0.3λ has arisen due to the genetic algorithm finding only points with IR in of less than 0.084 in this area, while the point from the two-function optimisation has an IR in of 0.24). Also shown in Fig. 2(d) is the intensity profile of three spots created by this simpler optimisation method. These have FWHM of 0.40λ, 0.33λ and 0.26λ with corresponding IR of 0.23, 0.88 and 21 and with IR in of 0, 0.08 and 0.2, respectively, for a FoV of 0.5λ.
The power of the CPSWF approach is shown in Fig. 3, where a range of 2D spots are shown with variation in FoV along the horizontal direction and spot size in the vertical direction. From this palette, the appropriate spot for any given experiment can be easily selected. These spots are generated using the two-function method, and so generation of the entire figure takes less than a second. Given this degree of control and speed, you can easily decide during an experiment that, say, you need a 20 % larger FoV and realise that spot at the click of a button.

Spot analysis and imaging simulations
Having shown above how to create optimised superoscillatory spots, we now demonstrate their performance in confocal imaging. Our numerical methods (described below) simulate a scattering confocal microscope, using two-function-optimised spots to illuminate the sample. We then analyse the images taken by five sample spots, as well as examining the trade-offs, over the whole parameter space, between efficiency and resolution inherent in superoscillatory imaging systems. Our analysis follows the thorough methodology of Kosmeier et al. [40], allowing easy comparison with previously described optimisation approaches.

Simulation methodology
In a standard laser scanning confocal microscope, a beam is focused on a sample through an objective, and scanned over the sample using a scanning mirror. Equivalently, as in our simulations, the sample may be scanned through the beam using a translation stage. The interaction of the beam with the sample is then imaged onto a pinhole in front of a detector. The role of the pinhole in a conventional confocal microscope is to block the out-of-focus light; in our system it also blocks light scattered from the sidebands by the sample, allowing imaging of objects larger than the FoV. The image is formed from the detected intensities at each point of the sample.
We consider a thin transmissive sample which is modelled as a binary distribution S(x, y) which equals 0 at those points where the sample is opaque and 1 at the points where it is transparent. For each position of the sample, the sample is multiplied (point-by-point) by the field of the focused beam E(x, y) resulting in a complex field A(x, y) = S(x, y)E(x, y). This field, A, is convolved with the amplitude point spread function PSF A of the detection objective -equivalent to filtering out those high spatial frequencies which are not collected by the objective. The modulus squared of this complex convolution P(x, y) = | A ⊗ PSF A | 2 is the intensity distribution that is imaged onto the pinhole. The detected intensity is the sum of those parts of P which pass through the pinhole. This process is repeated for each position of the sample to build up a scanned confocal image. In the simulations below, we use an effective pinhole diameter of 20 nm. This small pinhole (0.07 Airy units) improves the resolution of both the conventional and the superoscillatory imaging by the same factor, but the relative resolution improvement, as discussed below, is not strongly dependent on pinhole size. In the simulations, we simply use the spots as designed: in an experiment they would most easily be created using a spatial light modulator in the pupil plane of the illuminating objective, to control the amplitude and phase of the illuminating beam (as in, for example [48]). The required amplitude and phase pattern on the SLM can easily be calculated from the (numerical) Fourier transform of the spot amplitude.
As the spot design is carried out in a normalised set of coordinates, it is necessary to scale them to a physical axis for realistic imaging simulations. The normalised axes are in terms of wavelength, and with a (spatial) bandwidth equivalent to an NA of 1. Thus, to scale to a physical system, we simply scale the spots such that one wavelength in normalised space corresponds to λ/NA in the physically scaled system.
All imaging simulations are carried out on objects consisting of two holes with varying separation (characterised by the centre-to-centre distance). To determine the resolution achieved by a particular spot, we vary the distance between the holes, imaging at each step, until we find the distance at which the points are just resolved. This is taken to be the distance at which the intensity in the centre of the image falls below 73.5 % of the intensity at the holes, equivalent to the intensity dip seen in two spots resolved by the Rayleigh criterion. In some cases, although a resolution number could be calculated, the intensity pattern was not the expected double-peaked structure. In these cases, the sidebands have leaked into the pinhole and prevented the formation of a good image. These spots are excluded from subsequent analysis, but are still reported in our results.
It is important to note that the confocal system used here suppresses the sidebands in the image (in contrast to widefield superoscillatory microscopes [33]) and therefore can image objects much larger than the FoV without distortion [10]. It would be possible, in principle, to design superoscillatory PSFs that would be optimal for imaging a particular sample: as a simple example, if a sample is known to have a limited extent, and the spot designed such that all sidebands fall outside that extent, then imaging is possible at extremely high resolution (up to λ/14 is demonstrated in principle in Ref. [11]). We aim to design spots for a general super-oscillatory microscope and so do not use any a priori knowledge of the sample. Our imaging model uses scalar simulations to reduce the computational complexity. It is well known that the intensity pattern of a spot at the focus of a high-NA lens is significantly different from that predicted by the scalar approach (for an example compare Fig. 4(a) and Fig. 4(b)). However, for the case of imaging only, it can be shown that the scalar behaviour is recovered, and hence the scalar approximation is valid. As we show below, although the total intensity Table 1. Comparison of two-function method optimised spots (1 R . . . 5 R ) with opticaleigenmode-optimised spots (Table 1 of [40]; 1 K . . . 3 K ). Columns headed I.F. contain the improvement factor that has been achieved with the two-function method. Spots 1 R , 2 R , 3 R are compared with spots 1 K , 2 K , 3 K respectively. Spots 4 R and 5 R are both compared with spot 3 K . of the spot is significantly different to the scalar approximation (as expected for a high NA superoscillatory system [49-51]), the transverse component (which is the important quantity in an imaging system) is not badly distorted and can, in fact, be corrected to exactly reproduce the scalar intensity profile. For imaging systems, therefore, the scalar approximation is valid for high NA systems, even though vector approaches must be used for total spot intensity calculations. The simulations in this section assume linearly-polarised input light, an NA of 1.4 and a wavelength of 797 nm. In Fig. 4(a) we show the spot as designed (implicitly scalar) with a FWHM of 151 nm (equal to the design criterion of 0.27 λ/NA) in both the x and y directions. We assume that these spots are created by modulation in the pupil plane of the lens, and so to create a vector simulation of the same spot we Fourier transform the field of the (scalar) spot, giving the pupil plane amplitude and phase E pupil = F −1 {E focus }. We then transform this back to the focal plane, assuming x polarised input E pupil x = E pupil , E pupil y = 0 using the vector method of Boruah and Neil [52] and plot the total intensity in Fig. 4(b). While the total intensity profile is significantly distorted, this is not the important quantity: it has been repeatedly shown [53][54][55][56][57][58] that the longitudinal component of polarization is not registered in an imaging system, and therefore for the purposes of imaging, only the transverse components of polarization should be considered. The transverse intensity profile (see Fig. 4(c)) is a good approximation to the scalar spot with a 4 % error in the FWHM (x: 157 nm, y: 145 nm). We can, however, correct for this error in the transverse intensity as follows (using the notation of Boruah and Neil [52]). For x-polarised input, the y component of intensity remains negligible after focusing, and therefore E focus As expected, this collapses to the scalar approximation in the case of low NA (k x,y,r k 0 , k z ≈ k 0 =⇒ G X ≈ 1). The corrected field at the pupil is given by E corr After applying this correction, the intensity of the transverse polarisation can be brought back exactly to that of the scalar approximation (in Fig. 4(d)). Hence, the use of our (implicitly scalar) spots is valid, even in the high-NA limit and in the extreme case of superoscillatory signals.
Another important factor is the sensitivity of our superoscillatory spots to noise in the imaging system, as superoscillations in general can be disrupted in a noisy system [47,59]. While a full noise analysis is beyond the scope of this work, sample calculations show that these spots are within the reach of modern imaging systems. Considering the proposed applications, random amplitude or phase noise is not appropriate: in imaging systems, the primary cause of distortion is optical aberrations in the lenses, commonly parametrised using the Zernike modes [60]. As an example, we investigate the sensitivity of the spot of Fig. 4 to varying quantities of aberration of each of four Zernike modes (spherical aberration, coma, defocus and astigmatism). We find that, in general, about 0.05 waves (peak-to-peak) of aberration are required to change the spot FWHM by 5 %. About 0.1 waves cause a similar distortion in a conventional Airy disc focus. Hence, the superoscillatory spots are more sensitive to aberration than conventional spots, but only by approximately a factor of two. A good confocal microscope, especially if equipped with adaptive optical correction, should be able to produce the practical superoscillatory spots analysed below.

Simulation results
Many different approaches are used in the literature to report spot parameters and performance.
Here we compare our spots with Kosmeier et al. [40], where the optical eigenmode method is used to generate the spots, as they give clear tabulation of all the relevant parameters of both the spot and its imaging performance. As we show below, the two-function method generates either: spots with the same imaging performance as the eigenmode method but with considerably more energy in the central spot; or with similar energies but improved imaging resolution.
To perform this analysis we calculate five commonly used metrics: • the Strehl ratio, S: the ratio of the peak intensities of the superoscillatory spot to the Airy disk; • the relative intensity, I rel : the ratio of the peak intensity of the superoscillatory spot to the peak intensity of the nearest intense sideband (defined as having peak intensity≥ 20% of the superoscillatory spot). Note that this is the inverse of the intensity ratio (IR) used above; • the relative resolution, R rel , in two-point imaging: the ratio of the lateral resolution in the superoscillatory image to the lateral resolution of a conventional image; • the relative spot size: the ratio of the FWHM of the superoscillatory spot (w spot ) to the FWHM of the Airy disk (w Airy ); • the relative isolation: d SB /w Airy , where d SB is the distance to the nearest intense sideband.
The Rayleigh criterion for two circular apertures in the case of a confocal microscope is given by R = 0.56λ/NA [61]. For our simulated parameters of λ = 797 nm and NA = 1.4 the theoretical (lateral) resolution is R = 318.8 nm. Our simulations give R sim = 319 nm which is in excellent agreement with the theoretical result.
We demonstrate the power of our method by comparing five specific spots with those in Table 1 of Kosmeier et al. [40]. This comparison is shown in Table 1, with our spots denoted 1 R . . . 5 R and those from Kosmeier et al. denoted 1 K . . . 3 K . Spots 1 R and 2 R (also illustrated in Fig. 5(b) and Fig. 5(c) respectively) have the same resolution improvement as spots 1 K and 2 K , but with significantly more energy in the spot: spot 2 R has a 2 times better I rel and 31 times better Strehl ratio. Spot 3 R (in Fig. 5(d)) has the same imaging resolution as spot 3 K and a 46 times improvement in I rel and 3340 times improvement in Strehl ratio. This large difference in intensity leads to a qualitative difference in practicality: spot 3 R (I rel = 0.38) is a realistic proposition for experiments whereas the respective spot in [40] with I rel = 0.0083 is unlikely to be realised in practice. To demonstrate improved resolution with similar intensity measures, spots 4 R and 5 R of Table 1 (shown in Fig. 5(e) and Fig. 5(f)) have up to 26 % higher relative resolution than spot 3 K , significantly better Strehl ratio (3 and 1 orders of magnitude respectively), and better or equal I rel . However, although the relative intensity of those spots is much better than 3 K , it remains low (I rel < 0.04 -i.e. the central spot has 25 times lower intensity than the sidebands), and so these spots are probably of limited use in real imaging. The improvement we obtain in efficiency and/or resolution is driven by two key factors: the optimisation of intensity ratio (rather than the energy measures used in other works); and independent control of spot size and FoV. In our confocal system, the main distortion of the image is light scattered from the sidebands "leaking" into the pinhole. This level of leakage is controlled primarily by the intensity of the nearest significant sideband, and hence controlling the intensity of that sideband is of paramount importance to imaging quality. The importance of independent control of spot size and FoV is demonstrated in Fig. 6, where we cover a much larger fraction of the spot size/isolation parameter space than is reachable using the optical eigenmode method (compare Fig. 6 with Fig. 2 of Kosmeier et al. [40]) or other similar single-parameter optimisation techniques. This increased coverage allows us to find the spots that work best for a particular imaging situation, and hence improve the obtained resolution. A full analysis of this follows in the next section.

Analysis of superoscillatory trade-offs
It is generally accepted that superoscillatory spot intensity is reduced as the spot size becomes smaller and/or the distance to the sidebands is increased. Although this qualitative rule holds in general, and the relationship has been quantified in various specific cases, there is no general quantitative rule for this effect. In the following we examine the form of this relationship for our two-function optimised spots.
In Fig. 6, the Strehl ratio (S), the relative intensity (I rel ), and the relative resolution (R rel ) have been mapped with coloured circles against the relative spot size (w spot /w Airy ) and the relative isolation (d SB /w Airy ). At the top right of all graphs in Fig. 6 there is a group of spots with large relative spot sizes (approx 0.7-1) and large relative isolations (more than ten times the width of the Airy disk). These spots are very similar in form to the standard Airy disk and therefore all have Strehl ratios and relative resolutions of 1, and large relative intensities (> 10). There is another group in the same range of relative spot sizes, but with low relative isolations (< 2) in the bottom right of each panel of Fig. 6 with similar features to the one at the top of the graphs. These two groups demonstrate a transition between the Airy disk and the superoscillatory spots in which the features of the Airy disk are still dominant.
In Fig. 6(c) there is a clear trend towards higher resolution as the FoV decreases. At small isolations, however, the resolution deteriorates due to the sidebands leaking into the pinhole of the simulated confocal microscope and distorting the image (marked with stars). Figure 6(b) shows that those spots with both a small relative spot size and large relative isolation (towards the top left of the graph) have low relative intensities, as expected given the nature of the trade-offs in superoscillations. Most of these spots (shown by open circles) still produce good images in simulations, although those with the smallest spot sizes do not (stars). However, it should be noted that relative intensities of less than about 0.05 may still be impractical for experimental imaging simply due to loss of signal. From Fig. 6(b) and Fig. 6(c) it can be seen that, provided the relative isolation is large enough, a good image can be formed even if the relative intensity is very low (starred points are mainly concentrated in the small FoV region). Interestingly, this means that the smallest spots do not necessarily give the best resolved image.

Conclusion
We have introduced a new method for constructing superoscillatory spots that allows independent control of the spot size and field of view, and optimises the intensity ratio, which is particularly important in imaging experiments. The two-function algorithm runs in 20 ms, and gives excellent spots, showing significant improvements in resolution (for the same efficiency of focusing) or efficiency (for the same resolution) over comparable methods in the literature. The independent control of spot size and field of view allows a spot to be designed, in real time, for the precise needs of the optical setup and experiment at hand. Superoscillatory microscopes based on these spots should allow superoscillatory imaging down to resolutions of around λ/4, with experimentally reasonable parameters, simply by shaping light in the pupil plane of the microscope.