© M. Moccia et al., Published by EDP Sciences, 2016

1 Introduction

Nonlocal light-matter interactions [1, 2] are becoming increasingly relevant in a broad variety of application scenarios including dispersion engineering [3], ultrafast nonlinear optical response [4], artificial magnetism [5], additional extraordinary waves [6], enhanced spontaneous emission [7, 8], Dirac-point conical dispersion [9], scattering suppression [10], topological transitions [11], and light-based analog signal processing [12].

In a series of recent investigations [13, 14], we have been concerned with possible extensions of the well-established transformation-optics (TO) framework [15, 16] so as to systematically engineer metamaterials exhibiting desired nonlocal effects. As opposed to the conventional TO formulation, which exploits the form-invariant properties of Helmholtz or Maxwell’s equations with respect to spatial (and, hence, inherently local) coordinate transformations, the basic idea underlying our approach is to apply the coordinate-transformation machinery in the wavevector-frequency phase-space, accessed via spatial Fourier transform. In such phase space, nonlocal effects naturally manifest themselves in terms of wavevector-dependent constitutive properties, and can be associated with nonlinear wavevector transformations. In a series of application examples, we illustrated the insightful correspondences between typical nonlocal effects and the geometrical/analytical characteristics of the associated transformations, which indicate that the powerful geometrically-driven design, typical of conventional TO, is retained by our approach, with the actual synthesis problem reduced to the suitable approximation of analytically-derived constitutive “blueprints”. For instance, we showed that (non)reciprocal effects and the appearance of additional extraordinary waves [6] are directly related to the (non-)centersymmetry and multivaluedness, respectively, of the transformations [13]. Moreover, we derived the analytical properties of classes of transformations that could induce stationary points in the dispersion diagram, and we also showed that frequency-dependent wavevector transformations enable a finer tailoring in the phase space, thereby opening up the possibility to engineer complex spatio-temporal dispersion effects such as Dirac-point conical singularities [14].

In this paper, as a further illustration of the potential of the above approach, we apply it to the engineering of degenerate band edge (DBE) effects. Such exotic dispersion effects, first studied by Figotin, and Vitebskiy [17–21], are eliciting a growing attention (see, e.g., [22–30] and references therein) in view of their potential relevance to diverse applications including slow light, solid-state lasers, quantum-cascade lasers, sensors, optical delay lines, traveling-wave tubes, distributed solid-state amplifiers, and switches. Here, inspired by our nonlocal TO approach, we investigate an alternative metamaterial-based design, which is amenable to a multilayered implementation.

Accordingly, the rest of the paper is structured as follows. In Section 2, we introduce the problem statement and our general design strategy. In Section 3, we outline the synthesis procedure, from the derivation of the ideal constitutive blueprints to the actual multilayered implementation. In Section 4, we illustrate some representative examples, and we validate them via full-wave numerical simulations, by ascertaining the emergence of typical DBE physical “footprints”. Finally, in Section 5, we draw some brief conclusions, and provide few hints for future research.

2 Problem statement and geometry

In a series of influential studies [17–21], Figotin and Vitebskiy investigated the theoretical implications of stationary points in a dispersion relationship. Assuming a time-harmonic exp(−iωt) dependence, and given a general dispersion law ω(β) relating the angular frequency ω and the propagation constant β, a ν-th order stationary point at ω = ω₀ is characterized by$$\frac{{\mathrm{\partial }}^n\omega }{\mathrm{\partial }{\beta }^n}\left({\beta }_0\right)=0, n=1,\dots ,\hspace{1em}\nu -1, \frac{{\mathrm{\partial }}^{\nu }\omega }{\mathrm{\partial }{\beta }^{\nu }}\left({\beta }_0\right)\ne 0,$$(1)with β₀ = β(ω₀). This implies that, within a neighborhood of such point, the dispersion law behaves as$$\omega -{\omega }_0\propto {\left(\beta -{\beta }_0\right)}^{\nu }.$$(2)

The most trivial example of stationary point is the so-called regular band edge, corresponding to ν = 2 in equations (1) and (2), and naturally occurring in periodic structures such as photonic crystals. Less common and more subtle are higher-order cases such as ν = 3 and ν = 4. The former (ν = 3) is associated with the so-called “frozen-modes”, which can occur, e.g., in suitably engineered photonic crystals containing anisotropic material constituents. The latter (ν = 4) constitutes the DBE case of specific interest here, and differs fundamentally from the regular band edge case above, mainly due to the critical contribution of degenerate evanescent modes. The reader is referred to [17–21] for a thorough discussion of the theoretical details and related physical implications.

