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

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.


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 [6][7][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].
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 surfacewave propagation with NPV, whereas the low-temperature regime cannot.
In the notation adopted, vectors are underlined once and dyadics twice. An exp(−iωt) dependence on time t is implicit, with i = √ −1 and angular frequency ω. The triad û x ,û y ,û z contains the Cartesian unit vectors, while the position vector is denoted by r = xû x + yû y + zû z . The free-space wavenumber is k 0 = ω √ ε 0 µ 0 , with ε 0 and µ 0 being the permittivity and permeability of free space, respectively; and λ 0 = 2π/k 0 is the free-space wavelength.

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 nonmagnetic and non-magnetoelectric [38,39]. As the theory underpinning this setting is amply described in earlier works [9,40], only the bare essentials need be presented here. The half-space z > 0 is occupied by a uniaxial material labeled A. Its relative permittivity dyadic is given by wherein I =û xûx +û yûy +û zûz is the identity dyadic [41]. The optic axis of material A, whose direction is given by the unit vectorû = cos ψû x + sin ψû y , lies in the xy plane, at an angle ψ relative to the x axis. The half-space z < 0 is occupied by material B whose relative permittivity dyadic is ε B 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û x is considered. The electric field phasor in the half-space z > 0 is given in terms of amplitude vectors E A1 and E A2 as with the wavevectors and the decay constants The electric field phasor in the half-space z < 0 is given in terms of an amplitude vector E B as with the wavevector and decay constant Crucially, the inequalities Re {α A } > 0, ∈ {1, 2}, and Re {α B } > 0 must be satisfied for surface-wave propagation.
The normalized propagation constant q is delivered by solving the corresponding dispersion relation [9] Numerical methods, such as the Newton-Raphson method [42], are generally needed for this task. The following symmetries exist: if q = q satisfies Eq. (9) at the angle ψ = ψ then q = q also satisfies Eq. (9) for ψ = −ψ and ψ = π ± ψ . Accordingly, in the proceeding presentation of numerical solutions to Eq. (9) (in Sec. 4), only the range 0 ≤ ψ ≤ π/2 need be considered.

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 temperaturesensitive semiconductor InSb. The relative permittivity of InSb in the terahertz regime is provided by the Drude model [11,43] as Herein the plasma frequency ω p = N q 2 e /0.015 ε 0 m e is determined by the electronic charge q e = −1.60 × 10 −19 C and mass m e = 9.11 × 10 −31 kg, while the high-frequency relative permittivity ε ∞ = 15.68 and the damping constant γ = π × 10 11 rad s −1 . The temperature T (in K) dependence of ε InSb arises via the intrinsic carrier density (in m −3 ) [44][45][46] which depends upon the bandgap E g = 0.26 eV and the Boltzmann constant k B = 8.62 × 10 −5 eV 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 where f A InSb ∈ (0, 1) 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 ε HDPE . In the long-wavelength regime, the material filling z < 0 may be regarded as an isotropic continuum whose relative permittivity is estimated as where f B InSb ∈ (0, 1) is the volume fraction of InSb in the half-space z < 0. For the purposes of our calculations, the volume fraction f A InSb = 0.5 was selected. The frequency was fixed at f = 2.0 THz (i.e., λ 0 = 0.15 mm). At this frequency, the relative permittivity of HDPE is approximately constant, i.e., ε HDPE ≈ 2.387 + 0.006i, over the temperature range T ∈ [260, 290] K [52]. Plots of the real and imaginary parts of the relative permittivity parameters of the partnering materials A and B versus temperature are provided in Fig. 2

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

