Surface waves with negative phase velocity supported by temperature-dependent hyperbolic materials

TL;DR

This study explores how temperature-dependent hyperbolic materials can support surface waves, including those with negative phase velocity, with propagation characteristics changing with temperature and material composition.

Contribution

It demonstrates the existence of negative phase velocity surface waves in temperature-sensitive hyperbolic materials and analyzes their propagation behavior.

Findings

Negative phase velocity surface waves can occur at high temperatures.

The angular range for NPV propagation varies with temperature and material composition.

Multiple surface wave modes can propagate at high temperatures.

Abstract

A numerical investigation was undertaken to elucidate the propagation of electromagnetic surface waves guided by the planar interface of two temperature-sensitive materials. One partnering material was chosen to be isotropic and the other to be anisotropic. Both partnering materials were engineered composite materials, based on the temperature-sensitive semiconductor InSb. At low temperatures the anisotropic partnering material is a non-hyperbolic uniaxial material; as the temperature is raised this material becomes a hyperbolic uniaxial material. At low temperatures, a solitary Dyakonov wave propagates along any specific direction in a range of directions parallel to the planar interface. At high temperatures, up to three different surface waves can propagate in certain directions parallel to the planar interface; one of these surface waves propagates with negative phase velocity (NPV). At a fixed temperature, the range of directions for NPV propagation decreases uniformly in extent as the volume fraction of InSb in the isotropic partnering material decreases. At a fixed volume fraction of InSb in the isotropic partnering material, the angular range for NPV propagation varies substantially as the temperature varies.

Equations

$${\underline{\underline{\epsilon }}}_{\mathcal{A}}={\epsilon }_{\mathcal{A}}^{\mathrm{s}}\underline{\underline{I}}+\left({\epsilon }_{\mathcal{A}}^{\mathrm{t}}-{\epsilon }_{\mathcal{A}}^{\mathrm{s}}\right)\hat{\underline{u}}\hat{\underline{u}},$$

$$\hat{\underline{u}}=\cos \psi {\hat{\underline{u}}}_x+\sin \psi {\hat{\underline{u}}}_y,$$

$${\underline{E}}_{\mathcal{A}}(\underline{r})={\underline{\mathcal{E}}}_{\mathcal{A}1}\exp \left(i{\underline{k}}_{\mathcal{A}1}\cdot \underline{r}\right)+{\underline{\mathcal{E}}}_{\mathcal{A}2}\exp \left(i{\underline{k}}_{\mathcal{A}2}\cdot \underline{r}\right),$$

$${\underline{k}}_{\mathcal{A}\ell }=k_0\left(q{\hat{\underline{u}}}_x+i{\alpha }_{\mathcal{A}\ell }{\hat{\underline{u}}}_z\right),\ell \in \left\{1,2\right\},$$

$$\left.\begin{array}[]{l}\alpha_{{\cal A}1}=\sqrt{q^{2}-\varepsilon_{\cal A}^{s}}\vspace{8pt}\\ \alpha_{{\cal A}2}=\displaystyle{\sqrt{\varepsilon^{t}_{\cal A}\left[q^{2}\left(\frac{\cos^{2}\psi}{\varepsilon^{s}_{\cal A}}+\frac{\sin^{2}\psi}{\varepsilon^{t}_{\cal A}}\right)-1\right]}}\end{array}\right\}\,.$$

$${\underline{E}}_{\mathcal{B}}(\underline{r})={\underline{\mathcal{E}}}_{\mathcal{B}}\exp \left(i{\underline{k}}_{\mathcal{B}}\cdot \underline{r}\right),$$

$${\underline{k}}_{\mathcal{B}}=k_0\left(q{\hat{\underline{u}}}_x-i{\alpha }_{\mathcal{B}}{\hat{\underline{u}}}_z\right),$$

$${\alpha }_{\mathcal{B}}=\sqrt{q^2-{\epsilon }_{\mathcal{B}}}.$$

$${\epsilon }_{\mathcal{A}}^s\left({\epsilon }_{\mathcal{B}}\sqrt{q^2-{\epsilon }_{\mathcal{A}}^s}+{\epsilon }_{\mathcal{A}}^s\sqrt{q^2-{\epsilon }_{\mathcal{B}}}\right)\left[\sqrt{q^2-{\epsilon }_{\mathcal{B}}}+\sqrt{q^2\left(\frac{{\epsilon }_{\mathcal{A}}^t{\cos }^2\psi }{{\epsilon }_{\mathcal{A}}^s}+{\sin }^2\psi \right)-{\epsilon }_{\mathcal{A}}^t}\right]{\tan }^2\psi$$

