Dispersion engineering via nonlocal transformation optics

MASSIMO MOCCIA, GIUSEPPE CASTALDI, VINCENZO GALDI,* ANDREA ALÙ, AND NADER ENGHETA Waves Group, Department of Engineering, University of Sannio, I-82100 Benevento, Italy Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, Texas 78712, USA Department of Electrical and Systems Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA *Corresponding author: vgaldi@unisannio.it


INTRODUCTION
Since the seminal papers by Leonhardt [1] and Pendry et al. [2], transformation optics (TO) has been a major catalyst for the formidable advances in the field of electromagnetic (EM) metamaterials, by providing a systematic and versatile tool for the conception and design of artificial media with given field-manipulation capabilities (see [3] for a recent review).In essence, by exploiting the form-invariance of Maxwell's equations with respect to coordinate transformations, TO allows us to systematically tailor the spatial constitutive profile of a material so as to precisely manipulate the field distribution and propagation according to a given coordinatedistorted reference frame.The implied paradigm shift is the separation of the conceptual design from the actual material synthesis, with the former essentially driven by geometrical intuition and considerations, and the latter reducing to a suitable approximation of given ideal constitutive "blueprints."More recently, this framework has been extended to other disciplines beyond EM [4], such as acoustics, elastodynamics, thermodynamics, and even multiphysics scenarios [5].
Maintaining the focus on EM scenarios, several TO extensions have been proposed to broaden the class of materials and effects originally achievable, so as to accommodate, e.g., bi-anisotropic, nonreciprocal, nonlinear, nonlocal, non-Hermitian, epsilon/ mu-negative, and time-varying scenarios (see, e.g., [6][7][8][9][10][11][12][13][14][15] for a sparse sampling).In particular, in [14], we proposed an extension that enables the manipulation and control of nonlocal lightmatter interactions.Here and henceforth, "nonlocality" is intended as spatial dispersion, i.e., the dependence of the electric field displacement (or magnetic induction) at one point on the field applied in neighbor points (and, in view of causality, time instants) [16,17].These effects are negligible (by comparison, e.g., with explicit temporal-dispersion effects) in most natural materials, but they are assuming increasing importance in electrodynamics and optics, especially in view of their relevance in a variety of metamaterial applications, ranging from dispersion engineering [18] to ultrafast nonlinear optics [19].
In our proposed nonlocal-TO (NLTO) framework, the conceptual design of the desired response is guided by the direct geometrical manipulation of the dispersion characteristics in the reciprocal space of wavenumbers, and the needed (ideal) constitutive blueprints are systematically derived in terms of wavevectordependent constitutive tensors, and subsequently approximated via parameter matching with nonlocal effective models of physical metamaterial structures.Along with the increasing availability of nonlocal effective models of metamaterial classes [20][21][22][23][24][25][26][27][28][29][30][31], we expect this approach to establish itself as an attractive option for the systematic and versatile design of artificial materials with broad field-manipulation capabilities.Among the possible applications, we have already addressed the engineering of additional extraordinary-wave phenomena as well as nonlocal signal processing [14].In particular, the latter may open up intriguing venues in the recently emerged field of "computational metamaterials" [32].
In this paper, we present some further extensions and applications of the NLTO approach, with specific focus on dispersion engineering and an essentially threefold scope.First, we consider wavenumber transformations with explicit frequency dependence, which bring about additional degrees of freedom and design flexibility by giving access to the full frequency-wavenumber phase space.Second, we highlight the key role of complex-valued transformations in controlling the density of states.Third, we explore noncentrosymmetric transformations in order to engineer nonreciprocal effects.For illustration, we apply the above concepts and tools to the engineering of exotic dispersion effects, including one-way propagation, stationary points (e.g., "frozen" modes), and Dirac-point conical singularities.In all examples, we start from the geometry-driven conceptual design, proceed with the systematic derivation of the ideal constitutive blueprints, and finally address their physical implementation in terms of multilayered metamaterials.
Accordingly, the rest the paper is organized as follows, with a body of technical details and ancillary results relegated to the supplementary material (see Supplement 1).Section 2 contains a brief summary of the extended NLTO framework and its geometrical interpretation.Section 3 illustrates three end-to-end application examples, including their numerical (full-wave) validations.Finally, Section 4 provides some brief conclusions and hints for future research.

