Polarizability tensor retrieval for magnetic and plasmonic antenna design

A key quantity in the design of plasmonic antennas and metasurfaces, as well as metamaterials, is the electrodynamic polarizability of a single scattering building block. In particular, in the current merging of plasmonics and metamaterials, subwavelength scatterers are judged by their ability to present a large, generally anisotropic electric and magnetic polarizability, as well as a bi-anisotropic magnetoelectric polarizability. This bi-anisotropic response, whereby a magnetic dipole is induced through electric driving, and vice versa, is strongly linked to the optical activity and chiral response of plasmonic metamolecules. We present two distinct methods to retrieve the polarizibility tensor from electrodynamic simulations. As a basis for both, we use the surface integral equation (SIE) method to solve for the scattering response of arbitrary objects exactly. In the first retrieval method, we project scattered fields onto vector spherical harmonics with the aid of an exact discrete spherical harmonic Fourier transform on the unit sphere. In the second, we take the effective current distributions generated by SIE as a basis to calculate dipole moments. We verify that the first approach holds for scatterers of any size, while the second is only approximately correct for small scatterers. We present benchmark calculations, revisiting the zero-forward scattering paradox of Kerker et al (1983 J. Opt. Soc. Am. 73 765–7) and Alù and Engheta (2010 J. Nanophoton. 4 041590), relevant in dielectric scattering cancelation and sensor cloaking designs. Finally, we report the polarizability tensor of split rings, and show that split rings will strongly influence the emission of dipolar single emitters. In the context of plasmon-enhanced emission, split rings can imbue their large magnetic dipole moment on the emission of simple electric dipole emitters. We present a split ring antenna array design that is capable of converting the emission of a single linear dipole emitter in forward and backward beams of directional emission of opposite handedness. This design can, for instance, find application in the spin angular momentum encoding of quantum information.

emitters. In the context of plasmon-enhanced emission, split rings can imbue their large magnetic dipole moment on the emission of simple electric dipole emitters. We present a split ring antenna array design that is capable of converting the emission of a single linear dipole emitter in forward and backward beams of directional emission of opposite handedness. This design can, for instance, find application in the spin angular momentum encoding of quantum information.

Introduction
Metallic and dielectric nanoscatterers currently enjoy a surge of interest in photonics, due to the unusual optical properties that may be obtained through a suitable choice of material and geometry. In plasmonics, it is well established that Ag and Au single nanospheres, rods, wires, pyramids, triangles, cubes, stars or core-shell particles, as well as oligomers and arrays of such objects have very distinct scattering resonances that can be used for optical sensing, improvement of LEDs and solar cells, as well as plasmon-enhanced spectroscopy on the basis of large field enhancements near metals at the plasmon resonance [3][4][5][6][7][8][9]. In a related development, the field of metamaterials uses metal split rings, loops as well as so-called cut-wire pairs to generate a strong collective magnetic response [10][11][12]. The effective magnetic permeability and electric permittivity that is achieved arises from the strong electric and magnetic polarization obtained in each building block. Recently, the fields of plasmonics and metamaterials have come together in so-called 'metasurfaces', where non-identical subwavelength resonant scatterers are organized in a plane at subwavelength distances, in order to achieve arbitrary phase and amplitude masks that control the transmission, reflection, refraction and diffraction of light [13][14][15]. In all these developments, the response is fundamentally quantified by the geometrical arrangement of scatterers on the one hand, and the electric and magnetic polarizability of each building block on the other hand.
Numerical methods in electrodynamics play an increasingly important role in the design and understanding of nanostructures in an electromagnetic field. In daily practice, finitedifference time-domain codes (FDTD), finite element simulations (FEM) and boundary element methods (BEM) are used to replicate experiments, and extract expected observables such as transmission and reflection coefficients, or the brightness and directivity of localized sources. Remarkably, it is not common practice to use simulations to extract the fundamental parameter, i.e. the electric and magnetic polarizability, as well as possibly higher-order multipoles, which quantify how a building block scatters. A first step to this goal is a recent paper by Mühlig et al [16] that shows a retrieval of the multipolar moments induced in various metamaterial scatterers for a particular incident field. Here we propose a rapid and accurate method to retrieve electric, magnetic and magneto-electric polarizabilities of meta-atoms that can be applied to any electromagnetic solver (FDTD [17], volume integral equation method (VIE) [18] and BEM [19]). In our specific implementation, this method is based on surface integral equation (SIE) calculations to solve Maxwell equations for electric and magnetic fields exactly, and to calculate the induced effective electric and magnetic surface currents that quantify the scatterer response. We show how one can extract polarizabilities both from the calculated scattered field and, as an alternative method, also from the induced currents. By applying this method to, among other objects, split rings, we retrieve exciting insights regarding the electric and magnetic response of split rings. In particular, several authors have argued that split rings should be viewed either as an electric plus a magnetic dipole each with a polarizability [20], or as a cross-coupled object that also has a magneto-electric response [21]. We show that split rings are strongly magneto-electric, implying that a large magnetic dipole moment is most easily induced by electric driving. The particular phase relation between the electric, magnetic and magneto-electric polarizabilities further implies record-high per-building block optical activity in extinction and scattering. We show that the insights gained from the polarizability tensor can be used to construct new types of plasmonic array antennas that have the directivity of Yagi-Uda antennas [22,23], but with unique polarization properties. In particular, we show how metamaterial antennas allow control over the magnetic dipole content of emission and over the handedness of emitted light. On the basis of this type of control over emission, we envision applications in control of magnetic dipole emitters [24], photon spin angular momentum encoding in single-photon sources and enantioselective spectroscopies that employ near-field enhancement of chirality. This paper is organized as follows. In section 2, we present the SIE method (2.1) and the retrieval of polarizability tensors (2.2). In sections 3 and 4, we benchmark the retrieval for magneto-electric spheres, illustrating the zero-forward scattering paradox of Kerker. In section 5, we discuss the polarizability of split rings. In sections 6 and 7, we demonstrate how, on the basis of the extracted polarizability, split rings can be used for the rational design of antennas for emission control.

