Design of large scale plasmonic nanoslit arrays for arbitrary mode conversion and demultiplexing

We present an iterative design method for the coupling and the mode conversion of arbitrary modes to focused surface plasmons using a large array of aperiodically randomly located slits in a thin metal film. As the distance between the slits is small and the number of slits is large, significant mutual coupling occurs between the slits which makes an accurate computation of the field scattered by the slits difficult. We use an accurate modal source radiator model to efficiently compute the fields in a significantly shorter time compared with three-dimensional (3D) full-field rigorous simulations, so that iterative optimization is efficiently achieved. Since our model accounts for mutual coupling between the slits, the scattering by the slits of both the source wave and the focused surface plasmon can be incorporated in the optimization scheme. We apply this method to the design of various types of couplers for arbitrary fiber modes and a mode demultiplexer that focuses three orthogonal fiber modes to three different foci. Finally, we validate our design results using fully vectorial 3D finite-difference time-domain (FDTD) simulations. © 2014 Optical Society of America OCIS codes: (060.4230) Multiplexing; (130.3120) Integrated optics devices; (240.6680) Surface plasmons; (350.4238) Nanophotonics and photonic crystals. References and links 1. L. Yin, V. K. Vlasko-Vlasov, J. Pearson, J. M. Hiller, J. Hua, U. Welp, D. E. Brown, and C. W. Kimball, “Subwavelength focusing and guiding of surface plasmons,” Nano Lett. 5, 1399–1402 (2005). 2. Z. Liu, J. M. Steele, W. Srituravanich, Y. Pikus, C. Sun, and X. Zhang, “Focusing surface plasmons with a plasmonic lens,” Nano Lett. 5, 1726–1729 (2005). 3. Y. Akimov, W. S. Koh, and K. Ostrikov, “Enhancement of optical absorption in thin-film solar cells through the excitation of higher-order nanoparticle plasmon modes,” Opt. Express 17, 10195–10205 (2009). 4. D. S. Ly-Gagnon, K. C. Balram, J. S. White, P. Wahl, M. L. Brongersma, and D. A. B. Miller, “Routing and photodetection in subwavelength plasmonic slot waveguides,” Nanophotonics 1, 9–16 (2012). 5. D. A. B. Miller, “All linear optical devices are mode converters,” Opt. Express 20, 23985–23993 (2012). 6. N. Riesen, J. D. Love, and J. W. Arkwright, “Few-mode elliptical-core fiber data transmission,” IEEE Photon. Tech. L. , IEEE 24, 344–346 (2012). #197391 $15.00 USD Received 10 Sep 2013; revised 26 Oct 2013; accepted 28 Oct 2013; published 6 Jan 2014 (C) 2014 OSA 13 January 2014 | Vol. 22, No. 1 | DOI:10.1364/OE.22.000646 | OPTICS EXPRESS 646 7. J. D. Love and N. Riesen, “Mode-selective couplers for few-mode optical fiber networks,” Opt. Lett. 37, 3990– 3992 (2012). 8. N. Riesen and J. D. Love, “Weakly-guiding mode-selective fiber couplers,” IEEE J. Quantum Electron. 48, 941–945 (2012). 9. T. Sakamoto, T. Mori, T. Yamamoto, N. Hanzawa, S. Tomita, F. Yamamoto, K. Saitoh, and M. Koshiba, “Modedivision multiplexing transmission system With DMD-Independent low complexity MIMO processing,” J. Lightwave Technol. 31, 2192–2199 (2013). 10. C. P. Tsekrekos and D. Syvridis, “All-fiber broadband mode converter for future wavelength and mode division multiplexing systems,” IEEE Photon. Tech. L. 24, 1638–1641 (2012). 11. A. M. Bratkovsky, J. B. Khurgin, E. Ponizovskaya, W. V. Sorin, and M. R. T. Tan, “Mode division multiplexed (MDM) waveguide link scheme with cascaded Y-junctions,” Opt. Commun. 309, 85–89 (2013). 12. G. Stepniak, L. Maksymiuk, and J. Siuzdak, “Binary-phase spatial light filters for mode-selective excitation of multimode fibers,” J. Lightwave Technol. 29, 1980–1987 (2011). 13. M. Salsi, C. Koebele, D. Sperti, P. Tran, H. Mardoyan, P. Brindel, S. Bigo, A. Boutin, F. Verluise, and P. Sillard, “Mode-division multiplexing of 2 x 100 Gb/s channels using an LCOS-based spatial modulator,” J. Lightwave Technol. 30, 618–623 (2012). 14. J. Carpenter and T. D. Wilkinson, “Characterization of multimode fiber by selective mode excitation,” J. Lightwave Technol. 30, 1386–1392 (2012). 15. J. A. Carpenter, B. C. Thomsen, and T. D. Wilkinson, “Optical vortex based mode division multiplexing over graded-index multimode fibre,” (Optical Society of America, 2013), OSA Technical Digest (online),OTh4G.3+. 16. A. E. Willner, J. Wang, and H. Huang, “A different angle on light communications,” Science 337, 655–656 (2012). 17. J. Leach, M. J. Padgett, S. M. Barnett, S. Franke-Arnold, and J. Courtial, “Measuring the orbital angular momentum of a single photon,” Phys. Rev. Lett. 88, 257901 (2002). 18. G. C. G. Berkhout, M. P. J. Lavery, J. Courtial, M. W. Beijersbergen, and M. J. Padgett, “Efficient sorting of orbital angular momentum states of light,” Phys. Rev. Lett. 105, 153601 (2010). 19. H. Bulow, “Optical-mode demultiplexing by optical MIMO filtering of spatial samples,” IEEE Photon. Tech. L. 24, 1045–1047 (2012). 20. H. Bulow, H. Al-Hashimi, and B. Schmauss, “Spatial mode multiplexers and MIMO processing,” in OptoElectronics and Communications Conference (OECC), 2012 17th, (IEEE, 2012), pp. 562–563. 21. K. H. Wagner, “Mode group demultiplexing and modal dispersion compensation using spatial-spectral holography,” in IEEE Photonics Society Summer Topical Meetings, Space Division Multiplexing for Optical Communications, (2013), pp. 89–90. 22. Y. Jiao, S. Fan, and D. A. B. Miller, “Demonstration of systematic photonic crystal device design and optimization by low-rank adjustments: an extremely compact mode separator,” Opt. Lett. 30, 141–143 (2005). 23. V. Liu, D. A. B. Miller, and S. Fan, “Ultra-compact photonic crystal waveguide spatial mode converter and its connection to the optical diode effect,” Opt. Express 20, 28388–28397 (2012). 24. D. A. B. Miller, “Self-aligning universal beam coupler,” Opt. Express 21, 6360–6370 (2013). 25. D. A. B. Miller, “Self-configuring universal linear optical component,” Photonics Research 1, 1–15 (2013). 26. D. Miller, “Establishing optimal wave communication channels automatically,” J. Lightwave Technol. 31, 3987– 3994 (2013). 27. D. A. B. Miller, “Reconfigurable add-drop multiplexer for spatial modes,” Opt. Express 21, 20220–20229 (2013). 28. T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and P. A. Wolff, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature 391, 667–669 (1998). 29. H. J. Lezec, A. Degiron, E. Devaux, R. A. Linke, L. Martin-Moreno, F. J. Garcia-Vidal, and T. W. Ebbesen, “Beaming light from a subwavelength aperture,” Science 297, 820–822 (2002). 30. J. Lin, J. P. B. Mueller, Q. Wang, G. Yuan, N. Antoniou, X.-C. Yuan, and F. Capasso, “Polarization-controlled tunable directional coupling of surface plasmon polaritons,” Science 340, 331–334 (2013). 31. F. Afshinmanesh, J. S. White, W. Cai, and M. L. Brongersma, “Measurement of the polarization state of light using an integrated plasmonic polarimeter,” Nanophotonics 1, 125–129 (2012). 32. S. Zhang, W. Fan, N. C. Panoiu, K. J. Malloy, R. M. Osgood, and S. R. J. Brueck, “Experimental demonstration of near-infrared negative-index metamaterials,” Phys. Rev. Lett. 95, 137404 (2005). 33. C. Zhao and J. Zhang, “Binary plasmonics: launching surface plasmon polaritons to a desired pattern,” Opt. Lett. 34, 2417–2419 (2009). 34. T. Tanemura, K. C. Balram, D. S. Ly-Gagnon, P. Wahl, J. S. White, M. L. Brongersma, and D. A. B. Miller, “Multiple-wavelength focusing of surface plasmons with a nonperiodic nanoslit coupler,” Nano Lett. 11, 2693– 2698 (2011). 35. L. Li, T. Li, S. Wang, S. Zhu, and X. Zhang, “Broad band focusing and demultiplexing of in-plane propagating surface plasmons,” Nano Lett. 11, 4357–4361 (2011). 36. S.-H. Chang, S. K. Gray, and G. C. Schatz, “Surface plasmon generation and light transmission by isolated nanoholes and arrays of nanoholes in thin metal films,” Opt. Express 13, 3150–3165 (2005). #197391 $15.00 USD Received 10 Sep 2013; revised 26 Oct 2013; accepted 28 Oct 2013; published 6 Jan 2014 (C) 2014 OSA 13 January 2014 | Vol. 22, No. 1 | DOI:10.1364/OE.22.000646 | OPTICS EXPRESS 647 37. H. Liu and P. Lalanne, “Microscopic theory of the extraordinary optical transmission,” Nature 452, 728–731 (2008). 38. H. Liu and P. Lalanne, “Comprehensive microscopic model of the extraordinary optical transmission,” J. Opt. Soc. Am. A 27, 2542–2550 (2010). 39. X. Huang and M. L. Brongersma, “Rapid computation of light scattering from aperiodic plasmonic structures,” Phys. Rev. B 84, 245120 (2011). 40. Z. Ruan and M. Qiu, “Enhanced transmission through periodic arrays of subwavelength holes: the role of localized waveguide resonances,” Phys. Rev. Lett. 96, 233901 (2006). 41. E. Popov, M. Neviere, S. Enoch, and R. Reinisch, “Theory of light transmission through subwavelength periodic hole arrays,” Phys. Rev. B 62, 16100 (2000). 42. M. Paulus, P. G. Balmaz, and O. J. Martin, “Accurate and efficient computation of the Green’s tensor for stratified media,” Phys. Rev. E 62, 5797 (2000). 43. M. Paulus and O. J. F. Martin, “Light propagation and scattering in stratified media: a Green’s tensor approach,” J. Opt. Soc. Am. A 18, 854–861 (2001). 44. T. Tanemura, P. Wahl, S. Fan, and D. A. B. Miller, “Modal source radiator model for arbitrary two-dimensional arrays of subwavelength apertures on metal films,” IEEE J. Sel. Top. Quant. 19, 4601110 (2013). 45. P. Wahl, D. S. Ly Gagnon, C. Debaes, J. Van Erps, N. Vermeulen, D. A. B. Miller, and H. Thienpont, “B-CALM: an open-Source multi-GPU-based 3D-FDTD with multi-pole dispersion for plasmonics,” Prog. Electromag. Res. 138, 467–478 (2013). 46. S. Kim, H. Kim, Y. Lim, and B. Lee, “Off-axis directional beaming of optical field diffracted by a single subwavelength metal slit with asymmetric dielectric surface gratings,” Appl. Phys. Lett. 90, 051113 (2007). 47. B. Lee, S. Kim, H. Kim, and Y. Lim, “The use of plasmonics in light beaming and focusing,” Prog. Quant. Electron. 34, 47–87 (2010). 48. J. S. White, G. Veronis, Z. Yu, E. S. Barnard, A. Chandran, S. Fan, and M. L. Brongersma, “Extr


