Surface-plasmon-polariton wave propagation supported by anisotropic materials: multiple modes and mixed exponential and linear localization characteristics

TL;DR

This paper investigates surface-plasmon-polariton (SPP) wave propagation at interfaces involving anisotropic materials, revealing multiple modes, unique decay characteristics, and conditions for their existence, expanding understanding of plasmonic wave behavior.

Contribution

It introduces the concept of SPP-V waves with mixed exponential and linear decay in anisotropic materials and provides explicit solutions for their dispersion relations.

Findings

Multiple SPP modes can propagate in anisotropic interfaces.

SPP-V waves exhibit combined exponential and linear decay.

Conditions for SPP-V wave existence are analytically established.

Abstract

The canonical boundary-value problem for surface-plasmon-polariton (SPP) waves guided by the planar interface of a dielectric material and a plasmonic material was solved for cases wherein either partnering material could be a uniaxial material with optic axis lying in the interface plane.Numerical studies revealed that two different SPP waves, with different phase speeds, propagation lengths, and penetration depths, can propagate in a given direction in the interface plane; in contrast, the planar interface of isotropic partnering materials supports only one SPP wave for each propagation direction. Also, for a unique propagation direction in each quadrant of the interface plane, it was demonstrated that a new type of SPP wave--called a surface-plasmon-polariton-Voigt (SPP-V) wave--can exist. The fields of these SPP-V waves decay as the product of a linear and an exponential function of the distance from the interface in the anisotropic partnering material; in contrast, the fields of conventional SPP waves decay only exponentially with distance from the interface. Explicit analytic solutions of the dispersion relation for SPP-V waves exist and help establish constraints on the constitutive-parameter regimes for the partnering materials that support SPP-V-wave propagation.

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}}}_x{\hat{\underline{u}}}_x.$$

$$\left.\begin{array}[]{l}\underline{E}(\underline{r})=\underline{e}(z)\,\exp\left[iq\left(x\cos\psi+y\sin\psi\right)\right]\\[4.0pt] \underline{H}(\underline{r})=\underline{h}(z)\,\exp\left[iq\left(x\cos\psi+y\sin\psi\right)\right]\end{array}\right\}\,$$

$$\left.\begin{array}[]{l}\underline{e}(z)=e_{x}(z)\hat{\underline{u}}_{x}+e_{y}(z)\hat{\underline{u}}_{y}+e_{z}(z)\hat{\underline{u}}_{z}\\[5.0pt] \underline{h}(z)=h_{x}(z)\hat{\underline{u}}_{x}+h_{y}(z)\hat{\underline{u}}_{y}+h_{z}(z)\hat{\underline{u}}_{z}\end{array}\right\}\,$$

$$\left.\begin{array}[]{l}\nabla\times\underline{H}(\underline{r},\omega)+i\omega\varepsilon_{\scriptscriptstyle 0}\underline{\underline{\varepsilon}}_{{\cal A}}\mbox{ \tiny{${}^{\bullet}$} }\underline{E}(\underline{r},\omega)=\underline{0}\\ \vspace{-6pt}\\ \nabla\times\underline{E}(\underline{r},\omega)-i\omega\mu_{\scriptscriptstyle 0}\underline{H}(\underline{r},\omega)=\underline{0}\end{array}\right\}\,,\qquad z>0\,.$$

$$\left.\begin{array}[]{l}\nabla\times\underline{H}(\underline{r},\omega)+i\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{{\cal B}}\underline{E}(\underline{r},\omega)=\underline{0}\\ \vspace{-6pt}\\ \nabla\times\underline{E}(\underline{r},\omega)-i\omega\mu_{\scriptscriptstyle 0}\underline{H}(\underline{r},\omega)=\underline{0}\end{array}\right\}\,,\qquad z<0\,.$$

$$\frac{d}{dz}\left[\underline{f}(z)\right]=i\left[\underline{\underline{P}}_{\cal A}\right]\mbox{ \tiny{${}^{\bullet}$} }\left[\underline{f}(z)\right]\,,\qquad z>0\,,$$

$$\frac{d}{dz}\left[\underline{f}(z)\right]=i\left[\underline{\underline{P}}_{\cal B}\right]\mbox{ \tiny{${}^{\bullet}$} }\left[\underline{f}(z)\right]\,,\qquad z<0\,$$

$$\left[\underline{f}(z)\right]=\left[\begin{array}[]{c}e_{x}(z)\\[5.0pt] e_{y}(z)\\[5.0pt] h_{x}(z)\\[5.0pt] h_{y}(z)\end{array}\right]\,$$

$$\displaystyle\left[\underline{\underline{P}}_{\cal A}\right]=\left[\begin{array}[]{cccc}0&0&\displaystyle{\frac{q^{2}\cos\psi\sin\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon^{s}_{\cal A}}}&\displaystyle{\frac{k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}-q^{2}\cos^{2}\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ 0&0&\displaystyle{\frac{-k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}+q^{2}\sin^{2}\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon^{s}_{\cal A}}}&\displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ \displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{\omega\mu_{\scriptscriptstyle 0}}}&\displaystyle{\frac{-k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}+q^{2}\cos^{2}\psi}{\omega\mu_{\scriptscriptstyle 0}}}&0&0\vspace{8pt}\\ \displaystyle{\frac{k_{\scriptscriptstyle 0}^{2}\varepsilon^{t}_{\cal A}-q^{2}\sin^{2}\psi}{\omega\mu_{\scriptscriptstyle 0}}}&\displaystyle{\frac{q^{2}\cos\psi\sin\psi}{\omega\mu_{\scriptscriptstyle 0}}}&0&0\end{array}\right],$$

$$\left.\begin{array}[]{l}e_{z}(z)=\displaystyle{\frac{q\left[h_{x}(z)\sin\psi-h_{y}(z)\cos\psi\right]}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ h_{z}(z)=\displaystyle{\frac{q\left[e_{y}(z)\cos\psi-e_{x}(z)\sin\psi\right]}{\omega\mu_{\scriptscriptstyle 0}}}\end{array}\right\}\,,\qquad z>0\,.$$

$$\left.\begin{array}[]{l}\alpha_{{\cal A}1}=i\sqrt{q^{2}-k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal A}^{s}}\vspace{8pt}\\ \alpha_{{\cal A}2}=\displaystyle{i\sqrt{\frac{q^{2}\beta-2k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}\varepsilon^{t}_{\cal A}}{2\varepsilon^{s}_{\cal A}}}}\end{array}\right\}\,,$$

