Two Dyakonov–Voigt surface waves guided by a biaxial–isotropic dielectric interface

Electromagnetic surface waves guided by the planar interface of an orthorhombic dielectric material and an isotropic dielectric material were analyzed theoretically and numerically. Both naturally occurring minerals (crocoite, tellurite, and cerussite) and engineered materials were considered as the orthorhombic partnering material. In addition to conventional Dyakonov surface waves, the analysis revealed that as many as two Dyakonov–Voigt surface waves can propagate in each quadrant of the interface plane, depending upon the birefringence of the orthorhombic partnering material. The coexistence of two Dyakonov–Voigt surface waves marks a fundamental departure from the corresponding case involving the planar interface of a uniaxial dielectric material and an isotropic dielectric material for which only one Dyakonov–Voigt surface wave is possible. The two Dyakonov–Voigt surface waves propagate in different directions in each quadrant of the interface plane, with different relative phase speeds and different penetration depths. Furthermore, the localization characteristics of the two Dyakonov–Voigt surface waves at the planar interface are quite different: the Dyakonov–Voigt surface wave with the higher relative phase speed is much less tightly localized at the interface in the isotropic dielectric partnering material.

Any electromagnetic surface wave is guided by the planar interface of two dissimilar partnering materials 1,2 . The characteristics of a surface wave are determined by the constitutive properties of the two partnering materials. If the constitutive properties of both partnering materials are independent of direction (i.e., if both partnering materials are isotropic), then the properties of the surface waves will also be independent of the propagation direction; on the other hand, if the constitutive properties of at least one of the partnering materials depend on direction (i.e., if at least one of the partnering materials is anisotropic), then the properties of the surface waves will also depend on the propagation direction. For example, surface-plasmon-polariton waves are guided by the planar interface of a plasmonic material and a dielectric material (such as air or quartz) 3,4 , the real part of the permittivity dyadic of the plasmonic material being negative definite and the real part of the permittivity dyadic of the dielectric material being positive definite 5 . Likewise, Dyakonov surface waves are guided by the planar interface of two dielectric materials with at least one being anisotropic [6][7][8][9] . Unlike for surface-plasmon-polariton waves, dissipation in both partnering materials for Dyakonov surface waves can be negligibly small, which makes these surface waves attractive candidates for applications involving long-range optical communications.
In addition to Dyakonov surface waves, it was recently discovered that the planar interface of a uniaxial dielectric material and an isotropic dielectric material can guide the propagation of Dyakonov-Voigt surface waves 10,11 . As compared to conventional surface waves such as surface-plasmon-polariton waves and Dyakonov surface waves, Dyakonov-Voigt surface waves have quite different localization characteristics. That is, the fields of conventional surface waves decay exponentially with distance from the interface whereas the decay of fields of a Dyakonov-Voigt surface wave is specified by the product of an exponential function and a linear function of distance from the interface in the anisotropic partnering material. Furthermore, for the planar interface of a uniaxial dielectric material and an isotropic dielectric material, Dyakonov-Voigt surface waves propagate in only one direction for each quadrant of the interface plane; in contrast, conventional surface waves such as surface-plasmon-polariton waves and Dyakonov surface waves propagate in a continuous range of directions in the interface plane. Dyakonov-Voigt surface waves arise as manifestations of exceptional points 12  www.nature.com/scientificreports/ matrix that governs propagation in the anisotropic partnering material 13 , that matrix being non-diagonalizable at these exceptional points. Dyakonov-Voigt surface waves have some similarities to Voigt plane waves that propagate in certain directions in certain unbounded anisotropic 14,15 and bianisotropic 16,17 materials. Notably, Voigt waves arise as manifestations of exceptional points of a propagation matrix, when the matrix is non-diagonalizable 18 . And the fields of Voigt waves decay as the product of a linear function and an exponential function of propagation distance. However, the material supporting Voigt-wave propagation must be dissipative [19][20][21] (or active 22 ), unlike the case for Dyakonov-Voigt surface waves.
In the following sections, Dyakonov-Voigt surface-wave propagation is explored for the planar interface of a biaxial dielectric material, labeled A , and an isotropic dielectric material, labeled B . Both partnering materials A and B are homogeneous and nondissipative. Dyakonov surface waves for interfaces involving biaxial materials have been studied previously [23][24][25] but Dyakonov-Voigt surface waves have not. Our analysis reveals that one or two Dyakonov-Voigt surface waves can be guided for each quadrant of the interface plane, depending upon the birefringence of partnering material A . The manifestation of two Dyakonov-Voigt surface waves (rather than just one 13 ) marks a fundamental departure from the case involving the planar interface of a uniaxial dielectric material and an isotropic dielectric material 10,11 .
As regards the notation adopted, the permittivity and permeability of free space are written as ε 0 and µ 0 , respectively. The free-space wavenumber, wavelength, and impedance are written as k 0 = ω √ ε 0 µ 0 , 0 = 2π/k 0 , and η 0 = √ µ 0 /ε 0 , respectively, with ω being the angular frequency. Single underlining signifies a 3-vector and û x ,û y ,û z is the triad of unit vectors aligned with the Cartesian axes. Matrixes and column vectors are enclosed by square brackets. The superscript T denotes the transpose. The operators Re { • } and Im { • } deliver the real and imaginary parts, respectively, of complex-valued quantities; the complex conjugate is denoted by an asterisk; and dependence on time t is achieved implicitly through exp(−iωt) with i = √ −1.
canonical boundary-value problem 4 × 4 matrix ordinary-differential-equation formalism. In the canonical boundary-value problem for Dyakonov-Voigt (and Dyakonov) surface-wave propagation, material A fills the half-space z > 0 and material B the half-space z < 0 , as illustrated schematically in Fig. 1. Dissipation is assumed to be negligibly small in both materials. Material A is a biaxial dielectric material specified by the relative permittivity dyadic 26,27 with ε Aa , ε Ab , and ε Ac being the principal relative permittivity scalars. Material A belongs to the orthorhombic crystallographic class 28 . In the arena of crystal optics 29 , ε Aa > 1 , ε Ab > 1 , and ε Ac > 1 . The birefringence of material A is the excess of the square root of the largest principal relative permittivity scalar over the square root of the smallest principal relative permittivity scalar 25 . We take material B as an isotropic dielectric material specified by the relative permittivity ε B > 1.
The electromagnetic field phasors for surface-wave propagation are expressed everywhere as 2 (1) ε A = ε Aaûzûz + ε Abûxûx + ε Acûyûy , www.nature.com/scientificreports/ with q being the surface wavenumber. The direction of propagation in the xy plane, relative to the x axis, is prescribed by the angle ψ ∈ [0, 2π) . Substitution of the phasor representations (2) in the source-free Maxwell curl equations yields the 4 × 4 matrix ordinary differential equations 30,31 wherein the column 4-vector and the 4 × 4 propagation matrixes [P A ] and [P B ] depend on ε A and ε B , respectively. The x-directed and y-directed components of the phasors are algebraically connected to their z-directed components 27,31 . fields in material A. The 4 × 4 propagation matrix contains the scalar quantities that depend on the propagation characteristics of the surface wave. The z-directed components of the field phasors can be determined as follows: Dyakonov surface wave. The four eigenvalues of [P A ] can be written as ±α A1 and ±α A2 . The two with positive imaginary parts are wherein the q-independent scalar quantities 4ε Aa Scientific RepoRtS | (2020) 10:12894 | https://doi.org/10.1038/s41598-020-69727-z www.nature.com/scientificreports/ The eigenvectors of [P A ] corresponding to the eigenvalues α A1 and α A2 can be stated as follows: Hence, the general solution of Eq. (3) 1 is given as for fields that decay as z → +∞ . The complex-valued constants C A1 and C A2 have to be determined by application of appropriate boundary conditions at z = 0.
has only two eigenvalues, namely ±α A , each with algebraic multiplicity 2 and geometric multiplicity 1 32 . Unlike the case when material A is uniaxially dielectric 10 , there are two possible values of q that give rise to Dyakonov-Voigt surface waves, namely q + and q − ; from Eqs. (8), these are given as Herein the sign parameter σ = +1 for ψ ∈ (0, π/2) and σ = −1 for ψ ∈ (π/2, π) . When q = q + the eigenvalue α A takes the value α + A , and likewise the eigenvalue α A takes the value α − A when q = q − , with The square root in Eq. (13) is selected in order to ensure that Im α ± A > 0 which is necessary for surfacewave propagation. The propagation directions prescribed by ψ = ℓπ/2 , ℓ ∈ {0, 1, 2, 3} , are irrelevant to Dyakonov-Voigt surface-wave propagation.
An eigenvector of the matrix [P A ] corresponding to the eigenvalue α ± A can be written as The corresponding generalized eigenvector 32 satisfying can be written as   27) is too unwieldy for reproduction here, but q can be numerically extracted from Eq. (27), by the Newton-Raphson method 33 for example. If ψ is replaced by −ψ or by π ± ψ then the dispersion equation (27) is unchanged. Accordingly, in the following numerical investigation of Dyakonov and Dyakonov-Voigt surface waves, attention is restricted to the quadrant 0 ≤ ψ ≤ π/2.