Non-hyperbolic material A
Let us begin our numerical investigation of surface waves with the temperature regime in which material A is not hyperbolic and both materials A and B are weakly dissipative dielectric materials with Re The real and imaginary parts of q, along with the real parts of α A1 , α A2 , and α B , are plotted against orientation angle ψ in Fig. 3 for T = 260 K with f B InSb = 1. At this temperature ε s A = 3.4509 + 0.1222i, ε t A = 1.7106 + 0.0103i, and ε B = 5.9018 + 0.2445i. Over the angular range 0 • < ψ < 90 • , only one solution exists for ψ ∈ (0 • , 60 • ). For this solution Re {q} > 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., ε s A = 3.4509, ε t A = 1.7106, and ε B = 5.9018) then no Dyakonov waves would exist, since Dyakonov waves can only be found for nondissipative partnering materials if ε t A > ε B > ε s A [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 |E (zû z ) • n|, ∈ {A, B}, versus z/λ 0 with n ∈ û x ,û y ,û z , computed for ψ = 30 • and E B •û 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 Re {α A1 } = 0.03062, Re {α A2 } = 0.21760, and Re {α B } = 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., ε B = 1). The corresponding results are presented in Fig. 4. Like the case presented in Fig. 3, there is only one Dyakonovwave solution. It exists for all angles 0 • < ψ < 90 • . For this solution Re {q} < 1 as well. The imaginary part of the propagation constant q and the real parts of the decay constants α A1 , α A2 , and α B , are all generally smaller for Fig. 4 than they are for Fig. 3.
The profiles of |E (zû z ) • n|, ∈ {A, B}, 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 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.

Hyperbolic material
Most conspicuously, the real part of the propagation constant q for the branch-1 solution is negative for ψ > 11.5 • (and positive for ψ < 11.5 • ). Thus, for ψ > 11.5 • the phase of the branch-1 surface wave propagates in the direction of −û x whereas the surface wave attenuates in the direction of +û x . Furthermore, since Re {q} varies uniformly as the angle ψ is varied, as ψ approaches 11.5 • the phase speed of the branch-1 surface wave becomes unbounded because Re {q} ≈ 0. In contrast, Re {q} > 0 for the branch-2 and branch-3 solutions for all values of ψ. The magnitudes of Re {q} 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, Im {q}, 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 B is simply air, i.e., ε B = 1, which are presented in Fig. 6. As in Fig. 5, there are three solution branches for Fig. 6: branch 1 exists for ψ ∈ (0 • , 55 • ), branch 2 for ψ ∈ (34 • , 77 • ), and branch 3 for ψ ∈ (83 • , 90 • ). 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, Re {q} < 0 does not exist in Fig. 6. The profiles of |E (zû z ) • n|, ∈ {A, B}, provided in Fig. 6 are for the branch-1 solution at ψ = 30 • . 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.

Negative phase velocity
Let us return to the issue of Re {q} < 0 which arises in Fig. 5 for branch 1 when ψ > 11.5 • . For further exploration, profiles of |E (zû z ) • n| and |H (zû z ) • n|, ∈ {A, B}, versus z/λ 0 for n ∈ û x ,û y ,û z are plotted in Fig. 7, using the same parameter values as for the branch-1 solution in Fig. 5 with ψ = 30 • and E B •û y = 1 V m −1 . In addition, plots are provided of the corresponding profiles of the Cartesian components P (zû z ) • n, ∈ {A, B} and n ∈ û x ,û y ,û z , of the time-averaged Poynting vector where the asterisk denotes the complex conjugate. The dot product Re {q}û x • P A,B (zû z ) is plotted versus z/λ 0 too. Most significantly, Re {q}û x • P A,B (zû z ) < 0 for much of the range of z considered. Satisfaction of the inequality Re {q}û x • P A,B (zû 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 ψ = 5 • are provided in Fig. 8. In this case, Re {q} > 0 and Re {q}û x • P A,B (zû z ) > 0 for the entire z range considered. Thus, there is no NPV when ψ = 5 • . 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 ∆ψ NPV . For the parameter values used in Fig. 5, ∆ψ NPV = 36.5 • . In Fig. 9, ∆ψ NPV plotted against plotted against T ∈ (260, 290) K for f B InSb ∈ {0.585, 0.59, 0.6, 0.62, 0.7, 0.8, 1}. At a fixed temperature, ∆ψ NPV increases uniformly as f B InSb increases. At a fixed volume fraction f B InSb , ∆ψ NPV varies substantially as T increases, depending upon the value of f B InSb . The lowest temperature at which NPV solutions are supported increases uniformly as f B InSb decreases.

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 A was hyperbolic and partnering material B was plasmonic. That is, the conditions are necessary conditions for surface-wave NPV, but not sufficient conditions.