Here, we limit ourselves to emphasize that DBEs may give rise to giant slow-wave resonances which, by comparison with the regular-band-edge counterparts, tend to exhibit much larger energy accumulation and, consequently, much stronger light-matter interactions. This may lead, for instance, to giant gain enhancement [30].

DBEs can be obtained in photonic crystals containing uniaxial media with suitably tilted optical axes, as an effect of the coupling between two modes with different polarizations [17–21, 30]. Alternatively, they can also be engineered in optical fibers with multiple gratings [22, 24, 27], coupled periodic waveguides [23, 25, 26, 28], and periodically-loaded circular waveguides [29].

Here, we explore a different configuration based on a nonlocal metamaterial. As schematically illustrated in Figure 1a, we consider a 2-D scenario (with geometry and field quantities independent of y) where a transverse-magnetic-polarized plane-wave (with y-directed magnetic field) propagates in a homogeneous space entirely filled by a nonlocal, anisotropic medium. Such medium is assumed as nonmagnetic, and is characterized by a relative-permittivity tensor which depends on the wavevector k. Here and henceforth, the tilde ~ identifies wavevector-dependent quantities.

Figure 1. (a) Problem schematic in the assumed 2-D (x, z) reference system (with geometry and field quantities independent of y). We consider a transversely-magnetic-polarized plane-wave (with wavevector k) propagating in a homogenous space filled up by a nonlocal, nonmagnetic, anisotropic medium characterized by a relative-permittivity constitutive tensor . (b) Qualitative example of ideal dispersion surface [cf. Eq. (3)] pertaining to the assumed constitutive blueprints [cf. Eqs. (4) and (6), with ]. (c) Corresponding dispersion diagram (kx = 0 cut), with the black-cross marker indicating the DBE-type stationary point.

By focusing on a simple uniaxial anisotropy, the arising plane-wave dispersion relationship can be simply written in terms of the relevant components ( and ) as$$\frac{k_x^2}{{\tilde{\epsilon }}_{\mathrm{zz}}\left(\mathit{k}\right)}+\frac{k_z^2}{{\tilde{\epsilon }}_{\mathrm{xx}}\left(\mathit{k}\right)}=\frac{{\omega }^2}{c^2},$$(3)where k_x and k_z denote the x- and z-domain wavenumbers, respectively, and c is the speed of light in vacuum.

In what follows, inspired by our previously introduced nonlocal TO approach, we derive the ideal blueprints for the nonlocal medium so as to induce a DBE stationary point in the dispersion relationship (3). Moreover, we also address the synthesis of a multilayered metamaterial implementation that suitably approximates (within a neighborhood of the DBE stationary point) such ideal blueprints. Though still based on a periodic multilayer, our implementation is substantially different from those in references [17–21].

3 Synthesis procedure

3.1 Ideal constitutive blueprints

In reference [14], within the framework of nonlocal TO, we derived a class of wavevector transformations that could induce stationary points (of arbitrary order) in the dispersion relationship. In particular, we explored in detail the case ν = 3 in equations (1) and (2), corresponding to the “frozen-mode” regime. A similar approach could be applied to deal with the DBE (ν = 4) case of interest here. For brevity, we do not repeat here the analytical developments, which can be found in reference [14]. In fact, as discussed in reference [14], while providing a systematic, constructive derivation of the required material blueprints, our approach is not the only possible one, as the problem does not admit a unique solution. Moreover, for the particular 2-D, coordinate-separable and non-magnetic structure of the arising constitutive blueprints, their interpretation is rather straightforward.

