Electromagnetic multipole theory for optical nanomaterials

Optical properties of natural or designed materials are determined by the electromagnetic multipole moments that light can excite in the constituent particles. In this work we present an approach to calculate the multipole excitations in arbitrary arrays of nanoscatterers in a dielectric host medium. We introduce a simple and illustrative multipole decomposition of the electric currents excited in the scatterers and link this decomposition to the classical multipole expansion of the scattered field. In particular, we find that completely different multipoles can produce identical scattered fields. The presented multipole theory can be used as a basis for the design and characterization of optical nanomaterials.


Introduction
The classical electromagnetic multipole expansion [1] is a powerful tool for analyzing the electric and magnetic fields created by spatially localized electric charges and currents. Irrespective of the complexity of the charge and current distributions, the fields produced by them can be represented as the superposition of the fields created by a corresponding set of point multipoles. This correspondence provides a common basis for characterizing the fields radiated by localized charge and current excitations in arbitrary configurations.
In optics, the multipole expansion is well suited for describing the scattering of optical fields by small objects. Usually, if the wavelength of the field is large compared to the size of the object, the scattering is described mainly by the lowest-order multipole, the electric dipole, while the contributions of all higher-order multipoles are considered as mere perturbations. In recent years, it has been shown that in specifically designed optical nanomaterials [2], such as metamaterials, the contributions of the magnetic dipole [3] and the electric quadrupole [4] excitations to the scattering by the material's constituents can be made significant, which will substantially affect the optical properties of the material and lead to extraordinary phenomena, such as negative refraction [5]. In certain materials, higher-order multipoles can even completely overshadow the electric dipole contribution [6]. Thus, it is clear that higher-order multipoles have to be taken into account when evaluating the macroscopic electromagnetic characteristics of such materials [7].
In order to create a material with the prescribed optical properties, one should select the elementary unit of the material (often called the meta-atom) and optimize its scattering characteristics through adjusting the design. For an individual particle this can be done by numerically solving the Maxwell equations for the scattered field and applying the multipole expansion to determine the multipoles contributing to the scattering [8]. However, in a material composed of a large number of such elementary units, each particle interacts with the fields scattered by the other particles, which can significantly modify the amplitudes and phases of the excited multipoles. Since, at each point in the material, the field is a superposition that contains the fields scattered by all the particles, the approach developed for individual particles [8] can no longer be directly used, and one should perform multipole decomposition on the excited electric currents in each of the particles individually. This decomposition has only previously been introduced for a single localized electric current distribution in vacuum [1] and, therefore, a multipole theory suitable for analysis of nanoscatterers in an array and in an arbitrary homogeneous dielectric host medium has been missing. In this paper, we introduce such a theory. Our multipole expansion approach is particularly suitable for direct numerical implementation.
The geometry of the scatterer determines the types of electric current modes that can be excited in it by light. For the description of these modes, we propose a set of orthogonal electric current multipoles, which are composed of elementary point currents in simple configurations. Each element of the resulting current multipole tensor reflects the strength of one of these configurations, which enables one to visualize the real electric current modes that will be excited in a scatterer. We complete our theory by deriving expressions that relate the elements of the proposed electric current multipole tensor to the classical multipole expansion coefficients. In particular, these expressions reveal (i) perfectly dark multipole modes that do not create any electromagnetic field, and (ii) electric dipole radiation produced by electric currents with zero net electric dipole moment. These findings provide us with additional freedom in the choice of the particle geometry.

3
In section 2, the classical multipole expansion is adjusted to describe the electromagnetic field scattered by individual nanoparticles and nanoparticle arrays embedded in a dielectric host medium. In section 3, we map the coefficients in the multipole expansion to the electromagnetic fields created by sets of sub-wavelength current elements. Explicit mapping relations are derived up to the orders of electric octupole and magnetic quadrupole, which both describe third-order excitations in the multipole hierarchy [9]. In section 4, we summarize our results.