Introduction
Over the last decade, surface plasmon polaritons have generated a substantial amount of research interest because of their unique property of being confined to a metallic surface and the associated ability of metallic structures to confine light to sub-wavelength volumes. Surface plasmons have been used in a broad area of applications such as subwavelength focusing [1,2], light-harvesting [3], and photodetection [4]. Among those applications, mode coupling and demultiplexing are very interesting as any linear optical component can be described as a converter between two sets of orthogonal optical modes [5]. This problem of mode splitting or conversion has recently become very topical because of the possibility of using multiple modes in optical fibers. Various techniques have been explored in the literature to separate modes optically, including fiber or waveguide coupler-based devices [6,7,8,9,10,11], phase plates or spatial light modulators with free-space optics [12,13,14,15,16,17,18], spatial sampling into planar light-wave circuits or silicon photonics with subsequent waveguide interferometers [19,20], and holographic approaches [21]. Previous work with dielectric structures has shown the possibility of arbitrary mode coupling in photonic-crystal-like structures with custom regions [22,23]. With the additional incorporation of detectors and local feedback loops together with phase shifters and interferometers, schemes have been proposed that can automatically convert between arbitrary modes [24,25,26,27].
Structures that couple and demultiplex surface plasmons have mostly been based on metallic films perforated with subwavelength apertures or grooves [see Fig. 1]. Those have been used to show a variety of physical phenomena such as extraordinary optical transmission (EOT) [28], beam steering [29], beam steering with polarization control [30,31] and negative refraction [32]. While most work on these arrays has made use of periodic arrays, aperiodic arrays of metallic slits have also been used to achieve binary surface plasmon polariton fo- cusing [33] and spectrally discrete wavelength demultiplexing [34,35], which is not readily possible with periodic approaches.
In an array of slits, mutual coupling between the slits occurs as surface plasmons generated by one slit can be scattered by another slit and couple back into the first one. This phenomenon is numerically tractable for periodic arrays through the use of periodic boundary conditions [36] and for 1-dimsional (1D) and quasi-1D slit arrays [37,38,39]. However, an accurate modeling of the mutual coupling in aperiodic arrays of slits has often required rigorous exact 3-D fully vectorial methods such as Finite Difference Time Domain (FDTD) [40], Finite Elements Method (FEM) [39], rigorous coupled wave analysis [41] or Green's Tensor approaches [42,43]. Consequently, previous work on aperiodic arrays was operated in a regime where the slits are located relatively sparsely, so that the mutual coupling between slits could be safely neglected [34].
Recently, we have developed a novel semi-analytic model, namely the modal source radiator model, which allows us to compute the mutual coupling among arbitrarily located apertures with a significantly reduced computational cost [44]. We present a design method based on the modal source radiator model that, in contrast to the schemes used in [34], can account for the scattering by the slits of both the direct excitation and the surface plasmons all the slits generate. This knowledge allows us to calculate the optimal position of slits that scatter the surface plasmon waves more than they scatter the excitation beam. In addition, we are able to detect when a slit is "blocking" a focused surface plasmon and is basically scattering the field one is trying to focus; a better device is achieved by removing this particular slit. This paper is structured as follows: first we provide a brief description of the modal source radiator model which allows us to calculate the scattered field in a quantitatively accurate way. Next we describe our design method which is built on the modal source radiator model. We conclude with some design examples: the first of these allows the focusing of a higher-order The field inside the slits is decomposed into TE 01 and TE 02 eigenmodes of the slits.
L 31 fiber mode to an arbitrary location and the second is a structure that demultiplexes three orthogonal fiber modes to three different foci as depicted in Fig. 1. The latter design is confirmed with full-field FDTD simulations to prove the validity of our calculations.