Paralleling [14], we consider the following class of coordinate-separable material blueprints (henceforth identified by the superscript “BP”):$${\tilde{\epsilon }}_{\mathrm{xx}}^{\left(\mathrm{BP}\right)}\left(k_z\right)=\frac{k_z^2}{\tilde{Q}\left(k_z\right)}, {\tilde{\epsilon }}_{\mathrm{zz}}^{\left(\mathrm{BP}\right)}\left(k_x\right)={\epsilon }_{\mathrm{zz}0},$$(4)where ε_zz0 is a constant, and is a polynomial function. It can be readily verified, from equation (3), that$$\frac{{\mathrm{\partial }}^n\omega }{\mathrm{\partial }{\left(k_z\right)}^n}=\sum _{m=1}^n{\tilde{\mathrm{\Omega }}}_m\left(k_x,k_z\right)\frac{{\mathrm{\partial }}^m\tilde{Q}\left(k_z\right)}{\mathrm{\partial }{\left(k_z\right)}^m},$$(5)with denoting spectral functions. Therefore, by assuming propagation along the z-direction (i.e., k_x = 0, k_z = β), the desired DBE condition in equation (1) (with ν = 4) can be attained by choosing$$\tilde{Q}\left(k_z\right)=\gamma {\left(k_z-{\beta }_0\right)}^4+\tilde{R}\left(k_z\right)+k_0^2,$$(6)where γ is a nonzero real constant, and k₀ = ω₀/c = 2π/λ₀ denotes the wavenumber in vacuum at the design angular frequency (with λ₀ denoting the corresponding wavelength). Moreover, is an arbitrary function that vanishes (together with its first four derivatives) at k_z = β₀,$$\tilde{R}\left(k_z\right)=\sum _{m>4}{r}_m{\left(k_z-{\beta }_0\right)}^m,$$(7)with r_m denoting arbitrary coefficients.

As anticipated, the particularly simple analytical structure of the constitutive blueprints (4) [together with Eq. (6)] makes it possible to readily verify, by simple inspection, that the desired DBE conditions are fulfilled. In particular, via straightforward algebra, the coefficient γ in equation (6) can be directly related to the nonzero derivative term in the dispersion equation [cf. Eq. (1)],$$\gamma =\frac{{\omega }_0}{12c^2}\frac{{\mathrm{\partial }}^4\omega }{\mathrm{\partial }{\beta }^4}\left({\beta }_0\right).$$(8)

Figure 1b qualitatively shows the dispersion surface [cf. Eq. (3)] pertaining to the above class of constitutive blueprints (with ). As it can be observed, the quartic nature of the dispersion relationship gives rise to two modal branches, which degenerate at the design angular frequency, thereby yielding the sought DBE-type stationary point in the dispersion diagram, as more clearly visible in the k_x = 0 cut shown in Figure 1c.

It is worth highlighting that the class of constitutive blueprints in equations (4) and (6) is idealized, and does not necessarily fulfill the physical feasibility conditions (e.g., causality) in the entire (k, ω) phase space (see also the discussion in Ref. [14]). Likewise, the dispersion characteristics in Figures 1b and 1c should only be intended as conceptual illustrations, since explicit temporal dispersion (i.e., explicit ω-dependence) is neglected.

3.2 Multilayered metamaterial synthesis

Having determined a class of ideal constitutive blueprints that yields the desired dispersion effects, as in reference [14], the synthesis problem can be formulated as finding a physical metamaterial structure whose (nonlocal) effective constitutive parameters suitably approximate the targeted blueprints. It is important to stress that, in our problem, such approximation is needed only within a limited region of the phase space, i.e., in the vicinity of the DBE stationary-point (k_x = 0, k_z = β₀, ω = ω₀). As discussed in reference [14], this generally renders the synthesis problem tractable, even though the assumed blueprints may appear as unphysical if extended to the entire phase space (see also the discussion above).