SUMMARY OF EXTENDED NLTO APPROACH A. Theoretical Framework and Proposed Extensions
Our extended NLTO framework is schematically outlined in Fig. 1.Throughout the paper, vector and (second-rank) tensor quantities are identified by boldface symbols and double underlines, respectively, whereas reciprocal-space (wavenumber) quantities are identified by a tilde ∼.As in conventional (i.e., spatial-domain) TO [2], we start considering an auxiliary space (Fig. 1, left panel) identified by primed coordinates r 0 ≡ x 0 ; y 0 ; z 0 and filled by a homogeneous medium (say vacuum), where a given distribution of electric (J 0 ) and magnetic (M 0 ) sources radiates a time-harmonic EM field (E 0 , H 0 ) with suppressed exp−iωt dependence.Next, unlike conventional TO, we access the associated reciprocal (wavenumber) space (Fig. 1, center panel) k 0 ≡ k 0 x ; k 0 y ; k 0 z , via the spatial Fourier transform We subsequently introduce a vector coordinate transformation e F to a new reciprocal space k ≡ k x ; k y ; k z (Fig. 1, right panel), where the recasting in terms of a metric tensor ̳Λ is a notationally expedient device whose convenience will become clearer hereafter, and T denotes the transpose.In spite of the formal analogies with the framework previously introduced by us in [14], we highlight the explicit dependence on the (radian) frequency ω that we have now introduced in the mapping (2).Moreover, we allow the mapping to be generally complex-valued.These represent key differences, which considerably broaden the range of achievable effects and functionalities.
Fig. 1.Schematic of the extended NLTO framework.Starting from an auxiliary space (left panel) r 0 ≡ x 0 ; y 0 ; z 0 filled by a homogeneous medium (vacuum) with a given distribution of electric and magnetic sources (J 0 and M 0 , respectively) radiating an EM field (E 0 , H 0 ), the associated reciprocal space (center panel) k 0 ≡ k 0 x ; k 0 y ; k 0 z is accessed via spatial Fourier transform.Via a frequency-dependent wavenumber transformation, this auxiliary reciprocal space is mapped onto a new reciprocal space k ≡ k x ; k y ; k z (right panel), with the transformed fields and sources related to the original ones via Eqs.(3).The implied nonlocal field-manipulation effects can be alternatively interpreted as pertaining to an undistorted reciprocal reference frame k associated with a physical space r filled with a homogeneous, anisotropic "transformation medium" with wavevector-dependent relative permittivity ̳ ε and permeability ̳ μ tensors given by Eq. ( 4).For notational compactness, the frequency dependence is not shown explicitly.
Paralleling conventional TO (see Supplement 1 for details), the transformed fields ( e E, e H) and sources ( e J, e M) in the mapped reciprocal space k [and associated, via Eq.( 1), spatial domain r ≡ x; y; z] can be related to the original ones (in the k 0 space) via The anticipated convenience in recasting the mapping (2) with the aid of the metric tensor e ̳ Λ appears clear now, which renders the relationships in Eqs. ( 3) and ( 4) formally analogous to those encountered in the conventional (spatial-domain) TO framework.Also retained in the analogy is the aforementioned separation between the conceptual design and the actual metamaterials synthesis, with Eqs.(3a) and (3b) allowing us to systematically design a desired nonlocal field-manipulation effect guided by intuitive geometrical considerations in a distorted reciprocal reference frame, and with Eq. ( 4) explicitly providing the corresponding constitutive blueprints of an ideal nonlocal medium yielding such a response.The synthesis problem can therefore be posed as approximating such ideal blueprints in terms of a physical metamaterial structure characterized by suitable nonlocal effective constitutive parameters ̳ ε eff k; ω; α and ̳ μ eff k; ω; α, with the parameter-array α embedding the physical and geometrical properties of the constituents.From the mathematical viewpoint, this can be formulated as a parameter-optimization problem, with the superscript (TO) tagging the ideal constitutive blueprints, and ‖ • ‖ 2 S indicating a suitable error metric defined over a region S of the k; ω phase space.
As is well known [16,17], basic physical properties (e.g., passivity, reciprocity, causality) dictate some restrictions in the structure of the constitutive tensors describing physically feasible nonlocal media, which translate more or less directly into restrictions on the metric tensor ̳Λ, and hence on the coordinate transformation e F (see Supplement 1 for details).For instance, it can be readily verified that the center-symmetry condition ̳Λk; ω ̳Λ−k; ω guarantees reciprocity.Clearly, the physical feasibility of the constitutive blueprints, as well as a reasonably good "structural" matching with available classes of nonlocal effective models, represent necessary (but not sufficient) conditions for the wellposedness of the synthesis problem in Eq. ( 5).Within this framework, it is also worth stressing that tailoring the nonlocal response over broad regions of the k; ω phase space remains a formidable task.Nevertheless, in many scenarios of practical relevance, several degrees of freedom may be allowed in the blueprints, and the tailoring may be required only within limited regions of the phase space, thereby considerably reducing the synthesis burden.
In what follows, for simplicity of illustration, we refer to an (x, z) two-dimensional (2D) scenario (with y-independent geometry and quantities) and transverse-magnetic (TM) polarization (i.e., y-directed magnetic field), and consider a rather general class of coordinate-separable transformations (with e P, e Q, and W denoting real-valued polynomial functions) which yields a reasonably good compromise between versatility and practical implementability.In particular, by exploiting the degrees of freedom available in the choice of the metric tensor ̳Λ (see [14] and Supplement 1 for details), we obtain a class of effectively nonmagnetic (e μ yy 1) transformation media characterized by relative-permittivity components assuming particularly simple variable-separated rational forms, viz., In Supplement 1, we show that such constitutive blueprints are structurally compatible with nonlocal effective models typical of multilayered metamaterials (along the lines of [14,20]), and reformulate accordingly the parameter-optimization problem in Eq. ( 5).