$$\left.\begin{array}[]{l}\underline{v}_{{\cal A}1}=\left[\begin{array}[]{c}0\vspace{8pt}\\ \displaystyle{\frac{k_{\scriptscriptstyle 0}\alpha_{{\cal A}1}}{q^{2}\sin\psi\cos\psi}\vspace{8pt}}\\ \displaystyle{\frac{\cot 2\psi}{\eta_{\scriptscriptstyle 0}}+\frac{\csc 2\psi}{\eta_{\scriptscriptstyle 0}}\left(1-\frac{2k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}}{q^{2}}\right)\vspace{8pt}}\\ \eta_{\scriptscriptstyle 0}^{-1}\end{array}\right]\vspace{8pt}\\ \underline{v}_{{\cal A}2}=\left[\begin{array}[]{c}\displaystyle{1-\frac{q^{2}\left(\cos 2\psi+1\right)}{2k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ \displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{k_{\scriptscriptstyle 0}^{2}\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ 0\vspace{8pt}\\ \displaystyle{\frac{\alpha_{{\cal A}2}}{\omega\mu_{\scriptscriptstyle 0}}}\end{array}\right]\end{array}\right\}.$$

$$\left[\underline{f}(z)\right]=C_{\mathcal{A}1}{\underline{v}}_{\mathcal{A}1}\exp \left(i{\alpha }_{\mathcal{A}1}z\right)+C_{\mathcal{A}2}{\underline{v}}_{\mathcal{A}2}\exp \left(i{\alpha }_{\mathcal{A}2}z\right)$$

$$q=\sigma \frac{k_0\sqrt{{\epsilon }_{\mathcal{A}}^s}}{\cos \psi },$$

$${\alpha }_{\mathcal{A}}=i\sigma k_0\sqrt{{\epsilon }_{\mathcal{A}}^s}\tan \psi ,$$

$$\underline{v}_{{\cal A}}=\left[\begin{array}[]{c}0\vspace{8pt}\\ \displaystyle{\frac{i\sigma}{\sqrt{\varepsilon^{s}_{\cal A}}}\vspace{8pt}}\\ 0\vspace{8pt}\\ \eta_{\scriptscriptstyle 0}^{-1}\end{array}\right].$$

$$\left(\left[\underline{\underline{P}}_{\cal A}\right]-\alpha_{{\cal A}}\underline{\underline{I}}\right)\mbox{ \tiny{${}^{\bullet}$} }\underline{w}_{{\cal A}}=\underline{v}_{{\cal A}}$$

$$\underline{w}_{{\cal A}}=\frac{1}{k_{\scriptscriptstyle 0}}\left[\begin{array}[]{c}\displaystyle{\frac{2}{\varepsilon^{t}_{\cal A}-\varepsilon^{s}_{\cal A}}\vspace{8pt}}\\ \displaystyle{\frac{\tan\psi}{\varepsilon^{s}_{\cal A}}\left(\cot^{2}\psi-2\frac{\varepsilon_{\cal A}^{s}-\varepsilon^{t}_{\cal A}\cot^{2}\psi}{\varepsilon_{\cal A}^{s}-\varepsilon^{t}_{\cal A}}\right)\vspace{8pt}}\\ \displaystyle{\frac{2i\sigma\sqrt{\varepsilon^{s}_{\cal A}}}{\eta_{\scriptscriptstyle 0}\left(\varepsilon^{t}_{\cal A}-\varepsilon^{s}_{\cal A}\right)}\vspace{8pt}}\\ 0\end{array}\right].$$

$$\left[\underline{f}(z)\right]=\left[C_{\mathcal{A}1}{\underline{v}}_{\mathcal{A}}+C_{\mathcal{A}2}\left(iz{\underline{v}}_{\mathcal{A}}+{\underline{w}}_{\mathcal{A}}\right)\right]\exp \left(i{\alpha }_{\mathcal{A}}z\right)$$

$$\displaystyle\left[\underline{\underline{P}}_{\cal B}\right]=\left[\begin{array}[]{cccc}0&0&\displaystyle{\frac{q^{2}\cos\psi\sin\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\cal B}}}&\displaystyle{\frac{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}-q^{2}\cos^{2}\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\cal B}}\vspace{8pt}}\\ 0&0&\displaystyle{\frac{-k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}+q^{2}\sin^{2}\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\cal B}}}&\displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\cal B}}\vspace{8pt}}\\ \displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{\omega\mu_{\scriptscriptstyle 0}}}&\displaystyle{\frac{-k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}+q^{2}\cos^{2}\psi}{\omega\mu_{\scriptscriptstyle 0}}}&0&0\vspace{8pt}\\ \displaystyle{\frac{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}-q^{2}\sin^{2}\psi}{\omega\mu_{\scriptscriptstyle 0}}}&\displaystyle{\frac{q^{2}\cos\psi\sin\psi}{\omega\mu_{\scriptscriptstyle 0}}}&0&0\end{array}\right],$$

$$\left.\begin{array}[]{l}e_{z}(z)=\displaystyle{\frac{q\left[h_{x}(z)\sin\psi-h_{y}(z)\cos\psi\right]}{\omega\varepsilon_{\scriptscriptstyle 0}\varepsilon_{\cal B}}\vspace{8pt}}\\ h_{z}(z)=\displaystyle{\frac{q\left[e_{y}(z)\cos\psi-e_{x}(z)\sin\psi\right]}{\omega\mu_{\scriptscriptstyle 0}}}\end{array}\right\}\,,\qquad z<0\,.$$

$${\alpha }_{\mathcal{B}}=-i\sqrt{q^2-k_0^2{\epsilon }_{\mathcal{B}}},$$

$$\left.\begin{array}[]{l}\underline{v}_{{\cal B}1}=\left[\begin{array}[]{c}\displaystyle{1-\frac{q^{2}\cos^{2}\psi}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}}\vspace{8pt}}\\ \displaystyle{-\frac{q^{2}\cos\psi\sin\psi}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}}\vspace{8pt}}\\ 0\vspace{8pt}\\ \displaystyle{\frac{\alpha_{\cal B}}{\omega\mu_{\scriptscriptstyle 0}}}\end{array}\right]\vspace{8pt}\\ \underline{v}_{{\cal B}2}=\left[\begin{array}[]{c}\displaystyle{\frac{q^{2}\cos\psi\sin\psi}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}}\vspace{8pt}}\\ \displaystyle{-1+\frac{q^{2}\sin^{2}\psi}{k_{\scriptscriptstyle 0}^{2}\varepsilon_{\cal B}}\vspace{8pt}}\\ \displaystyle{\frac{\alpha_{\cal B}}{\omega\mu_{\scriptscriptstyle 0}}\vspace{8pt}}\\ 0\end{array}\right]\end{array}\right\}\,.$$

$$\left[\underline{f}(z)\right]=\left(C_{\mathcal{B}1}{\underline{v}}_{\mathcal{B}1}+C_{\mathcal{B}2}{\underline{v}}_{\mathcal{B}2}\right)\exp \left(i{\alpha }_{\mathcal{B}}z\right)$$

$$\left[\underline{f}(0^{+})\right]=\left[\underline{f}(0^{-})\right].$$

$$\left[\underline{\underline{M}}\right]\mbox{ \tiny{${}^{\bullet}$} }\left[\begin{array}[]{c}C_{{\cal A}1}\\ C_{{\cal A}2}\\ C_{{\cal B}1}\\ C_{{\cal B}2}\end{array}\right]=\left[\begin{array}[]{c}0\\ 0\\ 0\\ 0\end{array}\right].$$

$$k_0^2{\epsilon }_{\mathcal{A}}^s\left({\epsilon }_{\mathcal{A}}^s{\alpha }_{\mathcal{B}}-{\epsilon }_{\mathcal{B}}{\alpha }_{\mathcal{A}1}\right)\left({\alpha }_{\mathcal{B}}-{\alpha }_{\mathcal{A}2}\right){\tan }^2\psi$$

$$={\alpha }_{\mathcal{A}1}\left({\alpha }_{\mathcal{B}}-{\alpha }_{\mathcal{A}1}\right)\left({\epsilon }_{\mathcal{A}}^s{\alpha }_{\mathcal{B}}{\alpha }_{\mathcal{A}2}-{\epsilon }_{\mathcal{B}}{\alpha }_{\mathcal{A}1}^2\right),$$

$$\left[\underline{\underline{N}}\right]\mbox{ \tiny{${}^{\bullet}$} }\left[\begin{array}[]{c}C_{{\cal A}1}\\ C_{{\cal A}2}\\ C_{{\cal B}1}\\ C_{{\cal B}2}\end{array}\right]=\left[\begin{array}[]{c}0\\ 0\\ 0\\ 0\end{array}\right]\,.$$

$$\left[2{\epsilon }_{\mathcal{A}}^s\left({\epsilon }_{\mathcal{B}}+{\epsilon }_{\mathcal{A}}^s\right)+\left({\epsilon }_{\mathcal{A}}^s-{\epsilon }_{\mathcal{B}}\right)\left({\epsilon }_{\mathcal{A}}^s+{\epsilon }_{\mathcal{A}}^t\right){\cot }^2\psi \right]$$

