Abstract

Influence of the short-range lateral disorder in the meta-atoms positioning on the effective parameters of the metamaterials is investigated theoretically using the multipole approach. Random variation of the near field quasi-static interaction between metaatoms in form of double wires is shown to be the reason for the effective permittivity and permeability changes. The obtained analytical results are compared with the known experimental ones.

1. Introduction

Metamaterials are artificial media that allow tailoring the macroscopic properties of the light propagation by a careful choice of a design for the microscopic unit cell (called the meta-atom). By controlling the geometrical shape and the material dispersion of the meta-atom, novel effects such as negative refraction [13], optical cloaking [49], as well as series of optical analogues to phenomena known from different disciplines in physics could be observed [1014]. Despite the possibility to rely on rigorous computations for describing the light propagation on the microscopic level, an enduring problem in metamaterial research is the question on how the effective material tensor looks like for a certain metamaterial. A simple and versatile analytical model describing propagation of electromagnetic waves in metamaterials has been recently developed following classical approach of Maxwell equation averaging procedure [15]. This transition from the microscopic to macroscopic system of Maxwell equations takes into account all peculiarities of carriers dynamics under the action of the resulted electromagnetic field through the introduction of multipole moments which are supposed to be represented as the functions of the macroscopic electric and magnetic fields [16]. One of the great advantages of this model is the ability to evaluate straightforwardly influence of the charge dynamics of the meta-atoms on the effective properties of the metamaterials. In fact, the multipole moments are calculated through the averaged charge dynamics in the meta-atoms. Any factors influencing the charge dynamics (e.g., interaction between the meta-atoms, extra coupling of the meta-atoms with the other objects, etc.) cause the changes in the multipole expressions, which in turn change the effective parameters. It is important for the analysis presented here that the interaction between the meta-atoms and hence its influence on the effective permittivity and permeability can be straightforwardly taken into account [17].

The interaction between the small particles, both dielectric and metallic, and propagation of an optical excitation in a regular chain of such particles have been intensively investigated [1823]. Interest to the chains of metallic nanoparticles is stipulated mainly by the request for subwavelength guiding structures for a new generation of the optoelectronic components for communication and information processing. Nevertheless, theoretical tools for the modeling of these chains (irrespective to the nature and sizes) remain invariant: the electromagnetic excitation in the particles are supposed to be described by taking into account all possible eigen modes [18, 20] and interaction between all particles in a chain. There are several approximations which are typically accepted in this kind of problems. First, depending on the size of the particles, the model can be restricted by consideration of dipole moment only (for metallic nanoparticles) [19, 23]; the higher moments can be taken into consideration as well in the case of investigation of magnetic response [22, 24]. Usually, for the problem of electromagnetic excitation propagation along the chain, the dipole approximation is enough [25], provided distance between particles is not less than about three times their sizes. 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 and in addition to the quasi-static term ( is the distance between dipoles). The problem possesses an exact solution for the infinite chain in the quasistatic case, while taking into consideration the retardation leads to known math difficulties and requires continuation into the lower half frequency plane [19]. Consideration of the finite chain is free from these excessive math problems but can be treated only numerically; the respective solutions for both longitudinal and transverse modes are presented in [19, 26].

Natural expansion of the developed models of the electromagnetic excitation transport on a chain of particles with randomly varying parameters revealed several interesting peculiarities. The problem of wave propagation through disordered systems attracts great attention in both quantum and classical physics [27]. 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 [28]. The existence of delocalized modes that can extend over the sample via multiple resonances and have a transmission close to 1 was found in [29, 30] and experimentally confirmed in [31, 32]. Disorder-induced change of the guiding properties in a chain of plasmonic nanoparticles under small random uncontrollable disorder was considered in [26], and analogy of the Anderson localization in a chain of such particles was theoretically investigated in [33]. In the present analysis, the effect of Anderson localization is not considered; nevertheless, it is believed, that the developed analytical tool in this paper turns out to be suitable for the treatment of the similar effect in metamaterials with different types of disorder.