B. Geometrical Interpretation
As previously mentioned, our NLTO approach retains from conventional (spatial-domain) TO the mindset based on geometrydriven conceptual design of a desired functionality.Instead of the geodesic path of light rays in the physical space, emphasis is now placed on the dispersion characteristics in the reciprocal phase space.In particular, starting from the standard plane-wave dispersion law in the auxiliary vacuum space r 0 (and associated reciprocal space k 0 ; cf.Fig. 1), denoting the corresponding speed of light in vacuum, we straightforwardly obtain [via Eq. ( 2)] the corresponding dispersion law in the transformed reciprocal space k, viz., For the simpler scenario involving frequency-independent wavenumber transformations in [14], we focused on the deformation of the equi-frequency contours (EFCs), relying on the wellknown relationships between their geometrical characteristics (e.g., presence and number of asymptotes, symmetries, inflection points, single/multi-valuedness) and the kinematical (wavevector and velocity) properties of the wave propagation and reflection/ refraction [33].In the present study, the explicit frequency dependence and possible complex-valued character introduced in the mapping function enable a finer control of spatiotemporal dispersion effects and density of states.
The interpretations and implications are schematically illustrated in Fig. 2, with reference to the 2D scenario of interest.In this case, the auxiliary vacuum space is characterized by a conical dispersion surface [Fig.2(a)], and two relevant cuts: a circular EFC [Fig.2(b)], and (assuming propagation along the z 0 direction, i.e., k 0 x 0) a linear, bi-directional dispersion diagram [cf.Fig. 2(c)].As is well known, this dispersion law dictates that wavenumbers jk 0 x j ≤ ω∕c yield propagating waves (i.e., real-valued k 0 z ), whereas jk 0 x j > ω∕c implies evanescence (i.e., purely imaginary k 0 z ).For a simple illustration, we first consider a linear, frequencyindependent transformation e F ≡ e F x ; e F z with a purely imaginary component, corresponding to e Pk x −k 2 x , e Qk z k 2 z , W ω ω in Eq. ( 7).Figures 2(d)-2(f ) show the transformed dispersion surface (now of hyperbolic shape) and related cuts.As an effect of the transformation, the circular EFCs [Fig.2(b)] are mapped onto imaginary values, with the exception of the points k 0 x 0; k 0 z ω∕c that are imaged onto themselves.Conversely, the region corresponding to jk 0 z j > ω∕c and purely imaginary values of k 0 x (i.e., evanescent waves along x 0 ) are now mapped onto real values of k x .This corresponds to the well-known scenario of hyperbolic (or indefinite) media, which allow the propagation of modes with large wavevectors, resulting in a very high (theoretically diverging) density of states (see [34] for a recent review).Though seemingly trivial, the above example provides some very useful insight into the manipulation of the density of states via complex-valued transformations.Starting from this simple configuration, further levels of sophistication can be introduced in order to engineer more complex responses, with the degree and the positive/ negative value of the polynomials e P and e Q affecting the nonlocality and density of states, and with W controlling the explicit temporal dispersion.For instance, we showed in [14] that multivalued transformations are naturally related to the appearance of extra branches in the dispersion surface, corresponding to additional extraordinary waves.
Figures 2(g)-2(i) qualitatively illustrate the mapping effects induced by a frequency-independent, noncentrosymmetric transformation with a purely imaginary component, chosen in such a way as to fold the resulting dispersion surface [Fig.2(g)] completely in the region k z ≥ k z0 > 0 of the transformed phase space.As is also evident from the corresponding mappings of the EFCs [Fig.2(h

APPLICATIONS TO DISPERSION ENGINEERING
To showcase the capabilities and potential of our proposed extensions, we now illustrate three end-to-end application examples, starting from the intuitive geometrical pictures in Fig. 2, and proceeding with the actual metamaterial synthesis.As previously mentioned, we focus on classes of multilayered metamaterials made of subwavelength, nonmagnetic material layers periodically stacked along the x direction.For these structures, in Supplement 1, we systematically derive (along the lines of [14,20]) nonlocal effective models that are structurally compatible with the NLTO blueprints in Eq. (8).
In what follows, for simplicity of illustration, we restrict ourselves to paraxial propagation along the z axis (i.e., k x ≈ 0), and therefore find it expedient to choose in all examples e Pk x p 2 k 2 x ; with p 2 denoting a constant real-valued parameter.This implies that the k x dependence in the NLTO blueprints ( 8) is lost, and we only need to deal with the k z and ω dependencies.
It is worth stressing that the emphasis at this stage is on the illustration of the methodology, rather than the specific applications.Although no particular optimization was carried out in this respect, the constitutive parameters of the material constituents are constrained within realistic bounds.However, for better illustration of the phenomena of interest, we neglect the effects of material losses and, whenever not instrumental, temporal dispersion.
Here, we apply our proposed NLTO framework to engineer one-way propagation effects, as previously discussed and schematically illustrated in Figs.2(g)-2(i).In its simplest terms, the problem can be mathematically formulated as finding a coordinate mapping so that, at a given radian frequency ω ω 0 and k x 0, the transformed dispersion relationship (10) allows propagation in only one direction along the z axis.Assuming the "forward" propagation as the allowed one, this implies that the propagation constant k z must be either real-positive or purely imaginary.As shown in Supplement 1, this can be attained, in the possibly simplest way, by choosing a frequency-independent, noncentrosymmetric transformation characterized by with k 0 ω 0 ∕c 2π∕λ 0 denoting the wavenumber in vacuum at the design frequency (and λ 0 the corresponding wavelength), k z0 > 0 denoting the desired propagation constant, and q 1 and q 2 being constant real-valued parameters subject to The mapping (13) identifies a broad class of transformation media that automatically satisfy the required one-way propagation condition.In principle, the degrees of freedom available could be exploited to enforce additional characteristics (e.g., group velocity).However, overconstrained NLTO blueprints may be very difficult (if not impossible) to synthesize over the broad spectral range (ideally, all real values of k z ) implied by the one-way propagation scenario of interest.In our example, in order to relax the synthesis procedure, we only prescribe the propagation constant k z0 4.5k 0 , thereby leaving two free parameters [including p 2 in Eq. ( 12)] and a partially constrained one [cf.Eq. ( 14)].
Figure 3 shows the results of our synthesis, whose details are given in Supplement 1.In order to approximate the associated NLTO blueprints (8), we consider a multilayered metamaterial based on a three-layer unit cell [shown as an inset in Fig. 3(b)] including a nonreciprocal gyrotropic material described by the relative-permittivity tensor The optimized material parameters and thicknesses (as well as transformation parameters) are given in the figure caption.Figures 3(a) and 3(b) compare the synthesized metamaterial (nonlocal) constitutive parameters with the prescribed NLTO blueprints.A generally good agreement is observed, especially in the vicinity of the nominal design point (k x 0, k z k z0 ). Figure 3(c) compares the NLTO-prescribed EFCs [cf.Eq. ( 10)] with those pertaining to the synthesized multilayered metamaterial (numerically computed via a rigorous transfer-matrixbased approach, as detailed in Supplement 1), from which the prescribed one-way character and the good agreement around the nominal design point are evident.Finally, Fig. 3(d) shows the finite-element-computed (see Supplement 1 for details) field distribution in the unit cell excited by a magnetic line-source, which provides a further confirmation and illustration of the one-way propagation phenomenon.

B. Stationary Points in the Dispersion Relationship
As a further application example, we consider the engineering of stationary points in the dispersion relation, which assume particular importance in slow-light scenarios (see, e.g., [35,36]).For propagation along the z direction (i.e., k x 0), at a given radian frequency ω ω 0 , we consider the following general conditions: which, depending on ν, are associated with different physical phenomena, including a regular band edge (ν 2), a "frozen" mode (ν 3), and a "degenerate" band edge (ν 4) [35,36].As shown in Supplement 1 [see also Figs. 2(j)-2(l) for a schematic illustration], these characteristics may be enforced via simple conditions on a frequency-independent coordinate mapping.For instance, we consider the possibly simplest form with q ν denoting a constant real-valued parameter.In essence, the condition in Eq. ( 17) enforces an νth-order zero in the dispersion relationship for k x 0. In what follows, we focus on the engineering of frozen-mode conditions (ν 3) with k z0 k 0 , and with the transformation parameters q 3 in Eq. ( 17) and p 2 in Eq. ( 12) left as degrees of freedom.In view of the inherent center-symmetry breaking (i.e., nonreciprocity), as for the previous example, our metamaterial synthesis relies on a three-layer unit cell featuring a gyrotropic material [cf.Eq. ( 15)], and enforces the parameter matching with the NLTO blueprints over a limited neighborhood (k x 0, 0.9k 0 < k z < 1.1k 0 ) of the nominal design point (see Supplement 1 for details).We remark that our chosen configuration is different from those in [35,36].In particular, our assumed direction of propagation along the z direction (i.e., parallel, rather than orthogonal, to the layer interfaces) allows us to maintain a 2D configuration and to avoid uniaxial material constituents with tilted optical axes.
The synthesis results are illustrated in Fig. 4.More specifically, Figs.4(a) and 4(b) show the good agreement between the NLTO blueprints and the synthesized nonlocal effective parameters over the relevant spectral ranges.The desired frozen-mode characteristics are confirmed by the presence of a cusp-like point in the EFCs [Fig.4(c)] and a stationary inflection point in the dispersion diagram [Fig.4(d)], in good agreement with the NLTO prescriptions.We note that, albeit nonreciprocal, the attained dispersion characteristics exhibit an additional backward-propagating (k z < 0) mode, which is not prescribed by the one-way NLTO template (17).This is not inconsistent, as the comparison is  8), ( 7), (12), and (17) with p 2 −20.349, ν 3, q 3 6.068k −1 0 , and k z0 k 0 .Blue-solid curves pertain to the nonlocal effective model of the synthesized multilayered metamaterial [three-layer unit cell shown in the inset of panel (b)] with ε 1 2.883, ̳ ε 2 given in Eq. (15) [with ε d −1.1 and ε g −0.03], ε 3 9.416, d 1 0.188λ 0 , d 2 0.049λ 0 , and d 3 0.193λ 0 .The parameter matching is enforced at k x 0 and within the (green-shaded) region 0.9k 0 < k z < 1.1k 0 .(c), (d) Prescribed NLTO EFCs and dispersion diagram (with a magnified view around the nominal design point) compared with actual synthesis results (black-solid).(e) Finite-element-computed magnetic-field (magnitude) map in false-color scale (normalized with respect to incident one) in the unit cell due to a plane-wave impinging from vacuum along the positive z direction.The blackdashed and thin gray lines indicate the vacuum-metamaterial and the layer interfaces, respectively.(f ) Corresponding magnetic density energy (normalized with respect to incident one, and averaged along the x direction) at z 5.42λ 0 as a function of the radian frequency (black-solid).Also shown (red-dashed), as a reference, is the fit in terms of the theoretically predicted behavior ∼jω − ω 0 j −2∕3 .meaningful only within the (green-shaded) localized phase-space region where the parameter matching is actually enforced.Figure 4(e) shows the finite-element-computed field map pertaining to a plane-wave excitation (at the design frequency, along the positive z direction) impinging from a vacuum half-space, from which the typical physical signatures of the frozen-mode regime can be observed [35,36].In particular, the wave penetrates the metamaterial half-space with small reflection, and it gets converted into an extended (frozen) mode whose amplitude gradually grows up to reaching a saturation value that is nearly two orders of magnitude larger than the incident one.As shown in Fig. 4(f ), the magnetic energy density W m (averaged across the unit cell) consistently peaks at the design frequency in agreement with the theoretical prediction ∼jω − ω 0 j −2∕3 [35,36].
Along the same lines, it is also possible to engineer other related effects, such as degenerate band edges [ν 4 in Eq. ( 16)].

C. Dirac-Point Conical Singularities
With the recent surge of graphene as a promising material for nanophotonics, there has been growing interest in reproducing some of its unique transport properties [37] in photonic structures.A particularly remarkable property of graphene is the presence of Dirac cones in its band structure, i.e., conical singularities at the six corners of the hexagonal Brillouin zone where the energy-momentum dispersion relationship is linear [37].Interestingly, EM analogues of such dispersion properties have been found, for instance, in 2D photonic crystals [45][46][47][48][49][50], in metamaterials featuring proper transitions and/or combinations of positive and negative refractive index [51,52], and in metallo-dielectric multilayers [53][54][55].In this latter case, of specific interest for our study, it is well known that Dirac points cannot exist at the center of the Brillouin zone (k 0).However, it is possible, in principle, to induce a conical singularity at a given radian frequency ω ω 0 and a nonzero wavevector (k x 0, k z k z0 ), viz., with γ x and γ z determining the slopes (and hence group velocities) of the linear dispersion in the principal planes.The dispersion relationship in Eq. ( 18) represents a very interesting and insightful testbed for our proposed NLTO approach as it involves a nontrivial interplay of spatial and (explicit) temporaldispersion effects.Indeed, it admits a very intuitive geometrical interpretation in terms of a simple phase-space shift (along k z and ω) of the conical dispersion surface pertaining to vacuum [cf.Eq. ( 9) and Fig. 2(a)].In principle, this can be achieved via the transformations (7), with Eq. ( 12) (assuming p 2 γ x ) and from which the instrumental role played by the explicit frequency dependence induced by the mapping function W is evident, at variance with the previous two examples.However, the transformation in Eq. ( 19) would yield NLTO blueprints [cf.Eq. ( 8)] characterized by the presence of odd powers of k z (i.e., nonreciprocity) and ω (which are physically feasible only in the presence of losses and/or gain).Nevertheless, we can easily derive a symmetric approximant (see Supplement 1 for details) The outcomes of our synthesis procedure, detailed in Supplement 1, are illustrated in Fig. 5.We consider a four-layer unit cell alternating dielectric and metallic layers.The latter are modeled via simple Drude dispersion laws with ε ∞1 ; ω p1 and ε ∞3 ; ω p3 included among the design parameters to be optimized.In order to relax the synthesis, we only prescribe γ x 0.25 in the k x plane, and leave as further optimization parameters the other parameter γ z as well as k z0 .Figures 5(a) and 5(b) compare the synthesized nonlocal constitutive parameters with the prescribed NLTO blueprints at the design radian frequency ω 0 , from which good agreement is observed around the nominal Dirac point (k x 0, k z k z0 0.778k 0 ).Qualitatively similar agreement is also observed at nearby frequencies.Figures 5(c) and 5(d) show the corresponding dispersion diagrams pertaining to the k z k z0 and k x 0 planes, respectively, from which the induced conical singularities can be observed, in good agreement with the NLTO prescriptions over the optimized frequency range.The slightly worse agreement in the k x 0 plane [Fig.5(d)] is attributable to a slightly worse approximation provided by the nonlocal effective model (see Supplement 1).As an independent check and validation, Fig. 5(e) shows the finite-element-computed response, at the design radian frequency ω 0 , due to a collimated Gaussian beam impinging along the positive z direction from a vacuum half-space.As can be observed, the beam transmitted in the metamaterial half-space undergoes a phenomenon that is typically associated with Dirac-point conical singularities [53][54][55].Such behavior is markedly different from what is observed at nearby frequencies, as exemplified by the response at ω 1.05ω 0 shown in Fig. 5(f).
The above example provides a clear illustration of the substantially enlarged capabilities granted by our proposed extension in comparison with the original approach in [14].The synthesis of such a complex spatiotemporal dispersion effect requires access to the full wavevector-frequency phase space, which is granted via the proposed extension, and it could not have been addressed by using frequency-independent transformations as in [14].

D. Some Remarks
It could be argued that a desired dispersion effect can be attained without necessarily resorting to a coordinate transformation from an auxiliary space, by directly tailoring the nonlocal constitutive tensors in the dispersion relationship that stems from solving the source-free, algebrized Maxwell's equations in the reciprocal space.The solution to this problem is not unique, in view of the fact that different materials can yield identical dispersion relationships.While not being the only possible approach, our transformation-based interpretation provides a systematic, constructive procedure to derive a possible realization of the required material parameters.For the simple 2D, coordinate-separable transformations and uniaxial, nonmagnetic materials considered in the examples, the relationship between the transformations and the material parameters turns out to be quite simple.Nevertheless, our approach can be applied to more general scenarios featuring more complicated (e.g., nonseparable) transformations and materials, where the relationship with the constitutive "blueprints" is much more complex.
Another concern could be the technological feasibility of the material constituents in the examples.As previously mentioned, although no specific emphasis was posed on this aspect, attention was paid so as to constrain the constitutive parameters within realistic bounds.Thus, the positive relative-permittivity values utilized in the various examples range from nearly 3 to about 9, whereas the models and parameters of the negative-permittivity and gyrotropic materials are consistent with certain configurations considered in the literature (see, e.g., [41,55]).In our proof-ofprinciple examples, we treated these parameters as continuous optimization variables, rather than specific discrete values in a given library.In more application-oriented implementations, one could instead consider a library featuring a discrete set of specific material-related, complex-valued, frequency-dependent parameters.This would essentially affect the optimization procedure, for which a genetic-algorithm-based approach [56] would probably be more apt.However, it would not substantially change the approach.

CONCLUSIONS AND PERSPECTIVES
In summary, the above examples of applications clearly illustrate the potential of our extended NLTO approach in dispersionengineering scenarios.In particular, retaining the attractive features of standard TO, our approach allows the conceptual design of the desired functionalities based on intuitive geometrical considerations in the reciprocal phase space, which gives natural access and control of important physical properties (e.g., nonreciprocity, density of states) as well as complex spatiotemporal dispersion effects.This is somehow separated from the subsequent metamaterial synthesis, which is posed as a parameter optimization of a physical structure described in terms of a nonlocal effective model so as to approximate the ideal blueprints over suitable regions of the frequency-wavenumber phase space.Within this framework, possible (geometrical) degrees of freedom in the conceptual design can be fruitfully exploited in order to relax the synthesis procedure.On the other hand, an educated inspection of the  8), ( 7), (12), and ( 20 conceptual blueprints can be very insightful in selecting among available "libraries" of nonlocal effective models (and associated physical structures) a problem-matched implementation platform.
We believe that, along with the increasing availability of nonlocal effective models, our NLTO approach may establish itself as a powerful and versatile framework for the engineering of complex dispersion effects.To this aim, we are currently working at further extensions and applications, including 3D scenarios, as well as more sophisticated (e.g., non-coordinate-separable) transformations and associated metamaterial classes.Another critical extension would be the possibility to relax the infinite-medium restriction, so as to apply the approach to limited regions.Within this framework, a simple spatial truncation of the infinitemedium synthesis would induce major distortion effects in the designed functionalities, in view of the unmodeled reflections and transmissions from the interface(s).We are currently working on a rigorous extension that should allow dealing with spatial discontinuities.The basic idea is to apply different phase-space transformations in different spatial regions and then derive the sought dispersion relationships by enforcing the suitable field-continuity conditions at the interface(s).This extension could enable the tailoring of surface-state dispersion and topological effects.
See Supplement 1 for supporting content.
)] and the dispersion diagram [Fig.2(i)], by comparison with the previous examples [Figs.2(a)-2(f )], waves propagating along the z axis are nonattenuated (real k z ) in the forward direction, and evanescent (purely imaginary k z ) in the backward direction.The response has thus acquired a markedly one-way (nonreciprocal) character.As a further example, Figs. 2(j)-2(l) illustrate the effect of a frequency-independent mapping exhibiting a stationary inflection point at k x 0, k z k z0 , and ω ω 0 .It can be observed that this translates into a cusp point in the relevant EFC [Fig.2(k)] and a stationary inflection point in the dispersion diagram [Fig.2(l)].From the physical viewpoint, this corresponds to a so-called "frozen" mode[35,36].Finally, Figs.2(m)-2(o) illustrate a very interesting example of a two-valued, frequency-dependent mapping that images a neighborhood (ω ≈ 0) of the original conical dispersion surface [Fig.2(a)]onto two cone-like branches at a nonzero radian frequency ω 0 and k x 0, k z k z0 [Fig.2(m)].As can be observed, at ω ω 0 , the EFCs degenerate in two points (k x 0, k z k z0 ) [Fig.2(n)], around which the dispersion law stays approximately linear [see, e.g., Fig.2(o)].Such Dirac-point conical singularities have recently elicited significant interest in view of their similarities with the band structure of graphene[37].

