1. Introduction

In the original sense, metamaterials are bespoke plasmonic, periodic structures designed to manipulate the propagation of light. Their success is rooted in their ability to generate otherwise unavailable electromagnetic material properties leading to real-life applications such as the perfect lens. The popularity of metamaterials is, however, just as much caused by their accessibility through effective material parameters that formally resemble those in natural materials, but with the potential to generate values that are otherwise unattainable [1–3] Perhaps the most notorious such optical characteristic is negative refraction, which has been demonstrated in the microwave regime [4–6] and subsequently for smaller wavelengths in the infrared [7]. However, the range of “unnatural” optical responses extends far beyond negative index materials, from science fiction cloaking to chiral metamaterials exhibiting polarization-dependent birefringence. In addition, they enable a broad range of applications in sensing, waveguiding, and imaging [8–14].

As the structuring at nanometer length-scales required for visible-light metamaterials is still very challenging, it is of the utmost importance to develop a strong theoretical understanding of how the observed optical responses are achieved from particular geometries. Thanks to the computing power that is available today, trial and error strategies are more feasible than ever, as brute-force numerical methods such as finite elements or finite difference discretization is readily available to facilitate the study of hypothetical metamaterials without their fabrication. Optical simulations are powerful tools to predict and understand experimental findings, and can even be employed for machine learning pathways towards bespoke metamaterial functionalities [15]. From a fundamental perspective, they, however, provide little to no insight into the physical mechanism behind a particular metamaterial response. An understanding of the underlying electrodynamic modes can instead offer guidance in the targeted template design for desired optical properties.

For optical metamaterials, theoretical models generally fall into the category of effective medium theories (EMTs), which either work on first principles, accessing limited geometries [16–18], or the EMT parameters are retrieved by analyzing either the metamaterial’s scattering behavior or its microscopic fields [19–21]. Even more complex EMTs cannot provide a full physical description as they, by definition, assume a homogeneous medium with effective material parameters. EMT parameters typically consist of an effective permittivity $\varepsilon_{\textrm{eff}}$ and a permeability $\mu_{\textrm{eff}}$ tensor [22,23]. In some advanced models, a non-local response (wave vector dependence) [24,25] or electromagnetic cross-coupling (chirality) is additionally considered [26–28]. Yet, even very crude EMTs of elliptical meta-atoms, such as the Maxwell-Garnett (MGA) or Bruggeman approximations, have proven to be surprisingly accurate even for metamaterials violating their underlying assumptions, which limit them to low metal fill fractions [29].

Here, we explore a more complete theoretical metamaterial description, inspired by the photonic crystal [30,31] and mechanical metamaterial [32] communities. Instead of a full homogenization, which naturally yields only one mode per polarization direction, this description employs a complex bandstructure, which is the set of complex-valued Floquet wave numbers found for a given (real) frequency and crystal inclination. These modes form a complete solution basis of Maxwell’s equations in a semi-infinite metamaterial domain and can be correlated with EMTs, and scattering and emission experiments alike. The complex bandstructure has been recognized as a powerful tool to understand and model the behavior of optical metamaterials [24,25]. While it can be obtained by parameter retrieval [21] and a scattering or transfer matrix approach [33,34], the former yields little physical insight and the latter is only efficiently applicable to simple layered geometries such as the metal-dielectric multilayer [35], 1D metal-dielectric multilayer [36, 37, 63], or the fishnet structure [38]. In this article, we combine a Greens and a Galerkin approach to transform Maxwell’s equations into a low-dimensional non-linear eigenproblem that is generally valid for any two-component metamaterial and acts on explicit currents in one material domain (generally the metal domain) only. We show that for small frequencies and wave vectors ($\omega,|\boldsymbol {k}|\,{\ll }\,2\pi /a$ with lattice constant $a$) and without explicit restrictions on the volume fill fraction, our method recovers the MGA formulae for cylinders[spheres] arranged on isotropic 2D[3D] lattices. We here define a square or hexagonal lattice in 2D, and a cubic lattice (with arbitrary centering) in 3D as isotropic. Formally, only these lattices have sufficient point symmetry such that no 2D/3D matrix (apart from a multiple of the identity matrix) is invariant under their point group operations [38].

