{'fraction_non_alphanumeric': 0.03144604972016139, 'fraction_numerical': 0.008980866848887154, 'mean_word_length': 4.8471841704718415, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fourier-Domain Electromagnetic Wave Theory for Layered Metamaterials of Finite Extent

The Floquet-Bloch theorem allows waves in infinite, lossless periodic media to be expressed as a sum of discrete Floquet-Bloch modes, but its validity is challenged under the realistic constraints of loss and finite extent. In this work, we mathematically reveal the existence of Floquet-Bloch modes in the electromagnetic fields sustained by lossy, finite periodic layered media using Maxwell's equations alone without invoking the Floquet-Bloch theorem. Starting with a transfer-matrix representation of the electromagnetic field in a generic layered medium, we apply Fourier transformation and a series of mathematical manipulations to isolate a term explicitly dependent on Floquet-Bloch modes. Fourier-domain representation of the electromagnetic field can be reduced into a product of the Floquet-Bloch term and two other matrix factors: one governed by reflections from the medium boundaries and another dependent on layer composition. Electromagnetic fields in any finite, lossy, layered structure can now be interpreted in the Fourier-domain by separable factors dependent on distinct physical features of the structure. The developed theory enables new methods for analyzing and communicating the electromagnetic properties of layered metamaterials.


Introduction
Recent contributions to the long tradition of inquiry into the electromagnetic properties of planar layered structures [1][2][3][4][5][6][7][8][9][10] have been sparked by the novel conceptualization of these structures -particularly those composed of sub-wavelength-thick layers of metal -as metamaterials. The metamaterial concept is used to describe a structure with sub-wavelength scale heterogeneity in terms of plane-wave parameters such as refractive index and impedance. Planar layered structures, which possess heterogeneity along just a single direction, are the simplest metamaterial form and provide an experimentally feasible template for metamaterial devices operating at visible frequencies and beyond due to the availability of thin film deposition techniques with layer thickness control on sub-nanometer scales. Exciting applications for layered metamaterials include flat lens imaging [11][12][13][14] and analog computation [15].
Classification of a heterogeneous structure as a metamaterial begins by seeking an analog homogeneous structure possessing a plane-wave solution that mimics the more intricate wave solution corresponding the original structure. This process, known as homogenization, yields familiar plane-wave parameters to approximate wave behavior in the heterogeneous structure. Each homogenization technique invokes a unique * Corresponding author: kenneth.chau@ubc.ca set of assumptions, which are not always justified, to arrive at its plane-wave parameters. Effective medium theory can be used to define an effective permittivity tensor through volumetric averaging of the local permittivity values [2,16,17], which, for a layered structure, simplifies to a thickness-weighted average of the layer permittivity values. Although effective medium theory is intuitive, it relies upon the electrostatic approximation which neglects time-derivative terms in Maxwell's equations. The scattering parameter method [18][19][20][21] is based on equating the reflection and transmission coefficients of a heterogeneous structure to those of an equivalent homogeneous structure. Drawbacks include non-uniqueness [22,23] and the absence of correlation to the fields inside the structure. It is possible to derive effective constitutive parameters by averaging local permittivity values weighted by the fields [24] or energy densities [25] inside a structure, although this method also suffers from non-uniqueness. If a structure is periodic, a unique set of homogeneous parameters can be defined using the Floquet-Bloch theorem [26,27]. The wave inside a structure is decomposed into a set of Floquet-Bloch modes k F B [28][29][30][31] and if one mode carries dominant power, it is assumed to approximate the entire wave and its plane-wave parameters are conferred to the structure [32][33][34][35][36][37][38]. The Floquet-Bloch modes of a periodic layered structure can be found by imposing translational invariance of the wave over a period within a multiplicative exponential factor [7], a procedure that implicitly assumes infinite extent and no loss. When the medium is lossy, the Floquet-Bloch modes are complexvalued [39][40][41], but no longer discrete [32]. When the medium is finite, translational invariance is altogether lost due to reflections from the end facets of the medium and a dependence on the excitation conditions in the external bounding media [29].
The validity of the metamaterial approximation depends upon consistency between plane-wave parameters derived by homogenization and complete wave solutions to Maxwell's equations. Direct comparison between the two can be challenging because plane-wave parameters are most naturally interpreted in the spatial-frequency domain, whereas wave solutions from numerical solvers and analytical techniques (such as those based on matrix formalism [3,7,9,10,42]) are commonly communicated in the spatial domain. The obvious resolution is to obtain Fourier-domain wave solutions, but the most common solution techniques based on Fourier expansion of a single field component (known as the plane-wave expansion method [43,44]) or Fourier expansion of the electric and magnetic field components into a series of plane waves [30,31,38] are limited to lossless, non-dispersive, and infinitely periodic media not applicable to finitesized, metal-dielectric metamaterials. Recently, Fourier transformation of the numerically-computed wave solution in a finite-sized, lossless periodic layered medium has revealed a correlation between the spatial frequency content of the wave and the Floquet-Bloch modes corresponding to an unbounded version of the periodic medium [22]. This correlation, however, has not been made mathematically explicit and has yet to be established for the case of loss, a condition vital for analysis of metamaterials which frequently incorporate metals.
The goal of this work is to analytically derive the Fourier-domain electromagnetic wave solution for a general layered configuration that accommodates periodic or aperiodic layer arrangements under the realistic constraints of loss and finite extent. The wave solution is distilled into a compact, factorized form consisting of a term due to reflections from the medium boundaries, another associated with the complex-valued Floquet-Bloch modes of the unit cell, and a third that is dependent on the layer thickness and composition. In contrast to plane-wave parameters obtained by homogenization, the factorized, Fourier-domain wave solution is fully consistent with Maxwell's equations and provides physical insight into the effects of loss, finite extent, and periodicity. We apply the wave solution to analyze the electromagnetic fields and power flow in several representative metal-dielectric layered metamaterials at optical frequencies and to validate plane-wave parameters obtained by homogenization techniques such as effective medium theory, scattering parameter retrieval, and Floquet-Bloch analysis. Of the three, Floquet-Bloch analysis provides the closest approximation to the complete wave solution, an observation that is consistent with previous reports [22] and can now be explained by the explicit appearance of Floquet-Bloch modes in the factorized, Fourier-domain wave solution.
A factorized Fourier-domain wave solution contributes to the resolution of the following research questions pertinent to layered implementations of metamaterials: How accurately do plane-wave solutions derived by homogenization approximate the complete wave solution? How do the conditions of material loss and finite extent, both of which are important for the analysis of any practical layered metamaterial device, affect the validity of the metamaterial concept? Is there a systematic method to understand the relationship between the internal electromagnetic fields and power flow in a layered structure and the physical parameters of the structure? Finally, are there better ways to describe the electromagnetic properties of layered metamaterials than through approximate plane-wave solutions provided by homogenization?

Generalized Description of a Layered Medium
We consider a generalized one-dimensionally periodic medium ( Figure 1) immersed in free space and composed of M repeated unit cells, each consisting of J layers, yielding a total of M J layers. The unit cells are referenced by the integer m = 0, ..., M − 1 and the layers within any unit cell are referenced by the integer j = 1, ..., J. Each layer in the medium is uniquely labeled by the integer = mJ + j = 1, ..., M J and, for sequential consistency, the free-space half-spaces to the left and right of the medium are labeled = 0 and = M J + 1, respectively. Layer has a thickness of d and its linear electromagnetic properties are generally specified by a complex-valued relative permittivity (the underline denotes a complex variable) and a complex-valued relative permeability µ , resulting in a complex-valued refractive index [45] which is permitted to have a real part that is positive (right-handed) or negative (left-handed). The total thickness of the layered medium is L = M J =1 d . Due to periodic repetition of the unit cell, the quantities d , , µ , and n , corresponding to layer j of an arbitrary unit cell m, can be equivalently denoted d j , j , µ j , and n j corresponding to layer j of unit cell m = 0. The plane of the layers is aligned parallel to the xy plane and we denote the location of the plane between layer and + 1 as z , setting the position of the plane of the first interface z 0 = 0 without loss of generality.
An electromagnetic plane wave is incident onto the medium from the half-space z < 0 inclined at an angle θ in the xz plane. Due to the independence of this configuration with respect to the y coordinate, any solution can be expressed as a linear combination of solutions obtained by assuming either transverse-electric (TE) polarization (electric field E aligned to the y axis) or transverse-magnetic (TM) polarization (magnetic field H aligned to the y axis). Here, we treat the case of TM polarization, noting that the transformations E → − H, H → E, and µ yield complimentary equations for TE polarization [9].

Wave Solution in the Spatial Domain
We derive a representation of the electromagnetic wave solution in the finite one-dimensional periodic layered structure for the case in which one of the end facets is subject to plane-wave illumination. An incident TM-polarized electromagnetic wave is given by H = H 0 e i(kx,0x+kz,0z)ŷ , where H 0 is the amplitude, k x,0 and k z,0 are the real-valued wave-vector components along the x-and z-axes, respectively, and k 0 = k x,0 + k z,0 is the free-space wave vector. The wave is time-harmonic where an e −iωt dependence is assumed but suppressed. Invoking field continuity across the interfaces, the magnetic field in an arbitrary layer can be written as a sum of two counter-propagating waves using matrix formalism where T denotes the transpose operator, A and B are the wave coefficients, and k z, is the wave-vector component in layer along the z-axis (note that k z, corresponding to layer j of an arbitrary unit cell m can be equivalently denoted k z,j corresponding to layer j of unit cell m = 0). The wave-vector component k z, is related to the layer refractive index by where, according to (1), n can have a real part that is either positive or negative, describing a right-or lefthanded medium, respectively. The wave solution can be solved by relating the wave coefficients A and B across the M J + 1 boundaries. The wave coefficients in an arbitrary layer can be related to the coefficients in an adjacent layer + 1 by where the propagation matrix P corresponding to layer is given by and the transmission matrix T corresponding to the interface between layer and + 1 is given by with p = ( +1 k z, )/( k z, +1 ). Assuming unidirectional wave excitation from the left half-space, the wave coefficients in layer 1 are related to the coefficients in the left half-space by where T 0 is the transmission matrix from the left halfspace into the first layer, r is the complex-valued reflection coefficient of the entire layered system, and the incident wave amplitude has been assumed to be unity. The wave coefficients in layer M J are related to the coefficients in the right half-space by where T M J is the transmission matrix from layer M J into the right half space, t is the complex-valued transmission coefficient. Relation of the wave across the M J + 1 boundaries yields 2M J + 2 linear equations, which is sufficient to solve for the 2M J + 2 unknowns (r and t, in addition to the 2M J wave coefficients in the M J layers).
Upon solving for the unknown quantities, H in each layer is completely specified and we can succinctly express the total field distribution in the spatial domain as where z c, is the location of the center of layer and the rect function is defined as Equation (9) offers a valid representation of the field as a piece-wise function subdivided into spatial intervals corresponding to the layer regions. Although this form is amenable to numerical routines for solving sets of linear equations, there are at least two disadvantages. First, compartmentalization of the wave solution into the individual layers does not afford physical insight into the collective behavior of the solution across repeated sets of layers. Second, representation of the solution in the spatial domain does not produce immediate connections to homogenization parameters, which are generally represented in the spatial-frequency domain. In the next section, we apply Fourier transformation to the piecewise wave solution and demonstrate the utility of this strategy for decomposing the wave solution into physically meaningful terms and for enabling comparisons with plane-wave parameters obtained by homogenization.

Factorized Fourier-Domain Wave Solution
We re-express the general wave solution given by (9) in the spatial-frequency domain by where κ x and κ z are the spatial frequency variables along the respective x and z directions. Substitution of (9) into (10) and development of the integrand using well-known Fourier relations and theorems yields where δ is the Dirac delta function.
In this section, we outline a series of mathematical manipulations that can be applied to (11) to produce a physically insightful factorized Fourier-domain wave solution. First, we re-write the single summation in (11) as a nested double summation over the number of layers in a unit cell and the number of unit cells by making the variable substitutions d = d j and k z, = k z,j and the index substitution = mJ + j, resulting in Carrying out the convolution operation in (12) and using the relation z c,mJ+j − z mJ+j−1 = d j /2 yields The unit cell summation in (13) can be simplified using the relationship where D = J j=1 d j is the thickness of the unit cell and z j−1 is the position of the interface between layer j − 1 and j within unit cell m = 0. Substitution of (14) into (13) gives Equation 15 expresses the wave solution in terms of wave coefficients A mJ+j and B mJ+j distributed throughout the medium, which provides little additional insight over the spatial-domain representation of the wave solution in (9). We can further simplify the solution in terms of the wave coefficients in just the first unit cell by using the matrix relationship between wave coefficients in different layers. The wave coefficients in layer are related to the coefficients in an arbitrary layer + s (where the integer s ≤ M J − ) within the layered medium by where the transfer matrix W +s, is determined from the transmission and propagation matrices by The coefficients A mJ+j and B mJ+j corresponding to layer j within an arbitrary unit cell m can be related to the coefficients A j and B j corresponding to layer j within unit cell m = 0 by where U j is the unit cell transfer matrix from layer j to j + J and can be expressed as We can now simplify (15) in terms of wave coefficients A j and B j distributed throughout a unit cell referenced from layer j The unit cell transfer matrix referenced from layer j, U j , can be related to the unit cell transfer matrix referenced from layer 1, U 1 , using the relation Eigenvalue decomposition of U 1 yields where Q is the eigenvector of U 1 and λ is the eigenvalue of U 1 whose diagonal elements are the corresponding eigenvalues λ + and λ − . Because the determinant of U 1 is unity, the eigenvalues are inverses of each other, λ − = 1/λ + . As a result, the eigenvalues can be related to the Floquet-Bloch mode, k F B , by To arrive at the final form of the magnetic field solution from (15), we apply eigenvalue decomposition of the unit cell transfer matrix and relate the wave coefficients in layer j to the wave coefficients in the left half space, resulting in where we have highlighted three distinctive matrix factors -a layer matrix L j dependent on the thickness and wave vector in the jth layer of the unit cell, a Floquet-Bloch matrix FB dependent on the eigenvalues of the unit cell, and a weighting matrix C dependent on the reflection coefficient. The magnetic field solution can be succinctly written as where the three matrix factors have the following general form In (26), we have distinguished "forward" and "backward" elements of each matrix factor using the superscript "+" and "−", respectively. Developing the matrix factors in (25) yields which, in contrast to the spatial-domain solution in (9), requires summation of only J terms corresponding to the layers of the unit cell and describes collective wave behavior across repeated sets of unit cells through the elements of the Floquet-Bloch matrix, whose behavior will be examined in the next section.
Given the vector spectral magnetic field H(κ x , κ z ) = H(κ x , κ z )ŷ and the vector spectral electric field E(κ x , κ z ) = E x (κ x , κ z )x + E z (κ x , κ z )ẑ (which, like the magnetic field, can be represented by a factorized form using steps similar to those described in this section), it is possible to define the spectral time-averaged Poynting vector This is similar to the spectral Poynting vector proposed in Ref. [30,31] and used in Ref. [38] to analyze energy propagation of discrete Floquet-Bloch modes in infinite, lossless dielectric photonic crystals, except now extended to accommodate a continuous range of Fourier field components in a finite, lossy periodic system. It should also be noted that the spectral time-averaged Poynting vector defined in (28) is not equivalent to the Fourier transform of the spatial time-averaged Poynting vector, which would involve the convolution of the spectral electric and magnetic fields. It does, however, enable the spatial frequency κ z present in the electric and magnetic fields to be envisioned as an electromagnetic plane wave having a well-defined time-averaged direction of power flow, which is useful, for instance, for the characterization of the wave in terms of forward-and backward-propagating wave components.

Understanding the Roles of Finite Extent and Loss
The elements of the Floquet-Bloch matrix are given by where the Dirac comb ∆ 2π describing an infinite comb of discrete spatial-frequency harmonics spaced by 2π/D. Because the Dirac comb is convolved with the sinc function in (29), the elements F B ± contain harmonics widened through the combined effects of finite extent (M < ∞) and material loss ( (k F B ) = 0). The forward element F B + has a principal harmonic centered at (k F B ) and the backward element F B − has a principal harmonic centered at − (k F B ).
For a medium of infinite extent (M → ∞), the elements of the Floquet-Bloch matrix approach In the absence of loss, F B ± are discrete spectra with peaks at harmonics of ± (k F B ). The resulting magnetic field solution consists of discrete forward and backward Floquet Bloch modes (equivalent to the classical Floquet-Bloch solution) with amplitudes that can be explicitly determined by (2π) 2 J j=1 L + j C + and (2π) 2 J j=1 L − j C − , respectively. In the presence of loss, F B ± are continuous spectra with peaks that are centered about harmonics of ± (k F B ) and broadened due to the effect of (k F B )

Analysis of Metal-Dielectric Layered Metamaterials
We consider the configuration in which a metal-dielectric layered system is subjected to normal-incidence, TMpolarized illumination (k x,0 = 0). If the system consists of a bi-layer unit cell, the factorized, Fourier-domain magnetic field solution takes the form where L ± 1 and L ± 2 are the forward and backward elements of layer matrices respectively associated with layers 1 and 2 of the unit cell. Let's take the case bi-layer unit cell comprising a 30-nm-thick Ag layer and a 30nm-thick TiO 2 layer. We will consider a system composed of two unit cells (M = 2) and another composed of ten unit cells (M = 10). The systems are bounded by free space and excited by a wave with a free-space wavelength λ 0 = 365 nm. We assume that Ag has a complex refractive index of 0.076 + 1.605i (interpolated from experimental data [46]) and that TiO 2 has a real refractive index of 2.80. Figure 2 highlights, for the M = 2 and M = 10 cases, the modulus squared of the magnetic field solution and the elements of its matrix factors, in addition to the z-component of the time-averaged spectral Poynting vector. The magnetic field spectrum results from the confluence of broad spectral envelopes defined by L ± 1 (Ag layer) and L ± 2 (TiO 2 layer) and finer spectral combs defined by F B + C + and F B − C − . The Floquet-Bloch mode of the unit cell is k F B = 32.5 + 0.4i µm −1 , corresponding to a Floquet-Bloch refractive index n F B = 1.89 + 0.02i. The combs defined by F B + C + and F B − C − are offset; the former with a principal peak located at (k F B ) and the latter with a principal peak located at − (k F B ). Increasing M from 2 to 10 narrows the peaks of the combs, which in turn narrows the peaks in the magnetic field spectrum. Due to the losses in the metallic layer, the spectral peaks defined by F B + C + and F B − C − have finite width even in the limit of an infinite system, underscoring the limitation of the classical Floquet-Bloch solution for describing the spectral composition of the wave in lossy periodic systems. For both the M = 2 and M = 10 cases, the magnitude of F B + C + is larger than that of F B − C − . As shown in Fig. 3, |C + | 2 generally exceeds |C − | 2 , with the latter gradually approaching zero as the number of repetitions increases. It is noteworthy that the negative spatial frequency components of the wave solution are characterized by a time-averaged spectral Poynting vector that is negative, suggesting that these wave components are actually forward-propagating waves (having parallel phase and energy velocities) that have been reflected in the system. This is counter to the widely held belief that the negative spatial frequency components of a Floquet-Bloch field decomposition are backwardpropagating waves [28].
The Floquet-Bloch refractive index can take on values with negative real parts, particularly for TM-polarized illumination of finely layered structures composed of metal at frequencies near the bulk plasma frequency of the metal [14,47]. The physical significance of the sign of the real part of the Floquet-Bloch refractive index, however, is not immediately clear. As seen here, the magnetic field spectrum is generally dispersed in peaks spanning positive and negative spatial frequencies and a number of factors conspire to determine the relative weights of these peaks. Moreover, in systems composed of right-handed materials, it has been consistently observed that the positive and negative spatial frequency components of the wave are both forward-propagating waves.

Comparison with Parameters Derived from Homogenization Techniques
Application of common homogenization techniques, such as effective medium theory, scattering parameter retrieval, and Floquet-Bloch analysis, to describe the electromagnetic properties of layered structures yield planewave solutions communicated in terms of a complexvalued wave vector along the z-direction, k z (perpendicular to the plane of the layers). These plane-wave solutions are naturally parameterized in terms of the wave vector along the x-direction, k x (parallel to the plane of the layers). Mapping k z as a function of k x over the range |k x | < k 0 yields a two-dimensional wave vector diagram, or band diagram, that describes planewave solutions that can potentially be used to approximate the layered structure as a homogeneous medium. (A more complete discussion on the application of effective medium theory, scattering parameter retrieval, and Floquet-Bloch analysis to describe layered systems is provided in the Appendix.) There are two uncertainties associated with homogenization that have been recognized in the literature: 1) it is unknown how accurate a plane-wave solution approximates the true wave solution in the medium and, 2) for homogenization methods that produce an infinite number of plane-wave solutions, such as scattering parameter retrieval and Floquet-Bloch analysis, it is unclear which plane-wave solution best describes the true wave solution [22,23]. The analytical Fourier-domain wave solution provides a useful tool to overcome these challenges. Concatenation of the field solutions obtained by systematic variation of κ x over the range |κ x | < k 0 yields a two-dimensional Fourier-domain field distribution mapped as a function of spatial-frequency coordinates κ x and κ z . The relative field amplitude at a particular (κ x , κ z ) describes the relative contribution to the total wave solution made by a plane wave with spatial frequencies κ x and κ z along the respective x and z directions. Thus, the Fourier-domain field distribution can be overlaid on the wave vector diagrams to determine 1) if the plane-wave solutions produced by homogenization methods are present in the wave solution and 2) if present, the relative contributions of each plane-wave solution to the wave solution.
As a test case, we examine a finite periodic layered system immersed in free space and consisting of two bi-layer unit cells, each composed of a 30-nm-thick Ag layer and a 30-nm-thick TiO 2 layer (again for TM polarization at a free-space wavelength λ 0 = 365 nm). The Fourierdomain wave solution is presented as the magnetic field amplitude spectrum mapped as a function of κ z and |κ x | < k 0 , and the wave-vector diagrams are presented as complex-valued k z mapped over |k x | < k 0 .
As shown in Fig. 4, homogenization yields plane-wave solutions described by wave-vector diagrams consisting of either a single discrete branch (effective medium theory) or a family of discrete branches (scattering parameter retrieval and Floquet-Bloch analysis), whereas the magnetic field spectrum consists of a series of broad bands. Effective medium theory and scattering parameter retrieval are particularly poor at describing the spatial-frequency features of the wave solution; the former yields too few solutions and the latter yields too many solutions. Floquet-Bloch analysis gives the closest approximation of the wave solution, with the two most prominent bands of the wave solution fairly well aligned to two branches in the wave vector diagram. This is not surprising given the explicit appearance of Floquet-Bloch modes in one of the terms of the factorized Fourier-domain wave solution. However, the more complex features of the wave solution -variations of the magnetic field amplitude both within a band and between the bands and the presence of narrower subsidiary bands with weaker amplitudes -are beyond the descriptive capabilities of wave-vector diagrams obtained by Floquet-Bloch analysis. It should be noted that scattering parameter retrieval applied to a single unit cell of a symmetric layered system with identical bounding spaces produces the same solutions as those obtained by Floquet-Bloch analysis due to the equivalence of the two procedures under these conditions [20].