Surface integral equation method and α-tensor retrieval
Any electromagnetic problem is completely specified by the Maxwell equations, together with a definition of the source, and the boundary condition. We use the equations in integral form [25]. We divide space into a region 1, defined as the embedding medium that we take to be a homogeneous dielectric with permittivity 1 and permeability µ 1 , and a region 2 that represents the volume occupied by the scattering material of dielectric constant 2 and permeability µ 2 . In the integral equation formalism, it is useful to solve for the electromagnetic response by first finding the effective auxiliary electric and magnetic surface current densities J and M on the interface between medium 1 and medium 2 that are used to satisfy the boundary conditions for continuity of tangential E and H, and normal B and D. Once these surface currents are solved for, they can be used to construct the electromagnetic field solution everywhere. Assuming harmonic time dependence (frequency ω), the currents are set by an electric field integral equation (EFIE), and a magnetic field integral equation (MFIE) that reads Here G(r, r ) is the electric dyadic Green function for each type of homogeneous medium 'i' (with i = 1, 2) and ∇ × G(r, r ) is the curl of the Green function. The integral runs over the surface S that contains the current densities. We used the method introduced by Kern et al [26], which is based on the method of moments [27], and coined the SIE method to solve these equations. In brief, in the SIE method any scatterer is represented by effective electric and magnetic surface current densities J and M that are discretized on finite elements over the surface of the scatterer with the help of the Rao, Wilton and Glisson (RWG) basis functions f n [28]. Consider the surface S meshed with triangles. We define n = 1, . . . , N nodes as the shared edges of the triangles. The basis function f n (r) is zero everywhere except on the triangle pair T ± that shares node n. Here the function is pyramid shaped, with where L n is the length of the shared node, A ± n is the area of the triangle pair and p ± n are the non-shared vertices of the triangles, as explained in [26]. The discretized strength and direction of the currents is accounted for through basis expansion coefficients α n and β n in the following way: By projecting the EFIE and MFIE equations onto the RWG basis functions, the integral equations transform into a set of linear equations for the α n and β n values: In this system of linear equations, the matrix M needs to be found only once and can be used for any incident field that is only contained in q. The matrix M is defined by where Note how M self-consistently contains the interactions between all the discretized current elements, as evident from the appearance of G i (r, r ). After this matrix is calculated, it can be inverted and multiplied by the vector q, which is the projection of the incident field that drives the scatterer into the RWG functions. Specifically, We order the variables so that the vector q has the projections of E inc on the N basis functions as the first N elements, and the projection of H inc on the basis functions as elements N + 1 to 2N .
As an important implementation note, one of the key features of this method is that the Green function, which is singular at r = r , is written as the sum of a smooth G(r, r 0 ) S and singular part G(r, r 0 ) NS as follows: where k i is the wave vector defined as k i = 2π/λ · √ µ i i and R = |r − r 0 |. The singular part of the integral is treated analytically. For a detailed explanation of this separation method, we refer to [26]. Without this separation, the matrix M would be highly inaccurate on its diagonal, as well as for elements D mn and K mn that correspond to close triangles. Moreover, subsequent retrieval of the scattered field from the calculated currents would be highly inaccurate close to the scatterers. Once found, the coefficients α n and β n completely specify the electric and magnetic surface currents that in turn allow one to find the scattered field and the total field by propagating the currents with the aid of the Green function in the following manner: We have implemented the described algorithm in MATLAB, using triangular surface gridding that defines the set of f n (r) as input that we generated using Gmsh [29]. The main objective of our paper is to discuss the retrieval of polarizabilities from the SIE calculations. On the basis of the current contributions found through SIE, two different approaches can be taken in order to find the polarizability tensor. On the one hand, the fields produced by the effective currents can be propagated with the aid of G and ∇ × G (equations (13) and (14)) onto a sphere that is centered around the structure under consideration. The projections of the fields on the sphere on vector spherical harmonics (VSH) directly define the multipole moments through the expansion coefficients a nm and b nm , as explained by Jackson [25, chapter 10] as well as by Mühlig et al [16]. Thus, for this first retrieval method we use two steps. First we use SIE to solve for the fields E and H and then we project these fields onto VSH to find the dipolar moments and hence the polarizability tensor. This means that E and H might as well be found by using any other full-wave calculation, for instance FDTD or FEM. As an alternative method, multipole moments can be defined directly from the current distributions, without calculating fields.
Here we first discuss the multipole expansion method and then the direct definition based on J and M.

