First Principles Homogenization of Periodic Metamaterials and Application to Wire Media

Here, we present an overview of a first principles homogenization theory of periodic metamaterials. It is shown that in a rather general context it is possible to formally introduce effective parameters that describe the time evolution of macroscopic (slowly-varying in space) initial states of the electromagnetic field using an effective medium formalism. The theory is applied to different types of"wire metamaterials"characterized by a strong spatial dispersion in the long wavelength limit. It is highlighted that the spatial dispersion may tailor in unique ways the wave phenomena in wire metamaterials leading to exotic tunneling effects and broadband lossless anomalous dispersion.


Introduction
The interactions between waves and matter play a fundamental role in most physical processes. It is usually rather challenging to characterize exactly the wave propagation in macroscopic systems formed by a large number of identical elements, e.g., in periodic or random composite materials, due to the complexity of the wave phenomena at the microscopic level. Fortunately, in many instances, the detailed microscopic behavior of a wave is of very limited practical interest. Instead, one can resort to effective medium theories that provide a simplified description of the wave phenomena in terms of a limited set of parameters. Effective medium theories are particularly successful when the wavelength is large with respect to the characteristic spatial period of the composite material. In this case, the material may be regarded as a continuum, and the homogenization formalism gives a simplified and insightful picture of the wave propagation.
Effective medium theories have a long history [1]. In the case of light waves, the concepts of "permittivity" and "permeability" of a material are as old as the electromagnetism itself. Similarly, in semiconductor theory the effects of a periodic electrostatic potential associated with the ionic lattice can be modeled by an effective electron mass [2]. In the last two decades, the interest in effective medium theories has been renewed by the emergence of the field of metamaterials . Metamaterials are composite media formed by properly shaped dielectric or metallic inclusions embedded in a host medium, which are designed to exhibit extraordinary behavior such as a negative index of refraction [29], subwavelength imaging [30,31] or other applications [32]. Usually, in metamaterials the radiation wavelength λ is only moderately larger than the lattice constant a, typically 5-10 times. This contrasts with natural media where the ratio, λ/a, is several orders of magnitude larger than that value, even at optical frequencies. This property imposes restrictions on the application of classical homogenization theories to artificial materials [5,7,18,19] due to the emergence of spatial dispersion.
In a spatially dispersive material the electric displacement vector in a given point of space cannot be written exclusively in terms of the macroscopic electric field at the same point, but ultimately may depend on the distribution of the electric field in a neighborhood that encompasses many unit cells [33]. This non-locality of the electromagnetic response has many important and nontrivial repercussions on the physical properties of a material [34].
The objective of this review article is to present an up to date comprehensive description of a general homogenization procedure first developed in the context of electromagnetic metamaterials [7] and later generalized to semiconductor superlattices [35]. The effective medium theory is applicable to a wide range of periodic physical systems and takes into account both spatial and frequency dispersion [10,36]. We illustrate the application of the formalism to "wire media". This class of metamaterials is particularly interesting, not only because it allows for an analytic treatment that describes almost exactly the actual microscopic response of the metamaterial, but also because of the richness of the wave phenomena it enables.
The review article is organized as follows: in Section 2 we describe the general homogenization scheme of Ref. [35] that uses as a starting point a time-domain perspective. In Section 3, we focus our analysis on nonmagnetic and periodic electromagnetic metamaterials and explain how to find the effective response in the frequency domain. The homogenization approach is applied to wire media in Section 4. The nonlocal effective models for different wire medium topologies are presented in Section 4.1. In Section 4.2, it is shown that the nonlocality of wire metamaterials emerges naturally from a quasi-static model with additional state variables that describe the internal degrees of freedom of the metamaterial. Some subtleties arising from the nonlocality of the electromagnetic response, such as the definition of the Poynting vector and the need for Additional Boundary Conditions (ABCs) are discussed in Sections 4.3 and 4.4, respectively. Finally, in Section 5 we describe some exotic wave phenomena due to the spatial dispersion in two distinct wire medium configurations.

Effective medium theory
In this section, we present the fundamentals of the homogenization method originally developed in Refs. [7,35,36]. We adopt the general perspective of Refs. [35,36] where the effective medium parameters are defined in such a way that they describe exactly the time-evolution of any macroscopic (slowly-varying in space) initial wave packet.