$${\epsilon }_{\text{InSb}}={\epsilon }_{\infty }-\frac{{\omega }_p^2}{{\omega }^2+i\gamma \omega }.$$

$$N=5.76\times 10^{20}{T}^{3/2}\exp \left(-\frac{{\mathsf{E}}_{\mathsf{g}}}{2k_{B}T}\right),$$

$$\left.\begin{array}[]{l}\varepsilon_{\cal A}^{s}=f^{\cal A}_{\text{InSb}}\varepsilon_{\text{InSb}}+\left(1-f^{\cal A}_{\text{InSb}}\right)\vspace{6pt}\\ \varepsilon_{\cal A}^{t}=\left[\frac{f^{\cal A}_{\text{InSb}}}{\varepsilon_{\text{InSb}}}+\left(1-f^{\cal A}_{\text{InSb}}\right)\right]^{-1}\end{array}\right\},$$

$${\epsilon }_{\mathcal{B}}=f_{\text{InSb}}^{\mathcal{B}}{\epsilon }_{\text{InSb}}+\left(1-f_{\text{InSb}}^{\mathcal{B}}\right){\epsilon }_{\text{HDPE}},$$

$${\underline{P}}_{\ell }(\underline{r})=\frac{1}{2}\text{mbox}Re\left[{\underline{E}}_{\ell }(\underline{r})\times {\underline{H}}_{\ell }^{\ast }(\underline{r})\right],\ell \in \left\{\mathcal{A},\mathcal{B}\right\},$$

$$\left.\begin{array}[]{l}\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{s}\right\}\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{t}\right\}<0\vspace{6pt}\\ \mbox{Re}\left\{\varepsilon_{\mathcal{B}}\right\}<0\end{array}\right\}$$

Full text

**Surface waves with negative phase velocity supported by temperature-dependent hyperbolic materials **

Tom G. Mackay***E–mail: T.Mackay@ed.ac.uk.

*School of Mathematics and Maxwell Institute for Mathematical Sciences

University of Edinburgh, Edinburgh EH9 3FD, UK

*and

*NanoMM — Nanoengineered Metamaterials Group

Department of Engineering Science and Mechanics

Pennsylvania State University, University Park, PA 16802–6812, USA

Akhlesh Lakhtakia†††E–mail: akhlesh@psu.edu

*NanoMM — Nanoengineered Metamaterials Group

Department of Engineering Science and Mechanics

Pennsylvania State University, University Park, PA 16802–6812, USA*

Abstract

A numerical investigation was undertaken to elucidate the propagation of electromagnetic surface waves guided by the planar interface of two temperature-sensitive materials. One partnering material was chosen to be isotropic and the other to be anisotropic. Both partnering materials were engineered composite materials, based on the temperature-sensitive semiconductor InSb. At low temperatures the anisotropic partnering material is a non-hyperbolic uniaxial material; as the temperature is raised this material becomes a hyperbolic uniaxial material. At low temperatures, a solitary Dyakonov wave propagates along any specific direction in a range of directions parallel to the planar interface. At high temperatures, up to three different surface waves can propagate in certain directions parallel to the planar interface; one of these surface waves propagates with negative phase velocity (NPV). At a fixed temperature, the range of directions for NPV propagation decreases uniformly in extent as the volume fraction of InSb in the isotropic partnering material decreases. At a fixed volume fraction of InSb in the isotropic partnering material, the angular range for NPV propagation varies substantially as the temperature varies.

Keywords: hyperbolic materials, Dyakonov waves, surface–plasmon–polariton waves, negative phase velocity, temperature control, dissipative materials

1 Introduction