Multipole moments based on the projection onto vector spherical harmonics
We use the VSH functions as defined by Mühlig et al [16], which are equivalent to the textbook definition of [30]. As proven in [30] the VSH form a complete and orthonormal set [30]. Therefore the field E(r, θ, φ) found from SIE has a unique expansion where the expansion coefficients a nm and b nm can simply be found by taking the inner product of the calculated E(r, θ, φ) with the VSH functions. Here N nm (r, θ, φ) and M nm (r, θ, φ) are the VSH functions, and the inner product is defined as the integration over the unit sphere. While in principle one could densely sample E(r, θ, φ) on the unit sphere to evaluate the inner product, it is particularly advantageous to use the fact that a discrete spherical harmonic transform is exact for sampling points and weights chosen as consistent with Legendre quadratures for Legendre polynomials of order N + 1. Thereby one can obtain a highly efficient and exact algorithm, that requires only very few field sampling points for multipole expansion coefficients up to order n = N [31], by carefully separating the VSH into ordinary spherical harmonics. As in any discrete Fourier transform, the only caveat for this exact method is that aliasing artifacts may occur if the radiated field contains significant contribution from multipoles of order higher than the truncation order of the transform. Therefore, we use a truncation order N = 5, corresponding to just 50 sampling points on the unit sphere, as we do not expect multipole moments beyond N = 2 to be significant throughout this paper. The coefficients a nm and b nm , and hence the retrieved moments finally, will depend on where the center of the sphere is chosen [16,25,32] but are independent of the radius of the sphere as long as the scatterers are fully enclosed. The dipolar moments p and m are calculated from the coefficients a nm and b nm using the procedure explained in [16].

Multipole moments based on effective currents
As a second method to obtain the induced dipole moments p and m (electric and magnetic dipole moments), we can directly use the effective magnetic and electric current densities M and J calculated as the intermediate solution step in SIE. In particular, where the integration is performed over the surface of the scatterer . In contrast to other brute force methods, SIE naturally provides the effective magnetic and electric currents as an essential part of its solution strategy. In standard implementations of, for instance, FDTD modeling, retrieving these currents with enough numerical accuracy would itself be a challenge. On a standard FDTD Yee grid, these inaccuracies arise from the approximation of curved boundaries into discretized Cartesian grids, as well as from the fact that, in general, the field components and their derivatives are not sampled right on the boundary. Consequently, right at object boundaries large inaccuracies of local fields, fluxes and currents are obtained unless one uses specially improved FDTD algorithms [33]. The definitions of the RWG basis functions imply that where r c− n and r c+ n are the centroid vectors of the two triangles that share node n, L n is the length of the shared line between the two triangles, and finally r 2 and r 3 are the vector positions of the shared vertices of the two triangles. Inserting these results into the discretized form of equations (16) and (17) allows one to retrieve p and m in terms of the coefficients α n and β n : Both the electric dipoles p arising from the effective magnetic currents and the magnetic dipoles m arising from the electric effective currents depend on the choice of origin. One of the potential advantages of the effective current approach over the VSH approach is that one can find the dipolar contributions of single scatterers in close proximity to other scatterers, for instance when examining the physics of multi-element plasmon antennas. Furthermore, one can even envision that one could calculate the multipole moments for structures close to an interface or inside lossy environments, while this is certainly not possible with the VSH approach.