numerical studies
Let us now present and discuss numerical solutions to the dispersion equation (27). Numerical values were chosen for the relative permittivity parameters for this purpose, keeping in mind that the inequality (14) must be satisfied for material A in order for Dyakonov-Voigt surface waves to exist. Candidates materials in the following two categories to function as material A were explored: (i) naturally occurring minerals with relatively modest birefringence and (ii) engineered materials with large birefringence.
naturally occurring material A. Three biaxial minerals were selected to function as material A for the first numerical study: crocoite, tellurite, and cerussite. The relative permittivity scalars of these three minerals, along with the values of their biaxiality Table 1. Both optic ray axes of these minerals are arranged to lie in the xy plane with the y axis as the bisector, the angle δ A being the half-angle between the two optic ray axes. Whereas crocoite has moderately positive biaxiality, tellurite has very weakly positive biaxiality, and cerussite has strongly negative biaxiality 34 . While all three minerals were taken to be orthorhombic 28 , we note parenthetically that crocoite also exists in monoclinic form 35 . Values of the relative phase speed of Dyakonov surface waves are plotted versus the propagation angle ψ ∈ (0, π/2) in Fig. 2a Fig. 2, denoted by a black star, indicates the existence of a Dyakonov-Voigt surface wave with surface wavenumber q = q − . The alternative value of q, namely q = q + , does not deliver a surface wave. Each Dyakonov-Voigt surface wave corresponds to an exceptional point of the matrix [P A ] . For each mineral, the propagation angle for each Dyakonov-Voigt surface wave gets closer to the lower bound of the angular existence domain for Dyakonov surface waves as ε B increases.  Fig. 2 are further delineated in Fig. 3 wherein v p , ε B , and the penetration depths in materials A and B , as given by respectively, are plotted against ψ . For all three minerals considered, the quantities v p , ε B , and B increase monotonically whereas A decreases monotonically, as ψ increases. For each mineral, the range of ψ for Dyakonov-Voigt surface waves in Fig. 3 is bounded below by the value of ψ for which Im {α A } = 0 and bounded above by the value of ψ for which Im {α B } = 0 . As ψ decreases towards its lower bound, the corresponding value of ε B approaches ε Aa .
The localization of the Dyakonov-Voigt surface when A is crocoite and ε B = 5.732 ( ψ = 32.2 • in Fig. 2a) can be observed in Fig. 4. Therein, the quantities |E(zû z ) • n| and |H(zû z ) • n| are plotted versus z/ 0 for n ∈ û x ,û y ,û z and C B3 = 1 V m −1 . Also plotted in Fig. 4 are the quantities P(zû z ) • n for n ∈ û x ,û y ,û z , with Figure 4. Spatial field profiles for a Dyakonov-Voigt surface wave: components of the quantities |E(zû z ) • n| , |H(zû z ) • n| , and P(zû z ) • n are plotted versus z/ 0 . Material A is crocoite, ε B = 5.732 , q = q − = 2.39639k 0 , ψ = 32.2 • , and C B3 = 1 V m −1 . Key: n =û x green solid curves; n =û y red dashed curves; n =û z blue brokendashed curves.  Fig. 4 is rather more tightly localized to the plane z = 0 in material B than in material A ; even so, in both half-spaces the surface wave is essentially contained within |z| < 4 0 . Also, there is no energy flow directed normally to the interface plane in either half-space. Generally, the spatial field profiles for the other eight Dyakonov-Voigt surface waves identified in Fig. 2 are similar to those presented in Fig. 4. Examination of those field profiles revealed the following trend: if the propagation angle of the Dyakonov-Voigt surface wave is closer to the upper bound of the angular existence domain of Dyakonov surface waves then the Dyakonov-Voigt surface wave is less tightly bound in the half-space z > 0 to the interface z = 0 , whereas if the propagation angle of the Dyakonov-Voigt surface wave is closer to the lower bound of the angular existence domain of Dyakonov surface waves then the Dyakonov-Voigt surface wave is less tightly bound in the half-space z < 0 to the interface z = 0. engineered material A. Now we turn to scenarios in which material A is an engineered material with high birefringence. Such engineered materials may be conceptualized as homogenized composite materials 36 , for example, arising from component materials made up from elongated particles or fibers. By judiciously embedding highly elongated component particles or fibers in a host material, homogenized composite materials that exhibit exceedingly high degrees of anisotropy may be realized 37 . The relative permittivity scalars and derivative constitutive parameters of three chosen engineered materials labeled I, II, and III are provided in Table 2, all with strongly positive biaxiality.
Values of the relative phase speed v p of Dyakonov surface waves are plotted versus the propagation angle ψ ∈ (0, π/2) in Fig. 5a-c. For these calculations, (a) Fig. 2 when material A has a small birefringence. Also, the angular existence domains for Dyakonov surface waves in Fig. 5 for the three engineered materials serving as material A are substantially larger than those in Fig. 2 for the three minerals in the same role. On each curve in Fig. 5 there are two stars. These represent the coexistence of two Dyakonov-Voigt surface waves: the black stars denote q = q − and the orange stars denote q = q + . Each pair of Dyakonov-Voigt surface waves corresponds to a pair of exceptional points of the matrix [P A ] . The separation between the propagation angles for q = q − and q = q + is the greatest for the engineered material with the largest difference between ε Aa and ε Ac .
Further insights into the characteristics of the Dyakonov-Voigt surface waves identified in Fig. 5 are gleaned from Fig. 6, wherein v p , ε B , A , and B are plotted against ψ . For all three engineered materials considered, A decreases monotonically and ε B increases monotonically as ψ increases. For q = q + , we have Im {α B } = 0 when ψ takes its smallest value and its largest value. The range of ψ for Dyakonov-Voigt surface waves in Fig. 6 for q = q − is bounded above by the value of ψ for which Im {α B } = 0 . As ψ decreases toward its lower bound, the corresponding value of ε B for q = q − approaches ε Aa .
Lastly, we turn to the localization of the pairs of Dyakonov-Voigt surface waves identified in Fig. 5. The spatial field profiles when material A is engineered material II and ε B = 2.5 are provided in Fig. 7 with C B3 = 1 V m −1 . The corresponding propagation angles for the Dyakonov-Voigt surface waves with q = q − and q + are ψ = 30.1847 • and 24.4068 • , respectively.
The spatial profiles for the Dyakonov-Voigt surface wave for q = q − in Fig. 7 resemble those provided in Fig. 4, except that the latter profiles are more symmetric about the interface z = 0 than the former. Examination of the spatial profiles for all nine Dyakonov-Voigt surface waves for q = q − identified in Fig. 5 revealed the following trend: if the propagation angle of the Dyakonov-Voigt surface wave is closer to the upper bound of the angular existence domain of Dyakonov surface waves then the Dyakonov-Voigt surface wave is less tightly bound in the half-space z > 0 to the interface z = 0 , whereas if the propagation angle of the Dyakonov-Voigt surface wave is closer to the lower bound of the angular existence domain of Dyakonov surface waves then the Dyakonov-Voigt surface wave is less tightly bound in the half-space z < 0 to the interface z = 0.
The spatial profiles for the Dyakonov-Voigt surface wave for q = q + in Fig. 7 are substantially different from those for the Dyakonov-Voigt surface wave for q = q − in the same figure: the surface wave for q = q + is much less tightly localized to the interface in the half-space z < 0 and slightly more tightly localized to the interface Figure 6. Plots of v p , ε B , A , and B versus ψ for Dyakonov-Voigt surface waves. Material A is: engineered material I with q = q + (green solid curves) and q = q − (dark green dashed curves); engineered material II with q = q + (purple solid curves) and q = q − (orange dashed curves); and engineered material III with q = q + (cyan solid curves) and q = q − (blue dashed curves).
Scientific RepoRtS | (2020) 10:12894 | https://doi.org/10.1038/s41598-020-69727-z www.nature.com/scientificreports/ in the half-space z > 0 . In addition, the magnitudes of the energy flows in the x and y directions associated with the surface wave for q = q + are much less than the corresponding energy flows for the surface wave for q = q − . The spatial field profiles for the other eight Dyakonov-Voigt surface waves identified in Fig. 5 for q = q + are qualitatively similar to those presented in Fig. 7. Examination of the spatial profiles for all nine Dyakonov-Voigt surface waves for q = q + identified in Fig. 5 revealed that, as ε B increases, the Dyakonov-Voigt surface wave becomes more tightly bound to the interface in the half-space z > 0 for fixed values of ε Aa , ε Ab , and ε Ac .