2. Method

2.1 Floquet modes in linear binary media

Consider an arbitrary binary periodic metamaterial with primitive unit cell $\Omega \,{=}\, \Omega_1 {\cup } \Omega_2\,{\in }\, \mathbb {R}^{d}$ (where $d\,{=}\, 2,3$ is the dimension of the lattice and $\Omega_1 {\cap } \Omega_2 \,{=}\, \emptyset$). The respective materials determine the wavelengthfrequency-dependent permittivity $\varepsilon (\boldsymbol {r},\omega ) = \varepsilon_i(\omega )$ for $\boldsymbol {r}\in \Omega_i$, $i=1,2$. Most commonly, such a metamaterial consists of a metal domain and a dielectric background. Using a monochromatic ansatz for a fixed angular frequency $\omega \,{\ne }\,0$, such that $\delta_t \,{\rightarrow }\, -\imath \omega$, the macroscopic Maxwell equations yield a Helmholtz-type wave equation [39] for the electric field $\boldsymbol {E}(\boldsymbol {r})$. Defining the vacuum wave number $k_0 \,{:=}\, \omega / c$ with the vacuum speed of light $c$, and the material wave number $k(\boldsymbol {r}) \,{=}\, k_i \,{:=}\, \sqrt {\varepsilon_i} k_0$ on $\Omega_i$ (with $\Omega_2$ the metal domain) the wave equation takes the form

By the Floquet theorem, the (non-dimensionalized) electric field is of the form

While there is a large number of solvers available that can calculate the real bandstructure and the bulk modes with complex eigenfrequencies, computing evanescent Floquet modes with complex wave numbers is much less straight-forward. We show in the following that Eq. (1) is an ideal starting point to compute evanescent Floquet modes in a computationally affordable and physically meaningful way.

2.2 Floquet modes from a low-dimensional eigenproblem

We begin with a plane-wave expansion of the electric field, which is obtained by expressing $\boldsymbol {u}$ by its $d$-dimensional Fourier series

As we show below, however, a very large number of reciprocal lattice vectors can be afforded numerically, as the computational cost scales linearly with the number of lattice vectors. These lattice vectors only contribute to a lattice sum and not to the dimension of the computationally expensive eigenproblem itself. Nevertheless, a more efficient expansion of $\boldsymbol {E}$ is conceivable, for example using finite elements. The driving currents in Eq. (1) are efficiently expanded through polynomial basis functions $P_\alpha$ with compact support on $\Omega_2$, for convenience scaled by the Bloch phase [42, 64]. In this work, the polynomial basis functions are chosen as monomials of the type $x^{m} y^{n} z^{l}$ with $m+n+l \leq 8$. This very low maximum degree is more than sufficient in the studied cases due to the exponential convergence behavior in the polynomial degrees [43, 44, 65]

Using the two expansions Eq. (2) and Eq. (3), we derive the Floquet eigenproblem in Supplementary Section 1. In summary, we first analytically invert the wave operater $\mathcal {H}$ via Galerkin testing of Eq. (1a) with the plane wave basis. Here lies the advantage of the plane-wave basis, in which $\mathcal {H}$ is diagonal and can be analytically inverted. Formally, the procedure is equivalent to using the spectral lattice Greens function [43]. In the second step, we substitute the result of the first step into Eq. (1b), and Galerkin test with the polynomial basis functions. Using $N$ polynomial basis functions $P_\alpha$, this results in a $3N$-dimensional homogeneous algebraic equation on the $N$-dimensional vector space of polynomial coefficients $c_\alpha$,

In order to understand Eq. (4) as an eigenproblem, we define a coordinate system such that the normal of a chosen crystal inclination $(hkl)$[$(hk)$] for a 3D[2D] lattice is the $z$-direction, which we also refer to as propagation direction. Equation (4) now forms a family of eigenproblems with the in-plane components of the Bloch wave vector $\boldsymbol {\kappa }_\parallel$ and the vacuum wave number $k_0$ as input parameters and $\kappa_z$ as (non-linear) eigenvalue. Generally, a complete solution set for the metamaterial half-space for a given $k_0$ is formed by the solutions for all $\boldsymbol {\kappa }_\parallel$ in the surface[edge] Brillouin zone for $(hkl)$[$(hk)$]. For a scattering problem with a plane-wave incident on a metamaterial slab with $z$ inclination, it evidently suffices to only consider the solutions with $\boldsymbol {\kappa }_\parallel \,{=}\, k_0\sin \theta {(\cos \varphi,\sin \varphi )}^{\phantom {}^{\intercal }}$, with the angle of incidence defined by the polar angle $\theta$ and the azimuth angle $\varphi$. The far-field radiation of a classical dipole embedded in the metamaterial can be obtained using the same idea, but with the complication that in oder to satisfy Maxwell’s equations on the dipole plane, the fields above and below this plane must be matched.