The planar interface of two dissimilar materials can support the propagation of a variety of types of surface wave, even when both partnering materials (on either side of the interface) are homogeneous, non-magnetic and non-magnetoelectric [1]. For examples, (i) the planar interface of a plasmonic material and a dielectric material can guide the propagation of surface-plasmon-polariton (SPP) waves [2]; and (ii) the planar interface of an isotropic dielectric material and an anisotropic dielectric material can guide the propagation of Dyakonov waves [3, 4, 5, 6]. Generally, SPP waves propagate for wide ranges of directions parallel to the interface plane, for wide ranges of constitutive-parameter values of the partnering materials. In contrast, in the case of nondissipative partnering materials, Dyakonov waves are restricted to much smaller ranges of propagation directions parallel to the interface plane, and only certain restrictive ranges of constitutive-parameter values of the partnering materials are necessary for their propagation [7, 6, 8]. In cases where the partnering materials are dissipative, even to a small extent, Dyakonov-wave propagation is possible for much wider ranges of propagation directions, for much wider ranges of constitutive-parameter values, and in such cases more than one Dyakonov wave can propagate in a given direction [9].

If the constitutive parameters of one (or both) of the partnering materials are sensitive to temperature, then the interface may guide surface waves of different types at different temperatures. In a recent numerical study [9], the transition from Dyakonov waves to SPP waves was mediated by varying the temperature for the interface of a temperature-sensitive isotropic partnering material and a temperature-insensitive anisotropic partnering material. The temperature-sensitive material chosen was the semiconductor InSb. Indeed, InSb may also be used to thermally tune metamaterials and metasurfaces operating in the terahertz spectral regime [10, 11]. This ability to control the transition from Dyakonov waves to SPP waves could potentially be harnessed for temperature sensing and thermal imaging applications [12, 13].

Hyperbolic materials [14], as exemplified by an anisotropic dielectric material whose permittivity dyadic has a real part with eigenvalues of opposite signs, are associated with exotic phenomenons such as negative refraction [15, 16, 17] and the closely related phenomenon of negative phase velocity (NPV) [18]. Such materials may be exploited in subwavelength imaging [19, 20, 21], for radiative thermal energy transfer [22, 23], as analogs of curved spacetime [24, 25], and for diffraction gratings capable of directing light into a large number of refraction channels [26], for examples. Also, surface-wave propagation supported by hyperbolic materials has attracted the recent attention of both theorists and experimentalists [27, 28, 36, 29, 30, 33, 34, 35, 37, 31, 32].

Parenthetically, the taxonomy of surface waves supported by hyperbolic materials is problematic: on the one hand such surface waves could be regarded as Dyakonov waves because of the anisotropy of the hyperbolic material, but on the other hand such surface waves could also be regarded as SPP waves because of the negative-valued eigenvalue(s) of the real part of the hyperbolic material’s permittivity dyadic [1].

In the following sections, we numerically investigate the propagation of surface waves guided by the interface of two temperature-sensitive partnering materials. Both partnering materials are non-magnetic, non-magnetoelectric, and engineered materials, one being isotropic and the other anisotropic. The anisotropic partnering material is not of the hyperbolic kind at low temperatures, but it becomes hyperbolic as the temperature rises. Our main result in this paper is that the high-temperature regime can support surface-wave propagation with NPV, whereas the low-temperature regime cannot.

In the notation adopted, vectors are underlined once and dyadics twice. An $\exp(-i\omega{t})$ dependence on time $t$ is implicit, with $i=\sqrt{-1}$ and angular frequency $\omega$. The triad $\left\{\hat{\underline{u}}_{x},\hat{\underline{u}}_{y},\hat{\underline{u}}_{z}\right\}$ contains the Cartesian unit vectors, while the position vector is denoted by $\underline{r}=x\hat{\underline{u}}_{x}+y\hat{\underline{u}}_{y}+z\hat{\underline{u}}_{z}$. The free-space wavenumber is $k_{\scriptscriptstyle 0}=\omega\sqrt{\varepsilon_{\scriptscriptstyle 0}\mu_{\scriptscriptstyle 0}}$, with $\varepsilon_{\scriptscriptstyle 0}$ and $\mu_{\scriptscriptstyle 0}$ being the permittivity and permeability of free space, respectively; and $\lambda_{\scriptscriptstyle 0}=2\pi/k_{\scriptscriptstyle 0}$ is the free-space wavelength.

2 Canonical boundary-value problem

The setting for our study is the canonical boundary-value problem for surface-wave propagation guided by the planar interface of a uniaxial material and an isotropic material [1], with both materials being non-magnetic and non-magnetoelectric [38, 39]. As the theory underpinning this setting is amply described in earlier works [40, 9], only the bare essentials need be presented here.

The half-space $z>0$ is occupied by a uniaxial material labeled $\mathcal{A}$. Its relative permittivity dyadic is given by

[TABLE]