Polarizability tensor retrieval
For both the VSH retrieval method and the direct current-base retrieval method, the electric dipole moment p and the magnetic dipole moment m are retrieved given a particular incident field. Motivated by recent developments in the field of metamaterials, we propose retrieving polarizability tensors that specify the response for any incident field, rather than induced moments for a particular incident field. We focus on objects with an electric and a magnetic dipole response, which we expect to be fully captured by a 6 × 6 polarizability tensor [21]: In this tensor the upper diagonal 3 × 3 blockᾱ E is the usual electric polarizability tensor, while the lower diagonal blockᾱ H is the magnetic polarizability tensor. The off-diagonal blocks represent magneto-electric response, i.e. the electric (magnetic) moment that might be induced through magnetic (electric) driving. This form of the polarizability tensor is commonly used in the field of bi-anisotropic and chiral media [34]. Evidently, one should choose six independent incident conditions, retrieve the induced moments and apply matrix inversions to obtain In our work, we use as incidence conditions standing waves constructed as plane waves incident from opposing Cartesian directions. To construct six independent conditions, we use the three Cartesian axes as incident directions, each with two orthogonal polarizations (also along the Cartesian axes). Due to the fact that SIE rigorously respects the linear superposition principle, this choice of incidence conditions is immaterial for the final result. Although this choice is entirely arbitrary, it has the esthetic appeal of corresponding exactly to each one of the six Cartesian basis vectors used for the driving fields. It is well known that although the choice of the basis vectors will affect the resulting α-tensor, the different retrieved tensors are related by a unitary rotation matrix consistent with the basis choice. As a final note on the retrieval protocol, we add that the definition of origin that is chosen to refer the dipole moments to, is a non-trivial matter, due to the well-known dependence of electric and magnetic dipoles on the choice of origin (more precisely, the contributions to the electric dipoles created by magnetic currents and to the magnetic dipoles created by electric currents depend on origin). We have made use of the Onsager relations that the polarizability tensor has to fulfill due to reciprocity [21]. Onsager relations in particular state that the upper diagonal and lower diagonal cross-polarizabilities are each other's negative transpose. Accordingly, we choose the origin for both retrieval algorithms as the position for which the sum α E H + α H E is minimum.

Benchmark of vector spherical harmonics and effective current density α-retrieval
In order to benchmark the retrieval of the polarizability tensor, we consider an entirely known object, i.e. a Mie sphere. We focus on a sphere that has both a dielectric and a magnetic response in order to benchmark both the electric and magnetic dipole retrieval. We compare with the rigorous theoretical values for electric and magnetic polarizabilities given by the Mie coefficients (labeled here as c TM 1 and c TE 1 ) [35,36]: For our benchmark, we fix µ = 4 and set equal to the dielectric constant of gold as fitted by Etchegoin et al [37]. While these values do not represent any currently physically realizable object, these values allow us to assess whether we can accurately separate electric and magnetic dipole moments. We use a fixed discretization of the sphere surface by 572 nodes composed of 1241 triangles. While the vertices of the mesh are exactly located on the assumed nanoparticle radius of 10 nm, we note that the triangulated surface is entirely located within the assumed sphere. For this reason, SIE simulations effectively underestimate the sphere size, to a degree that reduces with increasing number of nodes. We quantify the effective radius by calculating the mean distance from the center of the sphere to the surface of the meshing triangles. For the particular meshing conditions used here, the effective radius is 9.96 nm, which we use in the Mie calculations to which we compare the SIE results. In figure 1, we plot the diagonal elements α E x x and α H yy of the polarizability tensor over the wavelength range from 100 to 4000 nm. For both the VSH projection and the equivalent current retrieval, there is excellent correspondence between the retrieved dipole, moments and the polarizabilities predicted by equations (24) and (25). The origin used for the retrieval was found to coincide with the center of the sphere, as expected based on symmetry. The error between the VSH retrieval procedure and the theoretical dipole moments is less than 0.1% throughout the whole wavelength range of the simulation. This agreement is only possible due to the precision of the discrete spherical harmonic transformation over the sphere (taken here to have a radius of 10 µm), and of course also to the extremely good convergence of the SIE method. Furthermore, this almost perfect agreement to the theoretical values spans all the way to wavelengths equal to the diameter of the sphere. This agreement is hence beyond what is needed for metamaterial analysis where wavelengths around five to ten times bigger than the structures are commonly used. When examining the current retrieval procedure, we find that the error between the effective currents retrieval and the theoretical electric dipole moments is less than 0.006% at 4000 nm and grows monotonically up to 1% at 100 nm wavelength. For the magnetic dipole, it is 0.03% at 4000 nm and 8.7% at 100 nm. The difference between the rigorous values and those extracted from the effective currents method is due to the fact that the current-to-dipole expressions used in the current retrieval procedure (equations (16) and (17)) are only valid for kr max 1 as explained by Jackson [25]. These two equations, which are the equations commonly used in the metamaterial field [38][39][40], derive from exact formulae (9.167 and 9.168 in [25]) by replacing the involved spherical Bessel functions with their small-argument asymptotes. Therefore, the effective current method is only expected to be accurate for r λ/2π (i.e. kr ∼ 1). The error is thus not a numerical error but an error due to a poorly met approximation. This error, and whether it is larger for p than for m or vice versa, depends not only on the size of the scatterer but also on the specific weighting given by the current distribution. In contrast, the VSH retrieval through fields is valid for arbitrary frequency and arbitrary size of the radiating object.
The VSH retrieval method can be downloaded from our webpage, to be used with the fields calculated with any full-wave maxwell solver (www.amolf.nl/research/resonantnanophotonics).

