Interaction of metamaterials with optical beams

We develop a general theoretical approach to describing the interaction of metamaterials with optical beams. The metamaterials are allowed to be anisotropic, chiral, noncentrosymmetric, and spatially dispersive. Unlike plane waves, beams can change their field distributions upon interaction with metamaterials, which can reveal new optical effects. Our method is based on a vector form of the angular spectrum representation and a technique to calculate the wave parameters for all required directions of wave propagation. Applying the method to various metamaterial designs, we discover a new optical phenomenon: the conversion of light polarization by spatial dispersion. Because of this phenomenon, the refractive index and impedance cannot be introduced for many metamaterial designs. In such cases, we propose an alternative approach to treating the beam–metamaterial interaction. This work takes a step forward in describing optical metamaterials by moving from unphysical plane waves to realistic optical beams.


Introduction
Modern computational methods and nanofabrication techniques facilitate the design, fabrication, and study of crystalline optical metamaterials with nanostructured unit cells that are smaller than visible-light wavelengths. Optical waves propagating in such materials are not split by diffraction, and therefore can be considered to be interacting with a homogeneous medium. Materials like these are called metamaterials. The nanostructured artificial molecules of these materials are designed to give the material extraordinary optical properties. For example, if the unit cells are designed to show optical magnetism, one can achieve negative refraction, which can for example be used for near-field focusing and 'perfect' imaging [1][2][3][4][5][6]. Other well-known applications demonstrated and proposed for these materials are elements enhancing the optical density of states and energy transfer [7,8], as well as optical coherence [9][10][11], elements allowing aberration-free propagation [12], and a variety of novel elements composed of temporally dispersive [13], nonlinear [14,15], and optically inhomogeneous and anisotropic materials [16,17].
The theoretical basis for the description and design of optical metamaterials is still incomplete. For example, the majority of functional optical metamaterials are spatially dispersive, which is connected to their relatively large unit cells (required, e.g., for obtaining optical magnetism [18]). For this reason, these materials are characterized by the so-called wave parameters, which are the effective refractive index n and impedance Z experienced by a plane wave and dependent on its propagation direction in the material [19]. These plane-wave parameters are in fact what metamaterial designers are usually interested in when creating new metamaterials. Several effective-medium approaches have been introduced to describe even anisotropic and spatially dispersive metamaterials in terms of their plane-wave response [20][21][22][23][24][25][26][27][28]. However, there is a substantial difference between optical beams-which we always deal with in experiments-and plane waves when they interact with such anisotropic and/or spatially dispersive metamaterials. For example, due to optical anisotropy and spatial dispersion, the energy flow and the propagation of the wavefronts, which are determined by the Poynting vector and the wave vector, respectively, can have different directions. Furthermore, spatial dispersion can affect the angular spectrum of a beam interacting with a metamaterial even in the absence of optical anisotropy. In spite of this significant difference, a general model for the interaction of metamaterials with optical beams has not previously been introduced.
In this work, we introduce for the first time a self-consistent theoretical basis for describing the interaction of optical beams with metamaterials that can be spatially dispersive and anisotropic. We propose treating the beams in terms of their plane-wave decomposition, known as the angular spectrum representation, written in the rigorous vector form. Each plane-wave component of the beam is considered separately using, for example, the values of n and Z corresponding to its propagation direction and polarization. Then, the beam profile at any transverse plane is obtained by summing the propagating plane-wave components over the whole angular spectrum of the beam. In this way, we obtain a complete description of the interaction of the beam with the metamaterial.
By applying our approach to various types of metamaterials, we have found that the presence of spatial dispersion often leads to an inevitable polarization conversion of plane waves propagating in the material, even if the material is not chiral and the wave polarization is such that optical anisotropy does not affect it. In such cases, it is impossible to introduce polarization modes for the material (the waves preserving the polarization) and the parameters n and Z associated with them. To solve this problem, we propose an alternative approach to treat the plane-wave components of optical beams. This approach is based on the transmission and reflection matrices, taking into account the polarization conversion. We also show that for metamaterials with negligible evanescent-wave coupling between metamolecular layers [20,21], arbitrarily thick metamaterial slabs can be fully characterized by the transmission and reflection matrices of a single metamolecular layer. The theory presented in sections 2 and 3 allows us to describe, design, and optimize general polarization-converting metamaterials and their interaction with optical beams. Section 4 introduces two particular examples of spatially dispersive materials, one of which illustrates the mentioned polarization conversion, and the other illustrates the spatial filtering of an optical beam by its reflection from a metamaterial slab. Section 5 presents our conclusions.

