Optical wave parameters for spatially dispersive and anisotropic nanomaterials

Spatial dispersion is an intriguing property of essentially all nanostructured optical media. In particular, it makes optical waves with equal frequencies and polarizations have different wavelengths, if they propagate in different directions. This can offer new approaches to control light radiation and propagation. Spatially dispersive nanomaterials, such as metamaterials, are often treated in terms of wave parameters, such as refractive index and impedance retrieved from reflection and transmission coefficients of the material at each incidence angle. Usually, however, the waves are approximated as transverse, which simplifies the description, but yields wrong results, if spatial dispersion or optical anisotropy is significant. In this work, we present a method to calculate the wave parameters of a general spatially dispersive and optically anisotropic medium without applying such an approximation. The method allows one to evaluate the true impedances and field vectors of the effective waves, obtaining thus the true light intensity and energy propagation direction in the medium. The equations are applied to several examples of spatially dispersive and anisotropic materials. The method introduces new insights into optics of nanostructured media and extends the design of such media towards optical phenomena involving significant spatial dispersion. c © 2017 Optical Society of America OCIS codes: (160.4236) Nanomaterials; (310.6805) Theory and design; (160.3918) Metamaterials; (350.4238) Nanophotonics and photonic crystals; (310.6628) Subwavelength structures, nanostructures. References and links 1. B. Gompf, B. Krausz, B. Frank and M. Dressel, “k-dependent optics of nanostructures: Spatial dispersion of metallic nanorings and split-ring resonators,” Phys. Rev. B 86, 075462 (2012). 2. C. Menzel, C. Rockstuhl, T. Paul and F. Lederer, “Retrieving effective parameters for metamaterials at oblique incidence,” Phys. Rev. B 77, 195328 (2008). 3. C. Menzel, T. Paul, C. Rockstuhl, T. Pertsch, S. Tretyakov and F. Lederer, “Validity of effective material parameters for optical fishnet metamaterials,” Phys. Rev. B 81, 035320 (2010). 4. T. Paul, C. Menzel, W. Smigaj, C. Rockstuhl, P. Lalanne and F. Lederer, “Reflection and transmission of light at periodic layered metamaterial films,” Phys. Rev. B 84, 115142 (2011). 5. P. Grahn, A. Shevchenko and M. Kaivola, “Theoretical description of bifacial optical nanomaterials,” Opt. Express 21, 23471–23485 (2013). 6. P. Grahn, A. Shevchenko and M. Kaivola, “Interferometric description of optical metamaterials,” New J. Phys. 15(11), 113044 (2013). 7. A. Shevchenko, P. Grahn and M. Kaivola, “Internally twisted spatially dispersive optical metamaterials,” J. Nanophot. 8, 083074 (2014). 8. A. Shevchenko, P. Grahn and M. Kaivola, “Spatially dispersive functional optical metamaterials,’ J. Nanophoton. 9, 093097 (2015). 9. V. Kivijärvi, M. Nyman, A. Karrila, P. Grahn, A. Shevchenko and M. Kaivola, “Interaction of metamaterials with optical beams,” New J. Phys. 17, 063019 (2015). 10. V. Kivijärvi, M. Nyman, A. Shevchenko and M. Kaivola, “An optical metamaterial with simultaneously suppressed optical diffraction and surface reflection,” J. Opt. 18, 035103 (2016). 11. V. Kivijärvi, M. Nyman, A. Shevchenko and M. Kaivola, “Optical-image transfer through a diffractioncompensating metamaterial,” Opt. Express 24, 9806 (2016). 12. L. D. Landau and E. M. Lifshitz Electrodynamics of Continuous Media (Pergamon, 1984). 13. A. Sihvola, Electromagnetic Mixing Formulas and Applications (IEEE, 1999). 14. C. R. Simovski, “Bloch material parameters of magneto-dielectric metamaterials and the concept of Bloch lattices,” Metamaterials 1, 62 (2007). Vol. 25, No. 8 | 17 Apr 2017 | OPTICS EXPRESS 8550 #286144 https://doi.org/10.1364/OE.25.008550 Journal © 2017 Received 3 Feb 2017; revised 10 Mar 2017; accepted 21 Mar 2017; published 3 Apr 2017 15. A. Fang, T. Koschny and C. M. Soukoulis, “Optical anisotropic metamaterials: negative refraction and focusing,” Phys. Rev. B 79, 245127 (2009). 16. J. Yao, Z. Liu, Y. Liu, Y. Wang, C. Sun, G. Bartal, A. M. Stacy and X. Zhang, “Optical negative refraction in bulk metamaterials of nanowires,” Science 321, 930 (2008). 17. Y. Liu, G. Bartal and X. Zhang, “All-angle negative refraction and imaging in a bulk medium made of metallic nanowires in the visible region,” Opt. Express 16, 15439 (2008). 18. T. Paul, C. Rockstuhl, C. Menzel and F. Lederer, “Anomalous refraction, diffraction, and imaging in metamaterials,” Phys. Rev. B 79, 115430 (2009). 19. E. B. Graham and R. E. Raab, “Multipole solution for the macroscopic electromagnetic boundary conditions at a vacuum-dielectric interface,” Proc. R. Soc. Lond. A 456, 1193 (2000). 20. E. B. Graham and R. E. Raab, “The role of the macroscopic surface density of bound electric dipole moment in reflection,” Proc. R. Soc. Lond. A 457, 471 (2001). 21. A. Shevchenko, V. Kivijärvi, P. Grahn, M. Kaivola and K. Lindfors, “Bifacial metasurface with quadrupole optical response,” Phys. Rev. Appl. 4, 014019 (2015). 22. J. Stratton, Electromagnetic Theory (McGraw-Hill, 1941) 23. M. G. Silveirinha, “Poynting vector, heating rate, and stored energy in structured materials: A first-principle derivation,” Phys. Rev. B 80, 235120 (2009). 24. A. Poddubny, I. Iorsh, P. Belov and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photon. 7, 958 (2013). 25. M. Y. Koledintseva, S. K. R. Chandra, R. E. DuBroff and R. W. Schwartz, “Modeling of dielectric mixtures containing conducting inclusions with statistically distributed aspect ratio,” Progr. Electromagn. Res. 66, 213 (2006). 26. H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Tamamura, T. Sato, and S. Kawakami, “Self-collimating phenomena in photonic crystals,” Appl. Phys. Lett. 74, 1212 (1999). 27. Z. Lu, S. Shi, J. Murakowski, G. Schneider, C. Schuetz, and D. Prather, “Experimental demonstration of selfcollimation inside a three-dimensional photonic crystal,” Phys. Rev. Lett. 96, 173902 (2006). 28. T. Paul, C. Menzel, C. Rockstuhl, and F. Lederer, “Advanced optical metamaterials,” Adv. Mater. 22, 2354 (2010). 29. B. Casse, W. Lu, Y. Huang, E. Gultepe, L. Menon, and S. Sridhar, “Super-resolution imaging using threedimensional metamaterials nanolens,” Appl. Phys. Lett. 96, 023114 (2010). 30. L. Novotny and B. Hecht, Principles of nano-optics (Cambridge University, 2006).


