Molding light with metasurfaces: from far-field to near-field interactions

2.1 Accuracy of metasurface homogenization

As mentioned in Section 1, the analysis of near-field interactions between a source and a metasurface requires to pay particular attention to the validity of the metasurface homogenization. In this section, we would like to provide the reader with qualitative insight into the accuracy of this modeling approach.

The “characteristic lengths” involved in the problem under consideration are the spatial periods of the metasurface elements p_x and p_y (see Figure 2), the wavelength in the surrounding medium λ₁, the distance of the source from the metasurface z₀, and the distance of an observation point z. In particular, the latter two parameters play a crucial role in the validity of the homogenization approach for a source/observation in the near field of the metasurface, whereas only the metasurface period relative to the incident wavelength is important in the case of a conventional metasurface problem with far-field illumination and observation point. The “first homogenization condition” requires that the incident field radiated by the point source should vary little over one spatial period of the metasurface along the x- and y-directions [42]. The field radiated by an ideal point source can be represented as a continuous plane-wave spectrum containing (i) propagating waves with transverse (in-plane) wavenumber kT=kx2+ky2≤|k1| and (ii) evanescent waves with kT=kx2+ky2>|k1|. While all propagating waves satisfy the first homogenization condition for a metasurface with p≪λ₁ where p=Max(p_x, p_y), the situation is a bit trickier for the evanescent portion of the spectrum. First, if a surface mode is excited on the metasurface with wavenumber k_SW, it is necessary that the surface-mode wavelength is larger than the metasurface period: p≪λ_SW=2π/k_SW. This limit, combined with the fact that surface waves with larger in-plane wavenumber k_SW decay faster in the out-of-plane direction, allows setting an approximate lower bound on the distance of the source from the surface. In fact, for a sufficiently distant point source, the evanescent waves with transverse wavenumber large enough to “see” the periodicity of the metasurface may completely decay before even reaching the metasurface. Based on numerical experiments, it was found in Ref. [42] that these considerations lead to an approximate lower bound for the source distance above the surface that reads: z₀>p.

The “second homogenization condition” requires that the evanescent field scattered by the metasurface should be negligible at the observation point. Assuming that the first condition is satisfied, the field scattered by the metasurface admits a Bloch-Floquet representation in which all higher-order space harmonics are evanescent. In this case, the second condition is satisfied if the space harmonic with smallest attenuation constant in the orthogonal z-direction has negligible amplitude at the observation point. In Ref. [42], it was found that the attenuation of this space harmonic can be approximated as α_z≃2π/p. Based on additional considerations motivated by numerical experiments, an approximate lower bound for the observation point can then be set as z>p [42]. If both homogenization conditions are met, an analytical model based on effective conductivity tensors and homogenized Green’s function is accurate and can be used to predict the near-field interaction of a point source with an arbitrary metasurface.

2.2 Metasurfaces as surface-impedance sheets, and limits on wavefront shaping

To complete our theoretical discussion, in this section, we briefly present two intuitive models of a metasurface in terms of surface-impedance sheets (see Figure 3) and discuss their applicability and usefulness in different problems. To make our discussion simpler, we assume here that the considered metasurface is purely electric, with no magnetic conductivity tensor and no magnetoelectric coupling.

Figure 3:

Modeling of metasurfaces in terms of surface-impedance sheets.

Left: Original problem of a metasurface (planar array of meta-atoms) in a host structure (in this example, a grounded slab). Right: Two models in terms of surface-impedance sheets based on different forms of boundary conditions, as discussed in the text: (A) one-sided surface-impedance model and (B) two-sided surface-impedance model, with the corresponding transmission-line representation shown on the right (thick lines represent transmission-line segments and black boxes indicate lumped impedances). A generic point sources is depicted above the metasurface. Adapted with permission from Ref. [46].