Fig. 2 .
Fig. 2. Geometrical interpretation of the NLTO approach.(a)-(c) Dispersion surface, EFC (at a given ω ω 0 ), and dispersion diagram (at k 0 x 0), respectively, pertaining to the auxiliary vacuum space [cf.Eq. (9)].(d)-(f) Corresponding maps in the transformed space for a frequency-independent transformation with a purely imaginary component [cf.Eq. (11)], yielding a hyperbolic dispersion law.(g)-(i) Same as above, but for a frequencyindependent, noncentrosymmetric transformation, yielding one-way propagation.(j)-(l) Same as above, but for a frequency-independent transformation yielding a cusp-like point in the EFC and an inflection point in the dispersion diagram (red dots), representative of a frozen mode.(m)-(o) Same as above, but for a two-valued, frequency-dependent transformation yielding Dirac-point conical singularities (red crosses).

Fig.
Fig. (a), (b) Constitutive parameters pertaining to a frozen-mode scenario at the design radian frequency ω ω 0 .Red-dashed curves represent the NLTO blueprints, obtained from Eqs. (8),(7),(12), and (17) with p 2 −20.349, ν 3, q 3 6.068k −1 0 , and k z0 k 0 .Blue-solid curves pertain to the nonlocal effective model of the synthesized multilayered metamaterial [three-layer unit cell shown in the inset of panel (b)] with ε 1 2.883, ̳ ε 2 given in Eq.(15) [with ε d −1.1 and ε g −0.03], ε 3 9.416, d 1 0.188λ 0 , d 2 0.049λ 0 , and d 3 0.193λ 0 .The parameter matching is enforced at k x 0 and within the (green-shaded) region 0.9k 0 < k z < 1.1k 0 .(c), (d) Prescribed NLTO EFCs and dispersion diagram (with a magnified view around the nominal design point) compared with actual synthesis results (black-solid).(e) Finite-element-computed magnetic-field (magnitude) map in false-color scale (normalized with respect to incident one) in the unit cell due to a plane-wave impinging from vacuum along the positive z direction.The blackdashed and thin gray lines indicate the vacuum-metamaterial and the layer interfaces, respectively.(f ) Corresponding magnetic density energy (normalized with respect to incident one, and averaged along the x direction) at z 5.42λ 0 as a function of the radian frequency (black-solid).Also shown (red-dashed), as a reference, is the fit in terms of the theoretically predicted behavior ∼jω − ω 0 j −2∕3 .

20 )
which contains only even powers of k z and ω and accurately reproduces the desired dispersion characteristics in the vicinity of the (symmetric) Dirac points (k x 0, k z k z0 ) at ω ω 0 [see Figs.2(m)-2(o) for a schematic illustration].This allows us to consider metallo-dielectric multilayers (as in [53-55]) as a possible implementation platform.