1. INTRODUCTION TO OPTICAL METAMATERIALS

It is now commonly accepted that the era of metamaterials (MMs) was ushered in with the publication of the widely known paper by Veselago [1], who first suggested that the basic principles of electrodynamics do not forbid the possibility of materials with negative values of the real parts of both the permittivity and permeability. One consequence of this suggestion is the existence of the so-called left-handed materials (i.e., the materials with simultaneously negative real parts of their permittivity and permeability), where the phase and group velocities are (in the simplest cases) opposite to each other. In fact, the idea of opposite directions of the phase and group velocities was first mentioned much earlier—credit has to be given to [2,3] back in 1904. In 1944, left-handed optical materials were mentioned in a series of lectures by Mandelshtam at Moscow State University [4], which unfortunately have not been published (as has happened many times before, for instance with Leibnitz and Newton, where the discoverer did not publish the respective results in time). It took another quarter of a century before Veselago formulated the fundamentals of electrodynamics with negative permittivity and permeability. At that time, a microscopic mechanism of achieving negative values for the permittivity and permeability was not even discussed, but one seminal work, predicting the crucial role of resonances in small particles and its influence on the anomalous values of the effective constants, had already been published [5]. The modern era of MMs began with the experimental verification of a negative refractive index, first in microwave [6] and later [7,8] in the optical domains. From this time onward, an explosive amount of publications appeared and continues to appear in scientific and popular papers regarding the fundamentals and applications of the MMs. However, it should be remembered that the physics of the MMs (as a branch of science) is still far from mature and there is much left to uncover. An optical MM is based on the nanophotonics [in order to provide meta-atoms (MAs), unit cells from which MMs are constructed; in the optical domain the MAs have to have nanosizes] and evolves in conjunction with the respective technological, experimental, and theoretical achievements. The development of this branch of modern science has already passed its embryonic stage and now seems to have crossed the invariable dip “nanochasm” in interest after the initial excessive expectation usually associated with all new ideas.

MMs are artificial media which tailor the macroscopic properties of light propagation by a careful choice of a MA from which they are constructed. By controlling the geometrical shape and material dispersion of the MA, novel effects such as negative refraction [9–11], optical cloaking [12–17], as well as a series of optical analogs to well-known physical phenomena from various disciplines within physics can be observed [18–22]. In addition to a bi-axial anisotropic (linear dichroism) material response [9–11,23–25], research was extended toward the exploration of MAs that affect the off-diagonal elements of the material tensors (elliptical dichroism), leading to, e.g., optical activity [26–30], bidirectional and asymmetric transmission [31,32], or chirality-induced negative refraction [33–35]. MMs can amplify evanescent waves and thus increase resolution of an optical device, forming a so-called perfect lens [36].

MM design in the optical domain is mainly carried out using rigorous Maxwell’s equation solvers like finite-difference time-domain simulations [37], finite-element methods [38], and Fourier modal methods (FMMs) [39]. Instead of these differential methods, integral ones, such as, e.g., the boundary-element method [40], offer an alternative choice of techniques. The discrete-dipole approximation [41] and the multipole method [42] are more physical ones, where the structure is represented by localized electric multipoles. Nevertheless, presently differential descriptions dominate in MM design.

In contrast to such numerical techniques, the analytical description of MMs is much less developed, and only a few fundamental papers have been published. Podolskiy et al. introduced the coupled-dipole equations, in order to approximate single and coupled metal wires [43,44]. The direct excitation of LC resonances (electrical circuit type resonances) with the magnetic field of the incident plane wave in a system of two coupled rods was proposed in [45] in order to explain the observed phenomena in terms of effective parameters. Following the “effective medium” theory, an investigation of dielectric and magnetic conducting inclusions was performed for spheroids [46]. In spite of the fact that the approach is limited to the realm of quasi-statics, the model was extended to describe dynamical problems. In order to simulate the current distribution in a coupled-wire structure and to calculate the permeability, the Green’s function technique was applied [47]. The application of LC circuit theory has been used to obtain the resonance frequencies and the quality factors of coupled split-ring resonators [48].

Several reviews regarding the state of the art in the physics of MMs have been published [49–53]. A new MA paradigm—hybridization with functional agents (like quantum dots, dye molecules, or biomolecules)—becomes one of the main points of the future development in this area. Taking into account that the functional agents are also supposed to be quantum systems, one can easily conclude the importance of bringing quantum concepts to bear on the physics of MMs, giving rise to the appearance of the new area of the quantum MM. The name quantum MM means that the internal dynamics of the MAs is described using quantum theory tools at least in part (for example, in the case of a spaser [54]) or in full (in the case of MM with superconducting MAs [55]).

The theoretical description of the MM was started from a priori introduced dispersive permittivity and permeability. The spectral dependence of both permittivity and permeability was supposed to follow analytical functions with a resonant denominator, which assumed a harmonic oscillator-type model for both effective parameters [1]. This model basically assumes local media response from small MAs, possessing not only the dielectric (like any dipole-like atom) but the magnetic response as well. The assumed smallness of the hypothetical MAs allows the commonly accepted homogenization theory to be invoked. The situation changed shortly after the first experimental realization of the MM, where the typical sizes of the MAs turned out to be not much smaller than the respective wavelengths and hence the usual homogenization model fails. In order to develop a homogenization procedure for MMs in general, more sophisticated models are required. Basically, the homogenization models for MMs appeared to be an extension of several approaches. One of them was developed in the optics of crystals [56] with the extension to an allowed magnetic response [57,58]; this extension follows a basic approach demonstrated by “A Course of Theoretical Physics” [59]. An extension of the exciton theory to MMs [60] has to be assigned to the same category. The second approach uses formalisms developed for photonic crystals with an appropriate extension over Bloch waves and the determination of the respective dispersion diagram in the form of Brillouin zones (see [61] and references therein) for a particular structure. The third approach appears to be a continuation of the methods developed in the theory of compound materials [62–64] and extended to the case of a magnetic response.

In contrast, a fourth approach, based on the approximation of the MAs as point-like multipoles, has not received much attention, but had been mentioned as one possible way of describing the behavior of MMs in Ref. [58]. This approach is not ab initio, and requires fitting with numerical or experimental data in order to fix the phenomenologically introduced parameters. Nevertheless, this approach allows us to develop a universal model for systematic consideration of practically all optical phenomena in MMs.

The necessity of taking into consideration spatial dispersion is true for all approaches. The problem appears to be in a gap between developed approaches and the basics of the homogenization, which have to be satisfied anyway. Ignoring of the basics of the homogenization can lead in some cases to the violation of the self-evident basic assumptions like causality and passivity (the review of these found in publications violations is given in Ref. [62]).

In some cases the different approaches like phenomenological [59] and multipole [65] appeared to be mixed [57] and different representations of the macroscopic Maxwell equations (the Landau–Lifshitz “L&L” representation [57] and Casimir representation [64]) sometimes appeared to be not clearly distinguished (actually, discussion about the different representations of the macroscopic Maxwell equations can rarely be found in publications at all). The connections between the different representations, the appearance of the spatial dispersion in different representations, and the mutual transformations between the different representations have not received enough attention in publications, to date.

The multipole approach developed originally in Ref. [65] assumes an expansion of the charge dynamics in atom/molecules, which in the case of the MM has to be replaced by the charge dynamics in the MA (note that in some publications under multipole expansion, an expansion of the fields over the wave vector is assumed, which is obviously a different math tool entirely; mutual relations between these expansions can be found in Ref. [66]). The original approach [65], where averaged parameters are elaborated based on the charge dynamics, gives us a unique avenue for the creation of a unified approach to all possible types of MMs: charge dynamics in the MA can be described in the frame of classic or quantum models, but the algorithm for the effective parameter calculation remains the same. It motivates us, in turn, to try to create a unified approach, which would unify classic, quantum, or semi-classic MMs, including toroidal/anapole structures [67,68].

It is worth considering the educational aspects of homogenization in connection with the creation of the model mentioned above. In modern courses of electrodynamics the homogenization (averaging) procedure sometimes does not receive enough attention, which can partially be explained by the fact that the commonly accepted approach (local frequency dispersive permittivity and permeability) gives in most cases pretty good correspondence with the experimental data (not for MM) and the necessity of the more sophisticated approaches is pretty low. The electromagnetic theory was developed a rather long time ago, but some basics of this theory (especially elaboration of macroscopic Maxwell equations—basics of the homogenization procedure) have not been revisited since the 1970s. Due to the appearance of the MM, the averaging procedure must now be considered more fully and has to become a part of the standard courses of electrodynamics in order to teach students in a self-consistent and unified manner.