Introduction
Essentially all optical nanomaterials are spatially dispersive, because their structural units are not negligibly small compared to optical wavelengths [1][2][3][4][5][6][7][8][9][10][11]. This makes the description of such materials in terms of electric permittivity and magnetic permeability tensors impractical, because optical plane waves propagating in the material in different directions experience different material parameters even if their polarizations are the same. This property, on the other hand, can offer new possibilities to control optical radiation. In general, spatial dispersion is caused by the excitation of higher-order multipole moments in the material's unit cells [5]. These moments differ at different propagation directions due to non-negligible phase delay of the wave across the unit cell and evanescent-wave coupling between the unit-cell parts. When spatial dispersion is significant, conventional effective medium approximations based on Clausius-Mossotti-or Maxwell-Garnett-type spatial averaging of dipole excitations in the medium [12,13] yield incorrect results [1,14,15]. Usually, spatially dispersive materials, such as metamaterials, are described in terms of plane-wave parameters, such as refractive index and impedance, that depend on both the polarization and the propagation direction of the wave [1][2][3][4][5][6][7][8][9][10][11]. The parameters should therefore be given as functions of the wave propagation angles. It is customary to divide the considered plane waves into transverse electric and transverse magnetic modes, for which either the electric or magnetic field vector is perpendicular to a certain symmetry plane of the field-matter system. However, when the wave parameters are retrieved from the transmission and reflection coefficients of a nanostructured material, it is common to assume that both the electric and magnetic fields in the material are transverse with respect to the wave vector [2,3,[5][6][7][8]. This approximation gives relatively accurate results, if the effects of spatial dispersion and optical anisotropy in the material are weak. Otherwise, the waves cannot be treated as transverse. As an example, in a transverse electromagnetic wave, the energy and wavefronts propagate in the same direction, while in reality these directions can differ significantly, such that a beam can show positive refraction in a material designed for negative refraction of wavefronts, and vice versa [16][17][18]. This difference is well known as a consequence of optical anisotropy, but not as a result of spatial dispersion. In fact, spatial dispersion can lead to a much larger deviation of the energy flow direction from that of the wave vector. In some cases, it can even make it impossible to introduce polarization modes and the associated wave parameters for certain propagation directions in the material [9].
In this paper, we derive effective wave parameters for spatially dispersive and optically anisotropic media in a systematic way without assuming that the waves are transverse in the electric and magnetic fields. We first introduce the general transmission and reflection coefficients at a boundary between two spatially dispersive and optically anisotropic media, assuming that the wave polarizations are preserved upon propagation. These coefficients are given in a polarization-independent form in terms of tangential wave impedances. Then, using these coefficients, we show how the effective refractive index, the directions of the electric and magnetic field vectors, and the true wave impedance can be retrieved from the transmission and reflection coefficients of a spatially dispersive material slab. Finally we apply the model to several examples of spatially dispersive and optically anisotropic media and show that the non-transverse character of waves in these media is important. In particular, we show that the average power flow in such a material is perpendicular to the evaluated effective electric and magnetic field vectors, the true impedance determines the intensity of optical beams and consequently also the radiation efficiency of optical sources inside such media, while the tangential impedances determine the power reflectance and transmittance of optical beams at material surfaces. These optical effects can be used to design novel optical materials in which spatial dispersion plays a significant role.