$$+2\sqrt{{\epsilon }_{\mathcal{A}}^s}\left({\epsilon }_{\mathcal{A}}^s+{\epsilon }_{\mathcal{B}}\right)\sqrt{{\epsilon }_{\mathcal{A}}^s+\left({\epsilon }_{\mathcal{A}}^s-{\epsilon }_{\mathcal{B}}\right){\cot }^2\psi }=0.$$

Full text

**Surface-plasmon-polariton wave propagation supported by anisotropic materials: multiple modes and mixed exponential and linear localization characteristics **

Chenzhang Zhou

*NanoMM — Nanoengineered Metamaterials Group

Department of Engineering Science and Mechanics

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

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

*NanoMM — Nanoengineered Metamaterials Group

Department of Engineering Science and Mechanics

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

Abstract

The canonical boundary-value problem for surface-plasmon-polariton (SPP) waves guided by the planar interface of a dielectric material and a plasmonic material was solved for cases wherein either partnering material could be a uniaxial material with optic axis lying in the interface plane. Numerical studies revealed that two different SPP waves, with different phase speeds, propagation lengths, and penetration depths, can propagate in a given direction in the interface plane; in contrast, the planar interface of isotropic partnering materials supports only one SPP wave for each propagation direction. Also, for a unique propagation direction in each quadrant of the interface plane, it was demonstrated that a new type of SPP wave — called a surface-plasmon-polariton–Voigt (SPP–V) wave — can exist. The fields of these SPP–V waves decay as the product of a linear and an exponential function of the distance from the interface in the anisotropic partnering material; in contrast, the fields of conventional SPP waves decay only exponentially with distance from the interface. Explicit analytic solutions of the dispersion relation for SPP–V waves exist and help establish constraints on the constitutive-parameter regimes for the partnering materials that support SPP–V-wave propagation.

1 Introduction

Surface-plasmon-polariton (SPP) waves are guided by the planar interface of a dielectric material and a plasmonic material [1]. While SPP waves cannot be excited by direct illumination, their excitation can be readily achieved indirectly via coupling with a prism [2, 3, 4] or surface-relief grating [5], for examples. SPP waves are of major technological importance: they have been widely exploited for optical sensing [6, 5, 7] and microscopy [8, 9], and applications for optical communications [10, 11, 12, 13, 14] and harvesting solar energy [15, 16, 17] are on the horizon.

The theory underpinning SPP-wave propagation is firmly established in the case where the two partnering materials are isotropic [1]. The case where an isotropic plasmonic material is partnered with an anisotropic dielectric material has also been considered previously [18, 19, 20]. However, SPP-wave propagation in the case where an anisotropic plasmonic material is partnered with an isotropic dielectric material has received scant attention from theorists, even though several experimental studies on this case have been reported recently [21, 22, 23, 24, 25, 26].

As we demonstrate in this paper, when anisotropic partnering materials are involved, some previously unreported SPP-wave characteristics emerge. Most notably, two different SPP waves, with different phase speeds, propagation lengths, and propagation depths, can propagate in a given direction in the interface plane. Analogously, this multiplicity of surface waves can also arise in the case of Dyakonov-wave propagation supported by dissipative anisotropic materials [27], and has also been reported in the case of SPP-wave propagation supported by periodically nonhomogeneous dielectric materials [28, 29].

Additionally, we demonstrate that when anisotropic partnering materials are involved, for a unique propagation direction in each quadrant of the interface plane, a new type of SPP wave—with mixed exponential and linear localization characteristics—can exist. We call this surface wave a surface-plasmon-polariton–Voigt (SPP–V) wave, because it is closely related to a singular form of plane waves called Voigt waves that can arise in certain unbounded anisotropic dielectric mediums [30, 31, 32].

A Voigt wave’s amplitude is governed by the product of an exponential function of the propagation distance and a linear function of the propagation distance, in stark contrast to conventional plane waves that propagate in unbounded anisotropic mediums [33, 34]. The existence of Voigt waves was established in early experimental and theoretical studies based on pleochroic crystals such as andalusite, iolite, and alexandrite [30, 31, 35]. But greater scope for Voigt-wave propagation is presented by more complex mediums [36, 37], such as bianisotropic [38] and nonhomogeneous mediums [39]. Furthermore, the directions of Voigt waves can be selected in carefully engineered materials [40, 41, 42, 43, 44] A host anisotropic medium that is either dissipative [35, 45] or active [46] is a prerequisite for Voigt-wave propagation. However, as established in the following, SPP–V-wave propagation is possible for an anisotropic plasmonic material partnered with a non-dissipative (and non-active) dielectric material.

In this paper, a unified theory of SPP-wave propagation and SPP–V-wave propagation is developed by formulating and solving a canonical boundary-value problem. The cases of

• (i)

an anisotropic dielectric material partnered with an isotropic plasmonic material, and

• (ii)

an isotropic dielectric material partnered with an anisotropic plasmonic material,

are considered, with emphasis on new combinations of partnering materials. Explicit analytic solutions of the dispersion relation for SPP–V waves are derived and used to establish constraints on the constitutive-parameter regimes for the partnering materials that allow SPP–V-wave propagation. Representative numerical results are presented to illustrate the theoretical results. And some closing remarks are provided at the end.

In the notation adopted, double underlining denotes 3$\times$3 dyadics while single underlining denotes 3-vectors; double underlining and square parenthesis denotes 4$\times$4 matrixes while single underlining and square parenthesis denotes column 4-vectors. The identity 3$\times$3 dyadic is written as $\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}$ [33], with the triad $\left\{\hat{\underline{u}}_{x},\hat{\underline{u}}_{y},\hat{\underline{u}}_{z}\right\}$ comprising the Cartesian basis vectors. The free-space wavenumber is denoted by $k_{\scriptscriptstyle 0}=\omega\sqrt{\varepsilon_{\scriptscriptstyle 0}\mu_{\scriptscriptstyle 0}}$, wherein $\omega$ is the angular frequency, and the permittivity and permeability of free space are given as $\varepsilon_{\scriptscriptstyle 0}=8.854\times 10^{-12}$ F m*-1* and $\mu_{\scriptscriptstyle 0}=4\pi\times 10^{-7}$ H m*-1*, respectively. The free-space wavelength and impedance are written as $\lambda_{\scriptscriptstyle 0}=2\pi/k_{\scriptscriptstyle 0}$ and $\eta_{\scriptscriptstyle 0}=\sqrt{\mu_{\scriptscriptstyle 0}/\varepsilon_{\scriptscriptstyle 0}}$, respectively. In addition, $i=\sqrt{-1}$.

2 Analysis of surface-wave propagation

2.1 Matrix ordinary differential equations

A general formalism for surface-wave propagation [47] is specialized to develop the canonical boundary-value problem for SPP-wave propagation guided by the planar interface of a uniaxial material, labeled ${\cal A}$, and an isotropic material, labeled ${\cal B}$. The two partnering materials ${\cal A}$ and ${\cal B}$ are both non-magnetic and non-magnetoelectric [34, 48]. Material $\mathcal{A}$ occupies the half-space $z>0$. As this material is a uniaxial dielectric material, it is characterized by an ordinary relative permittivity $\varepsilon_{\mathcal{A}}^{\rm s}$ and an extraordinary relative permittivity $\varepsilon_{\mathcal{A}}^{\rm t}$. With the unit vector $\hat{\underline{u}}_{x}$ pointing in the direction of the optic axis, the relative permittivity dyadic for material ${\cal A}$ is written as [33, 34]

[TABLE]