closing remarks
Our theoretical and numerical investigations have revealed that the planar interface of an orthorhombic dielectric material and an isotropic dielectric material can guide two Dyakonov-Voigt surface waves in each quadrant of the interface plane, provided that the birefringence of the orthorhombic partnering material is sufficiently large. The two Dyakonov-Voigt surface waves propagate at different phase speeds, in different directions, and with different penetration depths. Also, the Dyakonov-Voigt surface wave with the lower relative phase speed is Figure 7. Spatial field profiles for two Dyakonov-Voigt surface waves: components of the quantities |E(zû z ) • n| , |H(zû z ) • n| , and P(zû z ) • n are plotted versus z/ 0 . Material A is engineered material II, ε B = 2.5 , and C B3 = 1 V m −1 . Left: q = q − = 1.602k 0 and ψ = 30.1847 • . Right: q = q + = 1.58123k 0 and ψ = 24.4068 • . Key: n =û x green solid curves; n =û y red dashed curves; n =û z blue broken-dashed curves.
Scientific RepoRtS | (2020) 10:12894 | https://doi.org/10.1038/s41598-020-69727-z www.nature.com/scientificreports/ much more tightly bound to the interface in the isotropic dielectric partnering material. In contrast, the planar interface of a uniaxial dielectric material and an isotropic dielectric material can guide only one Dyakonov-Voigt surface wave in each quadrant of the interface plane. The canonical boundary-value problem investigated herein has yielded insights into the fundamental characteristics of Dyakonov-Voigt surface-wave propagation. Issues concerning the excitation of these exceptional surface waves in experimental configurations, such as the prism-coupled configuration 2 , shall be pursued in future studies. The numerical results presented in "Naturally occurring material A" are based on naturally occurring partnering materials while the numerical results in "Engineered material A" are based on realistic engineered partnering materials. Accordingly, in principle, the theoretical analysis presented in "Canonical boundary-value problem" and numerical results presented in "Numerical studies" may be readily confirmed by experimental means.
Lastly, let us emphasize that the coexistence of two Dyakonov-Voigt surface waves guided by the planar interface of a biaxial dielectric material and an isotropic dielectric material is wholly independent of the coexistence of two different plane waves, both extraordinary 26 and with different phase speeds, that can propagate inside a biaxial dielectric material 26,27 . Plane waves that propagate in the bulk of a homogeneous material are entirely different from surface waves guided by the planar interface of two dissimilar homogeneous materials.