Angular-spectrum representation of optical beams interacting with metamaterial slabs
The optical properties of metamaterials are often described in terms of effective plane-wave parameters, such as refractive index n and impedance Z, which depend not only on the wave frequency and polarization, but also on the direction of the wave vector k. These parameters are connected to and, as a rule, are also evaluated from the plane-wave transmission and reflection coefficients k t ( )and k r ( )of a slab of the metamaterial in question [29][30][31][32]. In the following, using the fact that optical beams are fully characterized by the properties of their planewave components, we introduce a general technique to treat the interaction of arbitrary optical beams with metamaterial slabs.
Our approach is based on the vector angular-spectrum representation [33]. Let us consider an optical beam that is incident on a metamaterial slab characterized by known plane-wave parameters k n ( )and k Z ( ), or equivalently, k t ( )and k r ( ). To obtain the field distributions in the reflected and transmitted beams, we apply the following three-step procedure. In the first step, the incident beam is decomposed into plane waves by calculating the angular spectrum of the beam at the entrance facet of the slab. Then, each plane wave is allowed to interact with the slab and split into the transmitted and reflected waves. The complex amplitudes of these waves are found by multiplying the incident-wave amplitudes with the corresponding coefficients, k t ( )and k r ( ). Finally, the transmitted and reflected beams are found by combining the plane-wave components obtained in the previous step. This Fourier optics based approach has not yet been applied to optical metamaterials.
The plane-wave decomposition of the incident beam can be written in the form of the following twodimensional inverse Fourier transform [33] ( ) The slab entrance and exit surfaces can be assumed to be located at z = 0 and z = D, respectively, in which case (1) describes the beam at z 0 < . The amplitudes E k k ẑ ( , ; ) x y i are found from the known electric field, E x y z ( , ; ) i , of the beam using the Fourier integral x y k x k y that reveals the angular spectrum. To treat the reflection and transmission of the obtained plane waves by the slab, it is necessary to split these waves into their s-and p-polarized components, respectively. Only then can their amplitudes be multiplied with the corresponding transmission and reflection coefficients. The vector Ê i can be written in terms of the s-and p-components using the following matrix equation . From now on, we denote an ordinary matrix product of â and b with abˆand assume all the vectors to be column vectors. Whenever we use notation a b · , we mean a dot product of vectors a and b. In general, the polarization of each s-or p-polarized plane-wave component can change upon reflection and transmission by the slab. It is therefore necessary to introduce transmission and reflection matrices, instead of scalar coefficients, such that x y x y x y x y x y r i = The reflection matrix, k k r ( , ) x y , in (5) includes the directional change of the p-polarization upon the reflection that takes place even if the material is isotropic.
Next, the fields E k k D ( , ; ) x y t and E k k ( , ; 0) x y r are inverse Fourier transformed, as in (1), which results in the final expressions for the transmitted and reflected beams, respectively: k r In these integrals, the vectors E k k tˆ( , ; 0) are written in terms of their Cartesian vector components. Equation (6) contains an exponential factor, ik D exp( ) z − , because the transmission coefficient t is defined in accordance with (4). Equation (6) is valid for z D > , while equation (7) is valid for z 0 < . These equations include both propagating and evanescent plane-wave components and allow one to study the action of a metamaterial slab on an optical beam, as long as the matrices t and r are known. In the following section, we show how these matrices can be evaluated.