In the literature, different methods have been used to model a metasurface as a surface impedance sheet depending on the specific problem at hand. These are generally based on either “one-sided” or “two-sided” surface-impedance conditions. In the one-sided surface-impedance modeling approach (Figure 3A), the metasurface, and the entire half-space behind it, which may include multiple layers, is modeled as a single impenetrable impedance sheet that acts as the load of a transmission line representing propagating/evanescent waves impinging on the metasurface (Figure 3A, right). Clearly, this model is suitable to study metasurfaces operating in reflection mode and, more generally, to problems that are concerned with the back-scattering properties of a multilayered system, as in the classical problem of calculating the modified emission rate (Purcell factor or PF) of a point dipole above a stratified medium [51]. In the one-sided surface-impedance approach, the boundary conditions in Eq. (1) can be reduced to z^×H1=σ_e⋅ET1, where a surface impedance tensor Z_s can be introduced as Z_s=σ_e−1, which is generally “spatially dispersive” (namely, it depends on the incident wavevector) [46]. As an example, for a metasurface composed of square periodic elements, the conductivity tensor is given by [54]

where K=1/(ZsTEZsTMkT2), and ZsTE and ZsTM indicate the parallel connection of the metasurface impedance Z_g and the input impedance of the entire region below the surface Z_d (namely, Z_s=Z_g||Z_d, where || indicates a parallel connection, as indicated in Figure 3A). In general, these impedances are different for TE and TM incidence. The expressions for the metasurface impedances, ZgTE and ZgTM, for typical metallization patterns can be found in Refs. [41], [42], and the transmission and reflection coefficients entering the Green’s function expansion in Eq. (3), as a function of surface impedance, are given in Ref. [46]. Moreover, as mentioned above, the one-sided surface-impedance tensor is generally spatially dispersive, which should be taken into account when calculating the power flow along the surface, as thoroughly discussed in Ref. [55] for a one-sided surface impedance model of a patterned metallic cladding over a grounded dielectric substrate.

The “two-sided surface-impedance method” instead simply models the metasurface as an impedance sheet Z_g, suspended in a host medium or host structure, and considers both the wave reflected and transmitted by the metasurface. An example is depicted in Figure 3B for the specific case of a metasurface on a grounded dielectric slab. This is clearly the method of choice for metasurfaces operating in transmission mode, for example, metasurfaces acting as flat lenses, designed to control the shape of the transmitted wavefront. For this modeling approach, the boundary conditions on the metasurface remain as in Eq. (1), with average fields on the right-hand side of the equations, Jes=σ_e⋅ETavg and Jms=σ_e⋅HTavg. The metasurface impedance tensor is again introduced as Z_s=σ_e−1. We also would like to note that, if, by applying this modeling approach, the obtained metasurface impedance inherently depends on the surrounding environment (e.g. the substrate), it is a clear sign that the employed metasurface homogenization process is inadequate (metamaterial or metasurface constitutive parameters are properly defined and meaningful only if they are an inherent property of the homogenized structure, independent of the surrounding environment [29], [56]).

In summary, when applied to the same geometry, as in Figure 3, either modeling approach, one-sided or two-sided, leads to identical results in terms of reflected fields. Therefore, to calculate the fields in the region that contains the source (backscattering or reflection problems), one can choose either the one-sided or the two-sided surface-impedance method. Instead, when calculating the transmitted fields through the metasurface, one needs to use the two-sided surface impedance approach and explicitly consider the presence of other layers behind the metasurface. In other words, the difference between these two modeling approaches originates from the two different forms of boundary conditions one may consider. In the case of one-sided surface impedance, there is no transmitted field, so the surface impedance we need to consider just connects the electric and magnetic fields above the impedance sheet, which therefore acts as the load of a transmission-line model of the structure under consideration (Figure 3A, right). In this case, the impedance of this load is composed of the actual metasurface impedance in parallel with the input impedance of the region below the metasurface. Conversely, in the two-sided surface-impedance method, the boundary conditions connect the electric and magnetic fields at the two sides of the metasurface through just the impedance of the metasurface itself, which therefore acts as a shunt impedance in a transmission-line model of the structure (Figure 3B, right). We refer the reader to Ref. [46] for additional details on these modeling approaches.