In view of the uniaxial character and analytical (rational) structure of the assumed blueprints in equations (4) and (6), multilayered metamaterials represent a particularly suited implementation. In the recent technical literature, several approaches have been proposed for the derivation of nonlocal effective models of metallo-dielectric multilayered metamaterials (see, e.g., [31–35]). Here, we utilize a generalized version of the approach proposed by Elser et al. [31], which determines the approximate nonlocal effective model in such a way the corresponding dispersion relationship matches the Taylor expansion of the exact one (around a desired phase-space point) up to a given order. First, we note that, in view of equation (6), the constitutive blueprints (4) that we need to approximate contain odd powers of k_z. The implied symmetry breaking (in the z-direction) could be physically induced in various ways, e.g., by considering uniaxially-anisotropic material constituents with tilted optical axes and/or gyrotropic (nonreciprocal) materials. Here, we consider this latter case, which allows us to directly exploit the results in reference [14]. Accordingly, we consider a generic periodic multilayer, stacked along the x-direction, and with a unit-cell composed of nonmagnetic, homogeneous, isotropic or gyrotropic layers. Hence, the p-th material constituent can be generally characterized by a relative permittivity tensor$${\underline{\underline{\epsilon }}}_p=\left[\begin{array}{lll}{\epsilon }_p& 0& i{\kappa }_p\\ 0& 1& 0\\ -i{\kappa }_p& 0& {\epsilon }_p\end{array}\right],$$(9)with κ_p = 0 yielding the isotropic case. It can be shown (see [14] for details) that, in the vicinity of the DBE stationary-point of interest k_x = 0, k_z = β₀, ω = ω₀, the above class of metamaterials can be described in terms of the nonlocal effective parameters$${\tilde{\epsilon }}_{\mathrm{xx}}^{\left(\mathrm{eff}\right)}\left(k_z;\underline{\alpha }\right)\approx \frac{2-{\chi }_0\left(\underline{\alpha }\right)}{k_0^2\stackrel{4}{\sum _{l=1}}{k}_z^{l-2}{d}^l{\chi }_l\left(\underline{\alpha }\right)}, {\tilde{\epsilon }}_{\mathrm{zz}}^{\left(\mathrm{eff}\right)}\left(k_x;\underline{\alpha }\right)\approx \frac{2-{\chi }_0\left(\underline{\alpha }\right)}{k_0^2{d}^2\left(1-\frac{k_x^2{d}^2}{12}\right)},$$(10)where d is the multilayer period (unit-cell thickness), the array symbolically denotes the set of geometrical and constitutive parameters of the multilayer (with the explicit dependence on ω₀ and β₀ omitted for notational simplicity), and the functions χ_l are chosen so that arising dispersion equation matches the Taylor series (about k_x = 0, k_z = β₀, ω = ω₀) of the exact one up to the second-order in k_x and fourth-order in k_z.

The synthesis problem can thus be cast as finding an “optimal” multilayer parameter set , so that the nonlocal effective model (10) approximately matches the desired constitutive blueprints [cf. Eqs. (4) and (6)] in the vicinity of the DBE stationary point. To further relax the synthesis problem, we only set and leave as degrees of freedom the parameters γ and β₀ in the blueprints, whereas, though not directly relevant in the assumed incidence conditions, the parameter ε_zz0 in equation (4) is obtained by simple matching at k_x = 0,$${\epsilon }_{\mathrm{zz}0}={\tilde{\epsilon }}_{\mathrm{zz}}^{\left(\mathrm{eff}\right)}\left(k_x=0;\underline{\alpha }\right)\approx \frac{2-{\chi }_0\left(\underline{\alpha }\right)}{k_0^2{d}^2}.$$(11)

We are thus led to an optimization problem, i.e., finding the multilayer parameter set as well as the blueprints degrees of freedom γ and β₀ that minimize a suitable error functional,$${\left[{\underline{\alpha }}^{\left(\mathrm{opt}\right)};{\gamma }^{\left(\mathrm{opt}\right)},{\beta }_0^{\left(\mathrm{opt}\right)}\right]=\underset{\underline{\alpha },\gamma ,{\beta }_0 }{\mathrm{argmin}}\sum _{j=1}^J{w}_j\left|{\tilde{\epsilon }}_{\mathrm{xx}}^{\left(\mathrm{BP}\right)}\left(k_{\mathrm{zj}};\gamma ,{\beta }_0\right)-{\tilde{\epsilon }}_{\mathrm{xx}}^{\left(\mathrm{eff}\right)}\left(k_{\mathrm{zj}};\underline{\alpha }\right)\right|}^2,$$(12)with k_zj and w_j, j = 1, …, J denoting discrete k_z-samples around k_z = β₀ and positive weight coefficients, respectively.

For the minimization of the error functional (12), we found that a combination of a standard differential-evolution algorithm and the Nelder-Meald method (see Appendix for details), while not necessarily guaranteeing the convergence to the global minimum, usually provided reasonably good results.

4 Representative results

As an illustrative example of application, we consider a multilayered metamaterial implementation with a four-layer unit-cell, containing three isotropic layers and a gyrotropic one (see the inset in Figure 2b). Since we are only interested in a very narrow frequency range (ω ≈ ω₀), we neglect the material dispersion of the constituents. Moreover, for a better illustration of the phenomena, we also neglect material losses.