Transmission and reflection matrices
Essentially all optical metamaterials demonstrated so far are both spatially dispersive and anisotropic. Hence, for them, the matrices t and r must be calculated at each angle of incidence separately, even if t and r can be expressed through the effective wave parameters, such as the refractive index and impedance. Sometimes, however, these parameters are worthless or impossible to introduce. For example, in the presence of a strong evanescent-wave coupling between metamolecular layers in a slab, the refractive index and impedance can depend on the slab thickness (D), so that their evaluation at each value of D does not make much sense. Also, the wave parameters cannot be introduced if the material does not support any polarization modes (e.g., if the polarization of an arbitrarily polarized plane wave will change upon the wave propagation).
If the interlayer evanescent-wave coupling is strong, the Cartesian elements of the matrices t and r defined by (4) and (5) and used in (6) and (7) must be evaluated directly for the whole metamaterial slab. The wave parameters, such as n and Z, should not be used because they will depend on the slab thickness. If, on the other hand, the evanescent-wave coupling is weak, the metamaterial can be considered to be homogenizable [21]. As a consequence, the matrices t and r can be written for an arbitrarily thick metamaterial slab in an analytical form in terms of the transmission and reflection matrices of a single metamolecular layer (τ andρ, respectively). In addition, if the wave parameters can be introduced, they will not depend on the thickness of the slab and will therefore characterize the material itself [21]. We also emphasize that usually both the near-and far-field interactions between the metamolecules in each layer are strong, but owing to this interaction, the evanescentwave decay length of the layer is short compared to that of a single metamolecule; it is, in fact, shorter than the lattice period [21], while for a single dipole the near-field decay length is on the order of λ.
To evaluate t and r for a given homogenizable metamaterial slab, we first find the monolayer matricesτ and ρ (e.g., numerically). These matrices are evaluated by considering the layer as an infinitely thin sheet located in the middle of the layer and surrounded by the host medium [20,21]. The monolayer matrices have the form , In (8), each matrix element, , ij τ is defined as the field ratio of the transmitted i-polarized plane wave to the incident j-polarized wave, and each ij ρ in (9) is analogously defined for the reflected waves. In numerical calculations, one must therefore consider both s-and p-polarizations. Having evaluatedτ andρ, one can proceed to the calculation of t and r. Below we present two alternative methods to do this. The first method (section 3.1) is based on a generalized transfer-matrix formalism, and the second method (appendix) makes use of recursive matrix relations and the Z-transform. In both of the methods, the following phase-shifted transmission and reflection matrices are used The quantity z Λ in these equations is the lattice constant of the metamaterial along the slab normal. Section 3.2 describes the conditions for the existence of the polarization modes and also presents a method for evaluating the effective wave parameters, n and Z, for these modes.

Transfer matrix method
We assume that the monolayer transmission matrices, fˆ+ and fˆ, − and the reflection matrices, ĝ + and ĝ , − for the plane waves propagating in the positive (+) and negative (−) z-directions are already known. The propagation angles of these waves are, respectively, θ and 180 -θ°, with θ evaluated with respect to the z-axis. Considering the waves propagating at these angles between some adjacent metamolecular layers n and n 1 + in the material, we can write the following relations for the Jones vectors J n ( ) is a 4 × 4 transfer matrix corresponding to layer n. This transfer matrix is a block matrix that is also valid for materials for which polarization modes do not exist, unlike transfer matrices conventionally introduced for photonic crystals and other periodic scattering structures [34][35][36]. The corresponding matrix of an arbitrarily thick N-layer metamaterial slab is then given by T Tˆ. (16) n N = Since T n is diagonalizable, the elements of T are easy to obtain analytically in terms of fˆ± and ĝ ± . Then, in analogy with (15), one can write the transfer matrix as This immediately yields the following set of equations for the slab transmission and reflection matrices: The obtained matrices can be used in (6) and (7) to treat the interaction of an optical beam with an arbitrary polarization-converting metamaterial slab.
The matrices t ± and r ± can be obtained in an alternative way with the help of a recursive matrix relation method, which is presented in the appendix. This alternative approach can in some cases be more convenient. Also, it can be used to check the validity of (18)-(21).