Polarizability retrieval applied to Kerker's paradox of zero-forward scattering spheres
As a more challenging benchmark, we also consider magnetoelectric spheres where material parameters are set to the very special condition that is the subject of Kerker's paradox raised in [1] and resolved in [2]. It was first noted by Kerker [1] that at a particular combination of and µ, spheres appear to have zero-forward scattering, yet non-zero extinction. This apparent paradox that occurs for very small spheres when = (4 − µ)/(2µ + 1) gained new interest in the framework of cloaking and invisibility [41,42]. Alù and Engheta [2] showed that these spheres indeed have very low yet non-zero forward scattering, thereby complying with the optical theorem. The almost zero-forward scattering results from destructive interference in the forward direction of the radiation of the generated electric dipole and magnetic dipole moments. Here we reproduce three of the examples studied by Alù [2], using the SIE method (see figure 2), and retrieve the polarizability tensor. First, in figure 2 the bistatic scattering cross section or differential scattering efficiency is plotted for the spheres treated in [2]. The spheres have different radii a = λ/100, λ/20 and λ/4. The permeability of the three spheres is µ = 3 while the permittivity is = 0.143, 0.121 and 0.315. It should be noted that for larger spheres, the condition of minimal forward scattering is shifted away from the criterion = (4 − µ)/(2µ + 1). The calculated efficiencies are in excellent quantitative agreement with the values reported by Alù [2]. It is evident that the forward scattering for the three spheres is close to zero. We report as a table in figure 2(b) the retrieved values of α expressed in units of a 3 for all the three cases, as extracted from the VSH method and expressed in units of the particle radius cubed. All offdiagonal elements are at least 10 5 times smaller than the diagonal elements, i.e. zero within numerical precision. The retrieved polarizabilities are isotropic to within 0.1%. We therefore only report the mean diagonals α E and α H . Evidently, for all three spheres the condition p = −m required for complete destructive interference in the forward direction is almost met, consistent with the conclusion derived in [2] and [1] that this is a necessary condition for zero-forward for a parallel and perpendicular scattered field from three different magnetoelectric spheres of radius 'a', as studied in [1]. The first sphere has =0.143 and µ = 3; the simulation is done at λ = 100a. The second sphere has = 0.121 and µ = 3; the simulation is done at λ = 20a. The third sphere has = 0.315 and µ = 3; the simulation is done at λ = 4a. The table shows the retrieved values of α expressed in units of a 3 for all the three spheres, as extracted from the VSH method and expressed in units of the particle radius cubed.
scattering. For increasing sphere size compared to the wavelength, the imaginary part of the diagonal elements of the tensor increases due to radiation damping. For the largest sphere, p deviates noticeably from −m, and forward scattering is noticeable.
This benchmark shows the usefulness of SIE to simulate magneto-electric scatterers with very high precision, and suggests that the retrieved α-tensor can be used on more complex systems to gain insight into the problem beyond that usually obtained from just brute force calculations.