Figure 2. (a), (b) Effective constitutive parameters at the design angular frequency ω = ω0. Red-dashed curves represent the optimized blueprints, obtained from equations (4) and (6), with , , and β0 = 0.46k0. Blue-solid curves pertain to the nonlocal effective model of the synthesized multilayered metamaterial [four-layer unit cell shown in the inset of panel (b)] with ε1 = 7.089, given in equation (9) (with ε2 = −1.104 and κ2 = 0.0365), ε3 = 3.35, ε4 = −1.957, d1 = 0.252λ0, d2 = 0.0815λ0, d3 = 0.102λ0, d4 = 0.0648λ0. The parameter matching is enforced at kx = 0 and within the (yellow-shaded) range 0.9β0 ≤ kz ≤ 1.1β0.

Figure 2 illustrates the synthesis results (in terms of nonlocal effective constitutive parameters) obtained by choosing, in the error functional (12), J = 21 k_z-samples equispaced within the interval [0.9β₀, 1.1β₀], with uniform weights (w_j = 1).

More specifically, for the optimized parameters given in the caption, Figure 2a compares the blueprint and synthesized components, over the relevant spectral range. A fairly good agreement is observed, especially within the (yellow-shaded) range [0.9β₀, 1.1β₀] where the parameter matching was directly enforced. Though not directly relevant for the assumed k_x = 0 propagation direction, Figure 2b shows the corresponding components, which are exactly matched only at k_x = 0 [cf. Eq. (11)], but do not substantially differ over a reasonably wide k_x-range.

From the optimized multilayer parameters given in Figure 2 caption, we note the presence of both positive and negative values of the permittivities. This is expectable, as strong nonlocal effects in multilayer configurations typically stem from the excitation and coupling of surface-plasmon polaritons (SPPs) supported at the interfaces separating the negative- and positive-permittivity layers [6–12].

For the optimized parameter configuration, Figure 3 shows the dispersion diagram for k_x = 0, numerically computed via a rigorous transfer-matrix-based approach (see [14] for details), in the vicinity of the design angular frequency ω₀, from which the attained DBE condition (at k_z = β₀ = 0.46k₀) is clearly visible. Also visible is the nonreciprocal character, manifested in the center-symmetry breaking; albeit not strictly required, this is a direct consequence of our design strategy relying on a gyrotropic material constituent (see the discussion in Sect. 3). The inset shows a magnified detail around the DBE stationary point, which highlights the good local agreement with the blueprint prescription.

Figure 3. Numerically-computed dispersion diagram (blue-solid curve) for the synthesized multilayered metamaterial (with parameters as in Figure 2 caption), for kx = 0 in the vicinity of the design angular frequency ω0. The inset shows a magnified view around the DBE stationary point, compared with the blueprint prediction (red-dashed curve).

Figure 4 shows instead a few representative (numerically-computed) equi-frequency contours pertaining to the synthesized multilayered metamaterial. It can be observed that, as the angular frequency approaches the design value ω₀, the equi-frequency contours tend to shrink and loose their center-symmetric character, and eventually degenerate to the DBE point.

Figure 4. Numerically-computed equi-frequency contours for the synthesized multilayered metamaterial (with parameters as in Figure 2 caption), at angular frequencies ω = 0.9ω0, ω = 0.99ω0, ω = 0.995ω0, ω = 0.999ω0 (red, green, blue, magenta curves, respectively). At ω = ω0, the equi-frequency contour degenerates to the DBE stationary point (indicated with a black-cross marker).

As a further, independent validation of the above synthesis, we now consider a half-space made of the synthesized multilayered metamaterial under transverse-magnetic plane-wave illumination (at the design angular frequency), normally-impinging along the positive z-direction from a vacuum half-space. Figure 5 shows a finite-element -computed (see Appendix for details) magnetic field map, over the multilayer unit-cell, which exhibits the typical physical “footprints” associated with the DBE regime [17–21]. As it can be observed, the field is transmitted (with small reflection) in the metamaterial half-space, where it gets converted into an extended mode with growing amplitude. At a distance of about ~5λ₀ from the vacuum-metamaterial interface, the amplitude reaches a maximum value that is over a factor 1000 larger than the incident one, and then it starts decreasing.