The last, but not least, comment is about the evident inclination in favor of numerical methods in order to describe the optical properties of MMs at the expense of physical intuition. There is no doubt that modern numerical algorithms and available computer facilities provide the main way to investigate more or less complicated problems. Nevertheless, the analytical approximate type models can provide a deeper understanding of the basic physical processes, stimulate discussion of new effects, and even provide a new paradigm for optimization of a particular design. The analytical models are complementary to the numerical ones, taking advantage of careful comparison with the results of rigorous numerical calculations, but at the same time remaining analytically treatable. In order to create this type of model, accurate approximations have to be made in order to simplify the respective consideration, at the same time keeping the main physical effects and interplay between them in the model. The analytical models allow us to show mutual self-consistency of the physical theory (in this particular case in application to the MM), which otherwise could be seen as a huge leap from independent experimental to theoretical facts.

One of the possible analytical models has been developed and presented in a set of publications; see [69]. In the presented review, analysis of different optical properties of MMs is given using the developed model as a universal tool. It is believed that it could help to create a more systematic approach for analysis of properties and potential applications of MMs.

2. HOMOGENIZATION OF MAXWELL EQUATIONS—PHENOMENOLOGICAL AND MULTIPOLE APPROACHES

A. Microscopic Maxwell Equations and Averaging Procedure

We consider as a starting point a system of microscopic Maxwell equations (MEs) in the following form:

Here $\overrightarrow{e}$ and $\overrightarrow{b}$ are the microscopic electric and magnetic fields, respectively; $\rho$ is the charge density; $q_i$, $\overrightarrow{p_i}$, $\overrightarrow{r_i}$, and $\overrightarrow{v_i}$ are the charges, momentum, coordinates, and velocities of charges; $\overrightarrow{j}$ is the microscopic current density; $\omega$ and $c$ are the frequency and the velocity of light in vacuum; and $[\overrightarrow{v_i}*\overrightarrow{h}]$ is the cross product. It is assumed that system (2.1) is strictly valid without any approximations [59]. Nevertheless, it has to be mentioned that the basic formulation of electrodynamics is still under discussion [70].

We consider propagation of an electromagnetic plane wave interacting with the medium in the case when the classical dynamics is supposed to be valid and the bulk material fills the whole space. System (2.1) in this case can be formally averaged over a physically small volume (or through statistic averaging), which results in the following:

The main problem here is to find the averaged current and charge distribution as functions of the averaged electric and magnetic fields:

Roughly speaking, the typical size of the averaging volume has to be in any case 5–10 times smaller than the wavelength; for MM it results in no more than in the best case three particles per averaged volume. Hence, the averaging concept for this case is required to be qualitatively different in comparison with the classical (for example, Lorenz–Lorentz) one. More details about applicability and limitations of the homogenization procedure can be found in Ref. [71].

There are two main approaches to the averaging, namely the statistical one [65,72,73], where averaging is taking over the ensemble of realization, and spatial, where the main questions are volume of the averaging and the averaging function, which set the minimal macroscopic scale at which the change of the macroscopic (averaged) functions is still significant [74].

A significant advantage of the statistical approach is the absence of characteristic scales, like the volume of the unit cell, because the averaging is performed over all possible realizations rather than over physical volume. A drawback of the statistical approach is the relatively complex math required for the elaboration of the model. The concept of the statistical averaging, being well developed for the electrostatics and magnetostatics, appears not to be completed for the microwave frequencies and optical spectra [75].

One more approach, called the scaling algorithm, has to be mentioned [76]. The method is based on specially introduced multipoles (different from the usually defined ones), which are assigned to the cells at each step of the averaging. The method (similar to [72]) is a universal one for diluted and dense composites, but is not constructive. A more comprehensive general review of different approaches to the problem of averaging can be found in Ref. [77].

B. System Under Consideration

It is methodologically appropriate to determine from the beginning the type of MMs which will be considered and keep this type in mind throughout the text. Here a medium consisting of artificial MAs embedded in a dielectric matrix will be considered (see Fig. 1, where only one layer of the considered material is presented).

 

Fig. 1. Artificial MAs (plasmonic nanoresonators) embedded in a dielectric matrix form a MM (only one layer is presented). Polarization of the electric and magnetic fields, and direction of the wave vector are shown.

Download Full Size | PDF

The MAs are assumed to be complex plasmonic structures, possessing so-called symmetric and antisymmetric eigenmodes (actually, in general for nonsymmetric structures the modes are called quasi-symmetric and quasi-antisymmetric)—see Fig. 2, where one possible structure (coupled nanowires, in general of different sizes) is shown.

 

Fig. 2. One of the possible shapes of MAs, possessing (a) symmetric and (b) antisymmetric modes. Electric field ${\overrightarrow{E}}_x$ of the incoming wave, propagating along the $y$-axis excites eigenmodes of the plasmonic MA.

Download Full Size | PDF

The structure consists of two nanowires with typical lengths (for the optical domain) of tens to hundreds of nanometers, placed one under another with the distance of several tens of nanometers, ensuring strong near-field interaction between both nanowires [78]. The structure possesses two fundamental eigenmodes, namely, (quasi)symmetric [Fig. 2(a)], which produces effective dielectric response, and (quasi)antisymmetric [Fig. 2(b)], which is responsible for the magnetic dipole and quadrupole moments and consequently for the magnetic response of the media.

The incoming electromagnetic wave interacts with the electrons of the plasmonic structure and effectively excites symmetric and asymmetric oscillations, provided that the frequency of the incoming wave is close to the respective eigenfrequency of the eigenmodes. In the case of symmetric oscillation [Fig. 2(a)], the MA formed by two wires of equal length does not produce any magnetic effects, but rather exhibits extra dipole moments leading finally to the change of the permittivity of the MM. The main interest to the MM is stipulated by the possibility to excite the antisymmetric modes [Fig. 2(b)]. In this case the structure presents to the first approximation a circle current, which provides (as it is well known from the school course of physics) a magnetic response. The fantastic peculiarity of the MMs is in the fact that the MMs provide magnetic response at optical frequencies, where no natural media has similar properties. This fundamentally distinguishes MMs from any natural materials, and makes such MMs (among others) extremely interesting objects for both fundamental research and various applications [49].

In spite of the wide range of different shapes of the nanoresonators considered in the publications (including extremely exotic ones [68]), the coupled nanowire was the first structure where a negative refractive index has been demonstrated [7,8]. Moreover, properties of this structure allow clear physical interpretation and relatively simple analytical treatment, which makes coupled nanowires a very good object for discussion of different physical models. One of the simple and at the same time rather good examples is in consideration of necessary conditions for the existence of antisymmetric modes. The antisymmetric mode can be excited due to the following reasons:

Both cases lead to rather different optical properties of the MM consisting of such MAs.

C. Different Representations of Material Equations

There is only one model (a multipole model [65]) where the averaging procedure for the $⟨\overrightarrow{j}⟩,⟨\rho ⟩$ as functions of microscopic dynamics is performed rigorously; all other models do not even try to make this step. The last equation in Eq. (2.1), $\frac{d\overrightarrow{p_i}}{dt}=q_i\overrightarrow{e}+\frac{q_i}{c}[\overrightarrow{v_i}*\overrightarrow{h}]$, is usually left out completely during the averaging, and the necessary equivalent information about charge dynamics is brought to the model phenomenologically.

The system of Eqs. (2.2) and (2.3) is rather useless in practice until we find analytical expressions for Eq. (2.3). Nevertheless, even without finding of an analytical form for Eq. (2.3), the averaged MEs can be analyzed and important conclusions can be made.

It is worth noticing that if we assume some analytical form for Eq. (2.3) (see, for example, [79] and references herein), then the averaging problem is basically fixed (or, better to say, bypassed). System (2.2) becomes self-consistent and can be solved for the electric and magnetic fields $\overrightarrow{E},\overrightarrow{B}$. Any further considerations (including introduction of $\overrightarrow{D},\overrightarrow{H}$ in different representations, as well as the permittivity and permeability) in this case are no more required. Thus, in what follows we assume that there are no explicit forms of Eq. (2.3) and it is necessary to elaborate Eq. (2.3) further in order to find some reasonable analytical expressions for the averaged charge and current densities. It has been shown [71] that the constitutive equations (2.3) can be written through one new function $\overrightarrow{P}$ with the accuracy of two more arbitrary functions ${\overrightarrow{F}}_1$ and ${\overrightarrow{F}}_2$:

The question about the physical meaning of the functions ${\overrightarrow{F}}_1$ and ${\overrightarrow{F}}_2$ in Eq. (2.4) remains opened. Both functions ${\overrightarrow{F}}_1$ and ${\overrightarrow{F}}_2$ are arbitrary and independent. There are different but countable numbers of choices for the possible representations of Eq. (2.4). The most general case is when both ${\overrightarrow{F}}_1$ and ${\overrightarrow{F}}_2$ are nonzero functions, namely,

In this case MEs include four functions $\overrightarrow{E},\overrightarrow{B},\overrightarrow{D},\overrightarrow{H}$:

Note that the case ${\overrightarrow{F}}_1=0$ and ${\overrightarrow{F}}_2=c*{\overrightarrow{M}}_C$ leads to the same form Eq. (2.6), where the rot part of the full polarizability is excluded (the physical meaning of this part—the presence of toroidal moment [67]—will be considered below).

Alternatively to Eq. (2.5),

