On the unambiguous determination of effective optical properties of periodic metamaterials: a one-dimensional case study - DTU Orbit (04/11/2017)

On the unambiguous determination of effective optical properties of periodic metamaterials: a one-dimensional case study We show how branch ambiguities in the extraction of effective parameters is arising as a direct consequence of the underlying Bloch state physics. The mutual importance of the different branches in general depends on the experimental context, and we show how the Fourier spectrum of the field inside the metamaterial can be used to access this. Different numerical examples illustrate how a predominant branch may be identified for l a while at higher frequency the power may be distributed over more branches. This is in particular true near bandedges and strong resonances. Extensions to two-and three-dimensional metamaterial designs are discussed.


INTRODUCTION
Metamaterials constitute a new paradigm in the optics and electromagnetic communities and the concept holds promises for new and extraordinary materials, reviving old proposals of e.g.negative index properties [1].Carefully engineered artificial sub-wavelength electromagnetic structures are envisioned to play a central role [2,3,4,5].In this context, homogenization theory plays an important role in understanding and interpreting experiments as well as for further developing and optimizing metamaterial properties [6].Common to the variety of homogenization approaches [7]- [18] is the desire to introduce a spatially homogenous material with effective material parameters causing slowly varying effective electromagnetic fields which mimic the envelopes of the rapidly varying fields in the real structure.Ultimately, a metamaterial would effectively behave like other ordinary homogeneous materials (but with a negative refractive index, possibly dispersive) and thus it could enter on an equal footing with existing materials in practical optical and electromagnetic designs.However, metamaterial properties have turned out to be far more complex.As an example, the effective optical properties are in general anisotropic and they may depend on the angle of incidence [19,20] posing serious questions to the concept of isotropic homogenized metamaterials [21].Metamaterial properties can hardly be defined without formulating the experimental context, i.e. how is the electromagnetic response to be probed.In some cases, thin slabs of metamaterial may even respond differently from the bulk properties [22].
The S-matrix approach developed by Smith and coworkers [8,23] has proved extremely useful as it allows an interpretation of metamaterial properties probed in a scattering configuration.Today, it is the workhorse for the extraction of effective parameters, taking either simulated or measured S-matrix parameters as input.However, the method poses (in some cases serious) pitfalls and involves ambiguities related to the extraction of the real part of the phase index.This is a topic of ongoing work in the metamaterial community.In this work we use a one-dimensional problem to illustrate the generality of this ambiguity and show how it necessarily follows as a consequence of the general Bloch properties of periodic structures.As with the coupling to photonic crystals [24], incident light will in general excite a linear combination of the Bloch states in the periodic metamaterial.A Fourier transform spectrum of the spatial field variations inside an excited Bragg stack (of finite length) is used to illustrate this.We identify a pronounced branch in the λ a limit, but at higher frequencies the power may be distributed over more branches, thus potentially jeopardizing the existence and identification of a pronounced branch.While the one-dimensional considerations help in simplifying our presentation of the problem, our general conclusion remains valid for also two and three-dimensional periodic structures.