Polarization modes
It is not always possible to introduce polarization modes for a given metamaterial. In this section, we obtain the conditions under which the modes can be introduced, and consequently, their effective refractive index, n, and impedance, Z, can be calculated. Introducing the modes provides physical insight into the interaction of the metamaterial with light, helps in designing metamaterials for prescribed applications, and allows one to analytically obtain the mode transmission and reflection coefficients, in terms of n and Z, valid for any media surrounding the material.
The polarization modes must preserve their polarization state and have the periodicity of the lattice. Hence, for such a mode we can write where z, γ ± is the effective propagation constant of the mode. Suppose that the slab has N metamolecular layers. Then, for the last layer, we can use (13) This equation implies that J n ( ) − must be an eigenvector of the transmission matrix, fˆ−. Analogously, one can use (22) and (12) and show that J n + is an eigenvector of fˆ+, (27) shows that J n ( ) + must also be an eigenvector of the double reflection matrix, g ĝ− + . Similarly, it can be shown that J n ( ) − is an eigenvector of the matrix g ĝ+ − . Physically this result expresses the fact that each mode reproduces its polarization after two reflections. To find the modes, one should therefore diagonalize the matrices fˆ± and g ĝ∓ ± simultaneously and find their common eigenvectors.
Once the eigenmodes are found, one can calculate the effective refractive index and impedance for them. Combining (22) and (A.1), we write In the following, we consider a symmetric metamaterial with ĝ ∓ and ĝ ± , which have common eigenvectors. For these metamaterials, eigenvectors of the double reflection matrix g ĝ∓ ± are also the eigenvectors of the single reflection matrices, ĝ ∓ . Hence, these matrices are diagonal, which leads to diagonal matricesα ± andβ ± obtained from (A.2) and (A.3), assuming that fˆ± is also diagonal. Equation (28) then yields where i, α ± and i, β ± are the diagonal elements of matricesα ± andβ ± and z e i z z , = γ Λ ± − ± . Subindices 1 and 2 refer to modes 1 and 2 of the material. In fact, (29) follows directly from (A.14). The solution of (29) is of the form m i ln 2 4 which is identical to the expression given in [19] and [21]; m is an integer. The effective refractive index is then found from The impedance is obtained as the ratio of the electric and magnetic fields. For example, the impedance of s-and p-polarized eigenmodes is given by the expression [19,21] ( ) ( ) where p is equal to 1 for s-polarization and −1 for p-polarization. Note that the metamaterial in question can be spatially dispersive, anisotropic, noncentrosymmetric, and internally twisted. For a flat boundary between two spatially dispersive metamaterials with known z, γ ± and Z ± , one can obtain the following generalized Fresnel coefficients [19,21]: where subindices i, t, and r refer to the incident, transmitted, and reflected wave, respectively. The coefficients above can be used to calculate the Fabry-Perot transmission and reflection coefficients for a metamaterial slab in a standard way, and by using (6) and (7), one can evaluate the field distributions in optical beams transmitted and reflected by the slab.

Examples of metamaterial slabs interacting with optical beams
In this section, we apply the developed theoretical methods to particular examples of interaction of spatially dispersive and polarization-converting metamaterials with optical beams. The first example demonstrates that even in the case of highly symmetric metamolecules such as silver discs, the metamaterial can exhibit significant polarization conversion, and the polarization modes for it can be impossible to introduce for all relevant incidence angles. This means that for many metamaterials, the beam-propagation characteristics can be obtained only by taking into account the polarization conversion of the beam plane-wave components. The calculations in such cases must be accomplished using the described matrix approach without resorting to the concepts of the effective refractive index and impedance. In the second example, a spatially dispersive metamaterial that acts as a reflective spatial filter for a focused laser beam is introduced. It consists of disc pairs aligned in the crystal lattice such that the polarization modes of the material are the s-and p-polarized plane waves. We evaluate the effective refractive index and impedance for these modes and, using them, we calculate the intensity profiles and polarizations of optical beams transmitted and reflected by a slab of this spatially filtering metamaterial.