Material $\mathcal{B}$ occupies the half-space $z<0$ and is characterized by the relative permittivity dyadic $\underline{\underline{\varepsilon}}_{\mathcal{B}}=\varepsilon_{\mathcal{B}}\underline{\underline{I}}$. For SPP waves to be guided by the interface of materials ${\cal A}$ and ${\cal B}$, one of the partnering materials must be a plasmonic material and the other partnering material must be a dielectric material. The canonical boundary-value problem is represented schematically in Fig. 1.

The electromagnetic field phasors that characterize a surface wave are expressed as [47]

[TABLE]

for all $z\in\left(-\infty,\infty\right)$. Herein the complex-valued scalar ${q}$ represents the surface wavenumber; the angle $\psi\in\left[0,2\pi\right)$ prescribes the direction of propagation in the $xy$ plane, relative to the $x$ axis; and the auxiliary phasors

[TABLE]

have complex-valued components, in general.

The field phasors (2) satisfy the source-free, frequency-domain Maxwell curl postulates [33]

[TABLE]

and

[TABLE]

Upon combining with the phasor representations (2), with the Maxwell curl postulates (4) and (5), respectively, the 4$\times$4 matrix ordinary differential equations [49]

[TABLE]

and

[TABLE]

emerge. Herein the column 4-vector

[TABLE]

contains the $x$-directed and $y$-directed components of the auxiliary phasors, which are algebraically connected to the $z$-directed components of the auxiliary phasors [34]. The forms of the 4$\times$4 propagation matrixes $\left[\underline{\underline{P}}_{\cal A}\right]$ in Eq. (6) and $\left[\underline{\underline{P}}_{\cal B}\right]$ in Eq. (7) are determined by the forms of $\underline{\underline{\varepsilon}}_{\cal A}$ and $\underline{\underline{\varepsilon}}_{\cal B}$, respectively.

2.2 Half-space $z>0$

The matrix on the right side of Eq. (6) is given as

[TABLE]

and the $z$-directed components of the auxiliary phasors emerge as

[TABLE]

2.2.1 Nonsingular case

The 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal A}\right]$ has four distinct eigenvalues, namely $\pm\alpha_{{\cal A}1}$ and $\pm\alpha_{{\cal A}2}$, in the nonsingular case. Each eigenvalue has algebraic multiplicity $1$ and geometric multiplicity $1$. The eigenvalues are given by

[TABLE]

wherein the parameter $\beta=\left(\varepsilon^{s}_{\cal A}+\varepsilon^{t}_{\cal A}\right)-\left(\varepsilon^{s}_{\cal A}-\varepsilon^{t}_{\cal A}\right)\cos 2\psi$. The signs of the square-root terms in Eqs. (15) must be selected such that $\mbox{Im}\left\{\alpha_{{\cal A}1}\right\}>0$ and $\mbox{Im}\left\{\alpha_{{\cal A}2}\right\}>0$, in order to ensure that fields decay as $z\to+\infty$. The following pair of eigenvectors of the 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal A}\right]$ match the eigenvalues $\alpha_{{\cal A}1}$ and $\alpha_{{\cal A}2}$, respectively:

[TABLE]

The general solution of the matrix differential eq. (6) for the nonsingular case is provided by

[TABLE]

for $z>0$. The constants $C_{{\cal A}1}$ and $C_{{\cal A}2}$ herein are determined by the boundary conditions at $z=0$.

2.2.2 Singular case

In the singular case, the 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal A}\right]$ has only two eigenvalues, namely $\pm\alpha_{{\cal A}}$. Each eigenvalue has algebraic multiplicity $2$ and geometric multiplicity $1$. This case arises when

[TABLE]

where the sign parameter $\sigma=+1$ for $\psi\in\left(0,\pi/2\right)$ and $\sigma=-1$ for $\psi\in\left(\pi/2,\pi\right)$. The eigenvalues are given by

[TABLE]

wherein the square-root term must be selected to have a positive real part in order to achieve $\mbox{Im}\left\{\alpha_{\cal A}\right\}>0$, which is required in order that fields decay as $z\to+\infty$ [47]. Accordingly, SPP–V-wave propagation is not possible for $\psi\in\left\{0,\pi\right\}$ because $\mbox{Im}\left\{\alpha_{\cal A}\right\}\leq 0$ for $\psi=0$ and $\pi$.

The following eigenvector of the 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal A}\right]$ matches the eigenvalue $\alpha_{{\cal A}}$:

[TABLE]

Furthermore, a generalized eigenvector that satisfies [50]

[TABLE]

is

[TABLE]

The general solution of the matrix differential eq. (6) for the singular case is provided as

[TABLE]

for $z>0$. The constants $C_{{\cal A}1}$ and $C_{{\cal A}2}$ herein are determined by the boundary conditions at $z=0$. Notice that the general solution (23) for the singular case contains a term that is linearly proportional to distance from the interface $z$, which is in stark contrast to the general solution (17) for the nonsingular case in Sec. 2.2.1.

2.3 Half-space $z<0$

The matrix on the right side of Eq. (7) is given as [33, 47]

[TABLE]

and the $z$-directed components of the auxiliary phasors emerge as

[TABLE]

The 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal B}\right]$ has two eigenvalues, namely $\pm\alpha_{{\cal B}}$. Each eigenvalue has algebraic multiplicity $2$ and geometric multiplicity $2$. The eigenvalues are given by

[TABLE]

wherein the sign of the square-root term must be selected such that $\mbox{Im}\left\{\alpha_{{\cal B}}\right\}<0$ to ensure that fields decay as $z\to-\infty$.

The following pair of independent eigenvectors of the 4$\times$4 matrix $\left[\underline{\underline{P}}_{\cal B}\right]$ match the eigenvalue $\alpha_{{\cal B}}$:

[TABLE]

The general solution of the matrix ordinary differential equation (7) is given as

[TABLE]

for $z<0$. Herein the constants $C_{{\cal B}1}$ and $C_{{\cal B}2}$ are determined by the boundary conditions at $z=0$.

2.4 Canonical boundary-value problem

2.4.1 SPP waves

The tangential components of the electric and magnetic field phasors across the interface $z=0$ must be continuous [33]. The four algebraic equations that consequently must be satisfied are compactly expressed as

[TABLE]

By combining Eqs. (17) and (32) with Eq. (33), the following equation emerges:

[TABLE]

The 4$\times$4 characteristic matrix $\left[\underline{\underline{M}}\right]$ herein must be singular for SPP-wave propagation [47]. The corresponding dispersion equation $\left|\left[\underline{\underline{M}}\right]\right|=0$ is equivalent to the equation

[TABLE]

from which the wavenumber $q$ can be numerically extracted, using the Newton–Raphson method [51], for example. From the symmetry of Eq. (2.4.1) it may be inferred that if a SPP wave propagates at the orientation specified by $\psi=\psi^{\star}$, then SPP-wave propagation is also possible for $\psi=-\psi^{\star}$ and $\psi=\pi\pm\psi^{\star}$.

2.4.2 SPP–Voigt waves

As discussed in Sec. 2.4.1, Eq. (33) follows from the continuity of tangential components of the electric and magnetic field phasors across the interface $z=0$ [33]. By combining Eqs. (23) and (32) with Eq. (33), the following equation emerges:

[TABLE]

The 4$\times$4 characteristic matrix $\left[\underline{\underline{N}}\right]$ herein must be singular for surface-wave propagation. And the corresponding dispersion equation $\left|\left[\underline{\underline{N}}\right]\right|=0$ reduces to

[TABLE]

The symmetries of Eq. (2.4.2) are analogous to those of Eq. (2.4.1). Hence, if a SPP–V wave propagates at the orientation specified by $\psi=\psi^{\star}$, then SPP–V waves can also propagate for the orientations $\psi=-\psi^{\star}$ and $\psi=\pi\pm\psi^{\star}$. Observe that Eq. (2.4.2) cannot be satisfied for $\varepsilon^{s}_{\cal A}=\varepsilon_{\cal B}$ unless $\varepsilon^{s}_{\cal A}=\varepsilon_{\cal B}=0$, but this eventuality may be dismissed as it is unphysical.