Figure 5. Finite-element-computed magnetic-field map, at the design angular frequency ω0, over a unit cell of the synthesized multilayered metamaterial (with parameters as in Figure 2 caption). A half-space configuration is assumed, with plane-wave excitation normally impinging from vacuum. The false-color-scale map indicates the field magnitude normalized with respect to the incident one. The field map is superimposed on the unit-cell (shown in trasparency).

Figure 6 shows instead the magnetic energy density (normalized with respect to the incident one, and averaged across the unit cell)$${\overline{W}}_m=\frac{1}{{\left|{\tilde{H}}_y^{\left(\mathrm{in}\right)}\right|}^{2}d}{\int }_0^d{\left|{\tilde{H}}_y\left(x,z_0\right)\right|}^2\mathrm{d}x,$$(13)as a function of the angular frequency, at a distance z₀ = 5.25λ₀ from the vacuum-metamaterial interface (where the resonant transmitted field is maximum, cf. Figure 5). As it can be observed, this quantity is strongly peaked at the design angular frequency ω₀, reaching a value that is over five orders of magnitude larger that the incident one, and then rapidly decreasing to negligible values for ω > ω₀ (i.e., in the bandgap). As evidenced in the magnified detail displayed in the inset, this is qualitatively consistent with the singular behavior predicted theoretically [17–21]$${\overline{W}}_m\propto {\left|\omega -{\omega }_0\right|}^{-1/2}, \omega \le {\omega }_0.$$(14)

Figure 6. Finite-element-computed normalized magnetic density energy averaged over a unit cell [cf. Eq. (13)] of the synthesized multilayered metamaterial (with parameters as in Figure 2), at z = 5.25λ0, as a function of the angular frequency (blue solid curves). The inset shows a magnified detail around the design angular frequency, and compares it with the theoretically predicted behavior [red-dashed curve, cf. Eq. (14)].

The above results clearly indicate that our designed multilayered metamaterial exhibits the typical physical footprints that characterize a DBE-type giant slow-wave resonance. Within this framework, a few remarks are in order. First, though still based on a multilayered structure containing anisotropic constituents, our proposed configuration differs substantially from those in references [17–21]. By comparison with these, our configuration is 2-D, as we do not rely on modes with different polarizations. Moreover, we assume a propagation direction that is parallel, rather than orthogonal, to the layer interfaces. Accordingly, the underlying physical mechanism is also different and, in our case, it essentially exploits the strong nonlocality arising from the excitation and coupling of SPP modes at the interfaces separating the negative- and positive-permittivity layers. Finally, although the main focus of this prototype study was on a simple proof-of-concept demonstration of the phenomenology, and we did not place special emphasis on application-oriented aspects, we did make sure that the material parameters were constrained within realistic bounds. In particular, the positive relative-permittivity values considered in our synthesis range from nearly 3 to about 7, while those of the negative-permittivity and gyrotropic materials are consistent with configurations already considered in the literature (see, e.g., [9, 36]).

5 Conclusions

In conclusion, this prototype study has shown the possibility to engineer DBE-type giant slow-wave resonances by harnessing the strong nonlocal effects that can occur in multilayered metamaterials featuring positive- and negative-permittivity constituents. Starting from ideal nonlocal constitutive blueprints (inspired by our previously-developed nonlocal TO extension), we have synthesized a multilayered metamaterial that (as verified via full-wave numerical simulations) exhibited the typical DBE physical footprints.

Our results may provide an alternative route for the engineering of DBE effects, and provide further evidence of the applicability and versatility of nonlocal-TO concepts.

Current and future investigations are aimed at exploring alternative implementations, not relying on gyrotropic (nonreciprocal) materials, and/or directly exploiting the material dispersion effects. Also of great interest is the use of gain materials, for loss compensation and/or possible applications to lasing.

References

Appendix

Details on the synthesis and numerical simulations