2.3 Numerical solution of the eigenproblem

To numerically solve the non-linear eigenproblem for any set of polynomial basis functions $P_\alpha$, we transform it into a matrix equation,

3. Results and discussion

General solutions of Eq. (4) require a numerical evaluation of Eq. (5) with cut-offs in both the reciprocal lattice vectors (through $N_G$) and the number $N$ of polynomial basis functions. We show in Section 3.1 that analytical solutions of Eq. (4) can be found for cylinder[sphere] packings for small wavenumbers (${\ll }2\pi /a$) and resulting slowly varying driving currents $\boldsymbol {C}$ ($N\,{=}\, 1$). The dispersion relation of our analytical solution matches the MGA prediction exactly. A comprehensive analysis of the full numerical solution for aligned cylinders is presented in Section 3.2. We thereby show that the above approximation produces reliable results for an electric field polarized along the cylinder axis over the whole optical spectrum and even into the near-ultraviolet region. For polarization perpendicular to the cylinder axis, the approximation predicts the mode with lowest imaginary part of the eigenvalue $\kappa_z'$ well for most of the spectrum, apart from the region close to the fundamental dipole Mie resonance, where two modes of the same symmetry classification with comparable $\kappa_z'$ exist.

3.1 Approximate eigenproblem for cylinders and spheres on a lattice

We here consider a metamaterial in the low wavelength limit, where we assume $k_0,\kappa \,{\ll }\,2\pi /a$. In lowest non-vanishing order, this immediately simplifies the inverse Helmholtz operator to $\mathcal {H}^{-1}_{\boldsymbol {\kappa }{+}\boldsymbol {G}}\overset {\boldsymbol {G}{\ne }0}{\approx }k_1^{-1}\hat {\boldsymbol {G}}\otimes \hat {\boldsymbol {G}}$, where $\hat {\boldsymbol {G}}\,{:=}\,\boldsymbol {G}/G$ is the direction of the reciprocal lattice vector. Given that the skin depth of gold and silver at optical frequencies is above 10 nm for optical frequencies and below [45], we further assume that the slowly varying field fully penetrates objects with less than ${\approx } {20}\;\textrm{nm}$ diameter. We therefore consider constant driving currents $\boldsymbol {C}(\boldsymbol {r})\,{=}\, \boldsymbol {c}_0\exp \{\imath \boldsymbol {\kappa }\cdot \boldsymbol {r}\}$ on $\Omega_2$, that is, with $P_0\,{=}\, 1$. With these approximations, the eigenproblem Eq. (4) simplifies to

The dyadic lattice sum $\mathbf {L}$ is a constant matrix, independent of $\kappa$ and $k_0$. While it generally needs to be computed numerically, we show in Supplementary Section 3 that it evaluates to

Let us first consider the conceptually simpler, fully isotropic 3D case, for which $\mathbb {1}_d$ is the 3D identity matrix. Equation (6) is evidently solved if $\boldsymbol {c}_0\,{=}\,\boldsymbol {\kappa }$ (longitudinal mode) or if $\boldsymbol {c}_0\cdot \boldsymbol {\kappa }\,{=}\, 0$ (transverse modes). The formal definition of longitudinal/transverse can be used in any case. Note that the longitudinal/transverse terminology, however, provides an exact description for real-valued $\boldsymbol {\kappa }$ only. Generally, the currents are elliptically polarized for $\kappa_z\,{\notin }\,\mathbb {R}$ and $\boldsymbol {\kappa }_\parallel \,{\neq }\,0$.

