Open Access
Add to BibSonomy
Abstract
A full-vectorial finite element method is developed to analyze the surface waves propagating at the interface between two media which could be dissipative particularly. The dissipative wave possessing a complex-valued propagation constant can be determined precisely for any given propagation direction and thus the property of losses could be thoroughly analyzed. Besides, by applying a special characteristic of the implicit circular block matrix, we reduce the computational consumptions in the analysis. By utilizing this method, the Dyakonov surface wave (DSW) at the interface between a dielectric and a metal-dielectric multilayered (MDM) structure which serves as a hyperbolic medium is discussed. Its propagation loss is smaller for larger period of the MDM structure but its field becomes less confined to the interface.
© 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Dyakonov surface waves (DSWs), first reported in 1988, are non-dissipative waves that propagate at an interface between a dielectric and a positive uniaxial material with its optical axis parallel to the interface [1,2]. In particular, DSWs propagate at a limited range of directions restricted by the birefringence of the involved anisotropic medium. However, the birefringence is tiny in real materials. Thus, DSWs were first observed with potassium titanyl phosphate (KTP) in experiment in 2009 [3], and people have tried to use some nanostructures such as photonic crystals to artificially create large effective birefringence [4–6], so as to benefit for DSWs to enlarge their ranges of propagation directions and field confinements. One of the simplest uniaxial nanostructures is the stacked layered structure with the effective optical axis being normal to the stacked interface [7]. However, if the nanolaminated structure is composed of two dielectric materials, that with the effective negative uniaxial medium cannot sustain DSWs. To construct an effective positive uniaxial material, metal-dielectric multilayered (MDM) structures were widely used [8–11], which behave as elliptic or hyperbolic media depending on the operating frequency since the metal is dispersive [9]. Recently, the application for the long range guiding was proposed with the experiment [12] and the incorporation of graphene in stratified structures was also demonstrated for tunable applications [13,14].
To analyze the DSWs based on MDM structures, for simplification, people used to replace the MDM structure by the effective medium with its principal permittivities estimated by the effective medium approximation (EMA) with the assumption of homogeneous field distribution. However, the surface plasmon polaritons (SPPs) [15,16] are excited at the metal-dielectric interfaces and cause the nonlocal effect, leading to the failure of the EMA [17,18]. Then, to take into account such nonlocal effect, some modified approximations were proposed [17,19,20] with the effective permittivity tensors as functions of the wave vector. Nevertheless, an additional extraordinary wave was revealed in the MDM structure [18, 21, 22] and particularly cannot be described by the effective medium. Thus, a more rigorous electromagnetic numerical approach like the finite element method (FEM) employed in this work is required to handle this kind of problems. Also, the inhomogeneous field distribution due to SPPs can be obtained in the meanwhile to acquire more interesting information.
Through the approach based on the two-dimensional (2D) FEM, the surface wave propagating along the xz-plane as depicted in Fig. 1(a) can be determined when this wave is non-dissipative. Its unknown propagation constant, k, is composed of the projected components in the x and z directions, kx and kz, respectively, and can be determined by solving the unknown kx with a guessed kz or vice versa. However, when the loss is present, k should be complex-valued and so are kx and kz. It suffers a serious puzzle since the guessed kx or kz is difficult to decide in the complex plane. Although some studies still simulate with a real-valued kx [9,22,23], it is not consistent with the condition in reality. Thus, the loss has been less analyzed precisely despite that this term is a critical characteristic of the propagating wave and some inaccuracy in the calculated k may appear when the ignored loss is enormous. To tackle this problem, we propose a method based on a 3D-FEM to directly calculate the complex-valued k of surface waves along arbitrary directions, and aim to investigate dissipative nanostructures thoroughly.
The rest of this paper is organized as follows. In Section II, the 3D-FEM for computing the oblique propagation is introduced and derived into a quadratic eigenvalue problem (QEP), in which a block circular matrix is observed and can be used to ease the computational consumptions. In Section III, the proposed method is validated and numerical results for the DSWs are given and discussed at the interface between a dielectric and a metal-dielectric multilayered (MDM) structure which serves as a hyperbolic medium. The conclusion is drawn in Section IV.
2. The numerical model
In this section, we describe the rigorous mode solver for the obliquely propagating situation to compute complex-valued k in detail. First, the standard vectorial 3D-FEM formulation is derived with the incorporation of the periodic boundary condition (PBC), turning into a QEP. Then, the utilized prism bases are given and the QEP is solved with the characteristics of the implicit circular block matrix. Finally, the computational complexity for the proposed method is discussed.
2.1. Governing equations
Under the source-free and time-harmonic conditions, the frequency-domain Maxwell’s curl equations can be expressed as
where ω is the angular frequency, μ0 and ∊0 are the permeability and permittivity of free space, respectively, and [μr] and [∊r] represent, respectively, the relative permeability and permittivity tensors of the general medium with possible nine non-zero elements, expressed as From Eqs. (1) and (2), the vectorial wave equation can be derived as where is the wavenumber in free space, the vector Φ is either the electric field E⃗ or the magnetic field H⃗, and the tensors [p] and [q] are, respectively, given by for Φ = E⃗ and for Φ = H ⃗. Through the application of the variational method, the solution of Eq. (4) can be obtained in 3D analysis by extremizing the functional [24] where the del operator ∇ can be split into the transverse and longitudinal parts as In the FEM, the computation domain is divided into several elements and the field distribution within an element can be expressed as where {U}, {V}, and {N} are the vectors for an element composed of the basis functions for the x, y, and z directions, respectively, is the vector composed of transverse coefficient unknowns for an element and is that composed of longitudinal coefficient unknowns, and T denotes transpose. By substituting Eqs. (8) and (9) into Eq. (7) and assembling all element matrices, we can obtain a general matrix equation for the vectors {ϕt} and {ϕz} for the entire domain as where the vectors, {ϕ} and {ψ}, and the matrix, [K], are given by with In Eq. (12), {ψt} and {ψz} are associated with the boundary conditions and will be cancelled at the interior interface of elements by intrinsic field continuity for the electromagnetic (EM) problems. As to the outer boundary of the computational domain, {ψt} and {ψz} are decided by the given boundary conditions (i.e., perfect electronic conductor (PEC), perfect magnetic conductor (PMC), PBC, scattering condition, etc.) and we will discuss the PBCs later.2.2. Periodic boundary conditions
Due to the periodicity, periodic structures can be entirely described by the concept of a unit cell along with the PBCs. Such unit cell as a cuboid is considered in this paper. Based on the Floquet theory, the field relations at the left (L), right (R), up (U), down (D), front (F), and back (B) boundaries for a unit cell are expressed as
where dx, dy, and dz are the lattice constants in the x, y, and z directions, respectively. For simplicity, the phase shift terms, e−jkxdx, e−jkydy, and e−jkzdz, are represented by φx, φy, and φz, respectively. To adequately incorporate these conditions into the numerical model like in [25,26], we can construct two matrices containing the above three phase shift terms. One is [Tr] which relates the coefficient unknowns {ϕ} with the reduced ones {ϕ′} (i.e., {ϕ} = [Tr]{ϕ′}), and the other is [Te] which eliminates the residual terms, {ψ} (i.e., [Te]{ψ} = {0}). By substituting these two matrices into Eq. (10) (i.e., [Te][K][Tr]{ϕ′} = [Te]{ψ} = {0}), a modified matrix equation from [K]{ϕ} = {ψ} is obtained as For a given frequency, the matrices [Ki] with i = 1, 2, .., 8 are constant, and Eq. (15) turns into a problem with three variables, φx, φy, and φz, obviously. Among them, φy is equal to 1 and can be disregarded. This is because we are interested in the waves propagating in the xz-plane as displayed in Fig. 1(a) and the fields of such waves decay along the y direction. Besides, if the structure is invariant in the z direction, dz can be chosen arbitrarily and thus φz can be equal to φx by choosing dz = tan θ · dx technically, where θ represents the angle between the propagation direction and the z-axis. As a consequence, Eq. (15) is derived to a QEP as At last, a full-vectorial FEM based eigenvalue problem is developed for analyzing the wave propagating along a given direction, θ, in the xz-plane of the z-invariant structures. The determined eigenvalue and eigenvector here are respectively the phase factor in the x direction, φx, and the vector, {ϕ′}, composed of the reduced coefficient unknowns of the field distribution. Through φx and θ, the propagation constant, k, and the effective index, Neff, of the determined wave can be obtained, where Neff is defined as k divided by k0.
Fig. 1 (a) Schematic of a planar interface between an anisotropic material and a metal. (b) Computational domain with properly placed PBCs and PMLs.
2.3. The prism elements
The structures we will consider are those with no structural variation along the z direction, and the prism element as depicted in Fig. 2(a) is an institutive type of elements for the 3D-FEM modeling. The functions of bases, {G}, are given by
where {W} and {O} represent the 2D and 1D basis functions, respectively. Following the earlier research [26,27], the basis functions for the curvilinear hybrid triangular element in 2D-FEM are utilized as Wi(x, y), where the edge element based on linear tangential and quadratic normal (LT/QN) vector basis functions and the quadratic nodal element are employed for the transverse and longitudinal field components, respectively. As for the 1D bases in the z direction, in order to satisfy the boundary conditions for the EM waves in the invariant z direction, which require the continuity of the field and its spatial derivative across the boundary, the cubic Hermite spline functions present a good candidate to satisfy these conditions. The mathematical expressions of such bases are expressed as which are illustrated in Fig. 2(b). In Fig. 2(b), it is observed with the characteristics that the nodal field and its derivatives at each terminal can be determined by the coefficients of these four bases. Thus, the continuity of fields and their derivatives along the z direction can be achieved easily. After two individual-base sets have been chosen, our utilized prism bases are obtained consequently.2.4. The quadratic eigenvalue problem solver
In numerical linear algebra, the implicitly restarted arnoldi method (IRAM) [28] is an iterative method typically used to cope with the large sparse linear eigenvalue problems. As to deal with the QEP, one approach was proposed [29] to transform the original n × n QEP
into a 2n × 2n generalized linear eigenvalue problem as where [C], [B], and [A] are n × n coefficient matrices, [I] is an n × n identity matrix, λ and {v} are the eigenvalue and the eigenvector of the oringinal QEP, respectively. The size of the generalized matrix will be twice that of the original one, and larger memory requirements and longer execution time are needed when it is solved by the IRAM directly. To avoid these difficulties, referring to [30], one can observe that solving the linear system (i.e., [T]{r} = {s}) is the most time-consuming portion in the IRAM, and can be simplified by the idea of the block matrix that the matrix equation gives where λ0 is the guessed eigenvalue and is an n × n matrix. With such strategy, we only need to handle a linear system with the original dimension, n × n, and some extra operations (i.e., multiplications), which cost much less consumptions compared with the n × n linear system. In consequence, by this strategy, the complexity of the QEP is almost the same as that of the generalized eigenvalue problem with the same dimension and the consumptions of memory and execution time are reduced.2.5. Linear system solver with circular matrix
For a large sparse linear system, direct methods and iterative methods are two general types of solving technique. Since in the IRAM, the linear system is dealt several times with the same coefficient matrix, we adopted the direct method and thus the factorization is needed only once, which is the most consuming part in the direct method. For better performance, the multifrontal method was utilized in the direct method and was in particular incorporated with parallel programming using the message passing interface (MPI) [24,31]. Moreover, in our FEM calculations, the computational domain can be recognized as a combination of several identical layers as shown in Fig. 3. During the assembly of elements, the temporary matrix for one layer, [Klayer], represents the accumulating contributions from all prism elements in one layer and is the same for each layer, except for the front (blue) and the back (yellow) layers modified by PBCs. Thus, we only stored the content of three matrices even though the number of layers is far beyond three and the consumption of memory here can be reduced, too. Furthermore, after merging [Klayer] of each layer, the matrix can be modified into a block toeplitz matrix as
As mentioned in Subsection 2.2, the eigenvalue indicates to the φx and its absolute value is close to unity when the determined wave is low-dissipative. Thus, it is a good choice for the guessed eigenvalue being unity (λ0 = 1). By this way, the block toeplitz matrix can be turned into a block circular matrix (BCM) with three submatrices [A0], [A1], and [A2]. As a consequence, the linear system can be solved more efficiently by utilizing the properties of the BCM [32].
2.6. Computational consumptions
The computational consumptions for solving a linear system with an ml × ml BCM are discussed in the following. In our proposed method, the l × l submatrices in such BCM are from [Klayer], and their dimension, l, is related to the division in the xy-plane. On the other hand, the order m corresponds to the number of layers. From [32], the BCM is inversed based on the inversion of m matrices with the same dimensions as submatrices. For these m matrices inversed with the multifrontal method, the computational consumption of memory is proportional to l, while the execution times for factorization and for solving the linear system are proportional to l3/2 and l, respectively. Thus, considering the linear system with the BCM, the complexities of computational consumptions are O(ml), O(ml3/2), and O(m2l) for the memory, the factorization time, and the solving time, respectively. Usually, the solving time can be omitted compared with the factorization time. Since the computational consumptions of the aforementioned linear system play the dominant roles in solving the QEP as discussed in Subsection 2.4, the complexities of computational consumptions in our proposed method are thus roughly O(ml) and O(ml3/2) for the memory and the execution time, respectively.
3. Numerical results
3.1. Validation
To validate our proposed method, a relative simple structure, i.e., a planar interface between an anisotropic medium and a metal [33] as depicted in Fig. 1(a), is analyzed first. In the analysis, the unit cell as the computational domain is shown in Fig. 1(b) with PBCs in x- and z- planes at the unit-cell boundaries and PECs in appropriate y- planes. Perfectly matched layers (PMLs) would replace these PECs when solving leaky waves with their fields leaking into the medium. Here, we adopt the anisotropic dielectric with the ordinary and extraordinary refractive indices no = 1.5 and ne = 3, respectively, and the metal with relative permittivity ∊m = −20 − 1.2632j at the wavelength λ = 0.694 μm. The optical axis is along the x direction and the propagation angle, θ, represents the angle between the propagation direction and the z- axis. We examine the method without the PMLs first and make sure the field-amplitudes of the determined waves approximate to 0 at the boundaries of the upside and the downside, so that the wave can be well acquired in the computational domain. The obtained Re[Neff] and loss versus θ are illustrated in Figs. 4(a) and 4(b), respectively. The results from this proposed 3D-FEM are presented by the symbol ‘×’ for comparing with those from the well-developed 2D-FEM [26] shown as the solid lines. Two different kinds of modes can be supported on this surface, one (green line) is a guided mode with its field confined near the interface in both media and the other (blue line) is a leaky mode with part of the wave leaking into the anisotropic substrate due to that its Neff is less than the θ-dependent extraordinary index (Ne(θ)) of the anisotropic material (the red dashed-dotted line). The θ-dependent ordinary index (No(θ)) is also plotted as the red dashed line. We can clearly see the agreement between the results of the two methods for both modes. In order to confirm the proposed method with PMLs, we disregard the metallic loss, so that the field of the leaky mode links into the anisotropic medium far away from the metal-anisotropic medium interface. The losses of the leaky mode are presented in Fig. 4(b) by the red line and crosses for the above two computing methods, and the results from the proposed 3D-FEM still get very good agreement with those from the 2D-FEM. Please note that the characteristic equations for these guided and leaky modes were given as Eqs. (1) and (2) in [33] and the 2D-FEM results were shown to agree excellently with the analytical solution results solved from these characteristic equations in Figs. 2 and 3 of [33]. The field profiles of the leaky mode versus the y position at θ = 20° for the lossless metal case are shown in Fig. 5, where the blue, green, and red lines represent Re[η0Hx], Re[Ex], and |Ex|, respectively. Here η0 is the free-space impedance. The field Hx which possesses only the ordinary wave component shows the great confinement around the interface, while the field Ex which possesses only the extraordinary wave component exhibits the leaky property with the radiating field growing along the y direction in the anisotropic medium and absorbed by the PML. As a consequence, our method provides good reliability for the surface waves with losses and also can handle the leaky waves with PMLs.
Fig. 4 Calculated Results from the 2D- and 3D-FEM for the structure displayed in Fig. 1(a). with no = 1.5, ne = 3, ∊m = −20 − 1.2632j, and λ = 694 nm. (a) Re[Neff] versus θ. (b) Loss in dB/μm versus θ.
Fig. 5 Ex and η0Hy versus y profiles for the wave at θ = 20° obtained from the 3D-FEM without the metallic loss.
3.2. The DSWs at the surface of the semi-infinite MDM structure
In addition, our FEM program is further applied to the analysis of DSWs supported at the interface of a realistic nanostructures comprised of an isotropic medium and an MDM structure [Fig. 6(a)]. The MDM structure is composed of alternatively stacked Ag and PMMA layers with its period, p, varied from 4 nm to 36 nm. The filling factor, f, defined as the fraction of Ag in one period is kept as a constant value of f = 0.25 for all the samples. The refractive index of PMMA is taken as 1.5, and the relative permittivity of Ag is ∊m = −11.7 − 0.83j at λ = 560 nm [34]. By the EMA, the effective parameters are calculated as ∊o = −1.2375 − 0.2075j and ∊e = 3.2043 − 0.0155j. Thus, these MDM structures serve as hyperbolic metamaterials. We first analyzed the characteristics of these MDM structures individually as shown in Fig. 6(b). Due to the uniformity along the yz-plane for such 1D multilayered structures, the projection of arbitrary k to the yz-plane can be simply represented by kz if one rotates the coordinate system properly. Therefore, we calculate the spatial dispersion of the MDM structure along the xz-plane for simplicity, and such analysis does not lose the generality. In addition, we take into account the metallic loss in our calculation and investigate the influence of metallic loss by comparing to the lossless cases (i.e, omitting ∊i,Ag in ∊Ag = ∊r,Ag + j∊i,Ag). As the dissipation of the DSW leads to a complex value of k, the prior related studies [9] obtaining the k with a real-valued kx actually cannot provide the precise results and thus our here-proposed method is necessarily applied to the later discussion of the DSW.
Considering the MDM structure, we ignored the metallic loss first for simplicity. Fig. 7(a) shows the spatial dispersion curves of the e-waves for different p’s, where Nx and Nz are the calculated Neff projected in the x and z directions, respectively. The angle, θ, represents the propagation direction with respect to the z-axis as defined previously. Similar to the previous studies using the exact transfer matrix method [9], our FEM results with different p’s (solid curves) exhibit hyperbolic-like dispersion and start to show differences from the EMA result (‘○’) as p increases. Such differences result from the nonlocal effect due to the excitation of SPPs at the metal-dielectric interfaces. From the perspective of the field distribution, Fig. 7(d) shows the |Hy| profiles of the e-waves propagating at θ = 20° and the degradation of field homogeneity can be observed especially in the metallic region, leading to the violation of the EMA assumption. As several factors such as the penetration depth of the SPP and the thickness of metal, wm, will influence the field distributions, the serious non-uniform field distributions for p = 36 nm leads to a larger difference between our calculation and the EMA result, corresponding to the serious nonlocal effect.
Fig. 7 Characteristics of the e-waves in the MDM structure. (a),(b) show the spatial dispersion contours in the xz-plane without and with the metallic loss, respectively. In particular, their differences, ΔNeff, are given relative to θ in the inset for p = 12, 20, and 36 nm. The point A represents the turning point of the curve for p = 4 nm. (c) Losses versus θ curves. (d) Mode-field profiles of |Hy| along the x-axis for the e-wave propagating at θ = 20° without considering the metallic loss.
Next, when the metallic loss is taken into account, the spatial dispersion curves of the e-waves are shown in Fig. 7(b). It can be seen that the dispersion for p = 4 nm with its curve bending back starting from the point A to the nearby of the origin is significantly different from others and converges to the lossy EMA result. Such agreement between the FEM calculation and the EMA result is consistent with the case ignoring ∊i,Ag due to the weak nonlocal effect for the small p cases. However, the significant difference between the lossless and lossy results for p = 4 nm indicates that the metallic loss should not be ignored when analyzing the dispersion curve for small p samples. In addition, the dispersions for p = 12, 20, and 36 nm show a similar trend to the lossless dispersions. For a close look, the inset of Fig. 7(b) shows the differences, ΔNeff, between the results with and without metallic loss for these three p’s, where ΔNeff = Neff,∊i,Ag =0 − Re[Neff]. One can find ΔNeff is smallest for p = 36 nm, which implies the influence of the metallic loss is reduced for larger p. This is because the modal fields inside the metallic region get smaller when p increases as can be seen in Fig. 7(d). Moreover, as seen in Fig. 7(c), this phenomenon also results in the less loss of the e-wave for the larger p since the loss is attributed to the metallic loss.
Then, the DSWs at the surface of the semi-infinite MDM structure depicted in Fig. 6(a) are analyzed. According to the above investigation on the MDM structure, considering the metallic loss, the spatial dispersion curve for p = 4 nm is close to that from the EMA result and is much different from those for larger p. Thus, in addition to the EMA result, the DSW is displayed for p = 4 nm in Fig. 8(a), where the solid, dashed, and dotted curves represent the dispersions of the DSW, the e-wave in the MDM structure, and the wave in the air, respectively. As we know, a propagating wave will be confined at the interface of two media if its Neff is larger than those of the waves propagating along the same direction in either medium. Therefore, comparing the dispersions of the DSW and the e-wave in the MDM structure, it can be seen that the DSW is getting more confined until θ ≈ 17° as shown by the point A in Fig. 8(a). For larger θ, the DSW turns to be loosely confined and even becomes a leaky wave with its field radiating into the MDM structure from the point B. However, such field radiating into the MDM structure is dissipative due to the metallic loss and cannot leak far from the surface. Thus, the DSWs can still be seen as confined waves until the field starts to leak into the air, which is associated with the point C.
Fig. 8 Spatial dispersion contours of the DSWs (solid curves) in the xz-plane, where the dotted and dashed curves represent the dispersions in the air and in the MDM structure, respectively. (a) The EMA result and p = 4 nm. (b) p = 12, 28, and 72 nm.
Fig. 8(b) shows the DSWs for larger p’s. As can be seen, Neff of the DSW decreases for larger p and once it becomes smaller than that of the e-wave in the MDM structure, the DSW should become a leaky wave which can be observed for p = 20 and 36 nm at smaller θ as shown in the inset of Fig. 8(b) even though they still behave as confined waves when metallic loss is involved. In addition, we can observe that the confinement of the DSW is better for the relatively smaller p, i.e, p = 12 nm. This phenomenon can be seen clearly in Fig. 9, where the magnetic-field profiles of the DSWs for p = 12 and 36 nm at θ = 20° (the points A and B in Fig. 8(b), respectively) are plotted over the same region along the y-axis. The penetration depths of Hy along the y direction in the MDM structure are about 55 and 88 nm for p = 12 and 36 nm, respectively. On the other hand, the losses of the DSWs are shown in Fig. 10(a) for different p’s and are reduced by making p larger. This phenomenon corresponds to the results of the waves in the MDM structure in Fig. 7(b), since the loss is mainly attributed to the components of DSWs in the MDM structure. Therefore, there is a trade-off between the propagation loss and the field confinement for the applications of DSWs. Finally, in Fig. 10(b), we compare the spatial dispersions of DSWs with and without the metallic losses represented by the solid and dashed curves, respectively. The difference between the two types of results is apparent when the DSW propagates along the z-axis (Nx = 0) and is reduced for the larger p or θ, where Neff of the e-wave much deviates from the EMA result.
Fig. 9 Magnetic-field profiles of the DSWs of the structure in Fig. 6(a) for the direction θ = 20°. (a) p = 12 nm. (b) p = 36 nm.
Fig. 10 (a) The losses of DSWs relative to θ. (b) Comparison between the spatial dispersions of DSWs with and without the metallic loss represented by the dashed and solid curves, respectively.
4. Conclusion
We have presented one full-vectorial 3D-FEM based on the eigenvalue algorithm for determining waves propagating at the interface of two media which could possess full anisotropy or be dissipative. The waves are obtained especially with complex propagation constants for a given propagation direction. Thus, the propagation loss which is the important property of the surface waves could be thoroughly analyzed while it was rarely discussed correctly restricted by the lack of proper simulation methods. In particular, we dealt with the proposed method with the help of an implicit block circular matrix, so that the computational consumptions for the memory and the computing time can be reduced and their complexity are O(ml) and O(ml3/2), respectively, where l corresponds to the divisions in the xy-plane and m is the number of the segmentation along the z direction for the applied prism elements in our FEM. It helps us to analyze denser or larger computation domains. Then, this method was carefully confirmed through the surface between a metal and an anisotropic material for both guided waves and leaky waves.
For the simulation results, the DSWs propagating at the surface of the semi-infinite MDM structure were analyzed for different periods, where such MDM structure is regarded as a hyperbolic medium by the EMA. We observed that the spatial dispersion of the DSW cannot be predicted without considering the metallic loss when the period is small. As for larger periods, the influence of the metallic loss on the spatial dispersion is reduced and the results may be approximated by ignoring the metallic loss. In particular, for these DSWs, we found the propagation loss can be reduced by making the periods larger while its field will become less confined to the interface. There is a trade-off between the propagation loss and the field confinement.
As a final remark, for dissipative surface waves, our proposed method shows a great benefit to the reliable results and makes the analysis of losses possible.
Funding
Ministry of Science and Technology (MOST) of the Republic of China (MOST 103-2221-E-002-048-MY2, MOST 105-2221-E-002-138-MY2); Excellent Research Projects of National Taiwan University (104R89081, 105R89081).
References and links
1. M. I. Dyakonov, “New type of electromagnetic wave propagating at an interface,” Sov. Phys. JETP 67, 714–716 (1988).
2. O. Takayama, L. C. Crasovan, S. K. Johansen, D. Mihalache, D. Artigas, and L. Torner, “Dyakonov surface waves: A review,” Electromagnetics 28, 126–145 (2008). [CrossRef]
3. O. Takayama, L. Crasovan, D. Artigas, and L. Torner, “Observation of dyakonov surface waves,” Phys. Rev. Lett. 102, 043903 (2009). [CrossRef] [PubMed]
4. S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, “Effective dielectric constant of periodic composite structures,” Phys. Rev. B 48, 14936–14943 (1993). [CrossRef]
5. F. Genereux, S. W. Leonard, H. M. van Driel, A. Birner, and U. Gösele, “Large birefringence in two-dimensional silicon photonic crystals,” Phys. Rev. B 63, 161101 (2001). [CrossRef]
6. D. Artigas and L. Torner, “Dyakonov surface waves in photonic metamaterials,” Phys. Rev. Lett. 94, 013901 (2005). [CrossRef] [PubMed]
7. A. Yariv and P. Yeh, “Electromagnetic propagation in periodic stratified media. ii. birefringence, phase matching, and x-ray lasers*,” J. Opt. Soc. Am. 67, 438–447 (1977). [CrossRef]
8. J. J. Miret, C. J. Zapata-Rodríguez, Z. Jaksić, S. Vuković, and M. R. Belić, “Substantial enlargement of angular existence range for dyakonov-like surface waves at semi-infinite metal-dielectric superlattice,” J. Nanophoton. 6, 063525 (2012). [CrossRef]
9. C. J. Zapata-Rodríguez, J. J. Miret, S. Vuković, and M. R. Belić, “Engineered surface waves in hyperbolic metamaterials,” Opt. Express 21, 19113–19127 (2013). [CrossRef] [PubMed]
10. C. J. Zapata-Rodríguez, J. J. Miret, J. A. Sorni, and S. Vuković, “Propagation of dyakonon wave-packets at the boundary of metallodielectric lattices,” IEEE J. Sel. Top. Quant. Electron. 19, 4601408 (2013). [CrossRef]
11. O. Takayama, D. Artigas, and L. Torner, “Practical dyakonons,” Opt. Lett. 37, 4311–4313 (2012). [CrossRef] [PubMed]
12. O. Takayama, D. Artigas, and L. Torner, “Lossless directional guiding of light in dielectric nanosheets using Dyakonov surface waves,” Nature Nanotech. 9, 419–424 (2014). [CrossRef]
13. Y. Xiang, J. Guo, X. Dai, S. Wen, and D. Tang, “Engineered surface bloch waves in graphene-based hyperbolic metamaterials,” Opt. Express 22, 3054–3062 (2014). [CrossRef] [PubMed]
14. A. G. Ardakani, M. Naserpour, and C. J. Zapata-Rodríguez, “Dyakonov-like surface waves in the THz regime,” Photon Nanostruct: Fundam. Appl. 20, 1–6 (2016). [CrossRef]
15. J. M. Pitarke, V. M. Silkin, E. V. Chulkov, and P. M. Echenique, “Theory of surface plasmons and surface-plasmon polaritons,” Rep. Prog. Phys. 70, 1 (2007). [CrossRef]
16. S. A. Maier, Plasmonics: Fundamentals and Applications (Springer, 2007).
17. J. Elser, V. A. Podolskiy, I. Salakhutdinov, and I. Avrutsky, “Nonlocal effects in effective-medium response of nanolayered metamaterials,” Appl. Phys. Lett. 90, 191109 (2007). [CrossRef]
18. A. A. Orlov, P. M. Voroshilov, P. A. Belov, and Y. S. Kivshar, “Engineered optical nonlocality in nanostructured metamaterials,” Phys. Rev. B 84, 045424 (2011). [CrossRef]
19. A. V. Chebykin, A. A. Orlov, A. V. Vozianova, S. I. Maslovski, Y. S. Kivshar, and P. A. Belov, “Nonlocal effective medium model for multilayered metal-dielectric metamaterials,” Phys. Rev. B 84, 115438 (2011). [CrossRef]
20. L. Sun, Z. Li, T. S. Luk, X. Yang, and J. Gao, “Nonlocal effective medium analysis in symmetric metal-dielectric multilayer metamaterials,” Phys. Rev. B 91, 195174 (2015). [CrossRef]
21. X. B. Kang, W. Tan, and Z. G. Wang, “Validity of effective medium theory for metal-dielectric lamellar gratings,” Opt. Commun. 284, 4237–4242 (2011). [CrossRef]
22. J. J. Miret, J. A. Sorní, M. Naserpour, A. G. Ardakani, and C. J. Zapata-Rodríguez, “Nonlocal dispersion anomalies of Dyakonov-like surface waves at hyperbolic media interfaces,” Photon Nanostruct: Fundam. Appl. 18, 16–22 (2016). [CrossRef]
23. J. A. Sorní, M. Naserpour, C. J. Zapata-Rodríguez, and J. J. Miret, “Dyakonov surface waves in lossy metamaterials,” Opt. Commun. 355, 251–255 (2015). [CrossRef]
24. J. Jin, The Finite Element Method in Electromagnetics (Wiley-IEEE, 2014).
25. M. Hofer, N. Finger, G. Kovacs, J. Schoberl, S. Zaglmayr, U. Langer, and R. Lerch, “Finite-element simulation of wave propagation in periodic piezoelectric saw structures,” IEEE Trans. Ultrason. Ferroelectr. Freq. Cont 53, 1192–1201 (2006). [CrossRef]
26. S. M. Hsu and H. C. Chang, “Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2d photonic crystals with arbitrary 3d anisotropy,” Opt. Express 15, 15797–15811 (2007). [CrossRef] [PubMed]
27. M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). [CrossRef]
28. R. B. L. Yang and D. C. Sorensen, and C., ARPACK Users Guide: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods (SIAM1997).
29. F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. 43, 235–286 (2001). [CrossRef]
30. A. Nicolet and C. Geuzaine, “Waveguide propagation modes and quadratic eigenvalue problems,” in “6th International Conference on Computational Electromagnetics,” (2006), pp. 1–3.
31. I. Duff, A. Erisman, and J. Reid, Direct methods for sparse matrices (Clarendon, 1986).
32. T. D. Mazancourt and D. Gerlic, “The inverse of a block-circulant matrix,” IEEE Trans. Antennas Propagat. 31, 808–810 (1983). [CrossRef]
33. H. H. Liu and H. C. Chang, “Leaky surface plasmon polariton modes at an interface between metal and uniaxially anisotropic materials,” IEEE Photonics J. 5, 4800806 (2013). [CrossRef]
34. E. D. Palik, G. Ghosh, and A. Press, The electronic handbook of optical constants of solids (Academic, 1999).