The minimization of the error functional in equation (12) was carried out by means of a Python-based implementation of the differential-evolution optimization algorithm available in the SciPy optimization library [37]. More specifically, we considered a population size of 100 elements (with random initialization), the “best1bin” strategy, a dithering of the mutation constant within the interval (0.5, 1), and a recombination constant equal to 0.6. We also found that slight improvements could be obtained by “polishing” the best element of the final population via the Nelder-Meald method. With reference to the specific example presented, the constitutive parameters of layer #2 (gyrotropic medium) were fixed to values already utilized in the literature [36]; for the remaining layers, the relative permittivities were assumed to vary with the realistic ranges 1 ≤ ε_1,3 ≤ 11 and −10 ≤ ε₄ ≤ −0.1. For the layer thicknesses, a variation range 0 < d_k ≤ 0.3λ₀, k = 1, …, 4 was assumed. Moreover, to avoid propagation of higher-order Bragg modes, the constraint d < λ₀/2 was enforced for the total unit-cell thickness; this was implemented by adding a penalty term in the error functional (12) (equal to 100 if the constraint was not satisfied, and zero otherwise). For the remaining two optimization parameters β₀ and γ, the variation ranges −2k₀ ≤ β₀ ≤ 2k₀ and were set.

The field map (Figure 5) and magnetic energy density (Figure 6) were computed by means of the finite-element-based commercial software package COMSOL Multiphysics [38]. More specifically, a single unit-cell was considered, extending 10λ₀ along the z-direction and terminated by periodic boundary conditions along x. The structure was excited by a plane wave impinging from a 2λ₀-long section of vacuum. In order to simulate a semi-infinite structure, the other side was terminated by a 60λ₀-long lossy section (not shown in Figure 5) with a linearly-increasing loss profile. The structure was discretized with a maximum mesh size of λ₀/100 (resulting into about 6.3 million degrees of freedom), and the MUMPS direct solver (with default parameters) was utilized.

All Figures

Figure 1. (a) Problem schematic in the assumed 2-D (x, z) reference system (with geometry and field quantities independent of y). We consider a transversely-magnetic-polarized plane-wave (with wavevector k) propagating in a homogenous space filled up by a nonlocal, nonmagnetic, anisotropic medium characterized by a relative-permittivity constitutive tensor . (b) Qualitative example of ideal dispersion surface [cf. Eq. (3)] pertaining to the assumed constitutive blueprints [cf. Eqs. (4) and (6), with ]. (c) Corresponding dispersion diagram (kx = 0 cut), with the black-cross marker indicating the DBE-type stationary point.

In the text

Figure 2. (a), (b) Effective constitutive parameters at the design angular frequency ω = ω0. Red-dashed curves represent the optimized blueprints, obtained from equations (4) and (6), with , , and β0 = 0.46k0. Blue-solid curves pertain to the nonlocal effective model of the synthesized multilayered metamaterial [four-layer unit cell shown in the inset of panel (b)] with ε1 = 7.089, given in equation (9) (with ε2 = −1.104 and κ2 = 0.0365), ε3 = 3.35, ε4 = −1.957, d1 = 0.252λ0, d2 = 0.0815λ0, d3 = 0.102λ0, d4 = 0.0648λ0. The parameter matching is enforced at kx = 0 and within the (yellow-shaded) range 0.9β0 ≤ kz ≤ 1.1β0.

In the text

	Figure 3. Numerically-computed dispersion diagram (blue-solid curve) for the synthesized multilayered metamaterial (with parameters as in Figure 2 caption), for kx = 0 in the vicinity of the design angular frequency ω0. The inset shows a magnified view around the DBE stationary point, compared with the blueprint prediction (red-dashed curve).
In the text

	Figure 4. Numerically-computed equi-frequency contours for the synthesized multilayered metamaterial (with parameters as in Figure 2 caption), at angular frequencies ω = 0.9ω0, ω = 0.99ω0, ω = 0.995ω0, ω = 0.999ω0 (red, green, blue, magenta curves, respectively). At ω = ω0, the equi-frequency contour degenerates to the DBE stationary point (indicated with a black-cross marker).
In the text

Figure 5. Finite-element-computed magnetic-field map, at the design angular frequency ω0, over a unit cell of the synthesized multilayered metamaterial (with parameters as in Figure 2 caption). A half-space configuration is assumed, with plane-wave excitation normally impinging from vacuum. The false-color-scale map indicates the field magnitude normalized with respect to the incident one. The field map is superimposed on the unit-cell (shown in trasparency).

In the text

Figure 6. Finite-element-computed normalized magnetic density energy averaged over a unit cell [cf. Eq. (13)] of the synthesized multilayered metamaterial (with parameters as in Figure 2), at z = 5.25λ0, as a function of the angular frequency (blue solid curves). The inset shows a magnified detail around the design angular frequency, and compares it with the theoretically predicted behavior [red-dashed curve, cf. Eq. (14)].

In the text