In this case MEs contain three functions $\overrightarrow{E},\overrightarrow{B},\overrightarrow{D}$:

Finally, we assume that the full polarizability contains only the rot part, namely:

In this case the system of MEs contains three functions $\overrightarrow{E},\overrightarrow{B},\overrightarrow{H}$ and reads

This set of equations can be used only in very special cases where the averaged charge density is zero. Thus, in general, the toroidal form cannot be used instead of the Casimir or Landau–Lifshitz forms. The physical object corresponding to such representation is a toroid [64,67], which is now of great interest in connection with the potential possibility of design of such structures at nanoscales for optical wavelength region application [81]. It is seen that the presence of toroidal moments is responsible for the function ${\overrightarrow{F}}_1$, and fixes the functions $⟨\rho ⟩,⟨\overrightarrow{j}⟩$ in general representation Eq. (2.4).

It is important to realize that there are no other choices for the constitutive equations. Any homogenization model has to start from the statement in which representation it will be developed; arbitrary mixing between several representations is not acceptable.

D. Serdyukov–Fedorov Transformations Between Different Representations

The Serdyukov–Fedorov transformations (SFTs) are relations between two sets of four vectors $\overrightarrow{E},\overrightarrow{B},\overrightarrow{D},\overrightarrow{H}$ and ${\overrightarrow{E}}^{\prime },{\overrightarrow{B}}^{\prime },{\overrightarrow{D}}^{\prime },{\overrightarrow{H}}^{\prime }$, where both sets satisfy Maxwell’s equations. The SFTs are usually written in the following form (for “C” representation):

“C” to “L&L” transformation: Starting from “C” representation, one can unambiguously reduce the MEs to the “L&L” form. It is important to emphasize that in general ($\overrightarrow{T_1}\ne 0$) both electric and magnetic fields are transformed and lose their initial physical meanings. The requirement of keeping the electric and magnetic fields the same in both representations is an additional one with respect to the SFTs.

“L&L” to “C” transformation: The reverse “C” to “L&L” transformation is in general undetermined. There are unlimited numbers of “C” forms which correspond to the same “L&L” form.

“C” to “T” transformation: The “C” to “T” transformation in general cannot be performed.

“T” to “C” transformation: The “T” to “C” transformation cannot be unambiguously determined.

“L&L” to “T” transformation: The transformation, similar to the case “C” to “T,” in general is not determined.

“T” to “L&L” transformation: The “T” to “L&L” transformation cannot be unambiguously determined.

3. MULTIPOLE EXPANSION

A. Multipole Approach

The multipole model was put forward in Ref. [65], and later developed in a similar form in Ref. [72]. The model is based on an averaging procedure using the probability distribution function (PDF) for the positions and velocities of all charges, included in the consideration—statistical averaging, which is supposed to be equivalent to the originally assumed averaging over volume. The model results in constructive expressions for ${\overrightarrow{P}}_C$ and ${\overrightarrow{M}}_C$ presented through the averaged dynamics of the charges in “C” form [82]:

The definitions clearly distinguish between microscopic ($\overrightarrow{r}$) and macroscopic ($\overrightarrow{R}$) coordinates [83], $q_k$ represents the charge, and $\eta$ is their density. The microscopic coordinates $\overrightarrow{r}$ designate the position vectors of the charges in a microscopic coordinate system, and $\overrightarrow{v}$ designate their velocities. $\overrightarrow{P}(\overrightarrow{R},\omega ),Q_{ij}(\overrightarrow{R},\omega ),\overrightarrow{M}(\overrightarrow{R},t)$ are the averaged polarizability, quadrupole, and magnetic moments, respectively, ${(\nabla •Q(\overrightarrow{R},\omega ))}_i\equiv \frac{\partial Q_{xi}}{\partial x}+\frac{\partial Q_{yi}}{\partial y}+\frac{\partial Q_{zi}}{\partial z}$. The center of the microscopic coordinate system is chosen to be the center of symmetry of the charge distribution (consideration of the dependence on the origin of the coordinate system will be given later). It is necessary to take into account both the electric quadrupole and the magnetic dipole terms, because they are of the same order in the multipole expansion series [74,82].

The multipole approach allows us to create a logical connection from the microscopic to macroscopic forms of the MEs without any methodological gaps. The fact that finally this approach results in “C” form serves as one more positive argument for the use of this model and its application to the problem of homogenization of MMs.

It has to be accepted that the basic conditions under which system (3.1) has been elaborated are met for typical MMs in the optical domain rather poorly, and the basic question about applicability of the multipole model to MMs remains open. In spite of the fundamental doubts about its applicability, one can easily bring several arguments in favor of the multipole model:

The multipole model created for MMs [69] contains parameters that can be tuned in order to compensate for the fundamentally stipulated discrepancies and finally fit the results of the model to the experimental and/or numerical data. It is believed that the combination of the multipole approach with final tuning of these coefficients makes this model an extremely simple and versatile tool for investigation of optical properties of MMs [84]. In Ref. [85] analytical expressions for the effective permittivity and permeability have been elaborated for MMs based on double-wire structures. The charge dynamics has been treated using two coupled harmonic oscillator equations, possessing symmetric and antisymmetric oscillation modes, excited by the electric field. Using the expressions for the symmetric and antisymmetric modes, the dipole, quadrupole, and magnetic dipole terms in Eq. (3.1) have been calculated as functions of the MAs and field parameters. Note that in [85] only symmetric structures have been considered; extension to the case of asymmetric structures was performed in Ref. [86].

Finally, the multipole model can be extended on higher-orders of expansion. Nevertheless, it is believed that the accepted here truncation by the second order allows us to combine the ability to introduce magnetic response of the metamaterials with analytically treatable description.

B. Dispersion Relation Elaboration

Assuming linear dependence of the multipole terms on the electric field (dependence on the magnetic field is negligible), one can write ([71], see Fig. 2)

The dispersion relation in the general form is

The constitutive equations are

Dispersion relation (3.3) basically solves the problem of propagation of plane waves in media with higher (up to second-order) multipoles.

C. Origin Dependence of the Multipole Moments

Analysis performed in Ref. [71] allows us to finally conclude that

These conclusions evidently contradict those accepted in the literature that the multipole model can be constructed in such a way that the requirement of the origin independency can be satisfied [82]. Note that in the approach developed in Ref. [82] initial expressions for the multipoles differ from one elaborated in Refs. [65] and [71]. The main difference is in following: in Ref. [82] the multipole expansion (i.e., expansion of the charge dynamics) is mixed with the expansion of the local field; in Ref. [71] the field expansion and consequent consideration of the local/averaged fields and the consideration of their mutual relations is not necessary for the conclusions presented above. The relation of the model developed in Ref. [71] and the one in Ref. [82] requires further investigations.

D. Impossibility of Unambiguous Effective Parameters Introducing for Bulk Materials

The scheme introduced above allows us to determine unambiguously the terms “effective parameter” and some other notations, which appeared in the literature in the contents of the homogenization procedure. The nonlocality (spatial dispersion) caused a discussion in the literature about whether it makes sense to call ${ϵ}_{\alpha \beta }(\overrightarrow{k},\omega )$ an effective material parameter. Following [79], the effective permittivity ${ϵ}_{\alpha \beta }(\overrightarrow{k},\omega )$ is an effective parameter, but is not an effective material parameter because it includes information not only about material, but also depends on the wave vector; in other words, it contains information about external electromagnetic fields. Similar considerations have been published in Ref. [62] where also two different notations—characteristic material parameters and effective material parameters—have been suggested.

In order to create a consistent terminological basis, it has to be mentioned that the introduced effective parameters in the case of strong spatial dispersion (for example, the permittivity in case of the “L&L” representation) solve the problem of homogenization, and to this extent do not require any other comments or discussions. If the functional form of ${ϵ}_{\mathrm{LL},\alpha \beta }(\overrightarrow{k},\omega )$ is found, then this function contains information about material properties (for example, eigenmodes of the electron oscillations in the metal nanoresonators) and the excitation conditions (wave-vector dependence). In general, both properties—eigenmodes and excitation conditions—are not separable. It has to be clearly realized that these coefficients contain information about charge dynamics in a particular excitation situation, and not only about the eigenmodes, i.e., material properties (see [71] for details).

The problem of the determination of effective parameters for bulk materials does not make much practical sense. From the theoretical point of view, the only problem which can be stated and (potentially) solved for the bulk materials is finding the dispersion relation, which does not assume even introduction of the effective parameters. The effective parameters can be introduced as some coefficients between polarizability and magnetization and the electric and magnetic fields, but it can NOT be done unambiguously—see SFTs, which give birth to an unlimited number of different effective parameters, each set of them nevertheless satisfy MEs.

4. MULTIPOLE APPROACH FOR HOMOGENIZATION OF METAMATERIALS: CLASSICAL METAMATERIALS

A. Charge Dynamics in Isolated Plasmonic Meta-Atoms: Antisymmetric Modes as a Source for Magnetization