An influence of various types of disorder on the effective properties of the metamaterials was intensively investigated as well. Light propagation and Anderson localization in superlattices were theoretically considered in [34, 35] using the model of multilayered system with phenomenological permittivity and permeability (positive and negative) in each of the layers. The effect of the statistical distribution of the sizes of the meta-atoms on the increase of losses in the operation frequency band was considered in [36] using generalized Clausius-Mossotti relation. 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 [37] using quasistatic expressions for the effective parameters. 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 [38]; the results were proven experimentally as well. Interaction in a chain of magnetic particles and its influence on the effective permeability were investigated in [24]. Using the introduced concept of “coherent” and “incoherent” metamaterials, authors of [39] showed that the influence of disorder on long-range correlated metamaterials is significantly more pronounced in comparison with the same effect in short-range ordered metamaterials. Random variation of interaction between meta-atoms was shown to be a main reason for the disappearance of the long-range correlation and consequently of the “coherent” state [39].

In the presented work, attention is primarily devoted to the extension of the multipole approach to describe in-plane disorder in metamaterials, which means the randomness in positioning of the meta-atoms in the plane of the substrate; see Figure 1. In [26], it was shown that the variations of the electromagnetic properties of the inclusions are less important than the disorder in their positions. Metamaterials formed by a self-organization display exactly this kind of disorder [4043]. Results of control experiments with the 2D metamaterials exhibiting such in-plane disorder [44] are used as a test of the model. The most notable discovery is the fact that although disorder has a deterrent effect on the permittivity, the permeability seems to remain practically unaffected. A theoretical model for such a class of random metamaterials should reproduce these observations.

The qualitative explanation of the influence of the spatial disorder on the effective parameters is in following. The positional disorder creates different conditions for the charged dynamics in the meta-atoms due to the interaction between them [21, 22]. This in turn leads to the changes of the averaged dipole, quadrupole, and magnetic dipole moments of the media and results in the changes of the effective parameters, which are expressed through these averaged multipole moments [15, 45].

This qualitative hypothesis requires further development of the existing theoretical multipole model; in particular, the interaction between the meta-atoms [17] has to be incorporated and adopted to the random character of this interaction. Let us assume that the charge dynamics in the microscopic multipole moments of the meta-atoms depends on the distance between them (see Figure 2). Following the approach of [15], it is necessary to average the resulted 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 meta-atoms, which mathematically is expressed as an integral over a probability distribution function ; namely, Here, is the microscopic multipole moment of the meta-atoms, and governs the distribution over all possible interseparation distances in a randomly arranged ensemble of the meta-atoms. In case of regular spatial distribution, each metaatom is affected by the same fields, is reduced to the delta function, and averaging (1) restores the microscopic multipole moments.

The quest to obtain such and the effort to incorporate the effect of disorder into the existing multipole model are discussed in details here. The paper is organized as follows. First, the probability model used to incorporate positional disorder into the multipole theory is described. As a test of principle and in order to create a systematic model, the approach is then applied to the simple case of randomly arranged dipoles. Then the treatment is extended to the case of randomly arranged quadrupoles. The probabilistic approach is applied to the specific case of randomly positioned meta-atoms, and the obtained results are compared with the experimental observations [44]; the mathematical procedures used to account for the other forms of disorder are highlighted.

The main difference of the approach here in comparison with the previous ones in the use of the multipole model is that the charge dynamics in meta-atoms is primarily considered and calculated taking into account the interaction between meta-atoms, which is expressed as a function of distance between them. Finally, averaging over all possible realization of the inter-metaatom distance gives the expression for the effective parameters. This paper is primarily devoted to the elaboration of the model and to the effective parameters calculation; further applications of the presented approach (disorder in propagation direction, transition “coherent-incoherent” states, influence of the Anderson localization on the effective parameters, etc.) will be done elsewhere. Interaction between meta-atoms is taken into account using a simplest way of dipole-dipole near field interaction in quasi-static limit; extrapolation of the model on the dynamic case is left for the future work. The interaction between quadrupoles is treated the same way, which makes the presented approach suitable for consideration of the magnetic properties of the metamaterials. In spite of the excessive simplification of the interaction, the presented model treats the effective parameters (especially magnetic response) in much more correct way than it was done before by just introduction of permeability and/or magnetic susceptibility, and is believed to provide a suitable platform for analytical or semianalytical treatment of the problems, appearing in the case of disordered metamaterials.

2. Modeling of Positional Disorder

The problem of the positional disorder modeling can be tackled in several ways. The most general formulation of the problem requires a Markovian treatment. As an illustration, the one-dimensional equivalent of the problem is considered here. Supposing that the meta-atoms are introduced one by one on a line of given length, the probability that a particle will take up a certain position on the line, and hence the probability of a particular inter separation distance, depends not only on the last particle, but also on the history and existing configuration. This is the essence of the Markovian approach. Standard techniques exist for tackling such problems, formulating a rigorous treatment. Nevertheless, due to its complex nature, an alternative simpler math treatment is developed in the present paper.