wherein $\underline{\underline{I}}=\hat{\underline{u}}_{x}\hat{\underline{u}}_{x}+\hat{\underline{u}}_{y}\hat{\underline{u}}_{y}+\hat{\underline{u}}_{z}\hat{\underline{u}}_{z}$ is the identity dyadic [41]. The optic axis of material ${\cal A}$, whose direction is given by the unit vector

[TABLE]

lies in the $xy$ plane, at an angle $\psi$ relative to the $x$ axis. The half-space $z<0$ is occupied by material $\mathcal{B}$ whose relative permittivity dyadic is $\varepsilon_{\mathcal{B}}\underline{\underline{I}}$. A schematic representation of the canonical boundary-value problem is provided in Fig. 1.

With no loss of generality, surface-wave propagation parallel to $\hat{\underline{u}}_{x}$ is considered. The electric field phasor in the half-space $z>0$ is given in terms of amplitude vectors $\underline{\mathcal{E}}_{\,{\cal A}1}$ and $\underline{\mathcal{E}}_{\,{\cal A}2}$ as

[TABLE]

with the wavevectors

[TABLE]

and the decay constants

[TABLE]

The electric field phasor in the half-space $z<0$ is given in terms of an amplitude vector $\underline{\mathcal{E}}_{\,{\cal B}}$ as

[TABLE]

with the wavevector

[TABLE]

and decay constant

[TABLE]

Crucially, the inequalities ${\rm Re}\left\{\alpha_{{\cal A}\ell}\right\}>0$, $\ell\in\left\{1,2\right\}$, and $\mbox{Re}\left\{\alpha_{\cal B}\right\}>0$ must be satisfied for surface-wave propagation.

The normalized propagation constant $q$ is delivered by solving the corresponding dispersion relation [9]

[TABLE]

Numerical methods, such as the Newton–Raphson method [42], are generally needed for this task. The following symmetries exist: if $q=q^{\star}$ satisfies Eq. (2) at the angle $\psi=\psi^{\star}$ then $q=q^{\star}$ also satisfies Eq. (2) for $\psi=-\psi^{\star}$ and $\psi=\pi\pm\psi^{\star}$. Accordingly, in the proceeding presentation of numerical solutions to Eq. (2) (in Sec. 4), only the range $0\leq\psi\leq\pi/2$ need be considered.

3 Temperature-controlled partnering materials

In order to control the permittivity parameters of the materials occupying $z>0$ and $z<0$ with temperature, both half-spaces are taken to be filled with homogenized composite materials containing the temperature-sensitive semiconductor InSb. The relative permittivity of InSb in the terahertz regime is provided by the Drude model [43, 11] as

[TABLE]

Herein the plasma frequency $\omega_{p}=\sqrt{Nq_{e}^{2}/0.015\,\varepsilon_{\scriptscriptstyle 0}\,m_{e}}$ is determined by the electronic charge $q_{e}=-1.60\times 10^{-19}$ C and mass $m_{e}=9.11\times 10^{-31}$ kg, while the high-frequency relative permittivity $\varepsilon_{\infty}=15.68$ and the damping constant $\gamma=\pi\times 10^{11}$ rad s*-1*. The temperature $T$ (in K) dependence of $\varepsilon_{\text{InSb}}$ arises via the intrinsic carrier density (in m*-3*) [44, 45, 46]

[TABLE]

which depends upon the bandgap ${\sf E_{g}}=0.26$ eV and the Boltzmann constant $k_{B}=8.62\times 10^{-5}$ eV $\mbox{K}^{-1}$.

The half-space $z>0$ is taken to be structured as a periodic multilayer, comprising alternating electrically thin sheets of InSb and air, with each sheet being parallel to the interface $z=0$. In the long-wavelength regime, the multilayer functions like a uniaxial continuum whose relative permittivity dyadic has the form given in Eq. (1). The relative permittivity parameters of this homogenized composite material occupying the half-space $z>0$ are estimated using the periodic-multilayer approximation [47, 48] as

[TABLE]

where $f^{\cal A}_{\text{InSb}}\in\left(0,1\right)$ is the volume fraction of the InSb layers occupying $z>0$. Parenthetically, a multilayer structure with alternate air layers can be fabricated by etching [49, 50], for example. Also, the engineered uniaxial material characterized by the relative permittivity dyadic in Eq. (1) could arise from the homogenization of identically oriented, electrically small, spheroidal particles dispersed randomly [18, 51]. Fabrication using standard technologies has also been reported [34, 35].