Radiator source model
In the modal source radiator model, the field inside the slits is decomposed into a superposition of a finite number of eigenmodes, such as TE 01 -like and TE 02 -like modes [see Fig. 2], which are the guided modes that would propagate in the z-direction inside the slit if the metal was infinitely thick [44]. Those modes can also be seen as the plasmonic equivalents of TE modes in rectangular waveguides made of a perfect conductor, where E z = 0. In addition, in the examples we consider below, we limited this decomposition to the TE 01 -like and TE 02 -like modes shown in Fig. 2. We can therefore write the |E and |H fields inside the slit using the Dirac notation as In Eq. (1), α runs over all the eigenmodes in all the slits, q α is the propagation constant of the slit mode, A α and B α are the complex amplitudes of the downward and the upward propagating modes respectively and |e ± α (|h ± α ) represents the mode profiles of the slits. The modal source radiator model is constructed based on the fact that the transverse magnetic field at the top interface |H I t can be expressed in two different ways using the following modal basis: where Eq. (2) is only valid for the area covered by the apertures of the slits. In Eqs. (2) and (3) |h tα is the transverse part of the magnetic field of the slit modes, |H I 0 is the field that would be present at the top interface in the absence of any slits and G I H |e ± α is the transverse magnetic field of the radiation pattern at the top metal interface, which is generated by the downward and upward propagating slit modes. Similar equations can be written for the bottom metal interface. Using Eqs. (2) and (3) at the top interface and analog expressions at the bottom interface, A α and B α can be calculated for all slits using a low-rank linear set of equations whose coefficients can be numerically calculated for arbitrary slit locations using a set of numerical simulations on a single slit [44]. The transverse magnetic and the normal electric field can be very efficiently calculated at any point on the interface using Eq. (3) and respectively, where G I E |e ± α can be calculated for every field component using one simulation on a single slit [44]. Figures 3 and 4 show the optimization algorithm we employ to design the mode demultiplexer, which couples multiple orthogonal optical modes to SPPs and focuses each mode towards a different position on either the top or bottom interface of the metal layer. Our iterative design method, which is illustrated in Fig. 3, uses the modal source radiator model at every iteration step. We define n s as the total number of slits and n f as the total number of modes and aim to focus the optical mode q onto the desired location r F q = (x F q , y F q ). First, in Step 1 of Fig. 3, we compute the full field of the slit array under illumination of each optical mode q, using the modal source radiator model, which we have described in the previous section. Consequently, for a certain slit configuration, A α and B α are known for every mode in every slit and the electric field in the normal direction |E I/II z can be calculated at both interfaces in a very computationally efficient way using Eq. (4).