The math treatment accepted in the presented work is stipulated by the real technological chains for the producing of the nanostructure. When performing control experiments, masks for random metamaterials are manufactured by e-beam lithography methods. The writing algorithms are modified so that instead of a periodic grid a randomized one is generated by the scanner. The extent of randomness can be specifically controlled, and statistically relevant parameters such as the mean period and the variance can be assigned to each mask. Translating the above approach to mathematics, one assumes that the meta-atoms are initially placed in well-defined slots of equal length (see Figure 1) and then are perturbed from their mean positions. The perturbation is described by . This PDF describes the extent of the perturbation and also ensures that the displacement is not beyond a certain space slot. However, is the positional disorder distribution function not the PDF for the metaatom interseparation ; the latter has to be found based on given , and can then be used in averaging procedure (1) to obtain the required material parameters. In this work, interseparation is obtained by employing a characteristic function approach. The characteristic function is given by a Fourier transformation of given : where satisfies the normalization condition: Alternatively, the characteristic function can be considered as an expectation value of function : The one-dimensional equivalent of the problem is formulated as follows. A periodic arrangement of meta-atoms on a given length is considered (see Figure 1). The spacing period for the slots is given by , and the location of the th metaatom can thus be given as . Now, the perturbation of the th metaatom from its mean position can be given using a random function , such that. So, the spacing between the meta-atoms is given by. The random functions and are completely independent from each other. For the present problem, the random function of interest is The mathematical form of has to be found. The form of the above probability function can be obtained by using characteristic functions. The characteristic function formed by a sum or difference of two or more PDFs is nothing but the product of the characteristic functions of the ingredient PDFs. That is, if one is interested in a probability distribution function for variable , given by where , and so forth are independent random functions and are their mutually independent characteristic functions, then the characteristic function of will be given by So, the form of can be simply obtained by an inverse Fourier transformation on the product of the characteristic functions.

Define the characteristic functions for and as Then making use of the above-stated method, the characteristic function of the required is. Then, required can be obtained by simply using the convolution theorem: Hence, required is the autocorrelation function of the positional disorder . The integral is taken over the displacement of the metaatom from its mean position and limited by the finite values of the slot length. The strength of the method is the fact that no explicit assumption has been made regarding the form of describing the positional disorder.

The mathematical procedure has to ensure that the perturbation does not become so large that the meta-atoms overlap each other. In the analysis, the particles are assumed to be placed in average in the center of the slots of a length equal to the mean spacing period. The particles can randomly move within their own slot, and the extent of the displacement from the center of the slot is given by . A consequence of such a restraint is that has to be restricted and normalized within this slot. Figure 5 shows how the autocorrelation function approaches a triangular function from its initial Gaussian form, as the position of the particle within the slot becomes completely random (i.e., takes a rectangular form). A simple algebraic form of the probability distribution function cannot be obtained due to this truncation. Hence, the following approach was adopted: the normalized versions of for the inter separation were obtained using numerical code, and they were subsequently used for numerical integration (according to (1)) for obtaining the effective material parameters.

3. Case of Randomly Positioned Dipoles

In this section, using the above mentioned principles, the effect of disorder in a chain of periodically placed dipoles is investigated. The geometry is given in Figure 3. The bold arrow shows the direction of propagation of the electromagnetic wave.

The system can be mathematically modeled as follows. Considering the coupling dynamics between two equal adjacent oscillators, one can write the equation describing their dynamics as The term on the right side is the same for both oscillators, and as the same field impinges on both of them. By substituting the time ansatz , the system can be easily solved for and as follows: Thus, Here, and are the coupling constant and the distance between the oscillator. It is assumed that the interaction between the oscillators is the near field one between the dipoles that stipulates the inverse cubic distance dependence in the second equation in (15).

The response of the system can thus be obtained by monitoring the susceptibility of the medium. The polarization of the system can be written as so that the effective susceptibility is To incorporate the effect of disorder, following (1), the averaged form of the above susceptibility can be obtained as or where is the inter separation PDF and here quantizes the amount of disorder presented in the system.

4. Case of Randomly Positioned Quadrupoles