2.5 Analytical solutions of the SPP–V dispersion equation

Unlike the SPP dispersion equation (2.4.1), the SPP–V dispersion equation (2.4.2) yields analytical solutions for the four variables $\varepsilon^{s}_{\cal A}$, $\varepsilon^{t}_{\cal A}$, $\varepsilon_{\cal B}$, and $\psi$, as follows.

• (i)

For fixed values of $\varepsilon^{t}_{\cal A}$, $\varepsilon_{\cal B}$, and $\psi\in(0,\pi/2)$, SPP–V-wave propagation is possible only when

[TABLE]

wherein the parameters

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

• (ii)

For fixed values of $\varepsilon^{s}_{\cal A}$, $\varepsilon_{\cal B}$, and $\psi\in(0,\pi/2)$, SPP–V-wave propagation is possible only when

[TABLE]

wherein the parameter

[TABLE]

• (iii)

For fixed values of $\varepsilon^{s}_{\cal A}$, $\varepsilon^{t}_{\cal A}$, and $\psi\in(0,\pi/2)$, SPP–V-wave propagation is possible only when

[TABLE]

wherein the parameters

[TABLE]

[TABLE]

and

[TABLE]

• (iv)

For fixed values of $\varepsilon^{s}_{\cal A}$, $\varepsilon^{t}_{\cal A}$, and $\varepsilon_{\cal B}$, SPP–V-wave propagation is possible only when

[TABLE]

In addition, an explicit formula for the surface wavenumber $q$ of a SPP–V wave is provided in Eq. (18).

2.6 Constraints on SPP–V-wave propagation

As well as the analytical solutions represented by Eqs. (38), (44), (46), and (50), constraints on the permittivity parameters of the partnering materials for SPP–V-wave propagation can be developed, as follows. We focus on the $\varepsilon^{t}_{\cal A}$ solution provided in Eq. (44). In the following the possibility of $\psi=\pi/2$ is discounted, because the only solution to emerge from the dispersion relation (2.4.2) for this propagation direction is $\varepsilon^{s}_{\cal A}+\varepsilon_{\cal B}=0$, which is impossible for dissipative materials.

2.6.1 Anisotropic plasmonic material ${\cal A}$ / isotropic dielectric material ${\cal B}$

Suppose that material ${\cal A}$ is plasmonic, i.e., $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}<0$, $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}<0$, $\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}>0$, and $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}>0$, while material ${\cal B}$ is a dielectric material that is generally dissipative, i.e., $\mbox{Re}\left\{\varepsilon_{\cal B}\right\}>0$ and $\mbox{Im}\left\{\varepsilon_{\cal B}\right\}>0$. We consider values of $\varepsilon^{t}_{\cal A}$ that support SPP–V-wave propagation in the direction specified by angle $\psi=\left(\pi/2\right)-\nu$, wherein $0<\nu\ll 1$. Since $\nu$ is taken to be a very small positive parameter, the approximations $\tan\left[\left(\pi/2\right)-\nu\right]\approx 1/\nu$ and $\sec\left[\left(\pi/2\right)-\nu\right]\approx 1$ are justified. Accordingly, Eq. (44) yields

[TABLE]

For fixed values of $\varepsilon^{s}_{\cal A}$ and $\varepsilon_{\cal B}$, the absolute value $|\varepsilon^{t}_{\cal A}|$ becomes increasing large as $\psi$ approaches $\pi/2$, since the possibility $\varepsilon^{s}_{\cal A}+\varepsilon_{\cal B}=0$ is forbidden. In order for $\varepsilon^{t}_{\cal A}$ to lie in the second quadrant of the complex plane, the following inequalities must hold:

[TABLE]

The inequalities (52) may be conveniently recast in terms of the real and imaginary parts of $\varepsilon^{s}_{\cal A}$ and $\varepsilon_{\cal B}$ as

[TABLE]

which are amenable to a geometrical interpretation. In the complex-$\varepsilon_{\cal B}$ plane, the inequality (53)1 prescribes the region outside a circle labeled U, of radius $R_{U}=\left|\varepsilon^{s}_{\cal A}\right|^{2}/\left|\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}\right|$ and centered at the point $C_{U}=\left[\left(\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}\right)^{2}/\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\},-\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}\right]$ in the third quadrant, while the inequality (53)2 prescribes the region outside a circle labeled V, of radius $R_{V}=\left|\varepsilon^{s}_{\cal A}\right|^{2}/\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}$ and centered at the point $C_{V}=\left[-\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\},\left(\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}\right)^{2}/\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}\right]$ in the first quadrant. The straight line connecting the circle centers $C_{U}$ and $C_{V}$ passes through the origin. The distances from the origin to $C_{U}$ and $C_{V}$ are

[TABLE]

respectively, while the distance between $C_{U}$ and $C_{V}$ is

[TABLE]

Some manipulation of the expressions (54) and (55) delivers the relations

[TABLE]

Consequently, it follows that the region in the first quadrant that lies outside the circle centered at $C_{V}$ is entirely outside the circle centered at $C_{U}$. Therefore, the constraint (53)2 is automatically satisfied provided that the constraint (53)1 is satisfied. The circles $U$ and $V$, and the parameter space of $\varepsilon_{\cal B}$ that supports SPP–V-wave propagation, are illustrated in Fig. 2, for a representative example.

In summary: for the case of a plasmonic material ${\cal A}$ and a dissipative dielectric material ${\cal B}$, the constraint (53)1 must be satisfied in order for SPP–V waves to propagate, in directions close to $\psi=\pi/2$. In the special case in which material ${\cal B}$ is nondissipative, i.e., $\mbox{Im}\left\{\varepsilon_{\cal B}\right\}=0$, the constraint (53)1 reduces to

[TABLE]

2.6.2 Anisotropic dielectric material ${\cal A}$ / isotropic plasmonic material ${\cal B}$

Suppose that material ${\cal A}$ is a dielectric material that is generally dissipative, i.e., $\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}>0$, $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}>0$, $\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}>0$, and $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}>0$, while material ${\cal B}$ is plasmonic, i.e., $\mbox{Re}\left\{\varepsilon_{\cal B}\right\}<0$ and $\mbox{Im}\left\{\varepsilon_{\cal B}\right\}>0$. As in Sec. 2.6.1, we consider values of $\varepsilon^{t}_{\cal A}$ that support SPP–V propagation in the direction specified by angle $\psi=\left(\pi/2\right)-\nu$, wherein $0<\nu\ll 1$. As the analysis follows in an analogous manner to that given in Sec. 2.6.1, the details need not be presented here. In the case of a dissipative dielectric material ${\cal A}$ and a plasmonic material ${\cal B}$, the constraint (53)2 must be satisfied in order for SPP–V waves to propagate, in directions close to $\psi=\pi/2$.

In the special case in which material ${\cal A}$ is nondissipative, i.e., $\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}=0$ and $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}=0$, a stronger result can be derived, as follows. Let us introduce the constant

[TABLE]

From the analytical solution (50), $K$ must be real valued and greater than zero for all values of $\psi\in\left(0,\pi/2\right)$. By equating real and imaginary parts, Eq. (58) gives rise to the pair of equations

[TABLE]

The inequality $K>0$ thus yields the twin inequalities

[TABLE]

which together imply the impossible result $-\left(\mbox{Re}\left\{\varepsilon_{\cal B}\right\}\right)^{2}>0$. Therefore, if material ${\cal A}$ is a nondissipative dielectric material and material ${\cal B}$ is plasmonic, then SPP–V-wave propagation is impossible for any value of $\psi$.

3 Numerical studies: SPP-wave propagation