Design method
Then, in Step 2, the complex contribution Ψ q,k of each slit to the field at each focal point is isolated by performing only a summation of all the modes α in slit k in Eq. (4) and evaluating at r F q as shown in Eq. (5).
The goal of the design method is to place the slits in such a way that for each focus q, all the n s contributions Ψ q,k would be in phase. Naturally, it is not possible to do this perfectly for a higher number of foci, or when strong coupling occurs between the slits, and we have to find a compromise using an optimization scheme that aims at minimizing the phase difference: where arg(Ψ q,k ) is the phase of Ψ q,k and is the target phase for each focus and is chosen to be fixed throughout the entire optimization process. The choice of the target phase is arbitrary in the example cases we study below, where we choose to optimize the intensity. Nevertheless, an arbitrary target phase can be set using our algorithm, should that be of importance for a particular application. Also, we found that for the examples we use in this paper, setting fixing the phase does affect the final position of the slits, but the overall performance remains the same. Using Eqs. (5)-(7), ∆ n Φ q,k is calculated for all slits at step 2.
In step 3, we calculate for each slit k the position r min k where ∆Φ q,k is as small as possible, compromising between all foci q. This is done by finding the minimum of a cost function CF k (r(x, y)) that is set up to be minimal when ∆ n Φ q,k is small for each focus in an area close to the slit and not in the immediate vicinity of other slits. In addition, CF k needs to take the field  4. Illustration of the design method. We assume that the surface plasmon propagates with phase velocity k spp . For slit S k , the dashed blue circles represent the locations where ∆Φ q,k = arg(Ψ q,k ) and the full black circles represent the locations where ∆Φ q,k = 2nπ. Also, δ 1 = ∆Φ q,k /k spp and δ 2 = (2π − ∆Φ q,k )/k spp . The design method has the goal to place the slit in a location where ∆Φ q,k ≈ 2nπ, like the area inside the red circle.
intensity in each focus into account so as to optimize towards equal intensities in each focus. We choose CF k as : where k spp denotes the surface plasmon phase velocity at the metallic interface and denotes the shortest distance between r and the circles near slit k where ∆Φ q,k = 2nπ, as depicted by full black circles in Fig. 4. Naturally, we cannot check this for all n in practice and only perform the minimization for |n|<2 throughout this work. In which represents a Gaussian distribution of width σ g around the slit k so that the minimum of CF k would be in the neighborhood of slit k located at r k . In addition, we introduce which is a sum over Gaussian distributions with width σ nn centered around all the slits l (except slit k); we do this to avoid having several slits evolve towards the same optimum, where they would collide.
To obtain a uniform intensity distribution along the n f foci we introduce a weight factor w q in Eq. (8) to give more weight to a focus with a relatively lower intensity. w q is initially defined as 1 for all q and is updated in step 4 according to the relation where max denotes the maximum over all foci and v is chosen as 0.1 throughout this work to avoid instability in the optimization [34].
In the examples we discuss later, we found that when ν << 0.1 we converged to a solution where one mode would be significantly favoured over the other modes, while in the case that ν >> 0.1, the algorithm would not converge at all as it would oscillate between solutions where one mode would be favoured over the others. Throughout the optimization, some slits will lie in the path of the focused surface plasmon generated by other slits and will consequently scatter the surface plasmon more than they scatter excitation source. This surface plasmon already has the correct phase and the scattering caused by this slit will not improve the merit function, as those slits are just "blocking" the propagation of the surface plasmon. Such slits can be detected as ∆ n Φ q,k does not tend to converge to a multiple of 2π. For this reason every n sc iterations we remove the slits for which where ℜ and ℑ denote the real and the imaginary part of a complex number respectively. P k is the sum over all foci of the projections in the complex plane of Ψ q,k , i.e., the contribution of each slit to focus q onto φ re f q , which is the target phase for each focus. The positions of the slits that remain will continue to be updated until the iteration number reaches the next multiple of n sc . Throughout this work we performed this operation (step 5) every thirty iterations (n sc = 30).
Finally, in step 6 , we move each slit k to a new position r new k given by where r min k = min(CF k (r)) is the minimum of the merit function for each slit k and w = 0.1 was used throughout this work. In the examples we discuss later, we found that when w << 0.1 we converged to a solution, albeit very slowly while in the case that w >> 0.1, the algorithm would not converge at all. During this operation it may happen that despite the introduction of G n k (r), two or more slits still converge to the same location when all potential local minima have been taken by other slits. To avoid that, after moving the slits to r new k , we check for slit pairs that are too close to each other and remove one of the two.