The half-space $z<0$ is taken to be composed of a random distribution of electrically small spheres made of InSb and a temperature-insensitive polymer, namely high-density polyethylene (HDPE), of relative permittivity $\varepsilon_{\text{HDPE}}$. In the long-wavelength regime, the material filling $z<0$ may be regarded as an isotropic continuum whose relative permittivity is estimated as

[TABLE]

where $f^{\cal B}_{\text{InSb}}\in\left(0,1\right)$ is the volume fraction of InSb in the half-space $z<0$.

For the purposes of our calculations, the volume fraction $f^{\cal A}_{\text{InSb}}=0.5$ was selected. The frequency was fixed at $f=2.0$ THz (i.e., $\lambda_{\scriptscriptstyle 0}=0.15$ mm). At this frequency, the relative permittivity of HDPE is approximately constant, i.e., $\varepsilon_{\text{HDPE}}\approx 2.387+0.006i$, over the temperature range $T\in\left[260,290\right]$ K [52]. Plots of the real and imaginary parts of the relative permittivity parameters of the partnering materials ${\cal A}$ and ${\cal B}$ versus temperature are provided in Fig. 2 for the range $T\in\left[260,290\right]$ K with $f^{\cal B}_{\text{InSb}}\in\left\{0.7,0.85,1\right\}$. The real part of $\varepsilon_{\cal B}$ decreases uniformly as $T$ increases, for the range of values of $f^{\cal B}_{\text{InSb}}$ considered. Specifically, for $f^{\cal B}_{\text{InSb}}=1$, $\mbox{Re}\left\{\varepsilon_{\cal B}\right\}$ is positive for $T<278$ K, negative for $T>278$ K, and close to null at $T=278$ K. In a similar vein, the real part of $\varepsilon^{s}_{\cal A}$ decreases uniformly as $T$ increases, with $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}$ crossing from positive to negative at $T=280.3$ K. In contrast, the real part of $\varepsilon^{t}_{\cal A}$ is positive for all values of $T$ considered, except that it is negative in the small interval $T\in\left(278.3,279.8\right)$ K. In particular, for the temperature range $T\in\left(260,278.3\right)\cup\left(279.8,280.3\right)$ K the material ${\cal A}$ is not hyperbolic since both $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}$ and $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}$ are positive, but for the temperature range $T\in\left(278.3,279.8\right)\cup\left(280.3,290\right)$ K the material ${\cal A}$ is hyperbolic since the product $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}<0$ for this range. The imaginary parts of $\varepsilon^{s}_{\cal A}$, $\varepsilon^{t}_{\cal A}$, and $\varepsilon_{\cal B}$ are generally small, i.e., less than 0.3, across the temperature range $T\in\left[260,290\right]$ K, with the exception of $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}$ which has a localized peak in the vicinity of $T\in\left(276,284\right)$ K.

4 Surface-wave solutions

In this section representative numerical results are presented of propagation constants, arising as solutions to Eq. (2), and decay constants, as derived from Eqs. (5) and (8), for surface waves guided by the interface $z=0$, across the temperature range $T\in\left[260,290\right]$ K.

4.1 Non-hyperbolic material ${\cal A}$