The idea of modeling a metasurface as a surface impedance suspended in free space has been used in Ref. [16] to derive “efficiency limits” on the ability of a generic passive metasurface to control the amplitude and phase of a propagating wave transmitted through it. In particular, for a single purely-electric metasurface, it was found that, (i) if the metasurface is isotropic, namely, with a scalar surface impedance, then the phase of the transmission coefficient t can be controlled only in a limited range [−90°, 90°] by varying the local properties of the metasurface, whereas, (ii) if the metasurface is anisotropic, with a tensorial impedance and polarization coupling, then the cross-polarized transmission phase ∠t_xy can be arbitrarily controlled, but its amplitude |t_xy| cannot exceed 0.25 for any transmission phase [16]. This fact sets a strict limit on the ability of inhomogeneous purely-electric metasurfaces to mold a propagating wavefront, namely, their ability to act as flat lenses. The efforts to go beyond these limits have led to more advanced metasurface designs operating in transmission mode: (i) multiple electric metasurface layers, inspired by radiofrequency transmit-arrays [16], [19], (ii) metasurfaces with magnetic, in addition to electric, response, realizing so-called “Huygens metasurfaces” that enable zero reflection (due to impedance matching to the surrounding medium) and ideal control of the transmission phase [18], [20], [21], [22], and (iii) metasurfaces with tailored bianisotropy and spatial dispersion for more general field transformations [24], [25], [26]. These ideas are currently at the basis of modern metasurfaces acting as flat optical/microwave components for far-field illumination.

3.1 Metasurfaces and ultrathin planar structures to control the emission rate of a localized emitter

One of the cornerstones of the field of nanophotonics is represented by the ability of controlling and tailoring the photonic emission of an excited atom/molecule or quantum dot (theoretically, a point source), by engineering the surrounding environment, a phenomenon known as the Purcell effect [51]. When placed in the vicinity of an optical/electromagnetic emitter, a metasurface, together with the surrounding medium, forms the environment that interacts with the source. Depending on the metasurface properties, the emission rate of the source can be modified owing to the inhibited or enhanced density of available states for photons, i.e. the photonic local density of states (LDOS), given by ρLDOS=(6/πω)Im[n^⋅G_(r, r, ω)⋅n^], where n^ represents the point-source polarization unit vector, and G(r, r, ω) is the total dyadic Green’s function at the source location [namely, the free-space Green’s function plus the scattered/reflected Green’s function given in Eq. (3)]. Based on the Green’s function analysis presented above, it is then relatively straightforward to investigate the effect of metasurfaces and any other stratified planar systems on the local density of photonic states, the spontaneous emission rate (i.e. the decay rate of a quantum emitter upon emission of photons due to interaction with the fluctuating electromagnetic vacuum field), and the so-called Purcell factor PF=(6π/k03)Im[n^⋅G_(r, r, ω)⋅n^] (i.e. the spontaneous emission rate with respect to the free-space scenario) [51].

Indeed, metasurfaces have recently attracted interest in the context of spontaneous emission control due to their guiding properties, optical resonances, and high field enhancements [57], [58], [59], [60]. Recent works have investigated the LDOS enhancement due to metasurfaces composed of lattices of discrete magnetoelectric scatterers [59] and the possibility of controlling the LDOS of a source on different graphene-based planar structures [60]. As an example, Figure 4A shows the tunability of the PF for a suspended graphene monolayer over a range of chemical potentials. It is clear that, in the low terahertz regime, the LDOS and PF can be tuned considerably by an external voltage bias. Figure 4B and C shows the PF for a graphene layer supported by a thin substrate and for a metasurface composed of two layers of graphene (a parallel plate-like waveguide). Finally, Figure 4D shows the PF as a function of source position and frequency.