A silver-disc metamaterial
Consider a metamaterial slab consisting of a three-dimensional array of silver nanodiscs embedded in glass. The discs form a tetragonal lattice with transverse periods 120 x y Λ Λ = = nm and a longitudinal period of 180 z Λ = nm. The discs have a 30 nm radius and are 20 nm thick. Their normals are aligned in the xy-plane at an angle of 45°to the x-axis, which makes the discs perpendicular to the entrance and exit facets of the slab. A 'monomolecular' layer of the material is shown in figure 1.
To find the polarization modes of this metamaterial and, for example, calculate the refractive index and impedance associated with these modes, we numerically evaluate the transmission and reflection matrices,τ andρ, of a single layer of such metamolecules in glass. We do this for various angles and planes of incidence at the wavelength of 633 nm of a HeNe laser. The calculations were done with the help of COMSOL Multiphysics [37] computer software. The optical characteristics of silver are taken from [38]. As an example, figures 2 and 3 show the polar plots of the matrix elements (absolute values and phases, respectively) of the calculatedτ andρ as functions of the incidence angle θ for the plane of incidence coinciding with the xz-plane. The full threedimensional surface plots for the elements ofτ andρ are also shown together with the corresponding cross sections. The basis vectors are the s-and p-polarization vectors, as explained in section 3. The off-diagonal elements of the matrices are not negligible, which implies that the s-and p-polarizations are not conserved upon their transmission and reflection by the layer. Therefore, the s-and p-polarized waves are not the modes of the material. For each plane of incidence, it can be found that the eigenvectors of the matricesτ andρ are the same only at θ = 0 and θ = 180°, in which case the modes have the electric-field vectors perpendicular and parallel to the discs. At other incidence angles, the modes are not possible to find, primarily because of the difference between the diagonal elements of the reflection matrix. This difference makes its eigenvectors differ from those  ofτ . We also want to emphasize that the off-diagonal elements ofτ are nearly equal, in both the amplitude and the phase, to those ofρ. This is as it should be, because the transmitted and reflected waves polarized orthogonally to the incident wave are the waves radiated in the forward and backward directions primarily by the excited dipole moments in the discs. The abrupt changes of the phase observed for the off-diagonal elements ofτ andρ in the three-dimensional plots of figures 3(b) and (c) are purely geometrical, and they are connected to the symmetry planes of the discs and the direction of the electric field vector of the incident s-polarized wave that ignores this symmetry. Returning to the polarization modes, we find that they can be introduced for all incidence angles corresponding to two particular planes of incidence, one parallel and the other normal to the discs. For all other planes of incidence, the polarization modes do not exist.
In spite of the fact that the modes and the related refractive index and impedance cannot be used to characterize the metamaterial, its interaction with optical beams can be comprehensively studied by using the transfer-matrix formalism of section 3 (or the recursive-relations based approach) and the angular spectrum representation of section 2. First, the plane-wave transmission and reflection matrices, t and r, of the slab are calculated for all relevant angles of incidence, using (18)- (21) or (A.4) and (A.5). We have verified that the two approaches give the same results. Then, using (2) and (3) the incident beam is split into its s-and p-polarized angular-spectrum components. Each of them is multiplied with the corresponding matrix to obtain the plane waves after their interaction with the slab [see (4) and (5)]. Finally, the obtained plane waves are superimposed by Fourier integration in (6) and (7) into the transmitted and reflected beams. Applying this procedure to a 15layer slab that interacts with a three-dimensional Gaussian beam, we obtain the intensity profiles shown in figure 4. The beam is focused onto the slab surface at normal incidence and has a divergence angle of 30°. It is linearly polarized, such that the tangential component of its electric-field vector is everywhere on the slab surface directed along the disc normals. In this case, the optical anisotropy of the material cannot change the beam polarization. Figure 4(a) shows the total intensity distribution of both the incident and transmitted beams, and  higher-order modes with an intensity minimum at the center. The orthogonal component originates from polarization conversion of plane waves that do not propagate in the symmetry planes of the discs. The conversion takes place independently of the initial polarization state of each of the waves and is caused by spatial dispersion. The total power converted to the orthogonal polarization comprises 7% of the incident-beam power. This new phenomenon would be impossible to reveal without our matrix approach taking into account the polarization conversion by spatial dispersion.
Our methods are also perfectly suitable for metamaterials exhibiting strong polarization conversion due to such well-known effects as optical chirality and anisotropy. For example, we have calculated that if the beam in the above example was polarized at an angle of 45°with respect to the disc-orientation direction-say, along the x-axis-the power conversion efficiency to the orthogonal y-polarization would be 9% as a result of optical anisotropy.