Discussion
Plane-wave illumination of a layered structure is an archetypal configuration that can be solved analytically or by numerical electromagnetic solvers such as those based on finite-difference techniques. Although these solutions are typically expressed in the spatial domain, it is straight-forward to apply numerical Fourier transformation to re-express them in the Fourier-domain [22]. Given the availability of these methods, is there need for deriving the analytical, factorized, Fourier-domain wave solution presented in Section 4?
In response, we must acknowledge that the value of a solution depends on both its consistency with fundamental laws (in this case, Maxwell's equations) and its ability to communicate physical meaning in a compact and efficient manner. In recent years, the metamaterial concept has been used to provide a succinct and intuitive way to describe the electromagnetic proper-ties of layered structures in terms of familiar plane-wave parameters. Although plane-wave solutions are easily communicated and understood, they generally sacrifice consistency with Maxwell's equations. On the other hand, analytical spatial-domain wave solutions in terms of wave coefficients dispersed throughout the layered medium are consistent with Maxwell's equations, but offer no insight into the collective behavior of the wave solution. Numerical wave solutions are also consistent with Maxwell's equations, but comprise data arrays that are valid only for a particular geometry and difficult to succinctly communicate. Thus, despite the availability of homogenization methods and analytical and numerical wave solutions for the layered geometry, the factorized, Fourier-domain wave solution developed here is important because it addresses a deficiency in methods to communicate the electromagnetic properties of layered metamaterials in a manner that is complete, consistent, and physically insightful.

Conclusion
Diverse classifications of layered metamaterials in the literature (using terms such as hyperbolic, anisotropic, and negative-index) stem from the use of different homogenization methods. Continued advancement of layered metamaterials will rely upon improved methods for communicating their electromagnetic properties. We have explored a new method to express the analytical wave solution for a lossy layered media of finite extent in a factorized form by analysis in the Fourier domain. Some interesting outcomes of this work include: 1. Fourier-domain representation allows the wave solution to be condensed into a product of three, physically-intuitive matrix factors.
2. Floquet-Bloch modes can be made to explicitly appear in the Fourier-domain solution, without a priori invoking the Floquet-Bloch theorem.
3. In the limit of no loss and infinite extent, the Floquet-Bloch matrix factor asymptotically approaches a discrete set of Floquet-Bloch modes.
The amplitude of the Floquet-Bloch modes can be explicitly determined from other matrix factors of the factorized, Fourier-domain solution.
4. The wave solution in periodic layered media of finite extent contains both forward (decay along +z) and backward (grow along +z) Floquet-Bloch modes.
5. For plane-wave illumination of layered structures composed exclusively of positive-index media, the positive and negative spatial-frequency components of the wave solution carry time-averaged power in the forward and backward directions, respectively. Thus, the negative spatial-frequency components are not backward-propagating waves with anti-parallel phase and energy velocities, but rather reflected forward-propagating waves with parallel phase and energy velocities.
6. Homogenization methods indicate the possibility of plane-wave solutions. A spatial-frequency representation of the wave solution is required to indicate if these plane-wave solutions are present and if so, their relative weights.
The theory developed here will help advance the science and technology of layered implementations of metamaterials by providing: 1) complete solutions to Maxwell's equations for these systems in a physically insightful and intuitive way and 2) a means to rigorously and directly validate current applications of plane-wave solutions critical to the very definition of a metamaterial.