The multipole approach can be equally used to describe optical properties of metamaterials consisting of any kind of meta-atoms. In order to demonstrate this method, the multipole expansion was applied here to describe the widely used double-wire geometry [8]. In what follows we assume again the geometry shown in Fig. 2 with the electric field ${\overrightarrow{E}}_x$ polarized along the long axis of the wires and propagation along the $y$-axis $(0,{\overrightarrow{k}}_y,0)$. In order to find the relations for the dipole, quadrupole, and magnetic dipole moments it is necessary to express charge dynamics in the MAs as the functions of the averaged fields. The microscopic interaction between charges and the electromagnetic wave is determined by the interaction with the electric field; the interaction with the magnetic field becomes significant only for relativistic velocities or extremely large magnetic fields [87]. A rigorous description of the charge dynamics in terms of eigenmodes can be rather straightforwardly done for the nanospheres [88] and can be easily adopted for the elliptical particles of various eccentricities [89], but the physical picture qualitatively remains unchanged—the dynamics is described to the first approximation by the harmonic oscillation equation(s) with the eigenfrequencies, damping rates, and mutual coupling coefficients. The coefficients depend on particular geometry of the MAs and can be found analytically or semi-analytically with some accuracy, which in most cases does not allow direct comparison with the experiments and a fitting with the experimental data is required. In Ref. [71] it was assumed that the MA consists of several “harmonic oscillators” coupled to each other—for example, for the double wires in Fig. 2 the system dynamics is supposed to be modeled by two coupled harmonic oscillators. The oscillator parameters (the eigenmodes, damping constants, and coupling coefficients) are introduced phenomenologically and are supposed to be found from the fitting with the experimental data or with the data from rigorous numerical simulation. This approach allows us to keep basic physics and at the same time does not overload the model with unreasonably complicated math. It is worth noticing that the same approach is adopted for the natural materials, for example, for the analytical description of the frequency dependence of the dielectric constant, and does not seem to be a significant drawback of the model. Under these restrictions the damped and driven harmonic oscillator equation, describing the dynamics of the charge $q_k$, takes the following form [90]:

In Eq. (4.1) ${\gamma }_k$ represents the damping constant, ${\sigma }_{ki}$ is the coupling between oscillators, and ${\omega }_k$ is the eigenfrequency of the charges in the microscopic coordinates ${\overrightarrow{r}}_k(t)$ of $k$-th oscillator. The oscillators are driven by local fields ${\overrightarrow{E}}_{k,\text{loc}}(\overrightarrow{R},t)$ at the point of the oscillator locations; the relation between the local and averaged fields will be discussed later. This set of equations of motion can be analytically solved and the system parameters ${\omega }_k^2$, ${\gamma }_k$, and ${\sigma }_{ki}$ can be evaluated in a phenomenological way by comparison with experimental or numerical data.

To apply this model to the double-wire geometry we propose the charge arrangement as shown in Fig. 3:

 

Fig. 3. Reference [85]: Double-wire MM geometry and corresponding suitable charge distributions that support electric dipole, electric quadrupole, and magnetic dipole moments. The dynamics including interactions between the top and the bottom wires is described by a coupled harmonic oscillator model, which is indicated by the red arrows.

Download Full Size | PDF

The electric dipole, the electric quadrupole, and the magnetic dipole moments can be calculated straightforwardly (in the frequency domain):

From Eq. (4.3) one can recognize that all second-order expansion moments (the electric quadrupole and the magnetic dipole moment) vanish for a symmetric mode $x_1=x_2$, while the polarization (electric dipole) is still nonzero. The symmetric system $x_1=x_2$ would consist of two classical dipoles with no influence of higher-order multipoles, and hence with no magnetization.

B. Dispersion Relations and Effective Parameters for Metamaterials: Symmetric Structures (Retarded Field)

For a polarization along the wire axis ($x$ direction) the complete equations describing the charge dynamics in the two wires are given by Eq. (4.1) and in the case of equal oscillators (symmetric structure) and not equal local fields are reduced to

The inhomogeneous solution of the system (4.4) can be obtained in the Fourier domain:

For the symmetric and asymmetric oscillations, entering the formalism:

The relations between local amplitudes on the electric field in the upper and bottom wires $E_{x1,\text{loc}}(y,\omega ),E_{x2,\text{loc}}(y,\omega )$ are assumed to be caused by the retardation and are taken into account by the following expressions:

Here $E_x(y,\omega )=E_x(\omega )\exp (i{k}_{y}y)$ is the macroscopic (averaged) field in Maxwell equations. The electric field propagating from the first to the second wire is determined by $k_y$, the complex wavenumber (which is supposed to be found later from the dispersion relation). In addition to this external electric field evolution, the excitation process is also governed by the near-field coupling of the two wires. This mechanism is taken into account by the empirical coupling constant $\sigma$ between the two wires and not by the additional electric fields on the right side of Eq. (4.4).

The interaction with the electric field can be expressed in terms of electric and magnetic multipole responses that consequently determine the dependence of all quantities on the electric field. The excitation of coupled-charge oscillations leads to a magnetic response that can be described by an effective magnetic permeability, which is again a consequence of the interaction with the electric (not the magnetic) field. It should be emphasized again that the physical picture of the magnetic response differs basically from that taking place in solid-state physics. The magnetic response in the latter case is caused by a magnetic field that induces or aligns existing magnetic moments of atoms or molecules (the free-electron magnetism effect again is caused by interaction with the magnetic components of the field). In the case of MMs the electric field excites localized or surface (like in case of fish-net structure) plasmon-polaritons, which contribute to both electric and magnetic responses, while the microscopic magnetic component does not participate in the light–matter interaction.

A quantitative comparison of results obtained by the outlined analytical approach and rigorous numerical calculations has been performed. In order to determine rigorously the dispersion relation of the geometry shown in Fig. 4, the Fourier modal method [91] has been used.

 

Fig. 4. Reference [85]: (a) Geometry of the simulated double-wire meta-atom is shown. (b) Three-dimensional (3D) bulk MM alignment to calculate the dispersion relation of a bulk MM, and (c) the slab arrangement which allows additionally the calculation of effective parameters.

Download Full Size | PDF

To describe the propagation of electromagnetic waves in 3D bulk media [see Fig. 4(b)] the calculation has been examined with periodic boundary conditions in all three space dimensions. The periods used in the $x$, $y$, and $z$ directions are ${\mathrm{\Lambda }}_x=600\mathrm{nm}$, ${\mathrm{\Lambda }}_y=500\mathrm{nm}$, and ${\mathrm{\Lambda }}_z=150\mathrm{nm}$, respectively.

The double wires are formed from gold with the sizes shown in Fig. 4(a), and are separated by a thin glass layer with $n=1.44$. The corresponding effective material parameters of the same geometry have been obtained by FMM calculations for one layer [see Fig. 4(c)]. The effective permeability and permittivity are calculated from the complex reflected and transmitted amplitudes, which are used in inverted expressions for the reflection and transmission of light at a homogeneous slab [92]; see Fig. 5. The details of the calculations can be found in Ref. [85].

C. Optical Activity in Metamaterials

Research was also extended toward the exploration of meta-atoms that affect off-diagonal elements of the effective material tensors ${\overline{\chi }}_{ij}(\omega )$ (elliptical dichroism), here $P_i={\overline{\chi }}_{ij}(\omega )E_j(\omega )$. It expands the number of observable optical phenomena, leading to, e.g., optical activity [26–29,63], bidirectional and asymmetric transmission [31,93,32], or chirality-induced negative refraction [33–35]. In general, investigating the geometry of the MM (the meta-atoms geometry and their arrangement) allows us to determine the form of the effective material tensors in the quasi-static limit, as extensively discussed in Ref. [27]. From such considerations it is possible to conclude on the symmetry of the plasmonic eigenmodes sustained by the MAs and on the polarization of the eigenmodes allowed to propagate in the effective medium [31]. But in order to determine the actual frequency dependence of the tensor elements, more extended models are needed which start in their description of the MA properties from scratch [32]. Such models are required to be universal, simple, and assumption free to the largest possible extent.

This model has been developed in Ref. [84] based on the same approach presented above, but for more complex resonator shapes. Specifically, it is shown that it is possible to directly infer that the model’s predictions are valid in terms of the Casimir–Onsager relations [94,95], the requirement for time reversal and reciprocity in linear media [96].

To reveal the versatile character of this approach the consideration is directly started by conceptually replacing the planar split ring resonator (SRR) geometry, shown in Fig. 6(b), by a set of coupled oscillators. Here only final results are presented; the detailed consideration can be found in Ref. [84].

 

Fig. 5. Reference [85]: (a) Dispersion relation obtained from the numerical calculations in an infinite 3D MM as shown in Fig. 4(b), and (b) the corresponding fitted analytical version.

Download Full Size | PDF

 

Fig. 6. Reference [84]: (a) The original SRR structure (left), the first modification, namely, the L (center), and the second modification, the S structure (right); (b) the SRR together with the carrier oscillators, marked by black dots, which are used to phenomenologically replace the SRR.

Download Full Size | PDF