DISPERSION PROPERTIES OF BLOCH STATES
We consider the simple problem of a one-dimensional Bragg stack.Though the unit-cell structure is of course much simpler than the more complex two and three dimensional unit cells of typical metamaterials, the present problem still serves the important purpose of illustrating homogenization properties of layered metamaterials [25], and in particular the physics issues of branch ambiguities.Obviously, the access to analytical solutions is a clear advantage.We consider a Bragg stack with two (i = 1, 2) alternating layers of, for simplicity, nonmagnetic materials (µ i = 1) with dielectric function ε i (possibly complex valued) and thickness a i .The periodically varying dielectric function, (z + a) = (z), supports Bloch waves.
The dispersion relation ω(κ) can be found with the aid of Bloch's theorem and it is governed by the exact analytical expression cos (κa with where κ = κ + iκ is the Bloch wave vector, a = a 1 + a 2 is the period, and c is the vacuum speed of light.We imagine a continuous wave (CW) excitation of the electromagnetic states so that the angular frequency ω = ck is real valued, with k = 2π/λ being the free-space wave vector and λ the free-space wavelength.From Bloch theory it is well-known that ω is periodic in κ , i.e. ω(κ) = ω(κ + mG) with m being an integer and G = 2π/a being the reciprocal lattice vector.This means that there is a multiple of real solutions for a given frequency, while there is an unambiguous solution for the imaginary part κ .This is illustrated in Figure 1 for the case of a dielectric Bragg stack with ε 1 = 1 and ε 2 = 5 while a 1 = a 2 = a/2.The left panel shows the real part of the Bloch wave vector κ (horizontal axis) versus ω (vertical axis), while the right panel shows the corresponding results for the imaginary part of the Bloch wave vector κ .Quite often, bands are only plotted in the reduced-zone picture, i.e. for 0 < κ a < π, but in the left panel we have deliberately chosen to plot the bands also outside the first Brillouin zone, to emphasize the multiple number of Bloch states at a given frequency.On the other hand, there is a unique value of κ associated with each value of ω.We note that in the particular case with 1 and 2 being real, κ is non-zero only inside the band gaps.
As is well known, and evident from the discussion of Figure 1 above, the phase and group indices of a periodic structure can be quite different, contrary to homogeneous materials (with a weak material dispersion).For the periodic structure, the phase index n p = n p + in p is given by and the Bloch properties necessarily lead to an ambiguity for the real part of the phase index, in accordance with that found for the S-matrix approach discussed above.Interestingly, for the present problem the ambiguity even arises in the same mathematical way as in the S-matrix approach, i.e. it is related to the branch-issue of the inverse cosine function.Clearly, since there is a multiple of solutions for κ for a given ω we get a real part given by with m being an integer.Notice that this expression merely tells that for a given choice of branch (and thereby κ ) the phase index of the other branches is obtained by adding an integer of λ/a.On the other hand, the imaginary part is unambiguously determined.We emphasize that the multiple branches of Bloch states are potentially all candidates for physical solutions and thus the integer m represents more than just a mathematical ambiguity associated with choosing a branch of the arccos function.Taylor-expanding Eq. ( 1) near the Γ-point, i.e. (ω, κ) = (0, 0) we get a linear dispersion relation with an effective permittivity eff = n 2 eff given by corresponding to a geometrically average of the dielectric function.In the context of discussing branch ambiguities, it is however convenient to have an asymptotic result which also applies at higher frequencies.The optical length of one unit cell is = √ 1 a 1 + √ 2 a 2 and from the Bragg condition k = κ a (= Mπ with M being an integer), we define the characteristic phase index which approximately accounts for the phase index n p both at low frequencies (first band) with the m = 0 branch, as well as for the related branches in higher frequency bands.In Figure 1, the dashed line shows the asymptotic result in Eq. ( 4).Branches of higher frequency bands with a comparable phase index have been highlighted in red.In the following we will measure the branch index m with respect to the above asymptotic result, i.e.
The group index is however unambiguously determined, which follows directly from the fact that is independent of m.Thus, the group index does not carry information of the branch index m as is also clearly seen from the example shown in Figure 1.Nevertheless, it is interesting to note that for the highlighted branches (shown in red), the phase index is comparable to the group index as we would expect for a homogeneous material (with no pronounced material dispersion).As we shall see, these are the modes where homogenization makes sense.
In the context of the homogenization hypothesis it is interesting to note that the effective wavelength also has an ambiguity.There is a consensus that λ a (or ka 1) is a natural condition for the structured material to effectively mimic a homogeneous material.However, the field variations inside the material are of course characterized by λ eff and thus one has the somewhat stronger requirement that also λ eff a (or κ a 1).The two criteria only differ by the phase index which could be ignored if n p was not too different from ñp .However, the above discussion shows that the situation can be more complicated for periodic structures.Clearly, since a we may have λ a while at the same time λ eff a for highorder branches with m > 0. A substantial population of highorder branches would clearly be conflicting with the underlying homogenization hypothesis, since the associated field variations would not be slow on the length scale of the lattice.
In the following we examine to which extend such branches are excited in a scattering experiment.

IMPEDANCE CONSIDERATIONS
We consider a plane-wave in air (z < 0) incident on a semiinfinite Bragg stack (z > 0) made from a non-magnetic material.The Fresnel reflection at the interface between two materials with different dispersion relations is In the case of air and a homogenous material with index n p this reduces to S 11 2 = (n p − 1) 2 /(n p + 1) 2 .Adapting this to a homogeneous material represented by the mth branch we would have If ñp is not too different from that of air, then the minimal reflection is associated with m = 0 when λ a.More generally, one would expect that the branch with a phase index closest to ñp would have the best impedance matching, while other branches would be expected less matched to the incident excitation since lim This suggest that the m = 0 branch will by far be the dominating one.However, from the above reasoning, one may not be able to fully neglect the excitation of other branches.If the air medium was replaced by a high-index material with n > ñp , then from an impedance point of view the excitation of other branches could be favored as well.As we shall see, this does to some extend occur, though the m = 0 is still largely responsible for the transmission of the main part of the power.
The predominant role of the m = 0 branch owes to the internal properties of the Bloch modes rather than being a consequence of impedance matching.