Let us begin our numerical investigation of surface waves with the temperature regime in which material ${\cal A}$ is not hyperbolic and both materials ${\cal A}$ and ${\cal B}$ are weakly dissipative dielectric materials with $\mbox{Re}\left\{\varepsilon_{\cal B}\right\}>\mbox{Re}\left\{\varepsilon_{\cal A}^{s}\right\}>\mbox{Re}\left\{\varepsilon_{\cal A}^{t}\right\}>0$. The real and imaginary parts of $q$, along with the real parts of $\alpha_{{\cal A}1}$, $\alpha_{{\cal A}2}$, and $\alpha_{\cal B}$, are plotted against orientation angle $\psi$ in Fig. 3 for $T=260$ K with $f^{\cal B}_{\text{InSb}}=1$. At this temperature $\varepsilon_{\mathcal{A}}^{s}=3.4509+0.1222i$, $\varepsilon_{\mathcal{A}}^{t}=1.7106+0.0103i$, and $\varepsilon_{\mathcal{B}}=5.9018+0.2445i$. Over the angular range $0^{\circ}<\psi<90^{\circ}$, only one solution exists for $\psi\in\left(0^{\circ},60^{\circ}\right)$. For this solution $\mbox{Re}\left\{q\right\}>1$, from which it is inferred that the phase speed of this surface wave is smaller than the phase speed of a plane wave in free space. Since both partnering materials are weakly dissipative dielectric materials at this temperature, the surface wave may be called a Dyakonov wave. Although the partnering materials are only weakly dissipative, the effects of dissipation are profound. If the corresponding nondissipative case were considered with real-valued permittivity parameters (i.e., $\varepsilon_{\mathcal{A}}^{s}=3.4509$, $\varepsilon_{\mathcal{A}}^{t}=1.7106$, and $\varepsilon_{\mathcal{B}}=5.9018$) then no Dyakonov waves would exist, since Dyakonov waves can only be found for nondissipative partnering materials if $\varepsilon_{\cal A}^{t}>\varepsilon_{\cal B}>\varepsilon_{\cal A}^{s}$ [4]. Furthermore, the angular existence domain of the solution represented in Fig. 3 is much larger than the angular existence domains typically associated with Dyakonov waves guided by the interfaces of nondissipative materials [6].

Also provided in Fig. 3 are profiles of |\underline{E}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}|, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$, versus $z/\lambda_{\scriptscriptstyle 0}$ with $\underline{n}\in\left\{\hat{\underline{u}}_{x},\hat{\underline{u}}_{y},\hat{\underline{u}}_{z}\right\}$, computed for $\psi=30^{\circ}$ and \underline{\mathcal{E}}_{\,{\cal B}}\mbox{ \tiny{{}^{\bullet}} }\hat{\underline{u}}_{y}=1 V m*-1*. These profiles indicate decay in an approximately exponential manner as the distance $|{z}|$ from the interface $z=0$ increases, in consonance with $\mbox{Re}\left\{\alpha_{{\cal A}1}\right\}=0.03062$, $\mbox{Re}\left\{\alpha_{{\cal A}2}\right\}=0.21760$, and $\mbox{Re}\left\{\alpha_{{\cal B}}\right\}=0.02050$ all being positive. This confirms that these solutions do indeed represent Dyakonov waves that are localized at the interface $z=0$.

By way of comparison and to highlight the effects of InSb in $z<0$, solutions were also computed for the analogous scenario in which the material occupying the half-space $z<0$ is simply air (i.e., $\varepsilon_{\cal B}=1$). The corresponding results are presented in Fig. 4. Like the case presented in Fig. 3, there is only one Dyakonov-wave solution. It exists for all angles $0^{\circ}<\psi<90^{\circ}$. For this solution $\mbox{Re}\left\{q\right\}<1$ as well. The imaginary part of the propagation constant $q$ and the real parts of the decay constants $\alpha_{{\cal A}1}$, $\alpha_{{\cal A}2}$, and $\alpha_{{\cal B}}$, are all generally smaller for Fig. 4 than they are for Fig. 3.

The profiles of |\underline{E}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}|, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$, provided in Fig. 4(f) reveal that the Dyakonov waves represented in Fig. 4 are less tightly localized to the interface $z=0$ than are the Dyakonov waves represented in Fig. 3. Thus, q=0.74302+{\color[rgb]{0,0,0}0.00067i}, $\mbox{Re}\left\{\alpha_{{\cal A}1}\right\}=0.03559$, $\mbox{Re}\left\{\alpha_{{\cal A}2}\right\}=0.00670$, and $\mbox{Re}\left\{\alpha_{{\cal B}}\right\}=0.00075$ for Fig. 4(f), whereas q=1.67076+{\color[rgb]{0,0,0}0.05150i}, $\mbox{Re}\left\{\alpha_{{\cal A}1}\right\}=0.03062$, $\mbox{Re}\left\{\alpha_{{\cal A}2}\right\}=0.21760$, and $\mbox{Re}\left\{\alpha_{{\cal B}}\right\}=0.02050$ for Fig. 3(f).

4.2 Hyperbolic material ${\cal A}$