Examples
In this section, we apply our optimization method to two different cases with arrays of slits with a width w s = 60 nm and a length l s = 300 nm in a silver layer of 200 nm thick on top of a silicon oxide layer with refractive index 1.45. We use a wavelength of 800 nm. First, we illustrate that our design method can be successfully used to design slit arrays that can focus an arbitrary mode onto an arbitrary spot on bottom side of the metallic sheet at the dielectric metal interface. To this end, we design two different slit arrays that focus an L 31 mode of a Corning ® SMF-28e ® optical fiber which has a diameter of 8.2 µm, a core index of 1.4585 and a cladding index of 1.4533 [ Fig. 5(a)] onto foci at the oxide-metal interface at r F 1 = (0, 0) µm and r F 1 = (12, 0) µm respectively. In both cases we start with a 40 x 34 array of slits that is initially centered around (0, 0) and with an initial slit separation of 400 nm and 500 nm in the x and y direction respectively. Also, in Eq. (8) we use σ g = 0.5 µm, σ nn = 0.05 µm ; k spp was computed using FDTD simulations and the mode profile itself was numerically calculated using COMSOL Multiphysics 3.5a. Our algorithm converged to a solution within 40 iterations and took about 50 seconds per iteration on a regular desktop computer with a Intel(R) Xeon(R) CPU X5650-2.67GHz processor and using 8 GB of RAM. In Figs. 5(b) and 6(a), the intensity of the normal electric field at the oxide metal interface is depicted for both cases. This intensity is normalized to the average   Fig. 7. |E z (r F 1 )| 2 at each iteration step for r F 1 = (0, 0) µm (a) and for r F 1 = (12, 0) µm (b). One clearly observes the improvement in field intensity obtained at iteration 30 when we remove the slits for which P k < 0. All the fields are normalized to the average field intensity in the fiber core.
intensity of the electric field inside the fiber core. The final positions of the slits are drawn as white lines in Figs. 5(b) and 6(a). The normalized field intensity at the focal point at each iteration step is depicted in Fig. 7. For both cases, it can be seen that the field intensity increases significantly throughout the optimization before saturating. The discontinuous improvement that occurs at the 30th iteration is caused by the removal of the slits for which P k < 0, which can be thought of as "blocking" the focused surface plasmon as explained earlier. This can also be seen in Fig. 5(b), where slits have been removed around the focus and Fig. 6(a) where slits have been removed in the area where the focused surface plasmon has already reached substantial intensity as it is being focused onto the focal point.
In a second example, we apply our optimization method to design a slit array aiming at demultiplexing three orthogonal fiber modes. In addition, to validate our method, we confirm the field patterns with full-field FDTD simulations. In order to keep the FDTD simulation size reasonable, we choose a smaller fiber core and a smaller array of slits. More specifically, we choose to demultiplex the L 01 , L 21 and L 22 [see Figs. 8(a)-8(c)] modes of a multimode fiber with a diameter of 2.2 µm, a core index of 2 and a cladding index of 1.6 onto three foci located at r F 1 = (7, 0) µm, r F 2 = (7, −3) µm and r F 3 = (7, 3) µm respectively. We start with a 13 x 11 array of slits that is initially centered around (0, 0) and with a initial slit separation of 400 nm and 500 nm in the x and y direction respectively. All other parameters are equal to those used in the previous example. As can be seen from Fig. 9, our design method converged to a solution containing 48 slits after 110 iterations and three clearly distinct intensity peaks can be observed at x = 7 µm corresponding to the yellow dotted line in Figs. 8(d)-8(f). To confirm the correct implementation of the modal source radiator model, we perform a full-field 3D-FDTD simulation of the area illuminated by the source. In order to limit the simulation size the focus plane was not included in the simulation area and only the area inside the yellow dotted box in Figs. 8(d)-8(f) was simulated with a cell size of 5 nm. This resulted in an FDTD grid of 1200x1200x400 cells. The FDTD simulation was performed using Belgium CAlifornia Light Machine [45] on 16 NVIDIA Tesla M2070 GPUs simultaneously. In this configuration, a simulation of 40000 time-steps completed in less than 140 minutes. As can be seen from Fig. 10, an almost perfect match is obtained between the modal source radiator model and the full-field FDTD, which confirms the validity of our model.