Figure 4:

Emission-rate manipulation using graphene-based ultrathin planar structures.

(A) PF for a suspended graphene layer in vacuum with a point source at 10 nm above the surface. (B) PF for graphene on a substrate with thickness d_s=10 nm, ϵ=4 and (C) for two layers of graphene separated by a 10 nm, ϵ=4 spacer. (D) PF as a function of position (y₀ is the point-source position above the surface) and frequency for the case in (b). Adapted with permission from Ref. [60].

3.2 Hyperbolic metasurfaces: extreme anisotropy on a surface

Simple isotropic metasurfaces (i.e. metasurfaces with scalar, not tensorial, conductivity) can support surface modes that are launched isotropically on the surface by a point source in the near field, resulting in in-plane omnidirectional surface-wave propagation (or quasi-omnidirectional depending on the orientation of the point dipolar source). In general, isotropic or weakly anisotropic metasurfaces do not allow realizing the idea illustrated in Figure 1 of launching and controlling electromagnetic energy from the source location toward specific targets. In drastic contrast to the isotropic case, strongly anisotropic metasurfaces enable the ability of guiding the radiated energy along determined propagation channels on the metasurface [61], [62], [63], [64], [65], as we will see in the following. In this context, hyperbolic metamaterials and metasurfaces are uniaxial structures with extreme anisotropy characterized by an effective material tensor (metamaterial) or conductivity/impedance tensor (metasurface) with components having opposite sign for orthogonal electric field polarizations [63], [64] (the term “uniaxial” refers to anisotropic systems characterized by a diagonal constitutive tensor in a suitable coordinate system, with one diagonal element different from the others). As further discussed in the following, hyperbolic metamaterials/metasurfaces exhibit hyperbolic dispersion for wave propagation, as opposed to the usual elliptic dispersion, essentially combining the response of transparent dielectrics and reflective metals in different directions [66]. In the context of 3D metamaterials, these exotic properties have led to new physical effects and novel optical devices for a wide range of applications, such as subwavelength imaging, nanolithography, light-emission engineering [66], negative-index waveguides [67], subdiffraction photonic funnels [68], and novel nanoscale resonators [69]. While the complexity of fabrication of 3D metamaterials, and their sensitivity to losses, have hindered the impact of these concepts, especially in the terahertz and optical regime [63], [64], planar metasurfaces may overcome many of these challenges [9], [15], offering a practical route to realize extreme hyperbolic dispersion.

Figure 5A shows a schematic of the geometry under consideration, which contains a hyperbolic metasurface with the following conductivity tensor:

Figure 5:

Anisotropic and hyperbolic metasurfaces.

(A) Anisotropic surface with conductivity tensor σ at the interface of two isotropic materials. (B) Equifrequency contours, ω(k)=const, for hyperbolic metasurfaces having σ_xx=0.003+0.25i mS and σ_zz=0.03–0.76i mS (blue hyperbola) and σ_xx=1.3+16.9i mS and σ_zz=0.4–9.2i mS (green hyperbola). For comparison, the isotropic case for σ_xx=σ_zz=0.03–0.76i mS (black circle) is also shown. The red dashed line at 45° with respect to the x-axis is plotted for visual guidance. (C) Polar plot of the electric field amplitude |E_y| excited by a y-directed dipole current above a graphene strip array, assuming chemical potential μ_c=0.45 eV, width strip W=56.1 nm, periodicity L=62.4 nm, which yields σ_xx=0.003+0.25i mS and σ_zz=0.03–0.76i mS [blue hyperbola in (B)]. (D) Similar to (C), but for a graphene-strip array with μ_c=1 eV, W=196 nm and L=200 nm, and σ_xx=1.3+16.9i mS and σ_zz=0.4–9.2i mS [green hyperbola in (B)]. Plots in (C) and (D) are calculated at frequency f=10 THz, ρ=0.2λ (radial distance from source) and y=0.005λ (vertical distance from the surface). Solid blue line is for numerical calculation of the Green’s function, and dashed red line for residue-branch-cut evaluation. Adapted with permission from Ref. [64].