The extension of the above model to metamaterials (i.e., taking into account the magnetic response) firstly requires that the interaction between the adjacent meta-atoms is taken into consideration. The system is taken to be similar as the one shown in Figure 3, but the dipoles are now replaced by quadrupoles (Figure 4). The long axis of the cut wires is oriented along the -axis. The cut wires forming the quadrupole are separated along the direction. The meta-atoms are arranged randomly (in terms of above described random positioning in the respective slots) along the direction.

Assuming that a plane electromagnetic wave now propagates through the ensemble along the direction, while its electric vector is polarized along the direction, the coupled dynamics of two meta-atoms can be modeled via four coupled oscillator equations: where and is the value of the coupling constant measured for the interseparation . The magnitude of the coupling constant varies inversely as the cube of the distance, and so its value can be obtained for other interseparations—here and the diagonal distance . The exponential phase factors in the right side take into account the retardation effect. It is clear that a change in the excitation conditions will affect the form of the right hand side of the above equations, while a change in the configuration of the meta-atoms can be accounted by a change in the form of the coupling coefficients. The procedure of determining the response of the medium then remains the same—one seeks to determine effective susceptibilities (corresponding to the symmetric and antisymmetric modes of oscillation), average them over all possible coupling configurations, and then use these values for ascertaining the effective material parameters.

The first step is to find the solution of the above set of equations. They can be transferred to the Fourier domain by using ansatz , , so that the system can be rewritten in a matrix form: The modes of oscillation of interest are given by . The system can be solved to obtain the values of and , and the modes of the system can be written as. and hence one can define the effective susceptibilities: where the and dependences are due to , , and . Due to this form of definition, the functional forms of the polarization, quadrupole moment, and magnetizaton remain the same:

The effect of disorder can then be taken into account by carrying out the extra averaging integration (1): or, more explicitly:

where is the mean period. The limits of the integration indicate that the autocorrelation procedure for has been already carried out. This integral can be solved numerically for a given value of frequency.

With the effective susceptibility as defined above, one may now consider a planar metamaterial, which is formed by the identical rows of the randomly positioned meta-atoms. The effect of randomness is taken into account by the averaging procedure, and hence the dispersion relation and the effective material parameters can be written in analogy to [16]. Adhering to the same conditions of geometry and excitation, the following expressions can be utilized:

The previous expressions can be easily carried over to a numerical code to obtain the material parameters of interest. The following section presents the results and compares them with the experimental observations.

5. Method of Numerical Implementation

For convenience, the integrations and other expressions have been converted to their normalized versions. The frequencies are normalized with respect to the resonant frequency of the independent cut wire, while the distances are normalized with respect to the cut-wire spacing . Specifically, the susceptibilities for the case of dipoles and quadrupoles are

where is assumed to be Gaussian: The excursion of the dipoles around their mean positions had to be limited within the interval ; recalculation into is given by (12). The integrals cannot be taken analytically and was done using the mathematical software MATLAB. Truncation of the positional PDF was achieved by coding. To obtain the autocorrelation of the PDF, a standard subroutine was used. The results of the operations are shown in Figure 5.

All the constants used in the analysis were taken from [16]. The spacing between the cut wires was taken to be 65 nm, and the resonant frequency of an isolated cut wire was taken as  rad s−1. The damping coefficient was taken to be  rad s−1. The mean periodic spacing (the mean spacing between the meta-atoms was taken to be 1.8 times the cut wire spacing ).

To verify the correct functioning of the code, the results for a very small disorder were compared with the result for a perfectly ordered system (with neighboring meta-atoms interacting with each other); see Figure 6.

6. Results

The results of the analysis for dipoles are presented in Figure 7, and the results of the analysis for quadrupoles are presented in Figures 8 and 9.

The analysis was carried out for two values of the spacing period , namely, (Figure 8) and (Figure 9). The positional PDF was taken for the four different values of the standard deviation , and consequently the inter separation was obtained using numerical coding in MATLAB. The effective susceptibilities were calculated by numerical implementation of the integration (29) and (31), and then the effective material parameters (28) were calculated.