Estimating the overall efficiency
An exact calculation of the overall efficiency is difficult, since the modes to which we couple are not well defined nor did we include a material in which the coupled light would be ab-  sorbed. Nevertheless, it is possible to estimate the efficiency from the field intensity |E z (r)| 2 as presumably most of the energy at the metal oxide interface is carried by surface plasmons. Therefore, we assume that the in-plane Poynting vector which denotes the integrated components of the Poynting vector that are parallel to the metal sheet over the z direction, is proportional to |E z (r)| 2 or In the Eq. (14), P(x, y, z) is the Poynting vector and z denotes a unit vector along the z direction. We calculate the proportionality factor τ in Eq. (15) numerically using a simulated surface plasmon propagating in the x direction and then estimate the efficiency η from η =´y 2 y 1 τ|E z (x F , y)| 2 dỹ P source z dxdy (16) where˜P source z dxdy denotes the integrated Poynting vector of the excitation source and y 1 and y 2 are chosen at 1 µm around at either side of y F . Using this method we estimated total efficiencies of between 0.8 percent for the designs depicted in Figs. 6 and 8 and 4 percent for the design depicted in Fig. 5. As is shown in Figs. 5 and 6, most of the surface plasmon energy is concentrated in the focus. However, for the demultiplexing design we observe multiple "peaks" of surface plasmon energy. While the main focus still contains most surface plasmon energy, we calculated that η p ≈ 40% , where η p is the ratio of the surface-plasmon energy that is actually concentrated in the focus, over the integrated surface plasmon energy around the slits. For this particular integration, we use a square with side length 2x f centered at (0,0). Also, η and η p remain stable when using a larger number of slits as the slits are already covering most of the illuminated area. In addition, η and η p remain stable when we permute the location of the foci in our design as shown in Fig. 11, where the intensity is plotted for a design where we optimize the slit positions for demultiplexing the L 01 , L 21 and L 22 fiber modes to r F 1 = (7, 3) µm, r F 2 = (7, 0) µm, r F 3 = (7, −3) µm, respectively. The behaviour of this design is completely analogous to the one shown in Fig. 8 and we obtained efficiencies η ≈ 0.8% and η p ≈ 40% . From Figs. 8 and 11, we observe that most of the non-focussed light is back scattered. As better unidirectional periodic grating couplers can be designed under illumination at an off-axis angle [46,47], it may be possible to improve η p and η by optimizing the design for illumination at an angle.
We emphasize that the goal of this paper is to explore an optimization scheme for mode couplers and demultiplexers that is based on the source radiator model for a certain slit. Consequently, the slits' size, metal thickness and wavelength are fixed and those parameters were somewhat arbitrarily chosen to illustrate the optimization of the positions of a number of given slits. In our structures, most of the light is still coupled into the oxide or reflected back and higher efficiencies may be obtained by optimizing the slit size and metal thickness for a specific wavelength [48]. Also we could consider using grooves instead of slits to prevent transmission when coupling at the top interface [49], or coating of the metal with a dielectric grating, trapping the light at the metal surface. The algorithms presented here could be used for the above-mentioned topologies, as those topologies could be modeled accurately with the source radiator model.