In order to use realistic relative permittivity parameters that can be conveniently varied, a homogenized composite material (HCM) is introduced to play the role of uniaxial material ${\cal A}$. The HCM arises from a mixture of identically oriented needle-shaped particles of two component materials labeled $a$ and $b$. Particles of both component materials are oriented with their long axes parallel to $\hat{\underline{u}}_{x}$. The volume fraction of component material $a$ is denoted by $f_{a}\in[0,1]$ whereas that of component material $b$ is $f_{b}=1-f_{a}$. For the numerical results presented here, $\varepsilon_{a}=-11.63+17.45i$ (cobalt at $\lambda_{\scriptscriptstyle 0}=600$ nm [52]) and $\varepsilon_{b}=5$ (a generic non-dissipative dielectric material).

Provided that the component particles are small in linear dimensions relative to the electromagnetic wavelength(s) involved, the composite material may be regarded as an effectively homogeneous uniaxial material, whose relative permittivity dyadic has the form given in Eq. (1) [53]. The relative permittivity parameters of material ${\cal A}$ are estimated by the Bruggeman homogenization formalism as [54]

[TABLE]

wherein the square-root term in the expression for $\varepsilon_{\mathcal{A}}^{\rm s}$ must be taken to have a positive-valued imaginary part.

The real and imaginary parts of $\varepsilon_{\mathcal{A}}^{\rm s}$ and $\varepsilon_{\mathcal{A}}^{\rm t}$, as provided by Eqs. (61), are plotted against $f_{a}\in\left[0,1\right]$ in Fig. 3. For $f_{a}<0.31$, material ${\cal A}$ is a dissipative dielectric material with $\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm s}\right\}>0$ and $\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm t}\right\}>0$. Specifically, $\varepsilon_{\mathcal{A}}^{\rm s}=6.30+0.94i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=3.34+1.75i$ for $f_{a}=0.1$; $\varepsilon_{\mathcal{A}}^{\rm s}=7.15+2.85i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=1.67+3.49i$ for $f_{a}=0.2$; and $\varepsilon_{\mathcal{A}}^{\rm s}=7.07+5.07i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=0.011+5.24i$ for $f_{a}=0.3$. For $f_{a}>0.71$, material ${\cal A}$ is a uniaxial plasmonic material since $\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm s}\right\}<0$ and $\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm t}\right\}<0$. Specifically, $\varepsilon_{\mathcal{A}}^{\rm s}=-0.17+12.31i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=-6.97+12.56i$ for $f_{a}=0.72$; $\varepsilon_{\mathcal{A}}^{\rm s}=-2.82+13.32i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=-8.30+13.96i$ for $f_{a}=0.80$; and $\varepsilon_{\mathcal{A}}^{\rm s}=-7.01+14.90i$ and $\varepsilon_{\mathcal{A}}^{\rm t}=-9.97+15.71i$ for $f_{a}=0.9$. In the regime $0.31<f_{a}<0.71$, material ${\cal A}$ is classified as a hyperbolic material since $\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm s}\right\}\mbox{Re}\left\{\varepsilon_{\mathcal{A}}^{\rm t}\right\}<0$ [55]. The hyperbolic regime is not considered in the following numerical studies, but this may be an interesting regime to investigate in the future — especially since hyperbolic partnering materials have recently been found to support surface waves with negative phase velocity [55].

3.1 Anisotropic plasmonic material ${\cal A}$ / isotropic dielectric material ${\cal B}$

Consider the case where material ${\cal A}$ is a plasmonic material, specified by the relative permittivity parameters (61) with $f_{a}>0.71$. Material ${\cal B}$ is taken to be a generic non-dissipative dielectric material with relative permittivity $\varepsilon_{\cal B}=5$.

Plots of the normalized phase speed

[TABLE]

and normalized propagation length

[TABLE]

as computed using values of $q$ extracted numerically from Eq. (2.4.1), versus $\psi\in\left(0,\pi/2\right)$ are provided in Fig. 4 for $f_{a}\in\left\{0.72,0.80,0.90\right\}$. Also provided in Fig. 4 are corresponding plots of the normalized penetration depths

[TABLE]

as calculated from Eqs. (15) and (30), respectively. In Fig. 4, the normalized phase speed, as well as the normalized penetration depths in both partnering materials, vary more as $\psi$ increases for smaller values of $f_{a}$. Also, the penetration depths in material ${\cal B}$ are substantially greater than the penetration depths in material ${\cal A}$. This observation is in line with what would be expected for a plasmonic/dielectric interface, regardless of anisotropy of the partnering material ${\cal A}$.

The nature of the SPP waves represented in Fig. 4 is further illuminated in Fig. 5 wherein $|{E_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$ and $|{H_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$ are plotted versus $z/\lambda_{\scriptscriptstyle 0}$ for the case $f_{a}=0.80$ with $\psi=40^{\circ}$. Also plotted are ${P_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}$ which represent the Cartesian components of the time-averaged Poynting vector

[TABLE]

where the asterisk denotes the complex conjugate. For these computations, we fixed $C_{\mathcal{B}1}=1$ V m*-1*. The localization of the SPP wave to the interface $z=0$ is clearly evident, with the degree of localization being substantially greater in the half-space $z>0$ than in the half-space $z<0$, as would be expected from the plots of the penetration depths in Fig. 4.

3.2 Anisotropic dielectric material ${\cal A}$ / isotropic plasmonic material ${\cal B}$

Consider the case where material ${\cal A}$ is a dissipative dielectric material, specified by the relative permittivity parameters (61) with $f_{a}<0.31$. Material ${\cal B}$ is taken to be a plasmonic material with relative permittivity $\varepsilon_{\cal B}=-11.63+17.45i$.

Plots of the normalized phase speed $v_{\text{p}}$ and normalized propagation length $\Delta_{\text{prop}}$, computed using values of $q$ extracted numerically from Eq. (2.4.1), versus $\psi\in\left(0,\pi/2\right)$ are provided in Fig. 6 for $f_{a}\in\left\{0.10,0.20,0.30\right\}$. Also provided in Fig. 6 are corresponding plots of the logarithms of the normalized penetration depths $\Delta_{{\cal A}1}$ and $\Delta_{{\cal A}2}$, and the normalized penetration depth $\Delta_{{\cal B}}$, as calculated from Eqs. (15) and (30), respectively. The SPP-wave solutions represented in Fig. 6 are both qualitatively and quantitatively different to those represented in Fig. 4. Most obviously, two solution branches exist for the case $f_{a}=0.10$: the first branch exists for all $\psi\in\left(0,\pi/2\right)$ whereas the second branch exists only for $0^{\circ}<\psi<8.07^{\circ}$. In contrast, only one solution arises for $f_{a}\in\left\{0.20,0.30\right\}$ and it exists for all values of $\psi\in\left(0,\pi/2\right)$. Also only one solution exists at each value of $\psi$ considered in Fig. 4. Furthermore, the SPP waves on the branch that exists for $0^{\circ}<\psi<8.07^{\circ}$ at $f_{a}=0.10$ penetrate much further into material ${\cal A}$ than do the SPP waves on the branch that exists for $0^{\circ}<\psi<90^{\circ}$ at $f_{a}=0.10$, but the penetration depths into material ${\cal B}$ solutions are similar for solutions on both branches. In addition, the penetration depths into material ${\cal A}$ for the SPP-wave solutions for $f_{a}\in\left\{0.20,0.30\right\}$ are much greater for $0^{\circ}<\psi\lessapprox 10^{\circ}$ than they are for $10^{\circ}\lessapprox\psi<90^{\circ}$, but this is not the case for the penetration depths into material ${\cal B}$.

Further light is shed onto the nature of the SPP waves represented in Fig. 6 by considering the field profiles in the direction normal to the interface $z=0$. In Fig. 7, $|{E_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, $|{H_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, and ${P_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}$ are plotted versus $z/\lambda_{\scriptscriptstyle 0}$ for the case $f_{a}=0.20$ with $\psi=40^{\circ}$, and $C_{\mathcal{B}1}=1$ V m*-1*. Unlike the case presented in Fig. 5, the degree of localization of the SPP wave to the interface $z=0$ in Fig. 7 is substantially greater in the half-space $z<0$ than it is in the half-space $z>0$.

4 Numerical studies: SPP–V-wave propagation

4.1 Anisotropic plasmonic material ${\cal A}$ / isotropic dielectric material ${\cal B}$

Next, material ${\cal A}$ is taken to be a plasmonic material while material ${\cal B}$ is taken to be a non-dissipative dielectric material.

Let $\varepsilon^{s}_{\cal A}=-1+0.1i$. The real and imaginary parts of $\varepsilon^{t}_{\cal A}$ that support SPP–V-wave propagation, as calculated from Eq. (44), are plotted versus $\psi\in\left(0,\pi/2\right)$ in Fig. 8 for $\varepsilon_{\cal B}\in\left\{2.5,5,10\right\}$. The plots of $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}$ converge to $-\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}$, and the plots of $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}$ converge to $-\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}$, as $\psi$ approaches zero. On the other hand, the plots of $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}$ and $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}$ both become unbounded as $\psi$ approaches $\pi/2$.