Considering the eigenmodes for $x$ polarization, it is observed that two parallel dipoles ($x_1=x_3$) are induced, while due to symmetry constraints no dipole is induced in the $y$ direction ($z_2=0$); see Fig. 7(c). By contrast, besides a dipole in the $z$ direction the $y$-polarized illumination induces oscillating dipoles in the $x$ direction in both arms [Fig. 7(d)].

 

Fig. 7. Reference [97]: The rigorously calculated (FMM) far-field transmission/reflection spectra compared with those obtained by the coupled dipole model (DM) for the two indicated polarization directions [(a) $x$ polarization and (b) $z$ polarization]. The stationary carrier elongation (normalized imaginary part of $x_{\mathrm{1,3}}$ and $z_2$) for the two corresponding polarizations [(c) and (d)]; (e) and (f) are the exactly retrieved parameters compared to calculations performed with the model.

Download Full Size | PDF

But due to the antisymmetric oscillation $x_3=-x_1$, the dipoles in the SRR arms [Fig. 6(b)] do not radiate into the far field because they oscillate out-of-phase and interfere destructively. Hence, no cross polarization is observed and the far-field polarization equals that of the illumination. As expected, the susceptibility tensor is diagonal. Due to the polarization-dependent carrier dynamics the SRR shows a linear dichroitic behavior.

The model correctly predicts the linear polarization eigenstates as required by the mirror symmetry with respect to the $x–y$ plane.

Results are shown in Figs. 7(a) and 7(b). The carrier oscillations are shown in the top of Fig. 7 and their horizontal position relates to the respective resonance frequency. Figures 7(e) and 7(f) show the rigorously retrieved effective material properties together with the ones of the model.

The derived effective material parameters are in excellent agreement with the rigorous results.

The next investigated low-symmetry structure is the L-MA [98,99]—see Fig. 6(a). The most significant change compared to the SRR is the appearance of off-diagonal elements in the susceptibility tensor ${\overline{\chi }}_{ij}(\omega )$, which is, however, symmetric leading also to ${ϵ}_{ij}(\omega )={ϵ}_{ji}(\omega )$. This symmetry relation is important because it is required for time reversal, known as the Casimir–Onsager relation [94,95]. As expected the tensor of the effective permittivity has the same form as that for planar optical active media [26,32] resulting in asymmetric transmission due to elliptical dichroism.

In order to check whether this simple description is valid and to reveal the relation between the SRR and the L-structure eigenmodes, numerical FMM simulations have been performed for the L-meta-atom and the results compared to the model in Fig. 8.

 

Fig. 8. Reference [97]: The far-field response of the L structure for (a) $x$ and (b) $z$ polarization. In addition to the numerical (FMM, circles) and the fitted data (DM, solid lines), the predicted spectra incorporating the SRR parameters (dashed–dotted lines) are plotted. (c),(d) In contrast to the SRR both eigenmodes can be excited for each polarization direction. The respective numerical cross-polarization contributions (FMM, circles) compared with the model-predicted (dashed–dotted lines) and the fitted (solid lines) values are shown in (e) and (f). Note that figures (e) and (f) are identical as required for such kind of effective media and are only shown for completeness.

Download Full Size | PDF

The associated eigenmodes for the carriers are shown in Figs. 8(c) and 8(d). The curves for both the copolarized [Fig. 8(a) and 8(b)] and the cross-polarized intensities [Figs. 8(e) and 8(f)] are in good qualitative agreement for the rigorous FMM results (dotted line) and the model (dashed line). Detailed consideration and discussion of the presented here L-type MA and S-type MA can be found in Ref. [84].

In order to verify the model, another modification of the SRR has been investigated, namely, the S structure [100,101]; see Fig. 6(a). The results are shown in Fig. 9. As expected, for both the S- and the L-meta-atoms the in-phase eigenmodes appear at smaller wavenumbers (larger wavelengths) compared to the out-of-phase ones. This is completely in agreement with arguments from plasmon hybridization theory [102]. Again the use of the oscillator parameters as optimized for the SRR structure [dashed–dotted lines in Figs. 9(a), 9(b), 9(g), and 9(h)] reveals the relationship between both structures, since the numerically (circles) calculated spectra agree with respect to the overall shape and the resonance positions [Figs. 9(d) and 9(e)] very well with the parametrical predictions (solid and dashed lines).

 

Fig. 9. Reference [97]: The far-field spectra of the S structure (a) $x$ polarization and (b) $z$ polarization obtained by numerical simulations (circles), predictions based on the SRR structure parameters (dashed–dotted lines), and adapting the parameters to fit the numerical values (solid lines). The carrier eigenmodes, i.e., an in line current over the entire structure and antiparallel currents (normalized imaginary part of $x_1(\omega ),z_2(\omega ),x_3(\omega )$) with respect to the center part of the S structure are observed for the two polarization directions (d) $x$ and (e) $z$. The comparison between the cross-polarization contributions $(T_{ij},R_{ij})$ for the numerical simulations (circles), the predicted lines from the SRR parameters (dashed–dotted lines), and the fitted parameter spectra (solid lines) for the found parameters deduced for fitting the copolarized response $(T_{ij},R_{ij})$. (c) The effective permittivity tensor, (f) diagonal, and (i) the off-diagonal components.

Download Full Size | PDF

A subsequent fitting again improves the results toward almost excellent agreement, which can also be observed for the cross-polarization observables [Figs. 9(g) and 9(h)]. Considering the tensor components of the effective permittivity [Figs. 9(c), 9(f), and 9(i)], one observes a difference between the diagonal entries due to the geometrical differences in the S-structure center and arms, while the antiresonance for the out-of-phase eigenmode is observed in the off-diagonal elements with the same origin as discussed for the L structure.

D. Metamaterial Analogy of Electromagnetically Induced Transparency

Electromagnetic-induced transparency is a quantum effect, which appears, for example, in a three-level system under appropriate conditions and special requirements for coherence of the fields involved. This effect can be explained by interference between different transition pathways in a three-level system under simultaneous action of two fields at different frequencies [103–105]. There are several classic systems exhibiting similar properties, namely coupled microresonators [106], resonators coupled with waveguides [107,108], metallic structures [21], and, in particular, plasmonic-induced transparency effect in MMs, which has been theoretically described [109] and investigated experimentally [22]. Another one similarity (with Fano resonances) has been both theoretically predicted and experimentally observed [110]. The results of calculation of the transmission and reflection spectra, obtained here based on the developed multipole model [84], have been compared with ones obtained in Ref. [109] from more rigorous calculations of the same structure [the structure is depicted in the insets in Figs. 10(a) and 10(b)]—see Figs. 10(c) and 10(d). It is clearly seen that the developed simple model, based on the coupled harmonic oscillator dynamics, gives very good results and can be undoubtedly used for estimations of the optical properties of MMs.

 

Fig. 10. Transmission/reflection spectra and respective permittivity of the MM with the MAs shown in insets in (a) and (b), calculated data and (c) and (d), data from [109], Fig. 3, are presented here for comparison.

Download Full Size | PDF

E. Toroidal/Anapole Metamaterials

In order to analyze an appearance of the magnetic response of the MAs, it is enough to extend the multipole expansion above the dipole level and take into account magnetic and quadrupole moments. It is commonly assumed that the next, third level of expansion containing magnetic quadrupole, octupole, and toroidal moments can be neglected if the previous second one moments (magnetic dipole and quadrupole) have nonzero contribution. It turns out that this statement can be violated, and third-order multipole contribution can interfere destructively with the first one—the toroidal moment (third level of multipole expansion) can compensate the dipole one (first level of multipole expansion), which results in a nonradiative structure, provided other moments are zeros. In order to understand this phenomena, one must refer to the general multipole expansion expressions. For the sake of generality, assume a nontrivial current distribution $\overrightarrow{j}(\overrightarrow{r},t)$ producing an electromagnetic field $E(\overrightarrow{r})$ [111]:

Here $B_{i\dots k}^{(l)}$ is a tensor of $l$-th rank. From these tensors various Cartesian multipoles in the form of sums of the unreducible tensors can be obtained. For example, $B_i^{(1)}\sim d_i^{(1)}$ determines the electric dipole moment $d_i^{(1)}$; $B_{ij}^{(2)}\sim Q_{ij}^{(2)}+{\mu }_i^{(1)}$ consists of electric quadrupole $Q_{ij}^{(2)}$ (symmetric) and magnetic dipole ${\mu }_i^{(1)}$ (antisymmetric) moments; and $B_{ijk}^{(3)}\sim Q_{ijk}^{(3)}+{\mu }_{ij}^{(2)}+T_i^{(1)}$ gives rise to electric octupole $O_{ijk}^{(3)}$, magnetic quadrupole ${\mu }_{ij}^{(2)}$, and toroidal dipole moments $T_i^{(1)}$. On the other hand, the radiation properties can be described by using the total scattering cross section in a canonical basis and can be written as a sum of intensities of spherical electric $a_E(l,m)$ and magnetic $a_E(l,m)$ scattering amplitudes [66]:

Thus, the total scattering can vanish if the spherical electric dipole scattering coefficient becomes zero, $a_E(1,\pm 1)\approx 0$, provided all higher-order scattering amplitudes are also close to zero. In order to zero spherical electric dipole, $a_E(1,\pm 1)=0$, the first-order (electric dipole) coefficients $B_i^{(1)}$ has to be compensated by the third-order coefficients $B_{ijk}^{(3)}$, which contains the toroidal dipole moments $T_i^{(1)}$. This simple consideration creates the basis for understanding of physics of the anapole mode, namely, mutual compensation of $B_i^{(1)}$ and $B_{ijk}^{(3)}$. The other two moments of the third order—octupole $O_{ijk}^{(3)}$ and magnetic quadrupole ${\mu }_{ij}^{(2)}$—are assumed to be negligible, which in fact takes place for the toroidal structure.

Typically, the toroidal moment is ignored as it appears in the third order of expansion, which is expected to be negligible. The results in Fig. 11 indicate that it is necessary to introduce toroidal moments for optically large particles to accurately describe the total scattered field, in particular for a dielectric nanodisk with high refractive index, such as the one of silicon. The total scattering cancellation is possible due to the fact that the radiation patterns of the electric and toroidal dipoles are equivalent.

 

Fig. 11. (Reference [81], Fig. 2): Decomposition of the contribution to the far-field scattering with electric dipole symmetry in terms of spherical and Cartesian multipoles. We consider scattering by a dielectric spherical particle inside as a function of diameter for refractive index $n=4$ and wavelength 550 nm: (a) scattering $|a_E(\mathrm{1,1})|$ and internal $|d_E(\mathrm{1,1})|$ Mie coefficients; (b) partial scattering cross section and energy density of the electric dipole; (c) calculated spherical electric dipole $|{\mathbf{P}}_{\text{sph}}|$ (black), Cartesian electric $|{\mathbf{P}}_{\text{car}}|$ (red), and toroidal $|{\mathbf{T}}_{\text{car}}|$ (green) dipole moments contributions to the partial scattering. These figures demonstrate that for small particles both contributions of the spherical and Cartesian electric dipoles are identical and the toroidal moment is negligible. For larger sizes, the contribution of the toroidal dipole moments to the total scattered field has to be taken into account. The anapole excitation is associated with the vanishing of the spherical electric dipole ${\mathbf{P}}_{\text{sph}}=0$, when the Cartesian electric and toroidal dipoles cancel each other.

Download Full Size | PDF

Anapoles have been introduced in the physics of elementary particles [67]. The electrodynamic analog of a stationary anapole is well-known toroid with a constant poloidal surface current, which is also associated with a toroidal dipole moment. It generates no field outside, but possible nonzero potential, which might lead to a violation of the reciprocity theorem and Aharonov–Bohm-like phenomena [112,113].

In the dynamic case, the oscillating toroidal dipole moment produces nonzero electromagnetic radiation with the pattern fully repeating the one from that of the electric dipole moment, but scaled by a factor of ${\omega }^2$. For the oscillating surface current, the radiationless properties can be kept by adding another dipole oscillating in antiphase with the toroid, resulting in complete destructive interference of their radiation due to similar far-field scattering patterns [114]. This radiationless electric and toroid dipole nontrivial current configuration has also been named “anapole” [114], which is Greek for “without poles.” Nevertheless, this compensation is not complete: the compensated toroidal dipole moment is a part of the third-order expansion, and all higher-order expansions remain radiative.

The anapole concept has attracted considerable attention in the metamaterial community as a possible realization of these radiationless objects [115]. The toroidal moment itself and the respective effects (including toroidal metamaterials [116]) have been deeply investigated theoretically [113,117]. Several experimental verifications in microwave [118] and optical [119] domains for toroidal moments confirm the theoretical conclusions.

Anapole mode in the optical domain using a simple silicon structure has been experimentally demonstrated for the first time in Ref. [81].

An anapole being a radiationless structure (in the frame of all above-mentioned limitations) is also obviously (with the same approximations) not sensitive to the external radiation. An anapole-like qubit design that is naturally insensitive to low-frequency noise and is well protected from other ambient noise sources, and therefore could be a good candidate for a superconducting qubit, has been proposed and theoretically investigated in Ref. [120].

It is interesting to compare the microscopically introduced toroidal moment with the phenomenologically elaborated equations for “T” form. Phenomenologically elaborated $⟨\overrightarrow{j}⟩=-i\omega \overrightarrow{P}-i\omega \mathrm{rot}[{\overrightarrow{F}}_1]+\mathrm{rot}[{\overrightarrow{F}}_2]$ allows two different allocations of term $i\omega \mathrm{rot}[{\overrightarrow{F}}_1]$; namely, it can be packed in polarizability $⟨\overrightarrow{j}⟩=-i\omega (\overrightarrow{P}+\mathrm{rot}[{\overrightarrow{F}}_1])+\mathrm{rot}[{\overrightarrow{F}}_2]$ or in magnetization $⟨\overrightarrow{j}⟩=-i\omega \overrightarrow{P}+\mathrm{rot}[-i\omega {\overrightarrow{F}}_1+{\overrightarrow{F}}_2]$. Final solutions of MEs of course do not depend on this relocation. This possibility of the different allocations corresponds to SFTs between the different representations.

Let us consider the case of material consisting of toroidal structures only (it means that the MAs have only toroidal moments and no other ones). The current in the frames of the developed phenomenological scheme is given by (taking into account that ${\overrightarrow{F}}_2$ is responsible for a nontoroidal magnetization and can be set to zero)

From the other side, referring to [111], the current is written through the toroidal moment (see (2.12) in Ref. [111]):

Alternatively, it has been shown [121] that the averaged current can be expressed in the following form (see (2.15) in Ref. [121] without nontoroidal contributions):

Here $\overrightarrow{T^e}$ and $\overrightarrow{T^m}$ are the contributions appearing due to the presence of structures with the toroidal moments in the frame of an alternative formulation of MEs with magnetic charges [122]. From the comparison of (3.82) with (3.83) and (3.84) in Ref. [122] it is seen that the introduced phenomenologically function ${\overrightarrow{F}}_1$ is fixed by the presence of the toroidal moments. From the other side, the exact functional form depends on the MEs initial formulation. In the case of traditional MEs it becomes

The toroidal moment determines function ${\overrightarrow{F}}_1$ introduced phenomenologically.

5. METAMATERIALS WITH INTERACTION BETWEEN META-ATOMS: PERIODIC AND RANDOM METAMATERIALS

A. Dispersion Relations for Periodic Chains of Coupled Dipoles and Quadrupoles

The interaction between the small particles (meta-atoms), either dielectric or metallic, and the propagation of an optical excitation in a regular chain of such particles has been extensively investigated [123–128]. Theoretical tools for the modeling of these chains (irrespective of the nature and sizes) remain invariant: the electromagnetic excitation in the particles is supposed to be described by taking into account all possible eigenmodes [123,125] and interactions between all particles in a chain. There are several approximations that are typically accepted in these kinds of problems. First, depending on the size of the particles, the model can be restricted by consideration of the dipole moment only (for metallic nanoparticles) [124,128]; the higher moments can be taken into consideration and similarly in the case of the investigation of magnetic response [127,129]. Usually, for the problem of only electromagnetic excitation propagation the chain of the dipole approximation is enough [130], provided the distance between particles is not less than about 3 times their dimensions. Second, the interaction between the particles in a chain can be considered in the frame of the quasi-static approximation, where no retardation between particles is retained; otherwise, interaction between dipoles contains terms proportional to the $1/r$ and $1/r^2$ in addition to the quasi-static term of $1/r^3$ ($r$ is the distance between dipoles). The problem possesses an exact solution for the infinite chain in the quasi-static limit, while taking into consideration the retardation leads to known mathematical difficulties [124]. Consideration of the finite chain is free from these excessive mathematical problems, but can be treated only numerically [124,131].

The multipole approach for the homogenization of the MM [69,85] has been extended in the case of regularly placed interacting MAs in the form of a double-wire structure. The MM consisting of the identical layers with the regularly spaced MAs in each is considered. The interaction between the MAs is assumed to be negligible in the longitudinal direction (wave propagation direction, perpendicular to the layer surfaces); in other words, the layers are assumed to be well separated from each other. The effect of interaction in the lateral direction (parallel to the layer surface) and its influence on the dispersion relation of the plane waves propagating in the MM have been calculated [132].

When there is a coupling between the MAs, the response of the medium is no longer truly local. As a result, depending upon the configuration of the MA, the medium responds differently to electromagnetic waves propagating in different directions. This phenomenon is called spatial dispersion (akin to the phenomenon of temporal dispersion, where the response of a medium at a given time depends on the history of its excitation). The direction dependence of the medium response is not just unique for the spatial dispersion—in the case of anisotropic media this effect appears as well. The qualitative markers for differentiation of these two effects have been considered in detail in Ref. [133] and [62] (see also references therein). One of the ways to describe spatial dispersion is to use a model of a chain of coupled harmonic oscillators. Such a model is adequate as a first approximation for the interactions between plasmonic nanoresonators. Eigenmodes of the response are obtained as wave solutions, giving the oscillations of the plasmonic charges in each nanoresonator (see Fig. 12).

 