at the interface of two different media with electric permittivity and magnetic permeability ϵ₁, μ₁ and ϵ₂, μ₂. Although here the conductivity tensor is assumed to be diagonal, hence representing a uniaxial metasurface, in a general scenario, off-diagonal conductivity tensor elements may arise if the subwavelength meta-atoms are nonsymmetrical with respect to the coordinate system [21], [28], [65], or because of the presence of magneto-optical effects [51], [65], or also nonlocal effects due to finite Fermi velocity of the electrons in the metasurface elements [65], [70]. In the interest of clarity, here we only focus on diagonal conductivity tensors. For a point source radiating above the metasurface, the complete electromagnetic response could be computed using the general formalism provided in the previous section. Although the numerical evaluation of the Green’s function (3) can be challenging [46], [62], advanced techniques based on complex analysis to efficiently perform the integration have been reported in the literature (see, e.g. Ref. [64]). In particular, it was shown in Ref. [64] that the 2D integral in Eq. (3) can be reduced to (i) a pole residue associated with the surface-wave contribution (representing the surface mode propagating on the metasurface) and (ii) a branch-cut integral associated with the so-called continuous-spectrum contribution (representing the radiation into the surrounding environment) [71], [72]. These two contributions represent well the electromagnetic response of the system: part of the energy radiated by the point dipole is coupled into the guided surface modes of the metasurface, whereas the rest is directly radiated into the environment.

To better understand the behavior of surface modes supported by an anisotropic hyperbolic metasurface with conductivity (5), it is insightful to inspect the surface-wave dispersion relation [64]:

where ϵ₁=ϵ₂=ϵ, μ₁=μ₂=μ, k=ω(ϵμ)^1/2, η=(μ/ϵ)^1/2, and a=(kx2+kz2−k2)1/2 is the vertical attenuation of the surface wave. One of the main differences between surface modes in isotropic and anisotropic metasurfaces is that, in the isotropic case, surface modes could be found in the form of pure TE or pure TM modes, whereas the modes of anisotropic metasurfaces are always a mixture of TE and TM field contributions. In addition, in an anisotropic structure, the direction of energy flow, defined by the group velocity vector [73], is not necessarily parallel to the direction of the wavevector k. As usually done, the direction of energy flow (equal to the group velocity in low-loss regions) can be determined by examining the “equifrequency surface” (EFS), ω(k)=const, of the surface mode [64]. Because the group velocity is defined as the gradient of the dispersion function ω(k) with respect to wavevector (i.e. ∇_kω(k)), the direction of energy flow is necessarily orthogonal to the EFS contours. Assuming that the conductivity is purely imaginary and lossless, σii=jσ″ii,i=x, z, and that k_x, k_z≫k, the zeros of Eq. (6) can be approximated as the solution of [64]

At a given frequency ω, this equation may describe either an ellipse or a hyperbola in the k_x−k_y plane depending on the relative sign of the two conductivity components. For example, “elliptic EFS topology” arises when σ″xx⋅σ″zz>0 corresponding to the black curve in Figure 5B. The length of the ellipse’s principal axes along k_x and k_z is proportional to σ″zz and σ″xx, respectively. In this case, the modes might be associated with quasi-TM (inductive surface Im[σ_xx]>0, Im[σ_xx]>0) or quasi-TE (capacitive surface Im[σ_xx]<0, Im[σ_xx]<0) polarization [65] [as mentioned above, the character of the modes can be purely TE or TM only if the metasurface is isotropic (σ_xx=σ_zz)]. It was shown in Ref. [65] that the quasi-TE surface modes present a free-space like dispersion, leading to weak light-matter interaction, whereas the quasi-TM modes exhibit more confinement near the surface and, therefore, stronger interaction.