Tangential reflection and transmission coefficients for spatially dispersive and/or optically anisotropic media
Let us consider two spatially dispersive and/or optically anisotropic media and an optical plane wave incident onto the boundary between them. We assume that the incident, reflected and transmitted waves are the polarization modes in the materials and, therefore, their polarizations do not change upon propagation and a refractive index and impedance can be introduced for them. As usually, two orthogonally polarized modes can be considered. We call them transverse electric (TE) and transverse magnetic (TM) modes, but allow the electric and magnetic fields, E and H, to be non-transverse with respect to the wave vector k. In fact, the electric displacement D and magnetic induction B, which according to Maxwell's equations are still perpendicular to k, can replace E and H in the definition of the TE and TM modes. The TE and TM polarizations are illustrated in Fig. 1. Note that the incidence and reflection angles are not necessarily equal even for an ordinary anisotropic medium.
To derive the reflection and transmission coefficients, we require that the tangential components E | | and H | | of E and H are continuous across the interface. This boundary condition is valid for most optical materials, excluding however such specific media in which certain electric quadrupole and higher-order multipoles can be excited to generate lower-order multipoles on the material surface [19,20]. We choose the directions of E and H such that for a zero phase shift, the sign of E | | does not change, but the sign of H | | changes upon reflection, as in Fig. 1. In this case, the following common continuity condition can be written for both TE and TM waves: Introducing the tangential reflection and transmission coefficients ρ E | | and τ E | | for the electric  Fig. 1. Reflection and transmission of (a) a TE and (b) a TM wave at a boundary between two spatially dispersive media. The electric and magnetic field vectors are not necessarily transverse with respect to the wave vector (shown by the black arrows) and the incidence and reflection angles can differ.
field and ρ H | | and τ H | | for the magnetic field as and a tangential impedance as we solve Eqs. (1) and (2) to obtain These four equations hold for both TE and TM waves. It can be seen that, if for a certain incidence angle the tangential impedances in the two materials are equal, the interface does not reflect the wave.