The multipole expansion of the scattered field
We consider a monochromatic plane wave, with the electric field amplitude E 0 , angular frequency ω and wave vector k, incident on a particle in an otherwise homogeneous lossless dielectric medium. In general, the scattered electromagnetic field can be written in spherical coordinates in the form of the following multipole expansion [1]: where X lm and h (1) l are the normalized vector spherical harmonics and the spherical Hankel functions of the first kind, respectively. The wave number k is taken to be that in the surrounding dielectric which has an impedance η. We have chosen the normalization of the multipole expansion such that the expressions for the scattering and extinction cross sections are compact. The vector functions in the multipole expansion form a complete basis for representing the electromagnetic field outside an arbitrary localized source [10].
Let us now assume that the scattered electric field E s is known in the region surrounding the scatterer (e.g. calculated numerically). Using the orthogonality properties of the vector spherical harmonics X lm and the scalar spherical harmonics Y lm , we can calculate the multipole coefficients from the distribution of the scattered electric field on any spherical surface enclosing the scatterer as The above expressions for these coefficients can also be written in terms of the scattered magnetic field H s . In this work, we use the magnetic coefficients expressed as For scattering by a symmetric particle, all the coefficients of higher order in l than some l max can be made equal to zero by properly selecting the origin of the coordinate system. The value of l max depends on both the size and the geometrical complexity of the particle. Next, we assume that the particle belongs to a large array of similar particles. The multipole coefficients must then be deduced from the distributions of the electric current density in the particles. For a single particle, such a derivation is presented in [1] under the assumption that the surrounding medium is vacuum. In order to allow for a dielectric surrounding, we define a quantity which we call the scattering current density. The electric field E = E i + j E s, j contains both the incident field E i (the field in the dielectric in the absence of any scatterers) and the field E s, j scattered by each particle j. In (6), r,d is the real-valued relative electric permittivity of the dielectric and r (r) the complex-valued relative electric permittivity at any coordinate r. We assume that each particle consists of a non-magnetic, isotropic and linear material. The assumption of discrete particles enables us to define a distinct scattering source current J S, j in each particle j, so that J S = j J S, j . Then, starting from the ordinary macroscopic Maxwell equations, we derive the following equations to hold for each particle j: In these equations, the incident field and the fields scattered by adjacent particles are implicitly present through J S, j . According to (7)-(10), the introduced J S, j describes the effective current density that creates the scattered field of the jth particle in the self-consistent solution of the Maxwell equations. Using (7)-(10), we derive the scalar wave equations [1] ( The solutions to (11) and (12) are inserted into (3) and (5) to obtain the multipole coefficients in the form (see [1] for comparison) where j l are the spherical Bessel functions. While the integrations in (13) and (14) are over the whole space, the integrands are equal to zero everywhere outside the particle in question. The spatial derivatives of J S, j make (13) and (14) cumbersome, especially for numerical calculations. We therefore use integration by parts to cast (13) and (14) in the form where l (kr ) = kr j l (kr ) are the Riccati-Bessel functions and l (kr ) and l (kr ) are their first and second derivatives with respect to the argument kr . The associated Legendre polynomials P m l are defined as in [1]. In (15) and (16) we have introduced the following functions and parameters: Equations (15) and (16) yield the same multipole coefficients as (3) and (4) for the light scattered by an isolated particle. However, equations (3) and (4) are not applicable to an array of scatterers. In contrast, since (15) and (16) only require knowledge of the total electric field inside the particles to calculate J S, j , they can be used to characterize the scattering of light by each particle in the array. The required total electric field can be calculated, e.g., by numerically solving the Maxwell equations. Thus, equations (15) and (16) make it possible to characterize the optical properties of nanomaterials, in which multipoles of arbitrarily high order can be excited.
For the scattering of light by a single particle, one can introduce the scattering cross section that describes the efficiency with which the particle removes energy from the incident plane wave into the scattered field. For our multipole coefficients, the scattering cross section can be derived in a similar way as in [11] as the following: The terms of the series in (20) allow one to determine the contribution of each multipole excitation to the overall scattering cross section of the particle. The extinction cross section can also be expressed in terms of the multipole coefficients. For an x-polarized incident wave of the form the extinction cross section is calculated as In contrast to (20), the expression for the extinction cross section depends on the choice of the polarization and propagation direction of the incident field. For example, if the incident wave is y-polarized, the extinction cross section is expressed as The extinction and scattering cross sections can also be calculated by using the optical theorem [1] and calculating the total power of the scattered light. The results can then be compared with those obtained by using (20), (22) and (23), which can serve as an additional check for computations.

The multipole expansion of the electric current density
The scattering of light by a particle can effectively be seen as a process where the incident light excites in the particle polarization and conduction currents that radiate. These currents can be decomposed into terms which we call current multipoles. To accomplish this decomposition, we choose to use the concept of point electric current elements introduced in [12]. A particle interacting with light can then be treated as a collection of such point elements. Before presenting the multipole expansion of the electric current density, we first consider a wire of length L that carries a time-harmonic electric current with a complex amplitude I. The wire is positioned at the origin of the coordinate system inside a homogeneous and isotropic dielectric. The oscillating current emits electromagnetic radiation of wavelength λ into the dielectric. By assuming that L λ, we can treat this current-carrying wire as a point current element with the complex amplitude of the current density given by This current element fully corresponds to an oscillating point electric dipole. Next we consider two point elements in which the currents oscillate in opposite directions. The element with a current +I is displaced from the origin in the positivex direction by a distance s/2, while the other element, with a current −I, is displaced in the opposite direction by the same amount. If s λ, we can treat this elementary current configuration as a secondorder current element. Considering s to be infinitesimally small, but such that the product I Ls stays finite, we obtain for the complex amplitude of the current density describing this secondorder element where the operator is defined as Similar displacements can be done in theŷ andẑ directions by applying the operators ς(ŷ) and ς(ẑ), respectively. The point current elements of the third and higher orders can be obtained by sequentially applying the operator in (26) to the current density of the lowest order point element of (24). In (26),û can be chosen asx,ŷ,ẑ or as any linear combination of them. The effect of the operator on a point current element of a certain order can be seen as follows. The operator makes a copy of the element and shifts the phase of the complex amplitude of this copy by π radians. The original element is then displaced in theû direction by a distance s/2 and the copy in the −û direction by the same amount (see figure 1). Finally, s is set to be infinitesimally small.
We describe the amplitudes of the point current elements by the following lth-order current multipole moments: The multipole expansion of the current density, obtained by repeatedly applying the operator to (24), can be written in the Cartesian coordinates as where M (l) (v, a, b), withv being equal tox,ŷ orẑ, are the elements of M (l) describing a multipole obtained by applying the operator ς(û) to av-oriented current element a times witĥ u =x, b times withû =ŷ and l − (a + b + 1) times withû =ẑ. The coefficients in (28) are chosen such that the elements of the multipole tensors are consistent with (27).
We wish to map the elements of M (l) onto the coefficients a E (l, m) and a M (l, m) in the multipole expansion of (1) and (2). For this we need to solve for the electromagnetic fields created by the current density distribution of (28). We start by defining the vector potential A through 8 In the Lorenz gauge, A satisfies the wave equation For a point electric current element, such as the one in (24), the solution for the wave equation where p = iI L/ω and = 0 r is the electric permittivity of the surrounding dielectric. The multipole expansion of the vector potential is obtained by applying the operator to (30) with J = J 1 and A = A 1 , until J on the right-hand side becomes that in (28). Since the spatial differential operators in the Cartesian coordinates commute with ∇ 2 , we obtain the multipole expansion for A from (30) as We found that the calculations can be significantly simplified by making use of a circular coordinate system (w, w * , z), where w = (x + iy)/ √ 2 and w * is its complex conjugate. The unit vectors in this system areŵ,ŵ * andẑ, whereŵ = (x − iŷ)/ √ 2. These circular coordinates are linearly related to the Cartesian ones and, therefore, the operator in (26) still commutes with itself. Thus, the vector potential can be expanded as where M (l) (v, a, b) are the elements of M (l) in the circular coordinate system. The magnetic field corresponding to (33) is given by (29). Outside the scatterer, the electric field is calculated as where (10) and (30) have been used. The multipole coefficients a E (l, m) and a M (l, m) are obtained by inserting ther components of E and H into (3) and (5). For each multipole, ther components of E and H can be written in terms of scalar spherical harmonics. The orthogonality of the spherical harmonics can then be used to solve for a E (l, m) and a M (l, m), without the need to manually calculate the integrals in (3) and (5). The elements of the current multipole tensors obtained in the circular coordinate system can be related to the tensor elements in the Cartesian coordinate system by using the tensor projections. In complex coordinates, the projection of the tensor M (l) onto the unit vectorŵ is M (l) ·ŵ * . After lengthy calculations, the following mapping relations are obtained (up to the electric octupole coefficients): where We have further verified the correctness of (35)-(48) numerically. We used the computer software COMSOL Multiphysics to numerically calculate the electromagnetic field created by electromagnetic point sources in subwavelength dipole, quadrupole and octupole configurations corresponding to each tensor element in (35)-(48). Using the obtained electric fields in (3) and (4), we evaluated the multipole coefficients for each tensor element in p, Q and O. The obtained coefficients are in full agreement with (35)-(48).
The presented theory can be used as follows. For an arbitrary nanoscatterer, including a nanoscatterer in an array, one first numerically evaluates the fields inside the scatterer and uses (6), (15) and (16) to obtain the multipole coefficients a E (l, m) and a M (l, m). Then (35)-(48) are used to find the essential current excitations in the scatterer. Each tensor element corresponds to a certain electric current mode in the scatterer. For example, the O x yz mode describes the current configuration obtained by operating with both ς(ŷ) and ς(ẑ) on an x-oriented current element. Following this recipe, we have, e.g., designed and characterized nanoscatterers in which incident light does not excite any electric dipole moment [6].
Let us consider some important properties of the derived current multipole tensors, starting with the quadrupole one. It can be seen that only eight coefficients [a E (2, ±2), a E (2, ±1), a E (2, 0), a M (1, ±1) and a M (1, 0)] are used in the multipole expansion, whereas there are nine elements in the Cartesian current quadrupole dyadic Q. This is explained by the fact that the spherically symmetric excitation with Q x x = Q yy = Q zz does not generate any electromagnetic field. This perfectly dark excitation corresponds to a radially oscillating positively charged spherical shell with an equal negative charge at the center. Thus, from the knowledge of the radiated electromagnetic field, the moments Q x x , Q yy and Q zz cannot be determined uniquely. However, the real currents should match the geometry of the scatterer, which allows one to make a unique, physically justified choice for the values of Q x x , Q yy and Q zz . This choice enables one to uniquely specify the excitation character. For the octupoles, there are 15 multipole expansion coefficients (7 electric with l = 3, 3 electric with l = 1 and 5 magnetic with l = 2), but 18 different elements in the Cartesian current octupole tensor. Similarly to spherically symmetric quadrupoles, the three symmetric octupole excitations 2O x x x = O yyx = O zzx , O x x y = 2O yyy = O zzy and O x x z = O yyz = 2O zzz are perfectly dark. One of the octupoles in each of these excitations can therefore be chosen arbitrarily.
From (35)-(48) we note that both current octupoles and dipoles contribute to the same three dipole coefficients [a E (1, ±1) and a E (1, 0)] in the multipole expansion. As a consequence, an octupole current distribution with zero dipole moment can create an electromagnetic field indistinguishable from the field of an electric dipole. For example, let us consider a current octupole with O = O(xxẑ +ŷŷẑ − 2ẑxx − 2ẑŷŷ), Q = 0 and p = 0. According to (35)-(48), this current distribution creates exactly the same electromagnetic field as a dipole with a momentp = 2Ok 2ẑ . This current octupole is illustrated in figure 2. The equivalence between the radiation patterns of the current octupole and dipole can also be seen by writing the vector potential in (32) for the octupole as where we have used the fact that ∇ 2 h (1) 0 (kr ) = −k 2 h (1) 0 (kr ). The second term is precisely the vector potential of a dipole with a momentp = 2Ok 2ẑ . The first term is a gradient of a scalar function that does not contribute to the radiated field (see (29)). This finding emphasizes the importance of the relations in (35)-(48) when using multipole coefficients to describe electromagnetic excitations. The fact that spatially orthogonal current excitations can create the same electromagnetic fields provides additional freedom to the choice of the particle geometry when designing functional optical nanomaterials.

Conclusions
In summary, we have introduced a theoretical approach to calculating the electromagnetic multipole excitations in a material consisting of localized nanostructures in a dielectric host medium. Propagation of light through an array of such nanostructures can be studied numerically by solving the Maxwell equations with appropriate boundary conditions. Then, using our theory, the multipole excitations in the structures can be revealed.
In order to obtain an intuitive picture of the real electric current excitations in the scatterers, we introduced a set of easily visualizable electric current multipoles such that any excitation can be represented by their linear superposition. The multipole expansion coefficients can be calculated numerically and utilized to find the real electric current modes in the nanostructure. The same equations can also be exploited in the reverse order to tailor the angular distribution and directionality of the scattered radiation. The theory presented in this work provides the reader with an exact recipe and all the necessary equations for the design and characterization of nanoscatterers and the optical materials composed of them.