Split ring polarizability
In the metamaterial community, the performance of a material is usually quantified through effective responses and µ. However, the fundamental parameter underlying the effective and µ is the α-tensor of the metamaterial building block, which is much less frequently studied. Here we use the benchmarked code to understand a single plasmonic magneto-electric structure and use this understanding to design a more complex antenna based on the unit cell characteristics.
Split rings are the archetypical structures used for metamaterials, and extensive literature has been devoted to explain the response of this structure [43][44][45] in terms of LC resonators. Here we report the full polarizability tensor of split rings resonant at 1.5 µm. We consider a gold split ring with dimensions of 30 nm height, 200 nm length and 200 nm width with a central hole of 140 nm by 80 nm. In our simulation, the split ring is placed in a homogeneous environment with = 1 and µ = 1 and the surface discretization used consists of 774 nodes.  The response of the split ring presents two resonant peaks in the wavelength range studied, i.e. from 400 to 1700 nm. The first resonance is centered around 1544 nm with a width of 124 nm and the second resonance is centered at 689 nm with a width of 44 nm, in excellent agreement with experiments and FDTD simulations [46]. The first resonant peak is also called the LC resonance as described in [46]. This resonance has a calculated maximum total scattering cross section of 0.13 µm 3 in very good agreement with measured data, which show an extinction cross section of 0.3 µm 3 with an albedo of 30% [46,47]. We focus on this fundamental resonance and its scattering characteristics, choosing the wavelength of 1544 nm for retrieval of the polarizability tensor. At the split ring resonance, the quadrupolar terms contribute less than 2% to the total extinction. Therefore, here we disregard any higher multipolar terms. The center used for the retrieval was found to be −25 nm from the geometrical center in the ydirection, i.e. closer to the base of the split ring. The SIE dipole polarizability retrieval procedure allows us to quantify both the diagonal values in the polarizability tensor and the crosscoupling between the magnetic and electric moments. We find the following values for the VSH and effective current retrieval procedures, where V denotes the geometrical volume of the split ring (V = 6.33 × 10 −4 µm 3 ).
All the other values in the polarizability tensor are 10 −3 below α E x . Both retrieval procedures indicate that the scatterer has magneto-electric nature given the values of the magnetic and cross polarizabilities. The difference in magnitude between both retrieval procedures is maximally 16%. Based on the MIE calculations, we already saw that the VSH retrieval procedure is accurate, whereas the common current-based definition is fundamentally limited and only valid for small objects of size r < λ/2π. This condition is not met for the split ring. Therefore, we will focus now on the values retrieved with the VSH procedure. The electric polarizability for the p x oriented dipole α E x is the largest polarizability in this structure, and is well in excess of the physical particle volume. Thereby the split ring is very much like a strongly plasmonic particle. That the retrieved α E x is mostly imaginary confirms that λ = 1544 nm corresponds to resonant driving. The magnetic polarizability for the m z oriented dipole α H z is 13 times smaller than the electric polarizability. The off-diagonal values α E H x z and α H E zx significantly exceed the magnetic polarizability. We note that the retrieval very well confirms fundamental constraints on the cross polarizabilities. In particular, the cross polarizabilities are the negative of each other to within 2%, as fundamentally expected from Onsager relations. Also the phase relations arg(α E H x z /α E x ) = π/2 and arg(α H z /α E x ) = 0 are satisfied to within 0.06 rad. These phase relations are consistent with the LC-circuit intuition that if a magnetic response arises through electric driving, i.e. through cross polarizability, it must lag by a quarter wave, as it is due to relaxation of the charge that accumulates in response to E across the capacitor. Finally, we note that α H E zx = i √ α H z α E x to within 0.4%. These data comply excellently with [45], which claims that for planar scatterers that can be described by a circuit model the α tensor must have 'maximally strong' cross coupling. The numerical values that we retrieve are in reasonable accordance with experimentally retrieved values [45], which were reported to be approximately |α E x | = 6.4 V, |α H z | = 0.9 V and |α H E zx | = 2.1 V. That the split ring in our model is comparatively even less magnetic than extracted in experiments is likely due to either one of two causes. Firstly, the split ring resonator (SRR) response depends sensitively on geometrical details such as the exact gap size and the rounding assumed for approximating the SRR shape. The SRR that we model is comparatively thin and rounded compared to the SRRs in experiments. Secondly, in the experiments the polarizability was retrieved rather indirectly, from a comparison of SRR array transmission to a lattice summation model. The fact that SRRs were located at an air-glass interface was disregarded. At a dielectric interface, polarizabilities can be significantly renormalized [48]. A strong response in the α E H cross-coupled polarizability elements implies that the structure possesses optical activity as explained in [45]. Indeed figure 3(b) shows strong optical activity for the split ring, as evidenced by the change of total scattering cross section for right-and left-handed circularly polarized light excitations at different incident angles. Thereby SIE confirms the optical activity expected for any magnetoelectric scatterer with a magneto-electric cross coupling term in its α-tensor [45]. The angle of maximum optical activity θ MAX is consistent with the right eigen-vectors found from the diagonalization of alpha. Indeed, we find that the maximum response of this structure occurs for an eigen-excitation with E = (0.964, 0, 0)E 0 and H = (0, 0, 0.0093 − 0.264i)E 0 /Z , and the minimum response occurs when E = (0.00396 − 0.261i, 0, 0)E 0 and H = (0, 0, 0.965)E 0 /Z (where Z is the characteristic impedance of free space). These excitations require a phase delay between E x and H z that can be generated by using an easily experimentally achievable circularly polarized plane wave under oblique incidence. Furthermore, it is straightforward to calculate that the maximum coupling of a circularly polarized plane wave with a split ring of polarizability as in equation (26) should occur at a polar angle ∼20 • whose full-wave SIE indeed shows as the angle of maximum scattering.