Appendix
For completeness, we detail here the application of effective medium theory, scattering parameter retrieval, and Floquet-Bloch analysis to derive homogenization parameters for planar layered structures.

10.A. Effective Medium Theory Applied to a Bilayered Structure
There are several forms of effective medium theory that can be used to obtain an effective permittivity tensor of a layered structure. The simplest and most widely used is the Maxwell-Garnett formulation, which can be derived for the case of a repeated planar bi-layer system (with layer thicknesses d 1 and d 2 and complex permittivity values 1 and 2 ) by considering the average electric and displacement fields of an incident plane wave for different propagation directions under the electrostatic approximation. For plane-wave propagation perpendicular to the plane of the layers, continuity of the electric field and volumetric averaging of the electric field over the bi-layer yields an effective permittivity component are the fractions of total volume occupied by the layers. For plane-wave propagation parallel to the plane of the layers, continuity of the displacement field and volumetric averaging of the displacement field over the bi-layer yields an effective permittivity component Thus, the bi-layer system acts like a uniaxial crystal with an optical axis perpendicular to the plane of the layer.

10.B. Scattering Parameter Retrieval Applied to a Finite Layered Structure
Scattering parameter retrieval is a general technique based on the conceptual replacement of a finite heterogeneous system with a hypothetical, finite homogeneous system, where the parameters of the homogeneous medium have been retrieved assuming reflection and transmission coefficients identical to that of the heterogeneous medium. Here, we detail its application to describe the electromagnetic properties of the finite, lossy, periodic layered system generally depicted in Fig. 1. The reflection and transmission coefficients of the medium can be related to each other using transfer matrix formalism by where the total transfer matrix W can be expressed as The reflection and transmission coefficients of the medium can be subsequently written as respectively. We next conceptually replace the medium with a fictitious, homogeneous, effective medium with an identical length L and immersed in free space. The homogeneous medium is characterized by an effective permittivity e , an effective permeability µ e , an effective refractive index n e , and an effective impedance p e . The reflection coefficient, r e , and transmission coefficient, t e , of the effective medium are related to each other by where the total transfer matrix of the effective medium can be expressed as Equating r e = r and t e = t, the parameters of the effective medium yielding reflection and transmission coefficients identical to that of the layered medium are given by where N is an integer. There are two notable limitations associated with these parameters. First, there are an infinite number of possible solutions for the effective refractive index. Second, there exist two sets of solutions for p e and n e , each having opposite sign (although for a system with loss only one set satisfies the physical constraints that (n e ) ≥ 0 and (p e ) ≥ 0). From the effective parameters given in (A.11), the remaining effective parameters can be obtained using the relations e = p e n e µ e = n 2 e e , (A. 12) for which there again exists two set of solutions, each having opposite sign. Care must be taken if solutions for e and µ e are then used to calculate n e using (1) in the main text, as only one set of solutions for e and µ e yield the correct value of n e . It should also be pointed out that if the bounding half-spaces are not identical or if the unit cell is not symmetric, scattering parameter retrieval yields different effective parameters in the forward and backward direction.