By raising the temperature, the regime is reached in which $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}>0>\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}>\mbox{Re}\left\{\varepsilon_{\cal B}\right\}$; i.e., material ${\cal B}$ is a plasmonic material while material ${\cal A}$ is a hyperbolic material (since the product $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}<0$). For example, at $T=290$ K with $f^{\cal B}_{\text{InSb}}=1$, the relative-permittivity parameters are $\varepsilon_{\mathcal{A}}^{s}=-2.1565+0.2624i$, $\varepsilon_{\mathcal{A}}^{t}=2.4569+0.0556i$, and $\varepsilon_{\mathcal{B}}=-5.3130+0.5248i$. Plots corresponding to those in Fig. 3 are provided in Fig. 5 for $T=290$ K. In this case there are three solution branches: branch 1 exists for $\psi\in\left(0^{\circ},48^{\circ}\right)$, branch 2 for $\psi\in\left(0^{\circ},44^{\circ}\right)$, and branch 3 for $\psi\in\left(0^{\circ},2^{\circ}\right)$. The profiles of |\underline{E}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}|, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$, provided in Fig. 5 are for the branch-2 solution at $\psi=30^{\circ}$. (The field profiles for the branch-1 solution are provided later in Figs. 7 and 8.) As in Figs. 3 and 4, the field profiles in Fig. 5 reveal an approximately exponential decay as the distance from the interface $z=0$ increases, in both the $+z$ and $-z$ directions, which confirms that these solutions do indeed represent surface waves that are localized to the interface $z=0$.

Most conspicuously, the real part of the propagation constant $q$ for the branch-1 solution is negative for $\psi>11.5^{\circ}$ (and positive for $\psi<11.5^{\circ}$). Thus, for $\psi>11.5^{\circ}$ the phase of the branch-1 surface wave propagates in the direction of $-\hat{\underline{u}}_{x}$ whereas the surface wave attenuates in the direction of $+\hat{\underline{u}}_{x}$. Furthermore, since $\mbox{Re}\left\{q\right\}$ varies uniformly as the angle $\psi$ is varied, as $\psi$ approaches $11.5^{\circ}$ the phase speed of the branch-1 surface wave becomes unbounded because $\mbox{Re}\left\{q\right\}\approx 0$. In contrast, $\mbox{Re}\left\{q\right\}>0$ for the branch-2 and branch-3 solutions for all values of $\psi$. The magnitudes of $\mbox{Re}\left\{q\right\}$ are less than unity on all branches, so the phase speeds are greater than the phase speed of a plane wave in free space for all solutions. In general, $\mbox{Im}\left\{q\right\}$, $\mbox{Re}\left\{\alpha_{{\cal A}1}\right\}$, $\mbox{Re}\left\{\alpha_{{\cal A}2}\right\}$, and $\mbox{Re}\left\{\alpha_{\cal B}\right\}$ in Fig. 5 are slightly larger than the corresponding values presented in Fig. 3.

The results in Fig. 5 can be compared with the corresponding results that arise when material ${\cal B}$ is simply air, i.e., $\varepsilon_{\cal B}=1$, which are presented in Fig. 6. As in Fig. 5, there are three solution branches for Fig. 6: branch 1 exists for $\psi\in\left(0^{\circ},55^{\circ}\right)$, branch 2 for $\psi\in\left(34^{\circ},77^{\circ}\right)$, and branch 3 for $\psi\in\left(83^{\circ},90^{\circ}\right)$. For all three solution branches, the phase speeds are small relative to the phase speed of a plane wave in free space. Unlike Fig. 5, $\mbox{Re}\left\{q\right\}<0$ does not exist in Fig. 6. The profiles of |\underline{E}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}|, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$, provided in Fig. 6 are for the branch-1 solution at $\psi=30^{\circ}$. As compared with the $T=260$ K case of Fig. 4, the surface waves represented in Fig. 6 decay much more rapidly in directions both parallel and normal to the interface $z=0$.

4.3 Negative phase velocity

Let us return to the issue of $\mbox{Re}\left\{q\right\}<0$ which arises in Fig. 5 for branch 1 when $\psi>11.5^{\circ}$. For further exploration, profiles of |\underline{E}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}| and |\underline{H}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}|, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$, versus $z/\lambda_{\scriptscriptstyle 0}$ for $\underline{n}\in\left\{\hat{\underline{u}}_{x},\hat{\underline{u}}_{y},\hat{\underline{u}}_{z}\right\}$ are plotted in Fig. 7, using the same parameter values as for the branch-1 solution in Fig. 5 with $\psi=30^{\circ}$ and \underline{\mathcal{E}}_{\,{\cal B}}\mbox{ \tiny{{}^{\bullet}} }\hat{\underline{u}}_{y}=1 V m*-1*. In addition, plots are provided of the corresponding profiles of the Cartesian components \underline{P}_{\,\ell}(z\hat{\underline{u}}_{z})\mbox{ \tiny{{}^{\bullet}} }\underline{n}, $\ell\in\left\{\mathcal{A},\mathcal{B}\right\}$ and $\underline{n}\in\left\{\hat{\underline{u}}_{x},\hat{\underline{u}}_{y},\hat{\underline{u}}_{z}\right\}$, of the time-averaged Poynting vector