Wave parameters from tangential reflection and transmission coefficients of a slab
We can now express the effective wave parameters in terms of the transmission and reflection coefficients of a slab of a spatially dispersive material. Let the slab be surrounded by an ordinary isotropic dielectric with a refractive index n s and a wave impedance η s = η 0 /n s , where η 0 is the wave impedance in vacuum (see Fig. 2, where medium 3 is the same as medium 1). The slab thickness is denoted by d. A plane wave incident on the slab at an arbitrary angle θ i forms two plane waves inside the slab, one propagating in the forward direction (see the blue arrow in the figure) and another in the backward direction (the red arrow). In general, these waves can experience different refractive indices and impedances [5,7]. We denote these parameters by n f and η f for the forward propagating wave and n b and η b for the backward propagating wave. The same pair of waves is produced in the slab, if the incident wave enters the slab at an angle of π − θ i (from the opposite side, as shown in Fig. 2). Since medium 1 is an ordinary isotropic dielectric, such as glass, the angles of incidence and reflection from the slab are equal. In addition, since media 1 and 3 are the same, the propagation angles of the incident and transmitted waves are also the same. Therefore, the slab reflection coefficient r and transmission coefficient t are equal to the corresponding tangential coefficients, i.e., This means that r and t can be written in terms of ρ | | and τ | | given in Eqs.
(3)- (6). All the necessary tangential reflection and transmission coefficients at the slab surfaces, ρ i j and τ i j , are shown in Fig. 2. They describe plane waves propagating from medium i to medium j. Furthermore, independently of the wave polarization, ρ i j and τ i j are given by Eqs. (8) and (9) for the electric component and by Eqs. (10) and (11) for the magnetic component of the field.
Reflection and transmission of plane waves by a slab of a spatially dispersive material. The arrows show the directions of the wave vectors. Parameters ρ i j and τ i j are the tangential reflection and transmission coefficients that for the electric field are given by Eqs. (3) and (4) and for the magnetic field by Eqs. (5) and (6). They describe the tangential components of the waves propagating from medium i to medium j at each surface of the slab independently of the polarization of the incidence wave.
Taking into account multiple reflections of plane waves inside the slab one can derive the following four equations for the slab transmission and reflection coefficients t and r describing the two incident waves depicted in Fig. 2: Here, k fz and k bz are the z-components of the wave vectors for the forward and backward propagating waves inside the material. They depend on the refractive indices n f and n b . The tangential component of the wave vector, k | | =xk x +ŷk y , is the same inside the material as outside of it, as required by the boundary conditions, and its length, k | | = n s k 0 sin θ i , is known.
Here k 0 is the wavenumber in vacuum. The parameters ρ i j and τ i j contain only two unknowns, η f | | and η b| | , as follows from Eqs. (3)- (6). Hence, the four equations (12)- (15) can be solved for the four unknowns that are k fz , η f | | , k bz and η b| | . Assuming that ρ i j and τ i j are the electric-field coefficients given by Eqs. (8) and (9), we find the solutions where ξ f | | = η f | | /η s| | and ξ b| | = η b| | /η s| | are the normalized tangential impedances, m is a properly chosen integer [2,5,7] (in particular such that the imaginary parts of k fz and k bz and the real parts of ξ f | | and ξ b| | are positive), and the parameters a and b are Equations (18) and (19) can also be obtained from Eqs. (16) and (17) by interchanging the subindices "f" and "b". The refractive indices n f and n b for the waves propagating in the material to the right (in the direction of the blue arrow in Fig. 2) and to the left (along the red arrow) are found from where the sign is chosen for the refractive index to have positive imaginary part. Note that, if the material is centrosymmetric and has a symmetry axis along the surface normal, we have t f = t b and r f = r b . For non-centrosymmetric materials, the latter equality does not hold even at normal incidence [5,7,21]. Similar equations have been derived previously by considering interference of plane waves reflected and transmitted by monomolecular layers in a metamaterial for TEand TM-polarized waves assumed to be transverse in both E and H [5,7]. The new Eqs. (16)-(23) allow for a non-transverse character of the waves and are of the form that is independent of their polarization. Using the above equations, the refractive index n(θ, φ) can be found for all possible propagation directions, with θ taken with respect to the z-axis and φ with respect to the xzplane. These polar and azimuthal propagation angles are obtained from the complex Snell's law [22]. For example, for the wave propagating to the right in Fig. 2, we have where n f = n(θ, 0). Note that also ξ f | | = ξ | | (θ, 0) given by Eq. (17) can be calculated for any propagation angle θ and used, e.g., to obtain the interface reflection and transmission coefficients in Eqs. (8)- (11). A surface r = Re{n(θ, φ)} plotted for a fixed frequency in spherical coordinates (r, θ, φ) forms an isofrequency surface for the real refractive index. In general, this surface depends on the wave polarization. Note that optical reciprocity requires that n(θ, φ) = n(π + θ, π + φ), owing to which it is enough to calculate n(θ, φ) only for waves with θ ∈ [−π/2, π/2], i.e, propagating in the forward direction.
For an absorbing material, n f is complex-valued, and therefore, also sin θ in Eq. (24) is complex. However, for optical nanomaterials of practical interest, the absorption must be weak, owing to which the right-hand side of Eq. (24) can be replaced with its real part, yielding a real-valued θ. The calculations that follow can be considered accurate, if the attenuation length of the field is long compared with the wavelength, in which case Re{n(θ, φ)} alone can be used to determine the directions of E and H. This condition is satisfied by a wide range of lossy nanomaterials, as will be discussed at the end of this section.
We are now about to calculate the field vectors E and H inside the material and the total impedances that similarly to n(θ, φ) can differ for different polarization modes. Let us first consider the transverse electric (TE) plane waves, with E perpendicular to the xz-plane in Fig.  2, assuming that they are the material's polarization modes. If neither H nor E is perpendicular to the plane of incidence, the wave polarization changes upon propagation and the wave parameters cannot be defined [9]. The magnetic field H must lie in the xz-plane and be parallel to the isofrequency contour [12] that is the intersection of the isofrequency surface with the