EXCITATION OF BLOCH STATES IN A FINITE STRUCTURE
In the following we imagine a Bragg stack of finite extension (L = Na with N being an integer), rather than being of infinite or semi-infinite length as above.We again consider the case where modes are excited by a plane wave incident from the left (z < 0), so that the field is given by while to the right of the structure (z > L), the transmitted field is given by Inside the structure (0 < z < L) we may expand the total field in the basis of the Bloch states, formally giving rise to an expansion of the form [26,27] E ω (z) with each term being a spatial harmonic.Although we have expanded the field as an infinite series of spatial harmonics, we emphasize that all the spatial harmonics must be simultaneously present in order for the total field to satisfy boundary conditions associated with the periodic variation of (z) [26].Thus, for general periodic dielectric structures a particular branch can in principle not be populated independently of the other branches and even though there may be a predominant branch, other branches will be at least weakly populated as well.However, the m = 0 branch with n p ∼ ñp will be the dominating one and the contributions from other branches become relatively weaker if the periodic function (z) does not have pronounced contributions from the higher spatial harmonics.In the appendix we show this explicitly.
In the above field expansion, S 11 and S 21 are the reflection and transmission amplitudes, respectively, while A m 2 and B m 2 carry information on the population of the forward and backward propagating Bloch branches inside the Bragg stack.Solving this set of linear equations (with the aid of the appropriate electromagnetic boundary conditions), in principle provides us with the unknown coefficients.Usually, main attention is paid to the scattering parameters S 11 and S 21 from which the phase index can be extracted [8,23], though with an ambiguity for the real part.In a way, by just calculating the scattering parameters S 11 and S 21 , all information on the population of the different branches m have been traced out, see the summa-tion in Eq. ( 14).On the other hand, the formal decomposition into Bloch states shows how we as a general principle may get access to the branch populations by Fourier transforming the total field inside the periodic structure.In the following we define which allows us to map out the excitation of different Bloch wave vector components and branches in a way closely related to the experimental Fourier-imaging technique [28].Fourier analysis have previously been applied to photonic crystals with emphasis on band-folding effects and negative refraction [29].
Figure 2 shows the results for a stack of N = 10 bi-layers probed in a scattering configuration in free space.The full solution for the electrical field is obtained by solving the wave scattering problem with the aid of a commercially available finite-element method (Comsol MultiPhysics).The integration in Eq. ( 15) is subsequently performed with a built-in integration routine.While being computationally time consuming, the method is powerful in highlighting the excitation and population of the different Bloch dispersion branches.As seen, for low frequencies only the m = 0 branch is excited.However, when going to higher frequencies more branches are moderately excited (difficult to visualize on a linear scale), though the predominant branch is in this case still the m = 0, i.e. the one closest to the asymptotic dispersion relation, Eq. ( 4).Note in particular how several branches become roughly equally populated for frequencies near band edges.
In Figure 3 we now consider the same structure, but with a surrounding environment different from vacuum.The structure has ñp 1.62 and we imbed it in a surrounding medium with = 8.This is equivalent to n p 2.83 corresponding to the linear dispersion relation indicated by the red curve in the left panel of Figure 3. From the simple discussions of impedance matching we should then expect matching to highorder branches (m > 0) at elevated frequencies, such as in the third band.Indeed, we do observe a weak population of the high-order branch.However, the results in Figure 3 still suggest that the branches with a phase index close to ñp are those predominantly being populated.Even though the incident wave may match reasonably with higher branches, those branches are themselves poorly represented in the Bloch solution.The population of other branches than the m = 0 is more clearly seen on a logarithmic intensity scale, as illustrated in the inset.The large reflections at the boundaries to the surrounding media with = 8 cause pronounced Fabry-Perot oscillations which are visible both in the intensity plot (left panel) as well as in the transmission spectrum (right panel).
Though periodic metamaterials have more complex unit cells leading to more complex band diagrams, they still share the same kind of Bloch physics.In particular, scattered fields may be decomposed in Bloch states inside the metamaterial with the excitation conditions (incidence angle, polarization etc.) determining the particular linear combination of Bloch states.The proper branch index may only be identified in situations where the expansion is dominated by a particular Bloch state.The effective negative magnetic response in arrays of highindex rods is particular interesting in this context.The resonance is associated with the first Mie resonance [9], but for the periodic arrangement of the rods the same Mie resonance is also responsible for the first bandgap, (e.g.see [11]).The utilization of the first Mie resonance for negative-µ physics will thus potentially involve excitation of more Bloch branches.