[TABLE]

where the asterisk denotes the complex conjugate. The dot product \mbox{Re}\left\{q\right\}\hat{\underline{u}}_{x}\mbox{ \tiny{{}^{\bullet}} }\underline{P}_{\,\mathcal{A},\mathcal{B}}(z\hat{\underline{u}}_{\,z}) is plotted versus $z/\lambda_{\scriptscriptstyle 0}$ too.

Most significantly, \mbox{Re}\left\{q\right\}\hat{\underline{u}}_{x}\mbox{ \tiny{{}^{\bullet}} }\underline{P}_{\,\mathcal{A},\mathcal{B}}(z\hat{\underline{u}}_{\,z})<0 for much of the range of $z$ considered. Satisfaction of the inequality \mbox{Re}\left\{q\right\}\hat{\underline{u}}_{x}\mbox{ \tiny{{}^{\bullet}} }\underline{P}_{\,\mathcal{A},\mathcal{B}}(z\hat{\underline{u}}_{\,z})<0 signifies NPV. The surface wave is localized to the interface $z=0$, albeit there is significant spreading of the fields into both the half-spaces $z>0$ and $z<0$, but the degree of localization is greater in the half-space $z<0$ than in the half-space $z>0$.

For comparison, the analogous plots for $\psi=5^{\circ}$ are provided in Fig. 8. In this case, $\mbox{Re}\left\{q\right\}>0$ and \mbox{Re}\left\{q\right\}\hat{\underline{u}}_{x}\mbox{ \tiny{{}^{\bullet}} }\underline{P}_{\,\mathcal{A},\mathcal{B}}(z\hat{\underline{u}}_{\,z})>0 for the entire $z$ range considered. Thus, there is no NPV when $\psi=5^{\circ}$. The degree of localization at the interface $z=0$ is similar for Figs. 7 and 8, for both $z<0$ and $z>0$.

Let us denote the extent of the continuous range of angular directions over which NPV is supported by the branch-1 solution in Fig. 5 by $\Delta\psi_{\text{NPV}}$. For the parameter values used in Fig. 5, $\Delta\psi_{\text{NPV}}=36.5^{\circ}$. In Fig. 9, $\Delta\psi_{\text{NPV}}$ plotted against plotted against $T\in\left(260,290\right)$ K for $f^{\cal B}_{\text{InSb}}\in\left\{0.585,0.59,0.6,0.62,0.7,0.8,1\right\}$. At a fixed temperature, $\Delta\psi_{\text{NPV}}$ increases uniformly as $f^{\cal B}_{\text{InSb}}$ increases. At a fixed volume fraction $f^{\cal B}_{\text{InSb}}$, $\Delta\psi_{\text{NPV}}$ varies substantially as $T$ increases, depending upon the value of $f^{\cal B}_{\text{InSb}}$. The lowest temperature at which NPV solutions are supported increases uniformly as $f^{\cal B}_{\text{InSb}}$ decreases.

5 Closing remarks

Hyperbolic materials are associated with a host of exotic electromagnetic phenomenons [14]. Our numerical studies have revealed that such materials, for certain constitutive-parameter ranges, can support the propagation of surface waves with NPV. Furthermore, hyperbolic materials can support a multiplicity of surface waves in a given propagation direction. One of these waves may exhibit NPV while the others exhibit positive phase velocity. In our numerical investigations, NPV was only found to be supported, for certain constitutive-parameter ranges, when partnering material ${\cal A}$ was hyperbolic and partnering material ${\cal B}$ was plasmonic. That is, the conditions

[TABLE]

are necessary conditions for surface-wave NPV, but not sufficient conditions.

Acknowledgments. TGM acknowledges the support of EPSRC grant EP/S00033X/1. AL thanks the Charles Godfrey Binder Endowment at the Pennsylvania State University for ongoing support of his research, with partial funding from US NSF grant number DMS-1619901.