If one of the two conductivity tensor components changes sign, such that σ″xx⋅σ″zz<0, then the topology of the equifrequency contour drastically changes and becomes “hyperbolic”. In this regime, the anisotropic metasurface behaves like a dielectric in one direction and a metal in the orthogonal direction. Two examples of hyperbolic EFS are shown in Figure 5B for two values of surface conductivity (blue lines: σ_xx=0.003+0.25i mS and σ_zz=0.03–0.76i mS and green lines: σ_xx=1. 3–16.9i mS and σ_zz=0.4–9.2i mS). The hyperbola asymptotes are defined by kz=±kx|σ″xx/σ″zz|, and the direction of energy flow of the surface mode is given by the angle [63], [64]

which is defined with respect to the z-axis (note that this is valid only in the hyperbolic regime in the asymptotic region, namely, for large enough wavenumber). A simple interpretation of the above equation is that the energy tends to flow mostly toward the direction with lower imaginary conductivity component [63].

Considering that the field radiated by a point dipole has, in theory, an infinite angular spectrum, containing all possible values of transverse wavenumber, this near-field illumination can excite a continuous set of surface modes, with different wavenumber, according to the EQS in Figure 5B (we also would like to note that this hyperbolic EQS has no cutoff at high values of spatial frequency, due to the “open” shape of the hyperbola, in drastic contrast with conventional elliptic dispersion; in practice, however, any realistic hyperbolic metasurface exhibits a high spatial-frequency cutoff, where the hyperbolic curve “closes” due to the effect of nonlocality, as the high spatial frequencies start “seeing” the discrete nature of the metasurface [70]). Because the vectors normal to the hyperbola point in roughly the same direction for a given sign of k_x and k_z, a dipole source will excite a narrow surface-wave beam confined on the surface of a hyperbolic metasurface [64]. For example, the asymptotes of the blue hyperbola in Figure 5B have an angle of 30° with respect to the z-axis, and thus the normal to the hyperbola (i.e. the group velocity vector) is pointing at 60° with respect to the z-axis.

As an example, an hyperbolic metasurface with the conductivity tensor corresponding to the blue hyperbola in Figure 5B can be implemented with an array of graphene strips, as shown in Figure 6A, with chemical potential μ_c=0.45 eV, strip width W=56.1 nm, and periodicity L=62.4 nm [64]. For this physical implementation, the electric field profile of the surface modes has been obtained in Ref. [64] using two methods: (i) numerically solving the Green’s function in Eq. (3), which is a computationally demanding operation, and (ii) using a faster residue-branch-cut evaluation, as discussed at the beginning of this section. The results are reported in Figure 5C, showing that the two methods are in good agreement, hence justifying the faster evaluation in terms of residue and branch-cut contributions. As can be seen, most of the surface-wave energy is flowing toward ϕ=60° from the z-axis (and mirror-symmetric directions with respect to the two axes), confirming the response predicted by the hyperbolic equifrequency contours. Figure 5D also shows a different case with energy flow toward an angle ϕ=36°, corresponding to the steeper green hyperbola in Figure 5B. This response can be realized with an array of graphene strips having μ_c=1 eV, W=196 nm, and L=200 nm [64].

Figure 6:

Graphene-based anisotropic and hyperbolic metasurfaces.

(A) Array of graphene strips of width W and period L. (B) Imaginary parts of σ_xx and σ_zz and (C) real parts of σ_xx and σ_zz (normalized by σ₀=e²/4ħ) for a graphene-strip array with τ=0.35 ps, μ_c=0.33 eV, W=59 nm and L=64 nm. Region 1 (light red) is hyperbolic and region 2 (light blue) is anisotropic with elliptical dispersion. The yellow/green transition region is characterized by high losses in the x-direction. Adapted with permission from Ref. [64].