For the longitudinal mode, all coefficients are proportional to $\varepsilon_1{-}\bar {\varepsilon }$, which generally has no solution. [47] For the transverse mode, Eq. (6) reproduces the well-known MGA result [16,48]

We now consider a wire metamaterial [49], consisting of aligned metal wires arranged on an isotropic 2D lattice. An example of such a metamaterial was fabricated in [50]. It consists of gold or silver cylinders with radius $R= {10}{nm}$ and lattice constant $a = {30}{nm}$, arranged on a hexagonal lattice with (10) inclination, as illustrated in Fig. 1. We define the wire direction as $x$. Let us first restrict the discussion to a vanishing ordinary (non-Floquet) wave number $\kappa_x\,{=}\, 0$. In this situation, the currents are polarized along the cylinder axis or perpendicular to it, so that a comparison with MGA is possible. Indeed, Eq. (6) is solved if the current modes are longitudinal ($\boldsymbol {c}_0\,{=}\,\boldsymbol {\kappa }$), transverse electric (TE, $\boldsymbol {c}_0\,{=}\,\boldsymbol {e}_x$), or transverse magnetic (TM, $\boldsymbol {c}_0\,{=}\,\boldsymbol {e}_x\,{\times }\,\boldsymbol {\kappa }$). The TE/TM symmetry classification is with respect to the $x\,{\mapsto }\,{-}x$ mirror and only possible for invariant $\boldsymbol {\kappa }$, that is, for $\kappa_x\,{=}\, 0$. As in the 3D case, the longitudinal mode does generally not exist. For the TE mode, Eq. (6) yields the MGA result for infinitely long ellipsoids with field polarization along their long axis [16]

 

Fig. 1. Schematic of a 2D-periodic metamaterial consisting of cylindrical metal wires stacked parallel to the boundary surface with (10) inclination on a hexagonal lattice.

Download Full Size | PDF

The two general analytical solutions of Eq. (6) for arbitrary [51] $\boldsymbol {\kappa }\,{=}\,\boldsymbol {\kappa }_\parallel +\kappa_x\boldsymbol {e}_x$ require considering the full three-dimensional problem. This, however, splits into a 1D problem and a 2D problem. The former yields a quasi-TM mode with $\boldsymbol {c}_0\cdot \boldsymbol {e}_x\,{=}\,\boldsymbol {c}_0\cdot \kappa \,{=}\, 0$ and the dispersion relation $\bar {\varepsilon }_{\textrm{TM}} k_0^{2}\,{=}\,\kappa ^{2}$. For non-zero $\kappa_x$ and $\boldsymbol {\kappa }_\parallel$, the 2D problem generally has only one solution. The corresponding (non-normalized) eigenvector is

3.2 Breakdown of the constant current approximation for TM modes

In this section, we concentrate on the hyperbolic wire metamaterial fabricated in [50], to investigate the accuracy of the approximate MGA predictions Eq. (8) and Eq. (9) above. This metamaterial consists of silver or gold cylinders (in a vacuum background), arranged on a hexagonal lattice as illustrated in Fig. 1. We here focus on the physically richer silver metamaterial with lattice constant $a = {30}\;\textrm{nm}$ and cylinder radius $R = {10}\;\textrm{nm}$ of [50]. More cases with gold wires and different radii are considered in Supplementary Section 4. As we have seen above, the solutions of the non-linear eigenproblem for the 2D wire metamaterial agree exactly with the MGA prediction in the low wavelength limit. If the latter restriction is relaxed, Eq. (4) remains three-dimensional, but the lattice sum becomes $k_0$ and $\boldsymbol {\kappa }$ dependent and needs to be evaluated numerically. For the metamaterial in question in the wavelength range of interest (${300}\;\textrm{nm}\,{<}\,\lambda \,{<}\, {800}\;\textrm{nm}$), we find that the analytical results above are reproduced with extremely good precision. The same holds for the constant current approximation in case of the TE modes, for which Eq. (8) is exact for any practical application. This is probably expected since the electric field is always tangential to the boundary surface between $\Omega_1$ and $\Omega_2$, so that Eq. (8) is simply the weighted volume average of the two permittivities. This formula, which is the same as for a binary Bragg reflector, can be immediately deduced from Fatou’s theorem [54], and approximates the dispersion relation well for frequencies below the fundamental bandgap (as long as $\lambda \gtrapprox 2\sqrt {\bar {\varepsilon }(\lambda )}a$). The resulting silver metamaterial TE bandstructure for $\boldsymbol {\kappa }_\parallel \,{=}\, 0$, for which the MGA and the exact solution are optically indistinguishable, is shown in Fig. 2.

 