Single split ring as a magnetic dipole converter
One of the most exciting features of plasmonic antennas is that since the plasmons are a combined oscillation of the optical fields and the free electrons in the metal, their resonances can be confined to very small modal volumes [3]. These small modal volumes make plasmonic antennas perfect candidates for coupling to single emitters since local fields and LDOS are enhanced [49]. Some of the functionalities that have been already experimentally proven for single emitters coupled to these antennas are change of polarization of the emitted field by using rod antennas [50] and directionality in the emission of the emitter through the use of Yagi-Uda antennas [51]. We present calculations of the interaction between a single emitter and an SRR.
A technical issue is that the field of a dipolar source driving the scatterer is singular at the position of the emitter. Therefore, unless a very fine discretization is used simulations are prone to large numerical errors. This holds for virtually any brute force method. In the case of SIE, this problem occurs when calculating the values of q (see equation (10)) for fields coming from a dipole in close proximity to the scatterer. However, since the field of the dipole source is given by the Green function of the environment, we can follow a similar procedure to equations (11 and 12), in which the integral over G(r, r 0 ) is separated into a smooth G(r, r 0 ) S and a singular part G(r, r 0 ) NS , leading to the following equation for q: where we used equation (10) and the fact that G(r, r ) T = G(r , r) and (∇ × G(r, r )) T = −∇ × G(r , r) for the free space Green function [30], and the identity ∇ × G(r 0 , r) = ∇G(r 0 , r) × 1 ([26, equation 28]). With µµ 0 we denote the permittivity of the environment where the emitting dipole is positioned. The smooth part can be calculated by a normal quadrature routine.
The singular part can be calculated by using the integration of the RWG function found in [52], i.e. by where K l j (T n ) are the integrals defined in [52] that are performed over the triangle T n or over the surface S n linked to the triangle with the same index. After having found q, we find the strength of the current densities J and M by finding α n and β n , as already explained. It is important to note that this same procedure can be used to find the scattered field at the source and thereby the local density of states (LDOS) [53] when using SIE, avoiding the common problems encountered when working with fields from a dipolar emitter close to scattering structures.
We performed simulations of an electric dipolar emitter located at different distances to the split ring and also at a fixed position in the middle of the split ring for different orientations of the emitter. Figure 4(a) shows the calculated electric and magnetic dipole moments found from just the scattered field of the split ring for different distances of the emitter to the center of the split ring. In other words, we calculate the induced dipole moments in the antenna. The induced electric dipole is given in units of the emitter dipole strength ( p 0 ), while the units of the magnetic dipole are given in terms of p 0 c. With this choice of units we can compare the magnetic and electric dipoles directly, since the magnitude of the radiated power produced by an electric dipole with strength p 0 is the same as that generated by a magnetic dipole with strength p 0 c. Figure 4(a) shows that when the dipole is far from the split ring, the induced dipole is fairly weak. Therefore, the total system emits only with an electric dipolar nature given by the emitter itself. As the emitter gets closer to the split ring to within 230 nm, the total electric dipole moment of the system increases to exceed that of just the emitter, as expected for a high local density of states position near a plasmonic structure. Furthermore, the nature of the lumped system starts to acquire magnetic character to the point that 30% of the emission is of magnetic nature. In the bottom part of the graph, we see that the phase of the driven electric dipole when the emitter is in close proximity to the split ring (50 nm from the geometrical center of the split ring) is delayed by π/2, as expected for a structure driven on resonance. On the other hand the magnetic dipole is in phase with the driving emitter, π/2 advanced with respect to the induced electric dipole. This is expected from [21] since α c = i √ α E α H in the polarizability tensor of a split ring. In a subsequent calculation, we have placed the electric dipolar emitter at the position of maximum radiative LDOS, i.e. 0.05 µm from the center of the split ring, and we varied the orientation of the dipole. From the total scattered plus emitted field of the lumped system, we calculate the effective dipole moments of the complete system. Figure 4(b) shows the electric and magnetic dipole moments as a function of the orientation of the electric emitter, where the emitter orientation is rotated around the z-axis that points through the split ring plane. It is evident that the coupling to the structure only occurs for the x-component of the dipole, i.e. when p x from the emitter couples to α E x of the split ring. The maximum total electric and magnetic dipole moment occurs when the dipolar emitter is aligned with the x-axis. For this alignment, p y < 0.005 p 0 is essentially zero, while p x = 15.1 p 0 and m z = 4.9 p 0 c. This result is commensurate with the relative magnitude of the purely electric and cross-coupled polarizability In (a) we show the calculated dipolar moment of the scattered field for different distances to the split ring center. In (b) we show the dipolar moment of the total field (scattered field plus single emitter field) for different orientation angles of the single emitter in a position 0.05 µm from the center of the split ring. The angles are rotated around the z-axis and therefore the x and y electric dipoles are shown as well as the z magnetic dipole. The other components of the electric as well as the magnetic dipoles are negligible in magnitude. In (c) we show the scattered field pattern of a split ring excited by an electric dipolar emitter at the position on maximum coupling. The electric field magnitude |E| 2 is calculated at a radius of 100 µm from the center of the split ring. The red continuous line shows the field in the plane 'x y' and the blue dashed line shows the field in the plane 'yz'. In (d) we show the calculation of the normalized total and radiative LDOS for different positions in a line along the y-axis through the center of the split ring.
of the split ring in equation (26), indicating that the dipolar scattering approximation of the split ring can be used for dipolar emitter excitations while still obtaining an agreement of 85% with the full-wave calculation. It is important to note that this agreement is dependent on the distance of the emitter to the split ring, since dipolar emitters in close proximity to a plasmonic structure (typically < 20 nm) can increasingly excite higher multipolar moments of the plasmonic structures due to the strong gradients in the exciting fields [54]. When the dipolar emitter is aligned with the y-axis, i.e. rotation angle π /2, then p y = 0.44 p 0 , p x = 0.12 p 0 and m z = 0.04 p 0 c. This result indicates, first, that p y hardly induces a magnetic dipole and, second, that this position has a local density of states for y-oriented dipoles lower than free space. Figure 4(c) shows a polar plot of the far-field intensity distribution of the scattered field for a split ring, excited with a dipole located at the position of maximum coupling and aligned along the gap of the split ring, i.e. along the x-direction. The |E| 2 distribution is evidently different from that of an electric dipole, since on the one hand the emission is asymmetric in the y-axis due to the front-to-back asymmetry of the split ring, and on the other hand the emission in the x-axis is different from zero, evidencing the partial magnetic nature of the scatterer. Finally in figure 4(d), we show the calculated total and radiative LDOS normalized to the vacuum LDOS. The calculations are done for different positions on the y-axis along a line that starts at the center of the split ring. The maximum total and radiative LDOS occurs at a position ∼50 nm away from the center of the split ring. While at this position the total LDOS for an x-oriented dipole is ∼75ρ 0 , the radiative LDOS for an x oriented dipole is ∼252.7ρ 0 . This radiative LDOS is consistent with the generated total electric dipolar moment of 15.1 p 0 and the magnetic moment of 4.9 p 0 c for the lumped system. These values for the dipole moments indicate a radiative LDOS enhancement of ( p 2 + m 2 )/ p 2 0 = 252. The fact that the total LDOS is three times the radiative LDOS is consistent with experimental measurements of the albedo of ∼30% measured for a single Au split ring by Husnik et al [47]. The relative magnitude of the total versus radiative LDOS indicates that the quantum efficiency of the system is η ∼ 33%.