DISCUSSION AND CONCLUSION
We have by explicit analytical solutions of a one-dimensional periodic system shown how the branch ambiguity in the extraction of effective parameters is arising as a direct consequence of the underlying Bloch state physics.Rather than being mathematical artifacts, the different branches are potentially all physically meaningful solutions to the electromagnetic wave problem.However, their mutual importance depends on the experimental context.Meaningful effective parameters can only be attributed in cases where a particular branch of the Bloch states with λ eff a dominates.The simple one-dimensional case illustrates how a scattering experiments predominantly excites m = 0 branches.The mutual population of branches is also influenced by the impedance matching to the incident field, but impedance considerations are by themselves not always sufficient in choosing the branch index.For one-dimensional structures the identification of a single predominant branch may be clear because only a single Bloch band is supported at a given frequency.For more complex two and three-dimensional unit-cell designs, the coexistence of more Bloch bands at the same frequency is not always prohibited, depending on the symmetry of the structure [22].In such cases the identification of a predominant mode and branch is not always meaningful.In addition, spatial dispersion has become an issue of major concern for the concept of homogenizing optical metamaterials [21].We emphasize that while the attenuation coefficient 2κ does not discriminate between branches of the same Bloch band, it could on the other hand serve to select one Bloch band over another.Such a discrimination would make thin slabs behave differently form thick slabs of the same metamaterial.
We believe that the present work should serve as an important input when turning homogenization theories into practise.Branch ambiguities of periodic structures, whether they appear in the S-matrix approach or in field-averaging techniques, should be addressed by for example, examining the mutual population of the different branches when excited in a given experimental configuration.

A HARMONIC EXPANSION
In the following we decompose the dielectric function into a Fourier series where in particular c 0 = eff , as given by Eq. ( 3).Likewise, we write the Bloch wave function as a harmonic expansion Substituting the two expansions into the wave equation we get an infinite set of coupled linear equations in the amplitudes A m .This system of linear equations may formally be written as the following matrix problem where we have used that in general c −m = c * m and c 0 = eff , as given by Eq. (3).In the case of a 1 = a 2 = a/2 the coefficients simplify to Introducing the phase index as defined in Eq. (2b), we have (κ − mG) 2 /k 2 = (n p − mλ/a) 2 .Obviously, c 1 c 2 c 3 . . .c m c m+1 and keeping only the c 1 terms we see that Eq. ( 19) mathematically corresponds to a nearest-neighbor tight-binding model with a parabolic on-site potential centered on the m = 0 site.The first eigenvalue, corresponding to the first Bloch band, is thus associated with an eigenvector which is strongly localized on the m = 0 site for κ in the first Brillouin zone, i.e. |A 0 | 2 |A ±1 | 2 |A ±2 | 2 . . .|A ±m | 2 |A ±(m+1) | 2 .A diagonalization supports that in general the m = 0 branch with n p ∼ ñp plays a pronounced role.This is also found in [29].

FIG. 1
FIG. 1 Bandstructure diagram ω(κ) for a Bragg stack with ε 1 = 1 and ε 2 = 5 while a 1 = a 2 = a/2.The left panel shows the real part of the Bloch wave vector κ = κ + iκ while the right panel shows the corresponding imaginary part.The dashed line in the left panel indicates the asymptotic result with a phase index ñ p , Eq. (4), while branches with a comparable phase index are highlighted by red solid curves.

FIG. 2 1 = a 2 = a/ 2 . 2 FIG. 3
FIG. 2 Plot of P(κ , ω)/ E 0 2 for the plane-wave excitation of a Bragg stack containing N = 10 bi-layers, sandwiched between homogenous materials with = 1 corresponding to air.The Bragg stack has layers with ε 1 = 1 and ε 2 = 5 while a 1 = a 2 = a/2.The superimposed black curves show the Bloch states of the corresponding infinite Bragg stack while the red linear curve shows the dispersion relation in surrounding material.