Another case of particular relevance is when the hyperbolic metasurface exhibits an extreme anisotropy in which the imaginary conductivity along a specific direction becomes very small compared to the one in the orthogonal direction, |σ″zz|≪|σ″xx|, which is sometimes called the “canalization regime”. In this extremely anisotropic hyperbolic metasurface, the hyperbolic equifrequency contour tends to become flat, which implies that any surface wave excited by a near-field source, with any transverse wavenumber k_x, would mostly propagate along the z-axis, yielding a narrow “diffraction-less” surface-wave beam that almost does not spread out as it propagates along the surface [74]. In practice, this scenario may again be implemented with a suitably engineered array of parallel graphene nanoribbons having complex conductivities. An analogous implementation has been studied in Ref. [74]: as shown in Figure 7A, the negative and positive conductivities are created electronically by gating a graphene layer over a triangularly corrugated plane. Figure 7C and D shows the electric field amplitude distribution for two point sources, separated by a subwavelength distance, on the plane z=0 over the modulated-graphene metasurface and over a homogenous graphene layer, respectively, clearly showing the surface-wave canalization in the modulated case. This behavior is of clear interest for imaging applications, as the surface waves avoid the usual diffraction expected on a homogeneous layer, allowing to distinguish the two point sources even at a large distance from the source plane.

Figure 7:

Surface-wave canalization on a modulated-graphene metasurface.

(A) A graphene monolayer with fixed gating voltage over a triangularly perturbed ground plane that realizes a tailored conductivity modulation. (B and C) Amplitude of the vertical component of the electric field for the surface plasmon polaritons launched by two closely located point sources over (B) modulated graphene and (C) homogeneous graphene. Adapted with permission from Ref. [74].

Further details on the use of graphene, and other 2D materials, to implement metasurfaces with extreme anisotropy are discussed in the next section.

3.3 2D Materials as anisotropic and hyperbolic metasurfaces

The most straightforward approach to realize uniaxial and hyperbolic metasurfaces is to reduce the profile of well-known bulk uniaxial materials resulting in thin layers that essentially act like 2D metasurfaces [65]. Even more interesting, the emerging class of 2D atomic crystals [75] represents the ultimate embodiment of metasurfaces in terms of thinness, combined with exciting performance in terms of tunability, flexibility, quality factor, etc. [64]. Some notable examples of 2D layered crystals include graphene and transition metal dichalcogenides and trichalcogenides [64]. Another crystal with similar properties is magnesium diboride, in which graphene-like layers of boron are alternated by densely packed layers of magnesium [65], [76]. Hexagonal boron nitride (hBN) should also be mentioned as another promising material in this category [65], especially due to its excellent compatibility with graphene optoelectronics [77]. Notably, hBN has allowed the experimental demonstration of low-loss hyperbolic phonon polaritons in the infrared [78], [79].

One of the most popular approaches to realize uniaxial and hyperbolic metasurfaces at terahertz and infrared frequencies is the use of dense arrays of graphene strips [63], [80], as shown in Figure 6A. The dispersion topology of this structure may range from elliptical to hyperbolic as a function of its geometrical and electrical parameters. The in-plane effective conductivity tensor of this graphene-based metasurface can be analytically calculated using an effective medium theory, obtaining [63]

where L and W are the period and width of the strips, respectively, G=L−W is the separation distance between two consecutive strips, σ is the conductivity of homogeneous graphene, and σc=jωε0Lπln(cscπG2L) is an equivalent conductivity associated with the near-field coupling between adjacent strips, obtained using an electrostatic approach [81]. These effective parameters are valid only when the homogenization condition L≪λ_SW is satisfied, where λ_SW is the wavelength of the guided surface wave in the in-plane direction perpendicular to the strips (x in this case), thus leading to a homogeneous 2D metasurface [64]. Figure 6B and C shows the real and imaginary parts of σ_xx and σ_zz in a wide range of frequencies in the terahertz range, for a graphene-strip array with scattering time τ=0.35 ps, chemical potential μ_c=0.33 eV, and geometrical parameters W=59 nm and L=64 nm. As can be seen in Figure 6, the effective conductivities are mostly imaginary, except in the resonant lossy window near 15 THz. Most importantly, at frequencies below the resonance, the imaginary conductivities have opposite signs in the two directions, whereas they are both negative above the resonance. As a result, this structure can exhibit both a hyperbolic response at low frequencies and a nonhyperbolic (but still anisotropic) response at higher frequencies, separated by a lossy region in which surface waves are strongly attenuated.