Fig. 2. Bandstructure of the 2D wire metamaterial: Real (solid) and imaginary (dotted) part of $\kappa_z$ of the TE mode (blue) and the two TM modes (red and yellow). The Maxwell-Garnett result for TE case coincides exactly with the calculated mode and is thus omitted, the TM case (black) coincides with either TM mode above and below the resonance frequency when the fields in the wires are approximately constant.

Download Full Size | PDF

The physics is naturally richer for the TM modes, where the fundamental dipole Mie resonance (the first pole of the scattering coefficients [55,56]) of the individual cylinders lies in the frequency range of interest. Within a small band close to this dipole resonance, the constant current approach is insufficient. To demonstrate this behavior, we solve the numerical problem of Eq. (5) in the spectral range of interest between the frequencies of $\nu_0 = {375}\;\textrm{THz}$ ($\lambda_0\,{\approx }\, {800}\;\textrm{nm}$) and $\nu_1 = {1000}\;\textrm{THz}$ ($\lambda_1\,{\approx }\, {300}\;\textrm{nm}$). As initial $\kappa_z$ for the Newton procedure, we use the MGA approximation Eq. (9) at $\nu_0$ and $\nu_1$, sufficiently far away from the dipole resonance, and iterate along the spectrum in both directions. Figure 2 shows the bandstructure for monomials up to degree $4$ in $y$ direction and degree $3$ in $z$ direction, denoted pd = (4,3), where the solutions converge sufficiently over the whole spectrum. The bandstructure is not as simple as MGA suggests and there are indeed two TM modes with reasonably small $\kappa_z'$, of which one (labeled TM1) converges to the MGA band for small frequencies, and the other (TM2) converges to the MGA band for large frequencies. The bandstructure for pd = (2,2) and a detailed convergence analysis is provided in Supplementary Section 5.

To better illustrate the breakdown of the constant current approximation, we calculate the well converged numerical currents for pd = (4,4) from the eigenvectors $\boldsymbol {c}$ using Eq. (3). The current intensities and field lines are shown in Fig. 3 for the two TM modes at three different frequencies well below (400 THz), above (900 THz) and at the dipole resonance (838 THz). In the long wavelength limit, the TM1 mode that is well approximated by Eq. (9) exhibits near-constant currents, as expected (Fig. 3(a)). The TM2 mode with relatively large $\kappa_z'$ is, on the other hand, strongly evanescent within the cylinder with strongly curved field lines (Fig. 3(d)), which explains why it is not contained in the constant current model. At the resonance, both modes exhibit varying currents to some degree (Fig. 3(b) and  3(e)). While they are clearly different, with concentrated currents at the left and the right of the cylinder in case of TM1, neither field is well resolved using a constant current approximation. The approximation Eq. (9) indeed moves from the TM1 branch to the TM2 branch between ${\sim }800$ and ${\sim } {850}\;\textrm{THz}$ (Fig. 2). At 900 THz, the TM1 mode is strongly evanescent within the cylinder (Fig. 3(c)) in agreement with the large imaginary part of $\kappa_z$ (Fig. 2). The field profile is generally relatively complex, explaining why the constant current approximation is not able to resolve this mode. The TM2 mode, on the other hand, exhibits a weakly varying field at 900 THz (Fig. 3(f)). The dispersion is consequently well approximated by Eq. (9) in this frequency regime.

 

Fig. 3. Heatmaps of normalized current intensities $|\boldsymbol {C}|(\boldsymbol {r})/\max \{|\boldsymbol {C}|(\boldsymbol {r})\}$ and field lines of the two TM modes TM1 (a-c) and TM2 (d-f) at 400 THz (a,d), 838 THz (b,e) and 900 THz (c,f).

Download Full Size | PDF