Discussion
We can compare the complexity of our designs here with previous analysis of the complexity required to make a particular optical component. As an example, we take our second design of a mode splitter above. Ref. [50] gives the complexity number N D of real numbers that need to be specified for a given component to be if we do not care about the absolute phase of the outputs [Eq. (11) of [50]]; in Eq. (17), M I is the number of input modes (the dimensionality of the input space), here 3, M C is the number of channels to be coupled through the device (here 3 also), and M O is the dimensionality of the output space. We can see from Fig. 9(b) that we are able to place our outputs along one line, with each output spot spaced approximately by one spot width from the neighboring spots. We might reasonably interpret this as saying we could place any one output at any of 5 distinct positions at the output along this line. In our device, with slits oriented in a "top-bottom" axis in Fig. (8), we can design to route light to the left and right faces, essentially. Light scattering into the "top-bottom" direction is forbidden by the polarization of the light and the symmetry of the slits. Since we could design this device to give its outputs along either the left or the right "faces" of the device, we could argue that our device has at least M O ∼ 10 different output modes it can choose from. Substituting these values in Eq. (17) gives a required minimum number of degrees of freedom of N D = 57 . In our design, we choose the x and y positions of 48 slits, corresponding to 96 degrees of freedom, just under a factor of 2 larger than N D . We found that approximately this number of slits was the smallest that we could use and still effectively separate the modes into the desired output spots. We should expect to use more degrees of freedom than the minimum number N D in a real design because that minimum number does not account for any control over the precise form of the output beams. To control that form, we should expect to be using a larger dimensional output space in practice. Our work here is suggesting a factor of ∼ 2 increase in output dimensionality is sufficient to allow the formation of convenient forms of output beams.

Conclusion
In this paper we successfully used a quantitatively accurate semi-analytic modal source radiator model for arbitrary two-dimensional non-periodic arrays of slits in order to design arbitrary mode couplers and demultiplexers. The use of the modal source radiator model allows us to calculate the contribution of each slit taking both the direct scattering of the source and the scattering of the surface plasmon into account in a quantitatively accurate and computationally efficient way. This ability enables the design of large arrays of closely-spaced slits. Our design method aims at focusing an arbitrary set of orthogonal modes onto different foci, by iteratively moving slits to minimize the phase difference between the contributions of each slit in each focus. Thanks to our model, we can detect the slits that are scattering the focusing surface plasmon more than the source and remove them. We successfully applied our method to design two large aperiodic slit arrays that focus a L 31 mode in the middle of the array itself or to a focal point outside they array. Also we used our design method to calculate aperiodic slit array that focuses the L 01 , L 21 and L 22 fiber modes onto 3 separated foci and confirmed the correctness of our calculations using full field FDTD simulations.