Microscopic theory
We consider a generic periodic in space physical system whose dynamics is characterized by a one-body Schrödinger-type equation of the form: Here,Ĥ is the operator that determines the time evolution of the system and ψ is the statevector that describes the state of the system. In general ψ is a multi-component vector (a spinor). Evidently, this type of formulation is suitable to characterize the propagation of electron waves in a bulk semiconductor or in semiconductor or graphene superlattices, and in such a contextĤ is the system Hamiltonian, ψ is the wave function and ħ is the reduced Planck constant [35,[37][38][39]. Importantly, the propagation of light can also be described using a similar formulation. Indeed, the Maxwell's equations can be written in a compact form as [35] where f = (e h) T is a six-element vector with components determined by the microscopic electric and magnetic fields and g = (d b) T is a six-element vector with components determined by the electric displacement and the magnetic induction fields. In electromagnetic metamaterials the f and g fields are related by a space-dependent material matrix M = M(r) through the constitutive relation g = M · f. In conventional isotropic media the material matrix is simply: where ε and µ are the permittivity and permeability, respectively. Hence, by definingĤ as: and identifying the state vector with the g field, ψ = g, the Maxwell's equations can be expressed as in (1). It should be noted that in the electromagnetic caseĤ is unrelated to the energy of the system, and should be simply regarded as an operator that describes the time evolution of the classical electromagnetic field. Moreover, in the previous discussion it is implicit that the relevant materials are nondispersive, i.e., the permittivity ε and the permeability µ are frequency independent. Yet, the formalism can be generalized to dispersive media, as it is always possible to get rid of the material dispersion with additional variables [40][41][42]. For lossy media, theĤ operator is non-Hermitian.

Spatial averaging and the envelope function
The envelope function is intuitively the slowly varying part, in space, of the state vector ψ. It is defined here as: where { } av is a linear operator that performs a spatial averaging. The averaging operator is completely determined by the response to plane waves, characterized by the function F (k) such that Thus, the action of the averaging operator on a generic plane wave with wave vector k yields another plane wave with the same wave vector, but with a different amplitude given by F (k). Because of the linearity of the operator { } av , its action on a generic function is determined by Fourier theory and is given by a spatial convolution. The envelope function can be explicitly written as: where N is the space dimension (e.g., N = 3 for a three-dimensional metamaterial). The weight function f is the inverse Fourier transform of F so that: Related ideas have been developed by Russakov in the context of macroscopic electromagnetism [43]. It is assumed that the averaging operator corresponds to an ideal low pass spatial filter such that: In this article the set B.Z. stands for the first Brillouin zone of the periodic lattice, but sometimes other choices can be relevant [36]. The envelope function Ψ(r, t ) has no significant spatial fluctuations on the scale of a unit cell, i.e., the microscopic fluctuations are filtered out by the averaging operator. Hence, Ψ(r, t ) determines the macroscopic state vector. In general, we say that a given state vector ψ is macroscopic when it stays invariant under the operation of spatial averaging: ψ(r) = ψ(r) av , (macroscopic state vector).
Importantly, a macroscopic state cannot be more localized in space than the characteristic period of the material.

The effective Hamiltonian
The effective Hamiltonian is the operator that describes the time evolution of the envelope function. Specifically, suppose that the initial state vector is macroscopic, so that ψ t =0 = Ψ t =0 .
In general, the time evolution of an initial macroscopic state does not yield a macroscopic state at a later time instant, i.e., ψ(r, t ) = Ψ(r, t ) for t > 0. We define the effective HamiltonianĤ ef such that Ψ(r, t ) calculated usingĤ ef is coincident with the spatially-averaged microscopic state vector {ψ(r, t )} av , where ψ(r, t ) is determined by the microscopic HamiltonianĤ [35,37]. These ideas are illustrated in the diagram of Figure 1. The time evolution of the macroscopic state vector is determined by a generalized Schrödinger equation: From the definition of the effective Hamiltonian it is clear that it must ensure that: Because of linearity, the action of the effective Hamiltonian on the wave function can be expressed as a convolution in space and in time [35]: Note that the kernel h ef is a function of r − r . We shall see below that this is possible because the spatial averaging operation essentially eliminates the spatial granularity of the system. In general, the kernel h ef (r, t ) is represented by a square matrix [h σ,σ ] because Ψ is a multi-component vector. In the photonic case the dimension of h ef is S = 6. Equation (13) shows that the effective Hamiltonian depends on the past history (0 < t < t ) and on the surroundings (r = r) of the observation point, rather than just on the instantaneous and local value of Ψ. It is convenient to introduce the Fourier transform of h ef (r, t ) defined as: The Fourier transform is bilateral in space and unilateral in time. The unilateral Fourier transform in time can also be regarded as a Laplace transform. In the Fourier domain, the action of the effective Hamiltonian reduces to a simple multiplication: Here, Ψ(k, ω) is the Fourier transform of the macroscopic state vector, and (Ĥ ef Ψ)(k, ω) is defined similarly. The integral in r is over all space. Note that for now the system is assumed to be periodic and unbounded, so that the effect of boundaries is disregarded. The convergence of the unilateral Fourier transform is guaranteed in the upper-half frequency plane, Im(ω) > 0.
The function H ef (k, ω) completely determines the effective Hamiltonian. Because of the properties of the spatial averaging operator, it is possible to enforce that: This property ensures that the effective Hamiltonian is a smoothened version of the microscopic Hamiltonian. In the following subsections, it is explained how H ef (k, ω) can be calculated for k ∈ B.Z.

Calculation of H ef (k, ω) with a time domain approach
Let us consider an initial macroscopic state of the form ψ t =0 ∼ e ik·r u l where the wavevector k can take any value in the B.Z. Here, (u l ) represents a basis of unit vectors that generates the S-dimensional vector space wherein ψ is defined. Because of the periodicity of the system, the microscopic time evolution of this initial state yields a state vector ψ(r, t ) with the Bloch property. In fact, ψ(r, t )e −ik·r is a periodic function in space for any fixed t . For the same reason,Ĥ ψ has also the Bloch property. Crucially, the operation of spatial averaging only retains spatial harmonics with wave vector inside the B.Z., and hence it follows that the dependence of {ψ} av and {Ĥ ψ} av on the spatial coordinates is of the form e ik·r for any time instant. In other words, within the effective medium approach the time evolution of a plane wave-type initial state yields another plane wave-type state, such that the homogenized structure behaves as a continuum. Moreover, for Bloch modes it is possible to write: with where Ω represents the unit cell and V cell is the respective volume. Taking now into account that Ψ = {ψ} av andĤ ef Ψ = {Ĥ ψ} av , and substituting (18) into (13), it is seen after straightforward manipulations that: In the above, ψ av (ω) and (Ĥ ψ) av (ω) stand for the unilateral Fourier (Laplace) transforms of the functions in (19). Hence, if we denote ψ (l ) , l = 1, . . . , S as the microscopic state vector determined by the time evolution of the initial state ψ (l ) t =0 = i/ħ e ik·r u l (the proportionality constant was fixed as i/ħ for convenience), it follows from the previous analysis that the effective Hamiltonian is given by: Here H ef and the two objects delimited by the square brackets are S ×S matrices. Thus, H ef (k, ω) can be written as the product of two matrices, whose columns are determined by the vectors ψ (l ) av (ω) and (Ĥ ψ (l ) ) av (ω).
In summary, for an arbitrary k ∈ B.Z. the effective Hamiltonian can be found by solving S microscopic time evolution problems associated with initial states of the form ψ (l ) t =0 = i/ħ e i k·r u l . The effective Hamiltonian is written in terms of the Fourier transforms in time of the functions (19).

Calculation of H ef (k, ω) with a frequency domain approach
The effective Hamiltonian may also be determined based on frequency domain calculations. To prove this we note that ψ av (ω) and (Ĥ ψ) av (ω) can be written explicitly as: where ψ(r, ω) is the unilateral Fourier transform of ψ(r, t ). Applying the unilateral Fourier (Laplace) transform to both members of the microscopic Schrödinger equation (1) and using the property , it follows that: Hence, ψ (l ) (r, ω) can be directly found by solving the above equation for −iħψ (l ) t =0 = e ik·r u l , with l = 1, . . . , S. Once ψ (l ) (r, ω) is known one can determine ψ (l ) av and (Ĥ ψ (l ) ) av using (22), and finally obtain the effective Hamiltonian from (21).
It is interesting to note that for −iħψ (l ) t =0 = e ik·r u l equation (23) implies that (Ĥ ψ (l ) ) av −ħωψ (l ) av = u l . Substituting this result into (21) one may also write the effective Hamiltonian as:

Stationary states
The spectrum of the effective Hamiltonian is exactly coincident with the spectrum of the microscopic Hamiltonian [35] (here, for simplicity it is assumed that there are no "dark states", for a discussion see [35]). The energy spectrum of the macroscopic Hamiltonian is determined by the nontrivial solutions of the stationary Schrödinger equation where E stands for the energy of a certain stationary state. For example, in the electromagnetic case the photonic band structure calculated with the effective Hamiltonian is coincident with the exact band structure obtained using a microscopic theory [7]. The enunciated result follows from the fact that in a time evolution problem (with no source excitation) the state vector can be written as a superposition of eigenmodes. The eigenmodes have a time variation of the form e −iω n t , being ω n = E n /ħ the relevant eigenfrequencies. Importantly, since the macroscopic and microscopic state vectors are related by the spatial-averaging operation (Ψ = {ψ} av ), both Ψ and ψ have the same-type of time oscillations. In other words, the averaging affects only the space coordinates, while the time coordinate is not averaged in any manner. As a consequence, the spectrum of the microscopic and macroscopic Hamiltonians must be the same. For a detailed mathematical proof of this property the reader is referred to Appendix C of Ref. [35].

The electromagnetic case
The formalism of the previous section when applied to the electromagnetic case (2) yields a 6 × 6 effective Hamiltonian of the form [36]: where M ef (k, ω) is the effective material matrix that links the averaged fields {f} av and {g} av of (2), [35,36]. For metamaterials made of non-magnetic particles the material matrix is of the form Thus, the homogenization problem reduces to the determination of the nonlocal effective permittivity ε ef (k, ω). The permittivity can be found using the source-driven homogenization theory developed for electromagnetic metamaterials [7]. As shown in [36] the effective response obtained with this theory is exactly coincident with the one obtained with the general theory of previous section. Below, we quickly review the main ideas of the source-driven homogenization, highlighting that the homogenization problem can be reduced to an integral equation [7]. We consider a generic nonmagnetic periodic metamaterial described by the periodic permittivity ε r (r, ω) = ε r (r + R, ω) with R a vector of the Bravais lattice. Assuming a time variation of the form e −iωt , the microscopic Maxwell equations in this system are where e, b are the microscopic electric and magnetic field, respectively and j e is an applied (macroscopic) electric current density that acts as a source of the electromagnetic fields. The applied current density is assumed to have the Bloch property and enforces a desired spatial variation within the unit cell. This means that the pair of parameters (ω, k) characterizing the time and space variations of the fields are independent of each other and do not need to be associated with an eigenmode. The applied current plays the same role as the initial state ψ t =0 in the formulation of last section. By applying the averaging operator (18a) to the microscopic Maxwell equations (28), one obtains the macroscopic Maxwell equations: where E av , B av and J e,av are the averaged e, b and j e , respectively, defined according to (22a). The averaged induced polarization P g is given by For system containing perfectly electric conducting (PEC) surfaces, the integration over the unit cell volume in the previous expression can be transformed into a surface integral, see [7,10] for more details. The nonlocal effective permittivity is defined through the relation between the averaged electric field and the averaged induced polarization: As shown in [7], for every pair (ω, k) the homogenization problem can be reduced to an integral equation. The unknown of the integral equation is the microscopic vector field p ind (r) = (ε r (r) − 1)e(r) and the excitation is the averaged electric field E av . A solution of the problem can be formally constructed using the Method of Moments (MoM). The unknown p ind is expanded as where the set of expansion functions w n,k has the Bloch property and is assumed to be a complete set in {r : ε r (r) − 1 = 0}.
For simplicity, next we focus on the case where the metamaterial inclusions can be modeled as impedance boundaries, characterized by some surface impedance Z s [44]. The surface impedance links the tangential electric field E tan at the boundary surface (∂D) with the current surface density J s =ν × H, as E tan = Z s J s [45]. Here,ν is the unit normal vector oriented toward the exterior of the inclusion. A PEC inclusion is described by the surface impedance Z s = 0. It can be shown that the effective permittivity is given by [7,10] In the above Φ p0 is the Green's function introduced in (35b) of [7], ∇ s stands for the surface divergence of a tangential vector field and the matrix [χ m,n ] is the inverse of [χ m,n ]. In the next section, we illustrate the application of the above formulas to the case of wire metamaterials. As shown in [7,10], Equation (33) can be generalized to the case of volumetric dielectric inclusions. The MoM formulation is particularly well suited to characterize the effective response of metamaterials made of metallic structures. Due to this reason, for dielectrics it is typically more practical to solve the homogenization problem with finite differences methods in the frequency [11] or in the time domain [14].

Application to wire metamaterials
Next, we apply the homogenization method to periodic arrays of thin metallic wires. Wire metamaterials are generically characterized by a strong spatial dispersion in the long wavelength limit.

Nonlocal effective models
In the following subsections we obtain the effective medium responses of three different wire metamaterials: the uniaxial wire medium, the double wire medium and the 3D connected wire mesh. In all cases, it will be assumed that the metallic wires are thin, R a, where R is the radius of the wires and a is the spatial period. The wires are modeled as impedance boundaries characterized by the surface impedance Z s = 2i/(ωε 0 (ε m − 1)R) where ε m is the metal relative permittivity. The wires are embedded in a host medium of permittivity ε h .

Uniaxial wire medium
The simplest example of a wire metamaterial is the so-called uniaxial wire medium. It consists of a square lattice of parallel and infinitely long metallic wires oriented along a fixed direction, here taken as theẑ direction as represented in Figure 2(a).
The study of such systems has a long history (dating back to the 1950s) that was renewed at the turn of this century after the discovery of negative index metamaterials [46][47][48][49][50][51][52].
As shown in [45], the application of the homogenization scheme of Section 3 to this wire metamaterial is particularly simple. Indeed, the current density induced on the metallic wires surface can be accurately modeled by a single expansion function: Note that the electric current density is proportional to p ind . Using w 1,k (r) in (33) and (34), it can be shown that the nonlocal effective permittivity reduces to [10,45,53] where f V = πR 2 /a 2 is the volume fraction of the wires and β p is the plasma wavenumber for an array of parallel PEC wires. The parameter β p depends solely on the system geometry (see the next subsection for the expression and Ref. [45] for further details). As seen, the effective permittivity of the uniaxial wire medium depends on the z component of the wavevector along the wires (k z ), which leads to a pole of the material response at low frequencies (for good conductors ε m → −∞ and the pole occurs for k z ≈ ω/c). Thereby, the spatial dispersion effects are rather strong. This feature has several nontrivial consequences, e.g., it implies that the medium may support two modes with the same polarization [50,53,54]. For a full discussion about the uniaxial wire medium modes the reader is referred to [53]. The uniaxial wire medium has interesting applications in subwavelength imaging when operated in the canalization regime [31,[55][56][57][58][59][60][61][62][63][64].

Double wire mesh
A more complex situation from the homogenization perspective occurs when a second array of parallel wires with a different orientation is inserted in between the first set of wires (see Figure 2(b)). Such structures are usually referred to as double wire meshes, and can have several interesting applications and rather exotic physics [65][66][67][68][69][70][71][72]. While the expression of the nonlocal effective permittivity of this metamaterial is well known [66], its direct derivation using the homogenization formalism of Section 3 was not reported previously in the literature. Since we believe that the proof is pedagogical we do so in the following.
The wire arrays are oriented along the generic directionsû 1 andû 2 . For simplicity we restrict our analysis to PEC wires (Z s = 0), orthogonal to each otherû m ·û n = δ m,n , with m, n = 1, 2, and consider a cubic lattice with period a. Similar to the case of the uniaxial wire medium, one expansion function per wire (two in total) is sufficient to obtain an approximate analytical expression for the effective permittivity. The expansion function that models the density of current induced in the nth wire oriented alongû n (here assumed parallel to one of the coordinate axes) is taken as w n,k (r) = e ik·r 2πRû n , n = 1, 2.
Substituting the above formula into (33), one finds that the effective permittivity can be written as: To obtain χ m,n , we substitute (37) into (34) and use the regularized lattice Green's function given by [7] Φ p0 (r|r ; ω, k) = 1 where k J = k+k 0 J with k 0 J = j 1 b 1 + j 2 b 2 + j 3 b 3 and the b i s are the reciprocal lattice primitive vectors. The second identity is valid in the long-wavelength limit, ω/c π/a and |k| π/a [45]. After straightforward calculations it is found that: where k i = k ·û i and β m,n is a quantity that depends only on the geometry of the system, and is given by where r n denotes the center of the nth wire in the unit cell (the nth wire translated by −r n is centered at the origin) and J 0 is the Bessel function of 1st kind and 0th order. For m = n, β m,m = β p is the plasma wavenumber for an array of parallel PEC wires mentioned in the last subsection [45]. For m = n the parameter β m,n is given by a simple series with an oscillating generic term due to the nonzero complex exponential coefficient. In contrast, for m = n the parameter β m,n is determined by a double series with the generic term of summation strictly positive. Due to this reason, one has |1/β 2 m,n | 1/β 2 p for m = n. The approximation is better for a larger physical distance between the two sub-lattices, as for a larger distance the complex exponential will oscillate faster. Thus, the off-diagonal terms of [χ m,n ] can be dropped, and with this approximation the inverse matrix elements are given by: Substituting this expression into (38) it is found that the dielectric function of the double wire medium is This result agrees with the nonlocal effective permittivity for perfect electric conducting wires derived in [73] using a slightly different approach. Similar to the uniaxial wire medium, the effective permittivity of the double wire mesh is strongly spatially dispersive. Remarkably, each wire array contributes independently to the permittivity function such that ε ef /ε 0 = I + n=1,2 where ε n is the permittivity of the nth wire array alone. The above derivation can be readily extended to plasmonic wires with a finite conductivity [66] and to triple non-connected wire arrays [73]. Furthermore, the proof can also be generalized to the case where the wire arrays are not perpendicular [70]. Also in this case, with similar approximations, one finds that each wire array contributes independently to the permittivity function.

3D connected wire mesh
The strong spatial dispersion characteristic of nonconnected wire arrays can be tamed by connecting the metallic wires, so that effectively the structure is formed by a single piece of conductor [73,74]. Here, we illustrate this by considering a 3D connected wire mesh formed by three orthogonal and connected sets of wires as represented in Figure 2(c).
In this system, because of the discontinuity of the induced current at the wire junctions, a single expansion function per wire is not enough to correctly homogenize the electromagnetic response. Instead, it can be shown that five expansion functions w n,k are required to obtain an approximate analytic expression of the effective permittivity [45,73]. Relying on an approach similar to that of the previous subsection (the details can be found in [45]), it can be shown that the effective permittivity of this metamaterial is where the transverse and longitudinal components are given respectively by In the above, l 0 = 3/(1 + 2β 2 p /β 2 1 ) and β 1 is a constant (with unities of wave number) that depends solely on the geometry of the structured material (see [45] for more details).
Remarkably, the 3D connected wire medium has a homogenized response equivalent to that of a plasma described by the hydrodynamic model [75]. In particular, the response to transverse waves (with electric field perpendicular to the wave vector) is described by the k-independent transverse permittivity ε t . However, the 3D connected wire medium remains spatially dispersive. The reason is that the response to longitudinal waves (with electric field parallel to the wave vector) is described by a k-dependent longitudinal permittivity ε l . The effects of spatial dispersion are several orders of magnitude stronger than in metal nanostructures at optics because the parameter l 0 is relatively small (l 0 ≈ 2). The effects of spatial dispersion can be further suppressed by loading the wires with metal plates, which leads to l 0 1 [74,76].
In general, the 3D connected wire mesh supports 3 electromagnetic modes: a longitudinal and two transverse plane waves. Propagation is only feasible above the effective plasma frequency. Thus, for long wavelengths the 3D connected wire mesh is completely opaque to radiation. For further details about the electrodynamics of the connected wire medium, the reader is referred to [45].

Quasi-static model
The nonlocal response of wire metamaterials can be explained by a quasi-static model developed in [76]. In the quasi-static model the macroscopic electromagnetic fields are coupled to the currents in the wires and to an "additional potential". The additional potential may be understood as the average voltage drop from a given wire to the boundary of the cell wherein it is contained [76]. Both the additional potential (ϕ) and the current are interpolated as continuous functions defined in all space. As reported in [77,78], the quasi-static model is particularly useful in problems involving interfaces, e.g., to obtain "additional boundary conditions", and to derive conservation laws [78].
For the case of the uniaxial wire medium (with wires oriented alongẑ) of Section 4.1.1 the quasi-static model is determined by: where E and H are the macroscopic electromagnetic fields (E = {e(r)} av and H = {b(r)} av /µ 0 ), E z = E ·ẑ, ε h is the permittivity of the host medium and C , L and Z w are the capacitance, inductance and self-impedance of a wire per unit of length, respectively, defined as in [76]. As seen, in this theory the macroscopic Maxwell equations are coupled to a set of differential equations governing the dynamics of the internal degrees of freedom of the medium (I z and ϕ).
The quasi-static model (47) fully describes the physical behavior of the uniaxial wire medium, as it can be transformed into the nonlocal model (36) by expressing I z and ϕ in terms of the macroscopic fields [76]. Importantly, the quasi-static model is l oc al as it corresponds to a standard partial-differential system. The differential operators act on the 8-component state vector (E, H, ϕ, I z ). The nonlocality of the electromagnetic response is a consequence of the fact that I z and ϕ are coupled to each other through a space differential operator (∂/∂z), different from conventional local media where the internal degrees of freedom are coupled through time differential operators (∂/∂t ).
Finally, it is worth mentioning that the quasi-static model is not restricted to the description of the uniaxial wire medium, as it can be extended to more complex connected and nonconnected wire medium topologies [76].

Poynting vector
In spatially-dispersive media, the energy density flux is not given by the standard textbook formula of the Poynting vector E×H [33,34,79,80]. For the case of lossless materials characterized by a nonlocal dielectric function the (time-averaged) Poynting vector must instead be calculated using: Here,l is a generic (real-valued) unit vector. It is implicit that the spatial dependence is of the form e ik·r with k real-valued and that the magnetic response is trivial. The formula can be generalized to a superposition of plane waves possibly associated with complex-valued wave vectors [79,81].
It was demonstrated in Refs. [79,80] that for a generic dielectric metamaterial, Equation (48) agrees precisely with the cell-averaged microscopic Poynting vector, provided the effective dielectric function is determined with the homogenization method of Section 3. Therefore, the macroscopic Poynting vector can be understood as a cell-averaged microscopic Poynting vector. Evidently, in wire metamaterials the Poynting vector can be determined using (48), using the relevant expression of the nonlocal permittivity in the formula. However, as previously mentioned, such formalism is only applicable to plane waves. A more general and useful expression for the Poynting vector can be obtained using the quasi-static model of Section 4.2. Indeed, based on (47) it is possible to derive a generalized Poynting theorem, which for the particular case of the uniaxial wire medium yields the following expression for the Poynting vector [78]: As seen, the Poynting vector is written in terms of the macroscopic electromagnetic fields and of the internal degrees of freedom (I z and ϕ) of the metamaterial. It can be verified that in the lossless case and for a spatial dependence of the form e ik·r with k real-valued the above expression reduces to (48). However, Equation (50) is more general than (48) as it can be applied to arbitrary electromagnetic field distributions. The stored energy in the wire metamaterial can also be expressed in terms of the state vector (E, H, ϕ, I z ), and for more details the reader is referred to [78].

Additional boundary conditions
One important consequence of spatial dispersion is that the usual Maxwellian boundary conditions, i.e., the continuity of the tangential E and H fields, are insufficient to solve wave propagation problems in the presence of interfaces [34,54,[81][82][83][84][85][86]. For example, consider a planar interface between two regions: a standard dielectric and a generic spatially dispersive material characterized by a nonlocal dielectric function. Suppose that a plane wave propagating in the dielectric illuminates the spatially-dispersive material half-space. The standard approach to find the scattered waves is to expand the electromagnetic fields into plane waves in the two regions and then to match the fields at the interfaces by imposing the standard Maxwellian boundary conditions. In standard dielectrics, there are exactly two plane-waves associated with an energy flow propagating away from the interface, i.e., there are only two polarization states per propagation direction. The potential problem is that in a nonlocal material the allowed number of polarization states per propagation direction may be greater than two, i.e., the medium may support "additional" waves. For example, a uniaxial wire medium typically supports three independent polarization states [85]. Consequently, it is generally impossible to solve a scattering problem relying only on the Maxwellian boundary conditions because the number of unknowns (number of waves that can be excited) is greater than the number of equations (number of boundary conditions). The problem is under-determined and additional boundary conditions (ABCs) are needed. The number of ABCs must be the same as the number of additional waves. For wire metamaterials, the ABC requirement is particularly clear from the quasi-static formulation of Section 4.2 where it is evident that in a scattering problem the boundary conditions for the internal degrees of freedom ϕ and I z must also be provided [77]. Thus, one needs to specify how the relevant internal variables behave at the interface. Unfortunately there is no systematic theory to find the ABCs, and their derivation must be based on the specific microscopic properties of the system under consideration. In particular, it is underlined that the ABCs (which are interface dependent) cannot be directly obtained from the nonlocal dielectric function, i.e., from the bulk response.
Here, we restrict our attention to an interface between a wire metamaterial and a standard dielectric. This situation covers the important case of an interface between wire media and air, which is of particular interest for scattering or imaging applications. Evidently, the microscopic electric currents in the metal wires are interrupted at the interface. Hence, for a system with N independent wires in the unit cell, it follows that at the dielectric interface J av ·û n = 0, n = 1, . . . , N where J av is the cell-averaged microscopic conduction current andû n is the unit vector oriented along the direction of the nth wire array [81]. The vector J av can typically be written in terms of the dielectric function of the medium [81].
In the particular case of a uniaxial wire medium, the ABC in the quasi-static model assumes the simple and intuitive form I z = 0. This ABC (together with the standard Maxwellian boundary conditions) can be expressed in terms of the electromagnetic fields as [85]: wheren is the unit vector normal to the interface, ε h is the host medium permittivity and ε d is the dielectric permittivity. Note that (52) is not equivalent to the continuity of the electric displacement vector, since the effective permittivity of the wire medium is different from ε h . Similar ideas are used to obtain the ABCs for the case of connected wire arrays [45], interlaced wire meshes [87] and for wires terminated with lumped loads [77,81].
It should be noted that wire metamaterials are amongst the very few examples of structured media for which there is a clear understanding of how to model the nonlocal effects near interfaces [77,81,85,88]. Another example, less well-developed, is the case of quadrupolar metamaterials characterized by weak spatial dispersion [89,90]. The general problem of characterizing the interface response of a generic nonlocal metamaterial is unsolved.

Anomalous refraction and light tunneling with wire metamaterials
To illustrate some of the unusual opportunities created by the spatial dispersion in wire metamaterials, we review in the next subsections the effects of anomalous light refraction and anomalous light tunneling.

Anomalous refraction in arrays of non-connected crossed metallic wires
As noticed in [68], a remarkable consequence of spatial dispersion is the possibility to achieve a low-loss and broadband regime of anomalous light refraction such that, contrary to what happens in a standard glass prism, longer wavelengths are more refracted than shorter wavelengths. This effect is forbidden by Kramers-Kronig relations in transparent and local materials. It may however occur in a prism made of a double-wire medium formed by nonconnected wires [68,72], see Figure 3. To understand the physical origin of this effect, we consider the wave propagation in an unbounded double wire medium made of perfectly conducting wires lying in the xoz plane and tilted by ±45 • with respect to the x-axis, as represented in Figure 3(b). For simplicity, we assume that the wave propagates along the z-direction (k = k zẑ ). For fields polarized along the x-direction the characteristic equation is where ε xx =x · ε ef ·x is the relevant component of the nonlocal effective permittivity for this polarization. According to the effective model (43), ε xx is given by Substituting ε xx into (53) and solving for k z , it is found that k z = (ω/c)n ef , where n ef is the effective refractive index of the double wire medium given by [68] n ef = 3 2 Remarkably, even though the metamaterial is lossless, the refractive index is a strictly decreasing function of the frequency. This unique property is only possible due to the spatial dispersion which makes the permittivity seen by the transverse field (polarized along x) dependent on a perpendicular wave vector component (here k z ). The same effect occurs for other propagation directions in the yoz plane. Due to the anomalous permittivity dispersion, a prism made of a crossed wire mesh can create a reverse rainbow as demonstrated theoretically in [68], and experimentally confirmed at microwave frequencies in [72]. In the experiment the prism is formed by a stack of dielectric slabs printed with the ±45 • -oriented metallic strips. A sample of the experimental results is presented in Figure 3(c). As seen, unlike conventional prisms, in the metamaterial prism the refracted beam comes out closer to the normal of the output interface when the frequency is increased. Materials with anomalous light dispersion may be useful for many applications, e.g., for the compression of light pulses or for the correction of achromatic aberrations [71].

Anomalous light tunneling in interlaced wire meshes
Here, we consider a metamaterial formed by two interlaced 3D connected wire meshes (mesh A and B ) separated by half-lattice period, see Figure 4(a) [87,91]. In what follows, we characterize the effective response of this "interlaced wire medium" and discuss a counter-intuitive tunneling effect rooted in the spatially dispersive response of the metamaterial.
Consider the general problem of homogenization of a metamaterial formed by two networks of inclusions A and B . In general, to find the effective response it is crucial to take into account the complex near-field interactions between the different types of scatterers. However, when the scatterers are physically distant in the unit cell it may be a good approximation to consider that each scatterer behaves as a "macroscopic source" from the point of view of the other scatterer. Essentially, this approximation is good when only the smooth (slowly varying) part of the fields radiated by one of the scatterers influences the currents on the other scatterer. It can be formally shown that in these conditions each component of the metamaterial contributes independently to the dielectric function such that [87,93]: where ε i ef with i = A, B is the nonlocal effective permittivity of the metamaterial formed only by the i th network of inclusions.
From the results of Section 4.1.2, it is readily recognized that a double-wired mesh of nonconnected wires provides a nontrivial example of a system in which the different types of scatterers interact with one another as "macroscopic sources". Interestingly, it turns out that the interlaced wire mesh of Figure 4(a) has the same property when the two 3D wire meshes are separated by the maximal possible distance (a/2) [87,92]. For the interlaced wire mesh ε i ef is the nonlocal effective permittivity of the (isolated) i th wire mesh given by (44).
Intuitively, the interlaced wire mesh should be opaque to radiation for frequencies below a certain effective plasma frequency. Surprisingly, that is not the case and it turns out that the metamaterial supports a longitudinal-type mode at arbitrary low-frequencies, as illustrated in the band diagram in Figure 4(b) [87,91,92]. This feature contrasts sharply with the properties of the individual 3D wire meshes, which do not support electromagnetic propagation for long wavelengths.
Remarkably, the low-frequency mode can originate a rather counter-intuitive tunneling effect. To illustrate this, we consider an interlaced wire mesh slab of finite length L (see the inset of Figure 5(a)). Using the effective permittivity model (56) and suitable ABCs, it is possible to find the transmission coefficient |T | of the slab [92]. Strikingly, as shown in Figure 5(a), provided the wire radii are different (r A = r B ) an incoming plane wave can tunnel through the metamaterial slab for large incidence angles. This transmission anomaly is due to a Fabry-Pérot resonance of the low-frequency longitudinal mode of the metamaterial. At the resonance the longitudinal wave vector satisfies the condition k z L = nπ, with n = 1, 2, . . . (see Figures 5(b) and (c)).
The physical origin of the tunneling anomaly is a Fano-type resonance [93] that occurs when r A = r B and enables the cancellation due to destructive interference of the scattering by the two subcomponents of the interlaced wire mesh. This metamaterial structure may be useful for angle-dependent filtering and sensing. For a detailed discussion of the physical properties of the interlaced wire mesh, the reader is referred to [92].

Conclusions
We presented an overview of a first principles homogenization approach based on an effective Hamiltonian that describes exactly the time evolution of the wave packet envelope when the initial state is less localized than the lattice period. The effective Hamiltonian determines completely the band diagram of the time-stationary states of the periodic system. The homogenization formalism can be applied to a wide range of physical systems. Its specific implementation for the case of nonmagnetic periodic electromagnetic metamaterials was described.
In particular, we focused our attention in the homogenization of wire metamaterials with diverse topologies. These structures are typically characterized by a strong nonlocal response in the long wavelength limit. In wire metamaterials formed by two or more non-connected networks, each metal network may contribute almost independently to the permittivity function. We underlined the nontrivial implications of the spatially dispersive response in different contexts, e.g., the emergence of additional waves, additional boundary conditions and a non-standard definition of the Poynting vector. Finally, we illustrated the richness of the physics of the wave propagation in wire medium, showing that it can lead to a counter-intuitive tunneling effect and anomalous frequency dispersion.