Generally, the MGA provides good results at long wavelengths if the radius of the cylinders is small compared to the lattice constant, as expected from the assumptions in its derivation. It is still a good approximation away from resonances for very high fill fractions, far beyond its theoretical limitations. Reasonable results are obtained even at a resonance, as we show in Supplementary Section 4. For gold, where the resonance is damped through losses, the MGA is valid over the whole visible spectrum up to $R/a\,{\sim }\,0.43$ ($\eta \,{\sim }\, {70}{\%}$), while for silver, the threshold is $R/a\,{\sim }\,0.26$ ($\eta \,{\sim }\, {25}{\%}$). At higher fill fractions above the critical radius, the MG prediction thus approaches the TM1 solution at low frequencies and the TM2 solution at high frequencies, and interpolates between the two solutions between 800 nm and 850 nm (Fig. 2). At lower fill fractions, on the other hand, the MG solution approximates a single exact solution. The two TM bands must therefore belong to one connected, topologically non-trivial solution manifold in the $R$-$\omega$ parameter space. We discuss the non-trivial topology, which stems from a band degeneracy with an associated winding number of $2$, in Supplementary Section 6.

4. Conclusion

We have developed a low-dimensional eigenproblem on the vector space of explicit metal currents to solve Maxwell’s Floquet problem in binary metamaterials. Our approach not only provides an intuitive and generalized method to solve for the complex bandstructure, but also justifies the use of the well-known Maxwell-Garnett approximation (MGA) for metal fill fractions far exceeding its original assumptions. We have indeed shown that the permittivities from the MGA are analytically reproduced by our modal solutions in the long wavelength limit for constant currents with no explicit restriction on the volume fill fraction. Our results rigorously apply to perfectly ordered isotropic arrangements, but it seems clear that the they can be extended to less ordered cases as long as a certain amount of order remains, regarding the distance between the individual spheres or cylinders, such as in hyperuniform packings [57].

At the same time, we revealed the necessity of two TM modes in the spectral vicinity of the single-wire dipole resonance of a 2D-periodic nanowire array, one of which approaches the MGA solution below the resonance and the other one above the resonance. We have shown that this is due to a topological transition that occurs in the bandstructure of the TM modes at a critical resonance strength (determined by the metal and dielectric background material and the volume fill fraction), above which the constant current approach fails to provide a sensible solution. The bandstructure of 3D plasmonic sphere packings is of course even richer. Perhaps unsurprisingly, even a dipole approximation, which coincides with with pd = (0,0) in the limit of small radii in our formulation, deviates from the isotropic MGA model close to the single sphere resonance [58]. However, our reults for the wire array suggests that the dipole approximation is hardly sufficient for sphere radii that are comparable to the lattice constant, particularly for silver spheres.

In this paper, we have restricted the discussion to cylinder and sphere packings, for which the inner Galerkin products can be evaluated analytically, but note that our approach is generally valid for an arbitrary connected metal domain $\Omega_2$. While many other sufficiently symmetric shapes of $\Omega_2$ allow an analytical evaluation of the sesquilinear products, an implementation for arbitrary geometries would be useful. This requires an efficient numerical algorithm to compute the associated integrals. Such an algorithm based on a simplex representation of the boundary of $\Omega_2$ in both 2D and 3D has been introduced in [59].

It is straightforward to further extend the method to non-connected domains, with $M$ connected sub-domains $\Omega_{2m}$ such that $\Omega_2\,{=}\,\cup_{m\,{=}\, 1}^{M} \Omega_{2m}$. These can for example be two separate spheres in the primitive unit cell. This simply requires monomial basis functions $P_{\alpha m}$ on each sub-domain. The sesquilinear products in Eq. (4) are then evaluated in the different sub-domains, while the dimension of the NLEVP becomes $3\times N\times M$. Compared to a naive use of a larger number of polynomials, it is expected to even proof more efficient to artificially subdivide topologically connected, but geometrically more complex metal domains, such as in 3D network-like metamaterials [60]. Finally, we emphasize that the sub-domains can even be filled with different materials, thereby generalizing our method to non-binary metamaterials. In Eq. (4), this only requires an $m$-dependent wave number $\delta k_m\,{:=}\, k_1{-}k_{2m}$, which comes at no additional conceptual complexity or computational cost.