TE TM
where n = Re{n(θ)} = Re{n f }. Note that H z is negative for the positive α shown in Fig. 3a.
The angle α is also the angle of energy flow with respect to the z-axis as long as this flow can be described by the Poynting vector S = Re{E × H * }/2 [23]. The total impedance η TE = E/H is therefore η TE = ξ | | η s| | | cos α| = ξ | | η s| | |n cos θ + (∂n /∂θ) sin θ| where and the modulus is used to keep the real part of impedance positive. Equations (25) and (26) give the direction of H and the total impedance for the TE case. The magnitude of H is then obtained as H = E/η TE . In the case of TM polarization, the magnetic field is perpendicular to the xz-plane and the electric field makes an angle β with the x-axis, as shown in Fig 3b. The angle β is found from where n(θ) = n f is calculated for the TM polarization. The total impedance is then given by where The magnitude of E is obtained as E = Hη TM . Now we have a complete set of equations for retrieval of the wave parameters and the field vectors E and H of polarization modes in a spatially dispersive and optically anisotropic material, in which optical absorption is not very high. A more elaborate condition for the allowable level of absorption can be derived by using frequency-domain Maxwell's equations k × E = ωB and k × H = −ωD and the fact that the waves attenuate exclusively in the z-direction. One can show that due to absorption the possible angular error Δγ for the direction of the Poynting vector determined by Eqs. (25) and (28) depends on its orientation angle γ = {α, β} and is approximately given by where γ is equal to α for TE polarization and β for TM polarization (see Figs. 3a and 3b). Hence, if Im{n} is small or the direction of the Poynting vector is close to the z-direction, the error is negligible. For relatively large Im{n} and large angles γ, the error is larger, but it can still be small even for highly absorbing materials. As an example, if n = 1.5+0.16i at λ = 1 μm, the material is essentially opaque, having a field attenuation length of 1 μm. However, the error Δγ would in this case be only 2.8 • at γ = 30 • and 4.8 • at γ = 60 • . Finally we emphasize that if for a certain propagation angle the polarization modes do not exist, the orientations of E and H are not defined and the wave parameters cannot be introduced. In such cases, other methods must be used to characterize the optical response of the medium [9].