Fig. 12. Reference [132]: Geometry of propagation: (a) the electric field is polarized along the long axis of the cut-wires, angular incidence gives rise to spatial modes in the ensemble; (b) the dipoles in an arbitrary triplet are labeled as $n$, $n+1$, and $n-1$. The positional coordinates of the charge clouds $x_n$ within these dipoles are also indexed with these labels.

Download Full Size | PDF

A similar approach has been used to study the effect of the interaction of the MAs on the bulk properties of the MM in the microwave frequency region [134–136]. The analysis provides a tool to analyze whether the coupling between MAs can enhance the available spatial spectrum of the propagating waves and hence increase resolution of the optical systems.

A chain of the periodically positioned dipoles (oriented with the long axis along the $x$ direction) is considered (see Fig. 12). For clarity, only one row is shown in the figure; it is assumed that rows of dipoles are placed along the $y$ direction.

The treatment of the problem remains the same, and coupling between neighboring rows is neglected. The arrangement of the dipoles is along the $z$ direction and the period of spacing is taken to be $z_0$.

The effect of coupling with adjacent oscillators is introduced via a coupling constant $\sigma$, which is a function of the distance between the oscillators. A well-known solution in the form of the transverse spatial modes can be straightforwardly obtained. Details of the respective calculations can be found in Ref. [137]. The above arguments have been extended to the case of the chain of the coupled quadrupoles in order to include magnetic response. The two double-wires forming the quadrupoles are assumed to be oriented with their long axes along the $x$ direction as before, while being separated by a small distance $2y_1$ in the $y$ direction. It is further assumed that the quadrupoles themselves are periodically spaced along the $z$ direction, with the spatial period $z_0$. The interaction of a single cut-wire with five of its nearest neighbors must be considered (see Fig. 13) for both cut-wires (that is, both $x_n$ and $x_n^{\prime }$), as they experience different excitation conditions due to the retardation of the wave propagating into the $y$ direction.

 

Fig. 13. Reference [137]: Nearest neighbor interactions—top view of the one-dimensional chain of the quadrupoles (two double-wires forming one from the three shown quadrupoles are surrounded by dashed frame). The dashed lines indicate the interactions that have to be taken into account for $x_n$. Point P indicates the center of the $n$-th quadrupole.

Download Full Size | PDF

The dispersion relations are presented in Fig. 14 for electromagnetic waves propagating in a medium with chains of the coupled dipoles and quadrupoles [137].

 

Fig. 14. Reference [137]: Electromagnetic dispersion curves for the system consisting of the one-dimensional chain of the coupled dipoles and quadrupoles for two spatial periods. (a), (b), (e) and (f) depict the dispersion relations for the dipole system, while (c), (d), (g), and (h) depict the dispersion relations for the quadrupole system. The first row depicts the real part of the normalized propagation vector $k_y{y}_1$, while the bottom row depicts the imaginary part. The values were obtained with the incident angle as parameter (blue, 0; green, $\pi /8$; cyan, $\pi /4$; red, $\pi /2$). Note the disappearance of the resonance associated with the quadrupole and magnetic dipole moments at the incident angle of $\pi /2$.

Download Full Size | PDF

The difference between the dispersion curves for dipoles and quadrupoles is in the presence of another set of spatial modes corresponding to quadrupoles, at a frequency lower than the dipole resonance frequency [see, e.g., Figs. 14(c) and 14(d)]. The response from the quadrupoles is much smaller and completely disappears for incidence at $\theta =\pi /2$, i.e., when light is propagating along the $z$ direction.

The angle of incidence has a clear and pronounced effect on the position of the symmetric response: the resonance peak shifts toward the blue part of the spectrum as the angle of incidence increases from 0 to $\pi /2$, while it is redshifted for the imaginary part [see, e.g., Figs. 14(g) and 14(h)].

More data for the dielectric permittivity and magnetic permeability have been published in Ref. [137].

B. Dispersion Relations for Randomly Positioned Dipoles and Quadrupoles

Random variation of the near-field quasi-static interaction between MAs in the form of double-wires was shown to be the reason for the effective permittivity and permeability changes in random metamaterials, as it was observed experimentally.

One of the great advantages of the developed model [69] is the ability to straightforwardly evaluate the influence of the charge dynamics of the MAs on the effective properties of the MMs. All factors influencing the charge dynamics (for example, interaction between the MAs, extra coupling of the MAs with the other objects, etc.) cause a change in the multipole expressions, which in turn changes the effective parameters. The interaction between the MAs and hence its influence on the effective permittivity and permeability can be, without difficulty, taken into account in the framework of the model presented here [137].

The problem of wave propagation through disordered systems attracts great attention in both quantum and classical physics [138]. In disordered chains of different dimensions, destructive interference between scattered waves gives rise to an existence of the localized modes, exponentially decaying in space—this effect has been originally found in solid-state physics and is known as Anderson localization [139]. The existence of delocalized modes that can extend over the sample via multiple resonances and have a transmission close to 1 was found in Refs. [140,141] and experimentally confirmed in Refs. [142,143]. Disorder-induced change of the guiding properties in a chain of plasmonic nanoparticles under small random uncontrollable disorder was considered in Ref. [131] and analogy of the Anderson localization in a chain of such particles was theoretically investigated in Ref. [144]. In the analysis presented here the effect of Anderson localization is not considered; nevertheless, it is believed that the analytical tool developed in Ref. [137] turns out to be suitable for the treatment of the similar effect in MMs with different types of disorder.

The influence of various types of disorder on the effective properties of the MMs has been thoroughly investigated as well. Light propagation and Anderson localization in superlattices was theoretically considered in Refs. [145,146]. The effect of the statistical distribution of the sizes of the MAs on the increase of losses in the operation frequency band was considered in Ref. [147]. A significant influence of a small (10%) deviation of the parameters of the microscopic resonances on the propagation wave in a wide frequency range was found in Ref. [148]. Averaging of the Lorenz-type expressions for the effective permittivity and permeability using a phenomenological probability distribution function showed that passband and negative refraction are still present under small positional disorder [149]; the results have been proven experimentally. Interaction in a chain of magnetic particles and its influence on the effective permeability was investigated in Ref. [129]. Using the introduced concept of “coherent” and “incoherent” MMs, the authors of [148] showed that the influence of disorder on long-range correlated MMs is significantly more pronounced in comparison with the same effect in short-range ordered MMs. Random variation of the interaction between MAs was shown to be the main reason for the disappearance of the long-range correlation and consequently of the “coherent” state [148]. In Ref. [131] it was shown that the difference in the electromagnetic properties of the inclusions itself is less important than the disorder in their positions. MMs formed by a self-organization display exactly this kind of disorder [150–153].

The qualitative explanation of the influence of the spatial disorder on the effective parameters follows presently. The positional disorder creates different conditions for the charged dynamics in the MAs due to the interaction between them [126,127]. This in turn leads to the changes of the averaged dipole, quadrupole, and magnetic dipole moments of the media and results in changes of the effective parameters, which are expressed through these averaged multipole moments [65,85].

Results of experiments with the 2D MMs exhibiting such in-plane disorder [154] were used as a test of the model developed in Ref. [132]. The most notable discovery is the fact that although disorder has a deterrent effect on the permittivity, the permeability seems to remain practically unaffected.

Let us assume that the charge dynamics in the microscopic multipole moments of the MAs depends on the distance ${\delta }_k$ between them (see Fig. 15). Following the approach of [65], it is necessary to average the resulting charge dynamics in the multipole moments over all possible representations. In other words, the microscopic multipole moments have to be additionally averaged over all possible distances between the MAs, which mathematically is expressed as an integral over the probability distribution function $\mathrm{PDF}({\delta }_k)$, namely,

 

Fig. 15. Reference [132]: Geometry of the MAs and their respective probability distributions; the spheres show MAs. The first (top) row shows a regular arrangement of the MAs, where each MA occupies the center of a slot with the length equal to the mean period. The second row depicts an arrangement of the MAs exhibiting random uncorrelated positional disorder (denoted by ${\rho }_k$), the extent of the disorder being governed by $\mathrm{PDF}({\rho }_k)$ as shown in the last (bottom) row.

Download Full Size | PDF

Here ${\chi }_{\text{micro}}({\delta }_k,\omega )$ is the microscopic multipole moment of the MAs, and $\mathrm{PDF}({\delta }_k)$ governs the distribution over all possible separation distances ${\delta }_k$ in a randomly arranged ensemble of the MAs.