The following features are clearly noted.(i)For the disordered dipole ensemble, the fall in the permittivity with increasing disorder is clearly visible (Figures 7(e), 7(f), 7(g), and 7(h)). The curve is symmetric for very small values of the variance (here, ). But as the disorder increases, the peak shifts towards lower frequencies, and the curves becomes broadened and asymmetric. The reason for these observed effects can be explained as follows. A disordered system can be thought to be made up of several different periodic systems. The resonant frequency for each periodic ensemble depends inversely on its spatial period. If the response of the disordered system is approximated by the sum of the responses of its constituent periodic systems, it becomes evident that the final curve will develop a tail approaching the blue end of the spectrum. The asymmetry can thus be attributed to an inverse power relationship between resonance frequency and interseparation. The effect of broadening is a consequence of particle conservation. On the other hand, the lowest frequency/largest wavelength of response is not a function of the periodicity, but is actually limited by the eigenfrequency of the independent oscillator. In fact, the resonance frequency approaches the eigenfrequency for a periodic assembly of dipoles, when the spatial period becomes large. Hence, as the disorder in the system increases, the curves become broadened and asymmetric, and the peak response shifts towards the eigenfrequency of the independent oscillator. (ii)In the case of the quadrupole ensemble, a decrease in the value of the electric permittivity is observed as is increased. This is in agreement with the experimental results. However, there is also a decrease in the value of the magnetic permeability. This decrease is more pronounced for as compared to . This is an unexpected result, as the magnetic response should remain almost unaffected. The reason for this discrepancy could lie in the simple form of the probabilistic model chosen to describe the randomness. (iii)Generally speaking, the final expressions for the permittivity and the permeability were derived under several approximations, associated with (1). The observed discrepancy could also be attributed to these approximations. Above all, the fundamental limitations of the multipole theory itself could affect the final results as well. These possibilities have to be investigated further.

In the light of the above arguments, 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 metamaterials with incorporated randomness.

7. Other Forms of Disorder

In the preceding analysis, the effect of positional disorder (arising due to aperiodicity) on the averaged material parameters was considered. In a random metamaterial other forms of disorder can exist as well. A particular case of interest is positional disorder along the cut-wire axis; see Figure 10.

If this form of positional disorder is taken into consideration along with the aperiodicity, the model would then be a step closer to emulate a true self-organized random metamaterial [46]. In the multipole model, the individual cut wires are replaced by dipoles. In case the quadrupoles are disarrayed, the coupling between them will also be a function of their relative angular positioning. This angular dependence can be introduced into the coupling terms of the dynamic equations. More specifically, the coupling constants and in the differential equations will include the angular dependence. All other mathematics remaining the same, the averaging procedure can now be carried out between angles .

The curves in Figures 12 and 13 summarize the results obtained by using the multipole approach.

Two forms of the distribution function were used in the analysis (Figure 11). The extent of disorder is correlated to the relative angular position of the dipoles. The first form of the angular PDF distribution function used was a rectangular function (see Figure 11(a)). The second form used was a Gaussian distribution (Figure 11(b)), the random variable being the relative angular position. The function is centered about 0 degrees, the extent of disorder being quantified by the standard deviation ; and are multiplied by the term to incorporate the angular dependence. Clearly then, when , there is no interaction between the cut wires. The Riemannian integration is limited between the values . The constants were again taken from the original reference [22], and the mean periodicity was taken to be . The results (Figures 12 and 13) show that both the effective permittivity and permeability are clearly affected by the angular disorder. In a similar fashion, in-plane and out-of-plane skew disorders of meta-atoms can also be accounted by the model, by making appropriate changes to the coupling terms.

8. Conclusions

In extending the multipole approach [16] to the case of random metamaterials, the effect of spatial distribution of the meta-atoms was taken into account by considering the near field coupling between neighboring meta-atoms. The effective susceptibility was expressed as a function of the inter-separation between meta-atoms, the ensemble averaged susceptibility was then obtained as an expectation value, weighed by the probability distribution function of all possible inter-separations. In the present work, the disorder was considered only along one direction (in-plane, perpendicular to cut-wire long axis), assuming that adjacent rows of meta-atoms do not interact with each other. Results obtained by the numerical implementation of the equations indeed confirm the experimental finding that increasing disorder has a more pronounced effect on the effective electrical permittivity than on the effective magnetic permeability. For smaller periodicities, however, the electrical permittivity and magnetic permeability are affected equally. This conflicting result could be caused by a coupling of the quadrupole moments of neighboring meta-atoms. This has not been explicitly considered in the present version of the model. Also, other factors such as the effect of incident polarization, or the coupling between adjacent rows, have not been considered in the present theory. The understanding gained from the study of this simple case can now be used to account for the above specific and more involved cases.