Examples
Let us first assume that the isofrequency surface of the material is known, being given by either an ellipsoid or a hyperboloid. The assumption of transverse electromagnetic waves used in [2,5,7] would lead to an approximate impedance when applied to such materials. The correct impedance can be obtained by multiplying the approximate one with a factor C(θ) = | cos θ/ cos β| = | cos θ| n 2 + (∂n/∂θ) 2 |n cos θ + (∂n/∂θ) sin θ| (32) for the TM-polarized waves, as follows from Eq. (29). For the TE waves the contour would be circular in the absence of spatial dispersion. Assuming that the z-axis is a symmetry axis of the material, we can write the function n(θ), in the xz-plane, as where "+" stands for an elliptic material and "−" for a hyperbolic material of a type with a twosheeted isofrequency hyperboloid [24]. If for the hyperbolic material the effective ordinary and extraordinary refractive indices n o and n e are chosen to be the same as for the elliptic one, the factor C(θ) is the same for the two materials. Figure 4a shows a polar plot of C(θ) calculated for n o = 1.5 and n e = 2. It can be seen that at θ = 0 the correction factor is equal to 1 (because the wave is in this case transverse) and increases monotonically with θ. Note that in the presence of spatial dispersion the impedance is not determined solely by the refractive index and can be rather arbitrary.
As another example, we consider a strongly spatially dispersive material designed to compensate for optical diffraction. It can be a photonic crystal [26,27] or a metamaterial [10,11,28,29]. Ideally, the refractive index of such a medium should be given by if the divergence-free propagation of light is chosen to be in the z-direction [10]. Here, n 0 is the refractive index at θ = 0. The factor C(θ) calculated for TM-polarized waves in this material is shown in Fig. 4b. It is clear that the waves must be considered as non-transverse in both cases of Fig. 4, as the mistake (shown by the red line) is significant already at small propagation angles. The final example treats a realistic metamaterial design that shows significant spatial dispersion and optical anisotropy. The material is composed of rectangular silver nanorods periodically distributed in glass (n glass = 1.5) with lattice periods of Λ x = Λ y = 130 nm and Λ z = 200 nm. The rods, 130 nm long and 30 nm thick, are tilted in the xz-plane by an angle of 45 • with respect to the z-axis. The structure can be seen in Fig. 5. Applying the derived equations to a slab of this material surrounded by glass, we have calculated the isofrequency contours of the refractive index and impedance in the xz-plane for TM-polarized waves at a vacuum wavelength λ 0 = 1.1 μm. The slab reflection and transmission coefficients have been calculated by using a commercial software COMSOL Multiphysics. The obtained isofrequency contours are illustrated in Figs. 6a and 6b. The blue contours show the real parts of the quantities and the red ones the imaginary parts multiplied by a factor of 100 to make them visible. When the imaginary part of η is negative, it is shown by the red dashed line. The impedance of glass is shown by the brown dashed line. Both impedances are normalized to the impedance of vacuum. The gray sectors indicate propagation angles θ that are not accessible for waves incident from the surrounding glass. It can be seen that at θ ≈ −45 • , i.e., for k parallel to the rods, the refractive index and impedance of the material are close to those of glass, which leads to reduced reflection at the glass-metamaterial interface. The real part of n is quite flat in this range, which leads to reduced optical diffraction. Due to spatial dispersion, the shapes of the isofrequency contours of Re{n(θ)} and Re{η(θ)} are not elliptical and cannot be determined from each other, in contrast to isofrequency contours of purely anisotropic materials. As an example, we apply an [13] and [25] (valid for anisotropic materials without spatial dispersion) to our structure and, by evaluating the elements i of the effective permittivity tensor, calculate the approximate refractive index and impedance. Within this approximation, the material is a uniaxial crystal, for which the principal values of the refractive index are given by n i = √ i / 0 and the relative impedance is η = 1/n. The resulting isofrequency contours are shown in Figs. 6c and 6d. Near the propagation angle of −45 • the longitudinal excitations of the rods are negligible and the approximate values of n and η are close to those of Figs. 6a and 6b. Far from this angle, the discrepancy is large and the approximation is not valid. We have verified that the discrepancy disappears only when the unit cells become about 10 times smaller than the wavelength, as a result of which spatial dispersion vanishes.
We first consider a wave propagating in the material at θ = 0. According to Eq. (28) and the retrieved Re{n(θ)}, the angle of energy transfer in the wave should be β = −25 • . To verify this result, we have numerically simulated, using COMSOL Multiphysics, the propagation of a wide optical beam with the wavefronts parallel to the xy-plane in the material. The beam radius in the x-direction is chosen to be 1 μm. The evaluated intensity, i.e., the magnitude of the time-averaged Poynting vector in the considered beam is shown in Fig. 7a. The beam propagates upwards, towards a horizontal boundary with glass. The blue dashed line shows the predicted energy transfer direction. In the material, it makes an angle of 25 • with the z-axis. It coincides perfectly with the axis of the numerically simulated beam, implying that the derived equations are correct. The phase distribution of the beam is shown in Fig. 7b. The wavefronts of the effective wave are clearly parallel to the xy-plane. From the phase distribution, one can estimate the effective wavelength to be about 560 nm that approximately corresponds to the retrieved effective refractive index of 2 at θ = 0. The beam crosses the interface with glass and propagates in the vertical direction, as it should.
In a second example, we generate an optical beam inside a metamaterial by introducing a continuous horizontal distribution of x-directed electric dipoles between two neighboring "metamolecular" layers of the material of the previous example. The horizontal black line in Figs. 8a and 8b marks the plane where the dipoles are positioned. The dipole amplitude has a Gaussian distribution along the x-direction with a 1/e radius of 1 μm. The phase of the dipoles is chosen to be proportional to k x x with k x being the wave-vector component of an effective wave propagating upwards at an angle θ = 19 • and a wave propagating downwards at an angle θ = 180 • − 27 • . The light source must therefore radiate two beams with wavefronts propagating   at these two angles. The energy transfer directions in the beams must be given by β that according to Eq. (28) must have the values of 0 and 180 • − 42 • for the upwards and downwards propagating beams, respectively. The intensity distributions of these two beams are shown in Fig. 8a. The slab of the material has 6 metamolecular layers and is surrounded by glass. According to our analytical calculations, the beams crossing the slab boundaries should propagate in glass at angles of 30 • and 180 • − 30 • with respect to the z-axis. In Fig. 8a, the predicted directions of light energy propagation in the beams are shown by the blue dashed lines. Clearly, the direct numerical calculations fully confirm the theory. The phase distributions of the beams are shown in Fig. 8b. While the wavefronts are not perfectly flat in the material, one can see that they propagate approximately at the predicted angles. The retrieved refractive indices for the upwards and downwards propagating waves are 2.30 and 1.65, respectively. The power of the beam radiated downwards is approximately the same as the power of the upwards propagating beam, but the average intensity in it is higher by a factor of 1.4. This factor is close to the ratio of the impedances of the two waves that is η(180 • − 42 • )/η(19 • ) ≈ 1.38. Hence, while the tangential impedance determines the overall power of a beam transmitted or reflected by the material surface, the total impedance given in Eqs. (26) and (29) determines the wave intensity inside the material and therefore reflects the density of states and possible enhancement or suppression of fluorescence into the modes of the material [30].

Conclusion
A set of equations for retrieval of wave parameters for spatially dispersive and optically anisotropic media is derived and verified numerically. The equations are valid for polarization modes in relatively weakly absorbing materials, such that the attenuation length of the mode is large compared to the wavelength. The theory predicts the refractive index and impedance for each propagation direction and polarization of the wave in the material, yielding also the directions of the effective electric and magnetic field vectors, average intensity, and the direction of energy transfer. The derived equations have been verified by several examples. The theory can be used to develop novel optical nanomaterials in which spatial dispersion plays a significant role. For example, a material can be designed to provide diffraction-free self-guidance of optical beams or make radiation of dipole sources highly directional and as such much more efficient.

Funding
Vilho, Yrjö and Kalle Väisälä Foundation of the Finnish Academy of Science.