10.C. Floquet-Bloch Analysis Applied to an Infinite Periodic Layered Structure
Floquet-Bloch analysis is commonly used to describe the properties of periodic layered media. Rigorously applicable to only infinite periodic media, it is based on the conceptual replacement of the heterogeneous unit cell of a periodic system with a hypothetical homogeneous unit cell, where the refractive index of the homogeneous unit cell has been selected to have a transfer matrix trace identical to that of the heterogeneous unit cell. To describe its salient features, we first assume an infinitely periodic layered medium having a unit cell with thickness D, where the fields at the boundaries of the unit cell are related by the unit cell transfer matrix U. We next replace the unit cell with a fictitious, homogeneous effective unit cell with an identical length D. Equating the traces of the transfer matrices for the unit cell and effective unit cell yields tr(U e ) = tr(U) 2 cos(k z,e D) = U 1,1 + U 2,2 . (A.14) Isolating the effective wave-vectorto one side of (A.14) yields k z,e D = ± cos −1 U 1,1 + U 2,2 2 ± 2πN, (A. 15) where, notably, the effective wave-vectoris non-unique. which describes a infinitely large family of discrete Floquet-Bloch modes. It is interesting to note that in contrast to scattering parameter retrieval, which generates two effective parameters based on two distinct equations (one for the transmission coefficient given by (A.9) and another for the reflection coefficient given by (A.10)), Floquet-Bloch mode analysis generates just one effective parameter based on one equation (equating the traces of the transfer matrices given in (A.14)).