Split ring array antenna
Having understood the split ring as a system composed of an electric and a magnetic coupled dipole moment whose maximal response to circularly polarized plane waves occurs for a certain polar angle θ MAX , and having studied the way in which electric dipoles couple to single split rings, we turn to the design of an array of split rings and to the study of the special properties that arise from it. In our design, we combine two of our earlier results. Firstly, figure 3 shows that the purest handed response is obtained at an off-normal incidence of 20 • . Secondly, excitation of a split ring array with a single molecule is most advantageous when placed 50 nm from the geometrical center and with a dipolar orientation along the gap. In addition, we know from [22,23,51,55] that one can attain directionality in the scattering of arrays of particles by placing them in a linear array with a pitch of ∼λ/3. Our design combines these three ideas in an array of five split rings tilted at θ MAX and excited by a dipolar emitter in the central element. Figure 5(a) presents the scattering pattern of the antenna which clearly shows directivity in its scattering, with scattered fields confined in a half angle < 40 • . By studying the complex fields obtained from the front and back scattering from the antenna, we can retrieve the polarization and plot it on the Poincare sphere, see figure 5(b). We find right-handed elliptically polarized light emanating from the front of the antenna (depicted by the blue point in figures 5(a) and (b)) and left-handed elliptically polarized light emanating from the back of the antenna. Both fields have an electric field seven times stronger in the x-direction than in the y-direction and the major axis of the ellipse is aligned with the x-axis. Thereby, metamaterial antennas allow new forms of control over emission compared to plasmon antennas. We foresee that with split rings with a stronger magnetic polarizability term, it would be possible to reach a totally circular polarized light regime. Reaching stronger magnetic polarizability currently seems easiest at mid-infrared and microwave frequencies [45]. We foresee interesting applications especially if one can reach this at optical frequencies. In this regime one could envision using split ring antennas to generate a single-photon source from a simple linear electric dipole emitter, or from a localized χ (2) nonlinear material that emits its photons in handed beams, or split into two narrow beams, where handedness and direction are entangled.

Conclusions
We have made use of the SIE method to retrieve the polarizability tensor of scatterers. This retrieval is performed in two different ways. The first method consists of a VSH projection of the scattered fields, which yields extremely good precision for any wavelength and size of the scatterer, thanks to the aid of a discrete harmonic transform on the sphere. With the second method, based on effective electric and magnetic surface currents, we can successfully retrieve the polarizabilities of small scatterers with the advantage that this retrieval can be done on non-isolated structures. We have used the lumped system of a dipolar emitter and a split ring to show how the radiation nature of the system changes drastically from a simple electrical dipole emitter. In the lumped system the emission is modified by the scattering of the split ring, which can imprint its magneto-electric nature on the emission. This realization further extends current research efforts that have shown how emission from a single electric dipole transition in a quantum dot can appear as if it originates from a multipole transition by strong coupling of the emitter to a plasmon antenna multipole resonance [56]. Also, such magnetic and magnetoelectric antennas may enhance the magnetic LDOS that magnetic transitions are sensitive to, as recently shown for rare earth ions near an interface [24]. Finally, we used our understanding of split rings to design an array antenna that splits emission from a point source into two beams of oppositely handed elliptical polarization. For ultimately strong magnetic scatterers, these findings might provide new ways of manipulating spins via light, and enhance enantioselective spectroscopies in the near field [57,58].