A silver-dimer metamaterial
In our second example, the metamaterial is composed of silver-dimer metamolecules arranged in a cubic lattice in glass. The dimers are composed of two discs that share a common axis and are separated by a gap. The discs have radii of 20 nm and 10 nm, and the gap is 40 nm wide. Both discs are 10 nm thick. The lattice constants are 175 The z-axis coincides with the normal of the exit surface of the slab, as in the previous example. The dimer axes are oriented normal to the slab surface, and such that the incident light faces the smaller disc first. A single layer of the metamaterial is illustrated in figure 5. These dimer metamolecules exhibit pronounced higher-order multipole resonances [39][40][41]. As a result, they can form a bifacial metamaterial [20], in which counterpropagating waves have different impedances.
Numerical calculation of the single-layer transmission and reflection matrices,τ andˆ, ρ reveals that the conversion between the s-and p-polarizations is absent for all planes and angles of incidence. Therefore, the matricesα ± andβ ± in (28) are diagonal in the basis of s-and p-polarizations, and we can use (30)-(32) to calculate the propagation constant, refractive index, and impedance for the s-and p-polarized eigenmodes. The results of these calculations are shown in figure 6. The polar plots of n and Z presented in this figure as functions of the incidence angle θ correspond to s-polarized waves; the curves for p-polarized waves are not shown, because they essentially overlap with the illustrated curves. The calculations were performed for λ = 525 nm. The refractive index is seen to be nearly equal to that of glass at all values of θ. It has a remarkably small imaginary part, especially at large angles of incidence, which implies that the absorption loss is low. Furthermore, the impedance of the material is close to that of glass at all values of θ larger than about 20°. This implies that plane waves incident onto the material from glass will be transmitted by the glass-material interface, rather than being reflected from it at these angles. At small angles of incidence, on the other hand, the wave impedance of the material significantly deviates from the impedance of glass and has a large imaginary part. As a result, the waves are expected to be efficiently reflected from the material.
The discussed spatially dispersive plane-wave reflection can be used to create a metamaterial-based spatial filter for optical beams. Assume that an optical beam is incident onto a slab of the described metamaterial. As the parameters n and Z have already been calculated, we can evaluate the Fresnel coefficients, τ and ρ, given by (33) and (34) and insert them into the following (Fabry-Perot) expressions for the transmission and reflection coefficients of the slab [20]   . Real (dotted lines) and imaginary (solid lines) parts of (a) refractive index n and (b) impedance Z as functions of incidence angle θ of plane waves at 525 λ = nm; n and Z of glass are shown with dashed lines. The impedance is normalized with respect to the impedance of glass. The imaginary part of the refractive index/impedance is multiplied/divided by 10 to make these quantities comparable with the real parts.

Conclusions
In this work, we introduce an approach to treating the interaction of optical beams with realistic, spatially dispersive metamaterials. The method allows us to describe the interaction of any continuous-wave beam with various crystalline optical metamaterials that can, for example, be optically anisotropic or chiral. For many metamaterials, spatial dispersion makes it impossible to introduce the polarization modes, which are the waves conserving the polarization state upon interaction with the material. In such cases, the wave parameters such as the refractive index and impedance, are also impossible to introduce. In order to still be able to describe the interaction of optical beams with such metamaterials, we have developed an approach based on transmission and reflection matrices that takes into account the possible polarization conversion of the beam's plane-wave components. For metamaterials with negligible evanescent-wave coupling between metamolecular layers, we have introduced an efficient semianalytical method to fully characterize the material in terms of the reflection and transmission matrices of a single metamolecular layer.
We have applied our methods to two particular examples of spatially dispersive metamaterials. The first example, dealing with highly symmetric metamolecules, underlines the importance of our new matrix-based approach that correctly describes all polarization-converting metamaterials. The second example illustrates the effect of spatial filtering imposed by a spatially dispersive metamaterial slab on a beam reflecting from the slab surface. The effect is completely insensitive to the coordinate of the focal spot of the incident beam, which cannot be achieved with conventional pinhole-based spatial filters.
The developments presented in this work take a significant step forward in describing the interaction of light with optical metamaterials by moving from plane waves to realistic optical beams and beam-like fields. Since the model presented in this paper is currently the only one that can properly treat this interaction, we anticipate that it will attract significant attention from the scientific community and lead to useful applications. In principle, it can also be adapted to radio-frequency beams and even acoustic beams interacting with anisotropic and spatially dispersive periodic structures. These expressions for matrices t ± and r ± can either be used to check the validity of (18)- (21) or directly inserted in (6) and (7) to calculate the fields of transmitted and reflected optical beams.