Also plotted in Fig. 8 are the corresponding plots of the normalized penetration depths in materials ${\cal A}$ and ${\cal B}$, namely $\Delta_{\cal A}$ and $\Delta_{\cal B}$, as defined in Eqs. (64) but with the symbol ${\cal A}\ell$ therein replaced by ${\cal A}$, and as calculated from Eqs. (19) and (30), respectively. Both penetration depths $\Delta_{\cal A}$ and $\Delta_{\cal B}$ converge to zero as $\psi$ approaches $\pi/2$. On the other hand, $\Delta_{\cal A}$ becomes unbounded as $\psi$ approaches zero whereas $\Delta_{\cal B}$ does not. Also, the plotted values of $\Delta_{\cal A}$ are almost independent of $\varepsilon_{\cal B}$ whereas the plotted values of $\Delta_{\cal B}$ are greater for larger values of $\varepsilon_{\cal B}$, especially so at smaller values of $\psi$.

The nature of the SPP–V waves represented in Fig. 8 is further illuminated in Fig. 9 wherein $|{E_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, $|{H_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, and ${P_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}$ are plotted versus $z/\lambda_{\scriptscriptstyle 0}$ for $\varepsilon_{\mathcal{B}}=5$, $\psi=40^{\circ}$ (which corresponds to $\varepsilon^{t}_{\cal A}=-2.85+0.015i$), and $C_{\mathcal{B}1}=1$ V m*-1*. The localization of the SPP–V wave to the interface $z=0$ is clearly demonstrated, with the degree of localization being substantially greater in the half-space $z<0$ than in the half-space $z>0$.

To allow a comparison between SPP and SPP–V waves in the same neighborhood of relative permittivity parameter values, in Fig. 10 the normalized phase speed $v_{\text{p}}$ and normalized propagation length $\Delta_{\text{prop}}$, computed using values of $q$ extracted numerically from Eq. (2.4.1), are plotted versus $\psi\in\left(0,\pi/2\right)$ using the same relative permittivity parameters as were used for Fig. 8, i.e., $\varepsilon^{s}_{\cal A}=-1+0.1i$ and $\varepsilon_{\mathcal{B}}=5$. The value $\varepsilon^{t}_{\cal A}=-2.85+0.015i$ was taken, which corresponds to $\psi=40^{\circ}$ in Fig. 8.

Also provided in Fig. 10 are corresponding plots of the normalized penetration depths $\Delta_{{\cal A}1}$, $\Delta_{{\cal A}2}$, and $\Delta_{{\cal B}}$, as calculated from Eqs. (15) and (30), respectively. The SPP-wave solution represented in Fig. 10 only exists for two disjoint $\psi$ intervals: $0^{\circ}<\psi<59.29^{\circ}$ and $68.30^{\circ}<\psi<90^{\circ}$. The penetration depths $\Delta_{{\cal A}1}$ and $\Delta_{{\cal B}}$ become unbounded as $\psi$ approaches $59.29^{\circ}$ from below and as $\psi$ approaches $68.30^{\circ}$ from above, whereas the penetration depth $\Delta_{{\cal A}1}$ remains bounded for all values of $\psi$. Notice that for Fig. 8, with $\varepsilon^{s}_{\cal A}=-1+0.1i$ and $\psi=40^{\circ}$, the corresponding value of $q/k_{\scriptscriptstyle 0}$ is $0.065+1.307i$, as delivered by Eq. (18), and this value agrees with the value of $q/k_{\scriptscriptstyle 0}$ plotted in Fig. 10 at $\psi=40^{\circ}$. Also, at $\psi=40^{\circ}$ the penetration depths $\Delta_{{\cal A}1}$ and $\Delta_{{\cal A}2}$ in Fig. 10 coincide and these depths agree with $\Delta_{{\cal A}}$ at $\psi=40^{\circ}$ in Fig. 8. And also at $\psi=40^{\circ}$ the penetration depths $\Delta_{{\cal B}}$ in Figs. 8 and 10 coincide. Thus, the solution presented in Fig. 10 represents a SPP wave for $\psi\neq 40^{\circ}$, but it represents a SPP–V wave at the singular orientation $\psi=40^{\circ}$.

4.2 Anisotropic dielectric material ${\cal A}$ / isotropic plasmonic material ${\cal B}$

As discussed in Sec. 2.6, SPP–V-wave propagation is not supported if material ${\cal A}$ is a nondissipative dielectric material and material ${\cal B}$ is a plasmonic material. Accordingly, material ${\cal A}$ is now taken to be a dissipative dielectric material while material ${\cal B}$ is taken to be a plasmonic material.

Let $\varepsilon_{\cal B}=-16.07+0.44i$ (silver at $\lambda_{\scriptscriptstyle 0}=600$ nm [56]). The real and imaginary parts of $\varepsilon^{t}_{\cal A}$ that support SPP–V waves are plotted versus $\psi\in\left(0,\pi/2\right)$ in Fig. 11 for $\varepsilon^{s}_{\cal A}\in\left\{2+0.1i,2+i,2+5i\right\}$. As in Fig. 8, the plots of $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}$ in Fig. 11 converge to $-\mbox{Re}\left\{\varepsilon^{s}_{\cal A}\right\}$, and the plots of $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}$ in Fig. 11 converge to $-\mbox{Im}\left\{\varepsilon^{s}_{\cal A}\right\}$, as $\psi$ approaches zero. As $\psi$ approaches $\pi/2$, the plots of $\mbox{Re}\left\{\varepsilon^{t}_{\cal A}\right\}$ and $\mbox{Im}\left\{\varepsilon^{t}_{\cal A}\right\}$ both become unbounded.

Also presented in Fig. 11 are the corresponding plots of the normalized penetration depths in materials ${\cal A}$ and ${\cal B}$, namely $\Delta_{\cal A}$ and $\Delta_{\cal B}$, as defined in Eqs. (64) but with the symbol ${\cal A}\ell$ therein replaced by ${\cal A}$, and as calculated from Eqs. (19) and (30), respectively. As in Fig. 8, both penetration depths $\Delta_{\cal A}$ and $\Delta_{\cal B}$ in Fig. 11 converge to zero as $\psi$ approaches $\pi/2$. As $\psi$ approaches zero, $\Delta_{\cal A}$ becomes unbounded whereas $\Delta_{\cal B}$ does not. Also, the plotted values of $\Delta_{\cal A}$ and $\Delta_{\cal B}$ are almost independent of $\varepsilon_{\cal B}$.