The inter-separation ${\delta }_k$ between the two subsequent MAs is a function of the random variables ${\rho }_k$ and ${\rho }_{k-1}$, and the analytic form of $\mathrm{PDF}({\delta }_k)$ can be obtained by the use of the statistical methods if the analytic form of $\mathrm{PDF}({\rho }_k)$ is given. The quest to obtain such $\mathrm{PDF}({\rho }_k)$ and the effort to incorporate the effect of disorder into the existing multipole model is of interest. The main discerning principle of the approach presented here is in the use of the multipole model: the charge dynamics in MAs is primarily considered and calculated taking into account the interaction between MAs, which is expressed as a function of distance between them. Finally, averaging over all possible realizations of the MA separation distances gives the expression for the effective parameters. Interaction between MAs is taken into account most simply by using dipole-dipole near-field interaction in the quasi-static limit.

The interaction between quadrupoles is treated the same way, which makes the approach suitable for consideration of the magnetic properties of the MMs. In spite of the excessive simplification of the model with regard to the interaction, the approach treats the effective parameters (especially magnetic response) in a much more correct way then was done before with just the introduction of permeability and/or magnetic susceptibility. Furthermore, it is believed to provide a suitable platform for analytical or semi-analytical treatment of the problems, appearing in the case of disordered MMs. The details of the model and comprehensive discussion have been published in Ref. [132].

The results of the analysis for dipoles and quadrupoles are presented in Figs. 16 and 17, respectively.

 

Fig. 16. Reference [132]: Effective material parameter curves for dipole ensembles exhibiting positional disorder. The effective permittivity and permeability curves for disordered dipole ensembles are presented for different values of disorder. The first column pertains to values obtained for a mean period of $z_n=1.2$, while the second column relates to those obtained for a mean period of $z_n=1.8$. For the respective periodicities: (a),(b) the positional disorder function; (c),(d) the respective inter-separation PDFs; (e),(f) scaled real part of the permittivity; (g),(h) scaled imaginary parts of the permittivity. Clearly, an increase in disorder brings about a fall in the maximums of the response of the system.

Download Full Size | PDF

 

Fig. 17. Reference [132]: Dispersion and effective material parameter curves for quadrupole ensemble with $z_n=1.8$. (a) Positional disorder function; (b) inter-separation PDF; (c),(d) real and imaginary parts of k-vector (normalized with the double-wire separation distance $a_0=2y_1$); (e),(f) real and imaginary parts of effective permittivity; (g),(h) real and imaginary parts of effective permeability.

Download Full Size | PDF

The following features are clearly noted:

It is concluded that as the observed positions of the resonances and the relative magnitudes of the parameters are within the limits of approximation, the analysis is valid, and can be used to roughly predict the properties of MMs with incorporated randomness.

6. NONLINEAR METAMATERIALS: MULTIPOLE (SECOND-ORDER) AND THIRD-ORDER NONLINEARITIES

A. Nonlinear Optical Properties: Second-Harmonic Generation

Magnetic response in correlation with the enhanced nonlinearity was found in metamaterials a rather long time ago [155]. The simultaneous consideration of the magnetic dipole and the electric quadrupole is required by nature since both occur in the same order of the multipole expansion. Though this has been extensively discussed in the literature [156–160], the quadrupole moment is frequently dropped. Besides the linear properties that can be covered by this expansion, the extension of the multipole description leads to the quadratic nonlinear optical regime. This procedure of introducing nonlinearity is known from the early works in nonlinear optics [161] and is supported by several papers that observed multipole-induced nonlinear optical effects in various plasmonic nanostructures [162–170].

The developed multipole formalism has been applied to the SRR structure (see Fig. 18), which was chosen because first experiments on the second-harmonic (SH) generation (SHG) were already reported in this structure, although in a configuration amenable for nanofabrication where the induced magnetic dipole is nonradiating [171,172]. The SRR constitutes a manageable carrier number, describing the fundamental plasmonic properties properly. An extension over other geometries is straightforward.

 

Fig. 18. Reference [97]: (a) SRR meta-atom and the intrinsic currents for the fundamental electric (black solid line) and magnetic (black dashed line) mode. (b) The associated auxiliary charge distribution (red points) with predefined degrees of freedom (black arrows).

Download Full Size | PDF

The wave equation incorporates multipoles up to the second order [85]:

Now in terms of the plasmonic eigenmodes of interest the respective meta-atom has to be mapped onto the point multipoles: electric dipoles and quadrupole $P_x(y,\omega )$ and magnetic dipoles $M_z(y,\omega )$. The quadrupole part of $P_x(y,\omega )$ being proportional to the multiplication of coordinates of the charges involved in dynamics [$Q_{ij}(\overrightarrow{R},\omega )=\frac{\eta }{2}⟨\sum _s^{\text{all charges}}{q}_s{r}_{i,s}{r}_{j,s}⟩$; see Eq. (3.1)] becomes proportional to the square of the electric field, provided each coordinate (dynamics in the $x$ and $y$ directions) is linearly proportional to the electric field. Thus, the nonlinearity appears even if both dynamics (along the $x$ and $y$ directions) are linear. This fact causes confusion: the combination of linear equations should not produce any nonlinear response. This paradox is resolved if we recall that the quadrupole also contains dynamics of electrons in the corners of the SRR (see 1 and 3 in Fig. 18), where the dynamics becomes nonharmonic due to the curved trajectories. It is worth noticing that it does not always take place; for example, it does not take place in a double wire system—the nonzero quadrupole moment does not guarantee the nonlinear response. This is nicely proven by comparison with the experimental results of [172].

The details of the calculations can be found in Ref. [97]. Interestingly in this model the nonlinear response of the magnetic dipole produces no nonlinear contributions, which is supported by rigorous simulations for a corresponding SRR configuration [173]. There the magnetic nonlinear contributions have been shown to be much smaller in comparison to a convective electric current [174], which is equivalent to the quadrupole contribution in the approach presented here.

At first the solution to the wave equations incorporating the slowly varying envelope approximation (SVEA) ansatz has been performed numerically; see Fig. 19.

 

Fig. 19. Reference [97]: Evolution of normalized electric field intensity for (a) the fundamental field (FF) and (b) the second harmonic (SH) as a function of the wavenumber of the fundamental. The red lines indicate the real (dashed) and the imaginary parts (solid) of the linear dispersion relation. (c),(d) The corresponding results for the undepleted pump approximation (UDPA).

Download Full Size | PDF

The fundamental field (FF) wave evolution $E_{\omega }$ is determined by the two eigenmode resonances (both indicated by the resonances in the red lined dispersion relation), where a strong damping is observed. For frequencies out of the spectral domain of these resonances the FF wave propagates without excessive losses, as expected. For the SH wave a strong contribution at the fundamental magnetic and electric resonance can be observed. In these calculations the SHG signal originating from the electric resonance around $6400{\mathrm{cm}}^{-1}$ seems to be much stronger than in the spectral vicinity of the magnetic resonance $4000{\mathrm{cm}}^{-1}$. This originates from the strong damping of the SHG wave for the magnetic resonance (which propagates at $8000{\mathrm{cm}}^{-1}$), because in this spectral domain an enhanced damping occurs due to the presence of the electric resonance. This changes dramatically for the SH wave induced by the electric mode, since at the SH frequency the imaginary part of $k(2\omega )$ is close to zero. Thus the second harmonic originating from the electric resonance propagates almost without damping.

The calculations presented in Ref. [97] provide a physical motivation, an estimated order of magnitude, as well as an expected dispersive dependence for second-order nonlinear effects occurring for such special types of nanostructures.

The nonlinearity caused by the quadrupole moment has been named “multipole nonlinearity.” The higher-order multipoles obviously could generate the respective orders of the nonlinear response as well, which is a topic for future research.

B. Third-Harmonic Generation

Plasmon-enhanced third-harmonic generation (THG) at the magnetic resonance of fishnet MMs was reported in Ref. [166].

It was shown that the THG spectra obey the principles of the local-field-enhanced nonlinear response. It was proposed that the wavelength dispersion of the THG efficiency is defined by the spectral line of the magnetic resonance cubed. The maximum of the THG at the angles of about 20° can neither be explained by means of dispersion of the local field factor at the fundamental frequency. Finally, the position of the maximum does not coincide with the angular position of the propagating diffraction order appearance. It was shown in Ref. [175] that, first, this feature is caused by retardation effects, and second, it is specific to the antisymmetric electric current structure of the magnetic resonance.

The phase difference between the oscillators in the two layers dictates whether the resonance is antisymmetric—currents in the two layers are antiparallel to each other, Fig. 19(a)—or symmetric—currents are parallel, Fig. 19(b). The sources of the radiation are oscillations in the gold layers at the third-harmonic frequency.

A model of coupled oscillators with a nonlinear extension was used [compare with Eq. (4.4)]:

Uncompensated charges are induced at the edges of the thick wires of the MM by the external electromagnetic field with a polarization along the thin wires, as shown in Fig. 20 [176].

 

Fig. 20. Reference [175]: Parameters of the model and uncompensated charge density distribution in the unit cell of the fishnet MM for (a) antisymmetric and (b) symmetric resonances and corresponding far-field radiation patterns. The blue area between the gold layers is shown for better understanding of the layout; no influence of the dielectric is assumed in the model.

Download Full Size | PDF