Another promising material platform for metasurfaces is represented by BP, which is also a layered material, similar to graphite, with each layer forming a “wrinkled” surface [64], as depicted in Figure 8A. BP is a thermodynamically stable phase of phosphorus at ambient temperature and pressure [82], and it has recently been exfoliated into its multilayers [83], [84]. BP possesses exciting optical and electrical properties, such as an intrinsic direct band gap that may range from ~2 eV in monolayers (phosphorene) to 0.3 eV in its bulk configuration; tunable electric response versus thickness, or via externally applied electric/magnetic fields and mechanical strain; and the possibility to support confined surface plasmons in the mid-infrared to near-infrared range [64], [65]. Similar to the case of graphene, the ultrathin nature of BP allows a simple electromagnetic characterization in terms of optical conductivity, which may be accurately derived applying the well-known Kubo formalism [85], [86]. Figure 8B and C shows the real and imaginary parts of the BP conductivity components.

Figure 8:

Naturally anisotropic 2D materials: the case of black phosphorus (BP).

(A) Lattice structure of monolayer BP. Different colors are used for visual clarity. (B) Imaginary part of BP conductivity components versus frequency for several values of chemical potential. (C) Real part of BP conductivity components versus frequency for several values of chemical potential. Solid, dashed, and dotted lines correspond to chemical potentials of 0.005, 0.05, and 0.1 eV, respectively. The considered BP thickness is 10 nm, its direct band gap is 0.485 eV, damping 5 meV, at a temperature of 300 K. Adapted with permission from Ref. [65].

Remarkably, a thin film of BP is an example of a natural hyperbolic metasurface. In fact, Figure 8B and C shows that, in the low-loss region where the conductivity is almost purely imaginary, BP can behave as either a hyperbolic or an elliptical anisotropic surface depending on the frequency and chemical potential. However, although BP exhibits a hyperbolic regime, the resulting values of conductivity are rather small; therefore, although a hyperbolic surface plasmon can be excited, this is generally not the dominant response, as discussed in Ref. [64]. Larger values of conductivities can be obtained in the nonhyperbolic (Drude) regime. For example, a 10-nm-thick BP film with electron doping level 10×10¹³/cm² has conductivity tensor components σ_xx=0.0008–0.2923i and σ_zz=0.0002–0.0658i mS at f=92.6 THz. As shown in Figure 9, this highly anisotropic, but nonhyperbolic, surface can support a directional surface wave mainly directed along one of the coordinate axes and tightly confined on the surface. In particular, Figure 9A shows the electric field profile of the surface plasmon polariton excited by a point source above the surface, calculated by numerically solving Maxwell’s equations using a finite-difference time-domain method (FDTD) provided by a commercial software (Lumerical [87]). Figure 9B shows the vertical variation of the field distribution calculated by Lumerical, showing strong light confinement near the surface.

Figure 9:

Directional surface plasmon polaritons on anisotropic black phosphorus (BP).

Field intensity distribution excited by a y-directed dipole current source on a 10-nm-thick BP layer with doping level 10×10¹³/cm² at f=92.6 THz. (A) In-plane field distribution and (B) vertical out-of-plane field distribution, dominated by the response of a surface plasmon polariton mainly directed along the x-axis and tightly confined on the surface. Both figures are calculated with a commercial FDTD solver. The black dot/arrow indicates a y-polarized dipole. Adapted with permission from Ref. [64].