The nature of the SPP–V waves represented in Fig. 11 is further illuminated in Fig. 12 wherein $|{E_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, $|{H_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}|$, and ${P_{\left\{x,y,z\right\}}(z\hat{\underline{u}}_{z})}$ are plotted versus $z/\lambda_{\scriptscriptstyle 0}$ for $\varepsilon^{s}_{\cal A}=2+i$ and $\psi=40^{\circ}$ (which corresponds to $\varepsilon^{t}_{\cal A}=8.93+0.94i$). For these computations, we fixed $C_{\mathcal{B}1}=1$ V m*-1*. The localization of the SPP–V wave to the interface $z=0$ is easy to see, with the degree of localization being substantially greater in the half-space $z<0$ than in the half-space $z>0$. Also, the SPP–V wave represented in Fig. 12 is localized to the interface $z=0$ to a substantially greater degree than the SPP–V wave represented in Fig. 9.

To allow a comparison between SPP and SPP–V waves in the same neighborhood of relative permittivity parameter values, in Fig. 13 the normalized phase speed $v_{\text{p}}$ and normalized propagation length $\Delta_{\text{prop}}$, computed using values of $q$ extracted numerically from Eq. (2.4.1), are plotted versus $\psi\in\left(0,\pi/2\right)$ using the same relative permittivity parameters as were used for Fig. 11, i.e., $\varepsilon_{\cal B}=-16.07+0.44i$ with $\varepsilon^{s}_{\cal A}=2+i$. The value $\varepsilon^{t}_{\cal A}=8.93+0.94i$ was taken, which corresponds to $\psi=40^{\circ}$ in Fig. 11.

Also provided in Fig. 13 are corresponding plots of the normalized penetration depths $\Delta_{{\cal A}1}$, $\Delta_{{\cal A}2}$, and $\Delta_{{\cal B}}$, as calculated from Eqs. (15) and (30), respectively. The SPP-wave solutions represented in Fig. 13 are organized in two overlapping branches: the first branch exists for $0^{\circ}<\psi<57.37^{\circ}$ while the second branch exists for $47.79^{\circ}<\psi<90^{\circ}$. At each orientation in the overlapping interval $47.79^{\circ}<\psi<57.37^{\circ}$, two SPP waves can exist. The penetration depths $\Delta_{{\cal A}2}$ and $\Delta_{{\cal B}}$ for the solution on the first branch become unbounded as $\psi$ approaches $57.37^{\circ}$ from below, and $\Delta_{{\cal A}2}$ for the solution on the second branch becomes unbounded as $\psi$ approaches $47.79^{\circ}$ from above. In contrast, the penetration depth $\Delta_{{\cal A}1}$ remains bounded for all values of $\psi$, for both solution branches. Notice that for Fig. 11, with $\varepsilon^{s}_{\cal A}=2+i$ and $\psi=40^{\circ}$, the corresponding value of $q/k_{\scriptscriptstyle 0}$ is $1.900+0.448i$, as delivered by Eq. (18), and this value agrees with the value of $q/k_{\scriptscriptstyle 0}$ plotted in Fig. 13 at $\psi=40^{\circ}$ (first-branch solution). Also, at $\psi=40^{\circ}$ the penetration depths $\Delta_{{\cal A}1}$ and $\Delta_{{\cal A}2}$ in Fig. 13 coincide and these depths agree with $\Delta_{{\cal A}}$ at $\psi=40^{\circ}$ in Fig. 11. And also at $\psi=40^{\circ}$ the penetration depths $\Delta_{{\cal B}}$ in Figs. 11 and 13 coincide. Thus, the solution branches presented in Fig. 13 represent SPP waves for $\psi\neq 40^{\circ}$, but the presented solution represents a SPP–V wave at the singular orientation $\psi=40^{\circ}$.

5 Closing discussion

The theoretical underpinnings of SPP-wave propagation supported by isotropic partnering materials are comprehensively described in the literature [1, 47]. And the case where an isotropic plasmonic material is partnered with an anisotropic dielectric material has also been considered previously [18, 19, 20]. However, the same is not true for SPP-wave propagation supported by anisotropic plasmonic materials.

The matter of anisotropic partnering materials is addressed in the preceding sections. In particular, our theoretical and numerical studies have revealed several characteristics of SPP-wave propagation that do not arise for isotropic partnering materials and are therefore attributable to the anisotropy of the partnering materials. These characteristics are as follows:

First, the phase speeds, propagation lengths, and penetration depths for SPP waves supported by anisotropic materials vary with propagation direction. Furthermore, SPP-wave propagation is not necessarily possible for all directions in the interface plane. That is, the angular existence domains of these SPP waves may be less than $360^{\circ}$, as is illustrated in Fig. 10.

Second, for certain relative permittivity-parameter regimes of the partnering materials and for certain propagation directions, the propagation of two distinct SPP waves is supported. As is illustrated in Fig. 13 for example, these two SPP waves have different phase speeds, propagation lengths, and penetration depths.

Third, for a unique direction in each quadrant of the interface plane and for certain relative permittivity-parameter regimes of the partnering materials, the propagation of SPP–V waves is possible. These SPP–V waves are fundamentally different from the conventional SPP waves insofar as the fields of SPP–V waves decay as the product of a linear and an exponential function of the distance from the interface in the anisotropic partnering material; in contrast, the fields of conventional SPP waves decay only exponentially with distance from the interface. A Voigt wave emerges in unbounded anisotropic materials when two planewave modes coalesce to form a singular wave whose amplitude varies with propagation distance [31, 35]. An analogous physical interpretation may be extended to the emergence of SPP–V waves.

The preceding numerical studies were based on realistic values for the relative permittivity parameters of the partnering materials ${\cal A}$ and ${\cal B}$. With these realistic values, the requirements for multiple SPP-wave propagation (i.e., a dielectric partnering material that is both anisotropic and dissipative) and constraints for SPP–V wave propagation (as established in Sec. IIF) could be satisfied. For the purposes of flexibility of presentation, the anisotropic partnering material was taken to be a homogenized composite material whose relative permittivity parameters could be conveniently varied. However, there is no reason to suspect that the requirements for multiple SPP-wave propagation and constraints for SPP–V wave propagation could not satisfied by an anisotropic partnering material with a simpler microstructure. Numerical results (not presented in this paper) that are qualitatively similar to those presented in Secs. 3 and 4 were obtained when the relative permittivity parameters were varied by modest amounts. Parenthetically, where an homogenized composite material is used as the anisotropic plasmonic partnering material, care must be exercised to choose a plasmonic component material that exhibits a moderately high degree of dissipation; we chose cobalt which has a relative permittivity of $-11.63+17.45i$ at $\lambda_{\scriptscriptstyle 0}=600$ nm [52]. This is because conventional homogenization formalisms that can be used to estimate the constitutive parameters of such homogenized composite materials, such as the Bruggeman formalism [53] adopted herein, can give unphysical estimates if the plasmonic component material is only weakly dissipative [57].

In closing, we note that the existence of multiple SPP waves (for homogeneous partnering materials), and the existence of SPP–V waves with mixed exponential and linear localization characteristics, have not been reported in previous SPP studies — and, in particular, these phenomenons have not been reported in previous studies of SPP waves involving anisotropic dielectric materials partnered with isotropic plasmonic materials [18, 19, 20]. The results reported herein have emerged from theoretical and numerical investigations of the corresponding canonical boundary-value problem. While the canonical boundary-value problem represents an idealization that does not take account of finite thicknesses of the partnering materials or the process(es) of excitation of the surface waves, it does yield useful information on the essential physics of surface-wave propagation. Further study is required to explore the excitation and propagation of multiple SPP waves and SPP–V waves for experimental scenarios.

Acknowledgments: This work was supported by EPSRC (grant number EP/S00033X/1) and US NSF (grant number DMS-1619901). AL thanks the Charles Godfrey Binder Endowment at the Pennsylvania State University for ongoing support of his research.