Unidirectional plasmonic edge modes on general two-dimensional materials

TL;DR

This paper explores the directional control of edge plasmons in 2D materials through polarization and dipole positioning, extending theoretical models to anisotropic and hyperbolic systems, and validating findings with simulations.

Contribution

It introduces a comprehensive analysis of unidirectional edge plasmon modes in general 2D materials, including anisotropic and hyperbolic systems, with validation against exact solutions.

Findings

Identification of conditions for unidirectional edge plasmons.

Extension of existing models to anisotropic and hyperbolic materials.

Validation of theoretical predictions with simulations.

Abstract

We investigate the field and spin-momentum coupling of edge plasmons hosted by general two-dimensional materials and identify sweet spots depending on the polarisation plane, ellipticity and the position of an electric dipole relative to the plane and edge. Exciting the dipole at these sweet spots by propagating light leads to uni-directional propagating edge plasmons or edge modes are totally suppressed. We also extent previous approximate treatments [A. Fetter Phys. Rev. B 32, 7676 (1985)] to include anisotropy and hyperbolic systems, elucidating its predictions for the existence of edge modes. A thorough assessment of the approximate description is carried out, comparing its spin-momentum coupling features in the near field with exact results from Wiener-Hopf techniques. Simulations are also performed confirming the overall picture. Our results shed new light on the quest of chiral plasmonics in 2D materials and should be relevant for future experiments.

Equations

$${\mathbf{A}}_m(x,z)=A_0\sum _{\mathbf{q}}{e}^{i\tau qx-q_m^{\prime }|z|}({\mathbf{e}}_x+sgn(z)\tau i\frac{q}{q_m^{\prime }}{\mathbf{e}}_z),$$

$$\psi =\frac{E_{\perp }}{i{E}_{\parallel }}=\frac{E_z}{i{E}_x}=sgn(z)\tau \frac{q}{q_m^{\prime }}$$

$$-{\nabla }^2\Phi (\mathbf{r},x)=\rho (x)\delta (\mathbf{r}),$$

$$g(\mathbf{r})=\frac{1}{2\pi }{\int }_0^{\infty }dk\frac{k{J}_0(kr)}{k^2+q^2}=\frac{1}{2\pi }{K}_0(qr)$$

$$\psi =\frac{E_{\perp }}{i{E}_{\parallel }}=\frac{E_r}{i{E}_x}=\tau \frac{K_0^{\prime }(qr)}{K_0(qr)},$$

$$\displaystyle\chi(y)=\chi f(y)\;,\;\rm{with}\;\chi=\left(\begin{array}[]{cc}\chi_{xx}&\chi_{xy}\\ \chi_{yx}&\chi_{yy}\end{array}\right)\;,$$

$$L_m^q(y,z)=\frac{1}{2}{\int }_{-\infty }^{\infty }\frac{d{q}_y}{2\pi }{e}^{i{q}_{y}y}{e}^{-q_m|z|}\frac{q_1^{\prime }{q}_2^{\prime }}{(q_1^{\prime }{)}^2{q}_2^{\prime }+(q_2^{\prime }{)}^2{q}_1^{\prime }}.$$

$${\mathbf{a}}_m(q,y,z)=$$

$$\times$$

$$L_m^q(y,z=0)\approx \frac{1}{\sqrt{8}}{e}^{-\sqrt{2}q|y|}$$

$$\displaystyle{\bf a}_{m,\pm}=\left(\begin{array}[]{c}\tau\\ \pm i\tilde{\kappa}_{\pm}\\ (-1)^{m}i\tilde{\kappa}_{z}^{\pm}\end{array}\right)$$

$${\tilde{\kappa }}_y\equiv {\tilde{\kappa }}_{-}=$$

$$2-{\tilde{\kappa }}_y^2=$$

$$\det (\chi )={\chi }_{xx}{\chi }_{yy}-{\chi }_{xy}{\chi }_{yx}>0,$$

$${\chi }_{yy}>0,$$

$${\omega }_e^F(q)=\sqrt{\frac{{\chi }_{xx}{\chi }_{yy}-{\chi }_{xy}^{\prime 2}-{\chi }_c^2}{2{\chi }_{yy}+{\chi }_{xx}-2\sqrt{2}\tau {\chi }_c}\frac{q}{{ϵ}_0ϵ}}$$

$$\det (\chi )={\chi }_{xx}{\chi }_{yy}-{\chi }_{xy}{\chi }_{yx}<0.$$

$$2\sqrt{1-{\tilde{\kappa }}_y^2}={\tilde{\chi }}_{xx}-{\tilde{\kappa }}_y^2{\tilde{\chi }}_{yy}-i\tau ({\tilde{\chi }}_{xy}+{\tilde{\chi }}_{yx}){\tilde{\kappa }}_y.$$

$${\psi }_{-}^F=sgn(z)\tau ,{\psi }_{+}^F=-\tau \sqrt{2}.$$

$$n_F(x,y)={\rho }_F(y)\frac{e^{i\tau qx}}{\sqrt{2\pi }},$$

$${\rho }_F(y)=|q|e^{|qy|/\sqrt{2}}\Theta (-y)+\sqrt{2}\delta (y),$$

$$\tan (\theta )=\frac{E_z}{E_y}=-sgn(z)\frac{{\kappa }_z}{{\kappa }_y}=-sgn(z)\frac{{\tilde{\kappa }}_z}{\sqrt{1-{\tilde{\kappa }}_z^2}}.$$

$${\tilde{\kappa }}_z={\left[\frac{{\omega }_e(q)}{{\omega }_b(q)}\right]}^2.$$

$$n(x,y)=\rho (\tilde{y})\frac{e^{iqx}}{\sqrt{2\pi }},$$

$$\tilde{\rho }(\tilde{k})={\tilde{\rho }}_{o}f(\tilde{k}),$$

$$f(\tilde{k})=\exp \left[-\frac{X_{+}(\tilde{k})}{2\pi i}\right],$$

$$\begin{gathered} X_{+}(\tilde{k})= \\ {\int }_{-\infty }^{+\infty }d{k}^{\prime }\left[\frac{1}{k^{\prime }-\tilde{k}-i{0}^{+}}-\frac{1}{k^{\prime }-i{0}^{+}}\right]\ln \left[{\tilde{\kappa }}_z^{-1}\sqrt{1+k^{\prime 2}}-1\right], \end{gathered}$$

$$f(\tilde{k})=\frac{\sqrt{1-i\tilde{k}}}{1-i\tilde{k}/{\tilde{\kappa }}_y}g(\tilde{k}),$$

$$g(\tilde{k})=\exp \left[-\frac{\tilde{k}}{i\pi }{\int }_1^{\infty }dx\frac{\pi /2-\arctan ({\tilde{\kappa }}_z^{-1}\sqrt{x^2-1})}{x(x-i\tilde{k})}\right].$$

$$\begin{gathered} {\rho }_{approx}=\sqrt{2\pi }{\tilde{\kappa }}_y\sqrt{1-{\tilde{\kappa }}_y}{e}^{{\tilde{\kappa }}_y\tilde{y}}\mathrm{erf}(\sqrt{(1-{\tilde{\kappa }}_y)|\tilde{y}|})\Theta (-\tilde{y}) \\ +\sqrt{2}{\tilde{\kappa }}_y\frac{e^{\tilde{y}}}{\sqrt{|\tilde{y}|}}\Theta (-\tilde{y}), \end{gathered}$$

$$\mathcal{E}(x,y,z)=\frac{e^{iqx}}{\sqrt{2\pi }}\mathbf{E}(\tilde{y},\tilde{z}),$$

Full text

I Introduction

Chiral nano-photonics has attracted increasing attention especially in view of new routes to information technology Novotny and Hecht (2012), among which plasmonic circuitries have the advantage of a reduced length scale, compatible with nano-electronics Ozbay (2006). Arguably, the most studied case of chiral plasmonics is the one of magneto-plasmons Fetter (1985); Volkov and Mikhailov (1988); Mikhailov (1995); Wang et al. (2011); Yan et al. (2012); Crassee et al. (2012); Kumada et al. (2013, 2014); Cohen and Goldstein (2018), leading to infrared topological plasmons in periodically patterned monolayer graphene Jin et al. (2017) and superlattices Pan et al. (2017). Non-reciprocal plasmons can further be obtained by an external source-drain bias Sabbaghi et al. (2015, 2018); Bliokh et al. (2018), leading to collimated plasmon beam in the case of large drift current Duppen et al. (2016). And a non-trivial Berry curvature of the band-structure can also lead to chiral propagation of the edge mode by appropriately breaking time-reversal symmetry via, e.g., circularly polarized light Kumar et al. (2016); Song and Rudner (2016). Exciting chiral (edge) plasmons at fixed frequency and momentum via a grating geometry would then lead to uni-directional (edge) plasmon propagation. Nevertheless, the opposite plasmons dispersions hardly differ from each other such that both propagation directions are usually excited. Also, plasmons are mostly excited by illuminating a point-dipole via propagating light that excites plasmons of all wave numbers and will thus also launch plasmons in the opposite direction.

Here, we will investigate a more fundamental route to chiral nano-plasmonics based on the inherent spin-momentum coupling Lee et al. (2012); Lin et al. (2013); Bliokh et al. (2015a); Mechelen and Jacob (2016); Lodahl et al. (2017), which (in some context) can be linked to the quantum spin Hall effect of light Bliokh et al. (2015b). From Maxwell’s equations it follows that for strongly confined plasmons, the electric field is either left or right circularly polarized depending on the direction of propagation. Providing a one-dimensional wave guide, circularly polarized light at oblique incidence can then only excite plasmons that propagate in one direction according to the spin-momentum coupling Lee et al. (2012). Anisotropic materials such as black phosphorousLow and Avouris (2014); Goncalves et al. (2017) or hyperbolic materialsMa et al. (2018) could thus provide a intrinsic propagation direction leading to plasmons with a high collimation of electric field energy Nemilentsau et al. (2019). The concept can also be extended to Janus and Huygens dipoles Picardi et al. (2018) and hyperbolic materials.

Nowadays, the most prominent 2D plasmonic material is graphene, mainly due to its high intrinsic mobility/life-times and easy tuneability of its plasmonic resonances Koppens et al. (2011); Grigorenko et al. (2012); Stauber (2014); Gonçalves and Peres (2016); Basov et al. (2016); Low et al. (2017). They can be excited and detected by means of infrared nano imaging Fei et al. (2012); Chen et al. (2012). Since graphene’s plasmons are strongly confined, a left circularly polarized dipole in the $xz$-plane on top of graphene placed at $z=0$ will excite plasmons propagating only in the positive $x$-direction and not in the opposite, negative $x$-direction. But graphene is an isotropic material and the circularly polarized dipole in the $xz$-plane is elliptically polarized in all other directions. Since elliptically polarized light can be decomposed into a left and right circularly polarized component, a dipole usually excites plasmons that propagate in both, forward and backward direction, albeit with different intensities. This naturally limits the usefulness of the spin-momentum locking in terms of collimation of the plasmonic energy current in only one direction.

We will, therefore, mainly consider plasmonic edge-modes of an isotropic medium such as graphene since they provide a natural wave guide at the border of the sample Nikitin et al. (2011). Recently, they have been revealed by infrared nano imaging, i.e., scattering-type scanning near-field optical microscopy (s-SNOM) Fei et al. (2015); Nikitin et al. (2016). It is generally argued that graphene edges should be separated by an effective working distance to avoid the overlapping of localized plasmon modes with the bulk modes, important for the design of graphene-based plasmonic circuits and devices.Xu et al. (2017) But here, we will focus on the edge modes as principle players of plasmonic circuits similar to the proposal of Refs. Pan et al. (2017); Bisharat and Sievenpiper (2017). In fact, our results will in general be valid for any two-dimensional system with an arbitrary conductivity tensor. For hyperbolic materials,Ma et al. (2018) i.e., for systems with a conducting and non-conduction principle axis, we can rule out the existence of edge plasmons within the Fetter approximation.Fetter (1985)

The paper is organised as follows. In Section II, we discuss the spin-momentum coupling for simple geometries modelling bulk and edge plasmonic modes. In Section III, the spin-momentum coupling for a two-dimensional plasmonic edge mode is discussed analytically within the FetterFetter (1985) approximation, extended to cover general elliptic and hyperbolic cases, studying the condition for the emergence of edge states. We also compare the approximate results for the near field with the exact Wiener-HopfVolkov and Mikhailov (1988) treatment. In Section IV, we then verify our conclusions by COMSOL simulations launching unidirectional plasmonic edge modes by circularly polarized dipoles. We close with conclusions and summary and add two appendices providing details of the analytical treatment.

II Spin-momentum coupling of simple geometries

The spin-momentum coupling is a consequence of the transverse nature of a confined gauge-field. But it can also be obtained from the static Poisson equation. In the following, we will discuss simple geometries related to bulk and edge plasmons and introduce basic concepts such as the ellipticity.

II.1 Plane interface

For an interface at $z=0$ and choosing in-plane momentum ${\bm{q}}=\tau q{\bm{e}}_{x}$ with $q>0$, unit vector along the $x$ axis ${\bm{e}}_{x}$ and propagation sense $\tau=\pm$, the vector potential in half-space $z<0$ ($m$=1) and half-space $z>0$ ($m$=2) can be written as

[TABLE]

where $q_{m}^{\prime}=\sqrt{q^{2}-\epsilon_{m}(\omega/c)^{2}}$ is the transverse wave number, $\epsilon_{m}$ the dielectric constants of medium $m$, and ${\bm{e}}_{z}$ the unit vector along the $z$ axis.

Due to the phase factor $i\frac{q}{q_{m}^{\prime}}$, there is a rotation of the polarisation vector in the $xz$-plane and the sense of rotation depends on the propagation direction. We can define the ellipticity of the electric fields ${\bf E}=-\partial_{t}{\bf A}$ as the ratio between the field perpendicular and parallel to the propagation divided by a phase $i$:

[TABLE]

Generally, the polarisation is elliptically polarized. In the case of graphene, however, the plasmons are strongly confined, such that $q^{\prime}\approx q$ and therefore $\psi\approx\pm 1$. This makes graphene’s plasmon mode * circularly* polarized. Also note that for bulk modes, the polarisation plane is always normal to the interface.

II.2 Infinite line charge

Let us consider the potential $\Phi$ of an infinitely long line charge along the $x$ axis described by the following Poisson equation:

[TABLE]

where ${\bm{r}}$ represents the perpendicular vector position. This could be a very elementary model for a plasmonic edge mode of a half-infinite plane as in Fig. 1, if we take a propagating wave in $x$-direction with wave number $q$, i.e., $\rho(x)=e^{i\tau qx}$. Writing $\Phi({\bm{r}},x)=e^{i\tau qx}g({\bm{r}})$, Eq. (3) then reduces to the problem of finding the Green’s function $g({\bm{r}})$ of a 2D massive particle:

[TABLE]

with Bessel function of first kind $J_{0}$, and modified Bessel function of second kind $K_{0}$. This model can be extended to a charged cylinder of radius $a$ without major changes.

Again, we are interested in the ratio of the absolute value of the electric fields with respect to the perpendicular and parallel direction with respect to the propagation:

[TABLE]

where $K_{0}^{\prime}$ means derivative. Note that this ratio depends on the radius $r$, and on the propagation sense $\tau$. For $r\gg 1/q$, this field becomes circularly polarized. Further, the polarisation plane is always that containing the line charge and the observation point.

This analysis suggests that the field of an edge mode outside the 2D material should be nearly circularly polarized and the polarisation plane is rotating around the wire. In the following section, we will discuss how this picture is modified when the extended, half-infinite, 2D material plane is taken into account.

III Spin-momentum coupling of two-dimensional plasmonic edge modes

Charged graphene, as any two-dimensional electron gas, hosts bulk density fluctuation excitations called plasmons. In addition, edge localized plasmons emerge when the geometry is restricted to a half-infinite plane. We choose the geometry depicted in Fig. 1, where the material half-plane is located at $z=0$ and occupies the region $y<0$, so the edge coincides with the $x$ axis. Although a complete characterization of this edge mode in the static limit can be achieved by means of Wiener-HopfVolkov and Mikhailov (1988) techniques, here we follow the simplified approach pioneered by FetterFetter (1985); Wang et al. (2011); Song and Kats (2016) to study its near field.

We are interested in the near fields produced by self-sustained currents inside a 2D conducting plane ${\bf j}=-\chi{\bf A}$ with local current-current response depending on the (say) $y$-component:

[TABLE]

where $f(y)$ denotes an arbitrary function. Note that we also introduce an arbitrary response matrix that could describe any possible physical system in 2D including hyperbolic materials. The response matrix $\chi$, related to the conductivity by $\chi=-i\omega\sigma$, can be termed the Drude matrix, as it reduces to the Drude weight in the isotropic case. The bulk plasmons propagating in $x$-direction between two dielectrics $\epsilon_{1}$ and $\epsilon_{2}$ have the dispersion $\omega_{b}^{2}=\tfrac{\chi_{xx}q}{2\varepsilon_{0}\epsilon}$ with $\epsilon=(\epsilon_{1}+\epsilon_{2})/2$.

We will solve the Maxwell equations with the constitutive equation of Eq. (8). In the Appendix A , we outline how retardation effects can be included in the basic equations characterised by the generating function

[TABLE]

With this, we can relate the vector potential ${\bf A}_{m}({\bm{r}})=A_{0}\sum_{q}e^{i\tau qx}{\bf a}_{m}(q,y,z)$ to the current:

[TABLE]

with $\Omega^{2}=\frac{\epsilon_{1}q_{2}^{\prime}+\epsilon_{2}q_{1}^{\prime}}{q_{1}^{\prime}+q_{2}^{\prime}}\frac{\omega^{2}}{c^{2}}$. We note that this equation is including retardation which has not been discussed so far. Retardation effects are important in the context of chirality as introduced by Tang and Cohen.Tang and Cohen (2010, 2011)

III.1 Analytical approximation

To proceed analytically, we will set $f(y)=\Theta(-y)$, neglect retardation effects and approximate the generating function as:

[TABLE]

a procedure first proposed by FetterFetter (1985) that amounts to replacing the non-local Kernel of Eq. 12 by one amenable to solution as a differential equation for the $y$ variable in the sheet plane. We then obtain ${\bf a}_{m}(q,\pm|y|,z)={\bf a}_{m,\pm}e^{-\kappa_{\pm}|y|}e^{-\kappa_{z}^{\pm}|z|}$:

[TABLE]

where we have introduced the dimensionless quantities $\tilde{\kappa}_{\pm}=\kappa_{\pm}/q$ and $\tilde{\chi}_{n,n^{\prime}}=\frac{q}{\epsilon_{0}\epsilon\omega^{2}}\chi_{n,n^{\prime}}$, with $\tilde{\kappa}_{+}=\sqrt{2}$. Defining to facilitate the notation in the material region $\tilde{\kappa}_{y}=\tilde{\kappa}_{-}$, then the putative edge mode emerges as the solution of the following two equations:

[TABLE]

as shown in the Appendix A.

III.2 Solutions

The solutions to the edge mode equations depend on the material Drude matrix, itself a function of frequency, a fact that apparently precludes the possibility of saying something general. Nevertheless, one can make progress by pretending that the Drude matrix at the target frequency remains constant for all frequencies. In this way, one can attain rather general conclusions on the existence of edge modes that are valid, at least, on the chosen $\omega=\text{constant}$ line in phase space, for instance, the horizontal lines of Fig. 2. The conclusions of such analysis for different material classes are described in what follows.

III.2.1 Elliptic (or ordinary) systems

We assume a non-absorbing material with, for the moment, time-reversal symmetry (TRS). The Drude matrix of Eq. 8 is then made of real entries with $\chi_{xy}=\chi_{yx}$ owing to TRS. The ordinary nature of the material is taken to mean a positive determinant for $\chi$:

[TABLE]

If, in addition, one assumes $\chi$ to be positive definite, a rotation to principal axes produces two positive eigenvalues which, in the isotropic case, reduces the matrix to a scalar, the Drude weight of an ordinary metal. Such system supports bulk plasmons with elliptic dispersion given by $\omega^{2}_{b}(q)=\frac{1}{2\epsilon_{0}\epsilon q}q_{n}\chi_{n,n^{\prime}}q_{n}^{\prime}$, leading in the isotropic case to the well-known expression $\omega^{2}_{b}(q)=\frac{Dq}{2\epsilon_{0}\epsilon}$, where $q$ is here the bulk wavevector $q=\sqrt{q_{x}^{2}+q_{y}^{2}}$, and $D$, the Drude weight.Stauber (2014)

As shown in the Appendix B, an edge state below the continuum of bulk plasmons for systems satisfying Eq. 21 will exist whenever

[TABLE]

which together with Eq. 21, amounts to a positive definite Drude matrix $\chi$, i.e., edge plamons will emerge whenever the bulk supports elliptic (ordinary or metallic) plasmons. This situation is schematically depicted on the (a) panel of Fig. 2. In the isotropic case considered by FetterFetter (1985), the edge mode is given by $\omega_{e}^{F}(q)=\sqrt{\tfrac{2}{3}}\,\omega_{b}(q)$, reasonably close to the exact solution of Volkov and MikhailovVolkov and Mikhailov (1988) $\omega_{e}^{VM}(q)\approx 0.9\,\omega_{b}(q)$. Therefore, our results generalize to the anisotropic case the existence of edge modes in ordinary systems.

In the absence of TRS such as, for instance, with a perpendicular magnetic field present, the off-diagonal entries of $\chi$ can acquire an imaginary component, $\chi_{xy}=\chi^{\prime}_{xy}+i\chi_{c}$, and $\chi_{yx}=\chi^{\prime}_{yx}-i\chi_{c}$, where the real part is again symmetric $\chi^{\prime}_{xy}=\chi^{\prime}_{yx}$, and the imagnary part, $\chi_{c},$ would be proportional to the cyclotron frequency. In that case, the analysis of Appendix B predicts, to lowest order in $\chi_{c}$, a linear splitting of the edge mode depending on the propagation direction, just as in the isotropic caseFetter (1985); Volkov and Mikhailov (1988).

Quite generally, if an edge plasmon exists, it is shown in the Appendix B that it satisfies the following dispersion relation within the present approximation:

[TABLE]

III.2.2 Hyperbolic systems

The Drude matrix of a hyperbolic, non-absorbing material with TRS is again real and symmetric, but now with a negative determinant:

[TABLE]

Such system also supports bulk plasmons with dispersion given by $\omega^{2}_{b}(q)=\frac{1}{2\epsilon_{0}\epsilon q}(q^{2}_{+}\chi_{+}+q^{2}_{-}\chi_{-})$, where $q_{\pm}$ are the wavevector components along principal axes, with corresponding Drude entries now of opposite sign: $\chi_{+}>0$ and $\chi_{-}<0$, therefore, a hyperbolic dispersion. Edge states, if present, would emerge above the bulk continuum, as illustrated in the (b) panel of Fig. 2. The analysis of the Appendix B shows that no such localized edge states exist in the hyperbolic case, irrespective of the degree of anisotropy.

III.3 Validity of the approximation

Let us estimate the quality of the Fetter approximation with respect to the boundary condition $\epsilon_{2}E_{z}(z=0^{+})-\epsilon_{1}E_{z}(z=0^{-})=-\omega\epsilon 2\tau\kappa_{y}=\rho/\epsilon_{0}$. For $y>0$, we have $\rho=0$, and the boundary condition is trivially satisfied. But for $z\to\infty$, we would expect a decay of the vector potential. Also, there should not be any discontinuity of the vector potential at $y=0$ for large $z$. Our solution can thus only be valid close to the edge and this restriction should be due to the approximation of Eq. (14).

For $y<0$, the boundary condition would imply $2\tilde{\kappa}_{y}=2-\tilde{\kappa}_{y}^{2}$ which is only fulfilled for $\tilde{\kappa}_{y}=\sqrt{3}-1$. We will thus only expect a solution to well describe the edge mode for decay length with $2/3\sim\tilde{\kappa}_{y}\sim 4/5$, allowing a relative error of up to $10$%.

We can proceed with an analytical discussion by noting that the in-plane potential of the edge-mode penetrates exponentially in the bulk, a feature expected to hold true also for the exact solution. This fact can be understood as an imaginary component for the $y$ component of momentum, $q_{y}=-i\kappa_{y}$. Making use of the sheet response and the continuity equation, the frequency of the edge plasmon must fulfill

[TABLE]

Comparing this with the dispersion obtained from the Fetter approximation, Eq. (20), we see that it is only consistent up to order $\mathcal{O}(\tilde{\kappa}_{y}^{4})$. The approximated dispersion is thus always red-shifted and for $2/3\sim\tilde{\kappa}_{y}\sim 4/5$, we have relative errors of $4-12$%. This analysis can be verified in the case of an isotropic response matrix where the exact solution of Volkov and Mikhailov gives $\omega_{e}^{VM}(q)\approx 0.9\,\omega_{b}(q)$ and the Fetter approximation $\tilde{\kappa}_{y}=1/\sqrt{2}$ and $\omega_{e}^{F}(q)=\sqrt{\tfrac{2}{3}}\,\omega_{b}(q)$, i.e., $\omega_{e}^{F}/\omega_{e}^{VM}=0.9$.

All results hold for edge modes in an arbitrary system as long as they exist. For the special case of an isotropic system (with or without breaking time-reversal symmetry), we can further relate the exact edge mode energy to the exact bulk mode energy as $\omega_{e}(q)=\omega_{b}(\kappa_{z}^{-})$. This result will be useful for the next subsection where we discuss the ellipticity and polarisation plane of the edge mode, and present a better assessment of the approximations.

III.4 Polarization plane and ellipticity

We are now in the position to discuss the polarisation plane and ellipticity for the edge mode of ordinary materials. Now, we have a field vector in all three directions and we have to group together the field components that have the same phase.

From Eq. (18) and for $\tilde{\kappa}_{\pm}<1$, Maxwell’s equations dictate $\tilde{\kappa}_{z}^{\pm}=\sqrt{1-\tilde{\kappa}_{\pm}^{2}}$; for $\tilde{\kappa}_{\pm}>1$, this would result in an imaginary decay length and we could take the average of the two branches and set $\tilde{\kappa}_{z}^{\pm}=0$. We would then have the following result for the ellipticity in the Fetter approximation, given for the geometry of Fig. 1 by $\psi_{-}^{F}=E_{yz}/iE_{x}$, where $E_{yz}$ and $E_{x}$ are the $yz$-plane and $x$-axis field components:

[TABLE]

The field is thus always circularly polarized for $y<0$ and elliptically polarized for $y>0$. Notice that the mentioned inconsistencies in the field are to be expected because the approximation replaces the exact Maxwell kernel for an approximate one.

A fairer assessment is provided by regarding the whole Fetter procedure as a means to get an approximate surface charge density for the edge mode. This approximate surface density profile is taken to be

[TABLE]

with $\rho(y)$ given up to a global scale by

[TABLE]

for the edge mode of the isotropic case, for simplicity. One can then obtain the near field associated with such charge using the exact Coulomb Kernel. The results of this calculation for the charge density and ellipticity in the isotropic case are shown in the (a) and (b) panels of Figs. 3, 4 and 5. They provide a rather satisfactory description of the near field of the edge mode when compared with the exact results, later obtained with far more involved Wiener-Hopt techniques.

A peculiarity of the near-field is that right on top and below the conducting sheet with isotropic response where $\tilde{\kappa}_{y}=1/\sqrt{2}$, the field orientation in the $yz$ plane, $\theta=\arctan(E_{z}/E_{y})$, reaches the asymptotic values of $\pm 135^{\circ}$, see (a) panels of Figs. 4 and 5. Being the charge localized near the edge, a naive expectation would suggest a field oriented along the $-{\bm{e}}_{y}$ direction deep into the graphene sheet, as for the line charge of section II.2. This expectation is ill founded, as can be shown with the following exact reasoning. From the general discussion of the exponential decay of the edge mode in the bulk, we obtain the following relation for the asymptotic orientation of the near-field just above and below the graphene sheet:

[TABLE]

Let us now use the relation $\omega_{e}(q)=\omega_{b}(\kappa_{z})$ valid for an isotropic system with or without a magnetic field. This leads to

[TABLE]

We believe this relation to hold true for the exact solution in the instantaneous approximation, as it only relies on the asymptotic exponential penetration of the edge excitation in bulk graphene, see discussion above. The Fetter solution only approximately fulfills this relation, owing to the lack of exact Coulomb self-consistency. Nevertheless, our value of $\theta=\pm 135^{\circ}$ in the case of a isotropic response compares favorably with $\theta\pm\approx 125^{\circ}$, a result obtained from Eq. (29) for the exactVolkov and Mikhailov (1988) $\omega_{e}(q)\approx 0.9\,\omega_{b}(q)$.

III.5 Wiener-Hopf (exact) solution

In this section we obtain the edge mode density and associated near field for an ordinary, non-absorbing, isotropic material using the Wiener-Hopf technique developed in Ref. Volkov and Mikhailov (1988). The treatment provides the exact solution, at least within the certainly valid instantaneous approximation. The comparison with the previous approximate calculation confirms the latter to be a very good approach.

We assume a half-infinite sheet in the plane $z=0$ with $y<0$ as in Fig. 1. The surface charge density for a propagating edge mode can be written as

[TABLE]

where $\tilde{y}=|q|y$. The density $\rho(\tilde{y})$ can be obtained from its Fourier transform, $\tilde{\rho}(\tilde{k})$, which the Wiener-Hopf techniques of Ref. Volkov and Mikhailov (1988) allow to write as

[TABLE]

with $f(\tilde{k})$ given by

[TABLE]

where

[TABLE]

with $\tilde{\kappa}^{-1}_{z}=(\tfrac{\omega_{b}(q)}{\omega_{e}(q)})^{2}=1.2183$, under the assumption of a perfect Drude metal without dissipation. $f(\tilde{k})$ can be obtained almost analytically as

[TABLE]

where $\tilde{\kappa}_{y}=\sqrt{1-\tilde{\kappa}^{2}_{z}}$, and the function $g(\tilde{k})$ is given by

[TABLE]

It is interesting to notice that $g(\tilde{k})$ is a smooth function and, if one ignores its $\tilde{k}-$dependence as $g(\tilde{k})\approx 1$, a simple analytical expression for the density can be obtained, with correct functional behavior in the limits $y\to 0$ and $y\to-\infty$. It is given by

[TABLE]

where $\operatorname{erf}$ and $\Theta$ stand for error and Heaviside functions, respectively.

In any case, the results to follow use the (numerically) exact density. Writing the electric field as

[TABLE]

the numerically obtained, exact Fourier transform of Eq. (32) allows its calculation by means of the following explicit formulas

[TABLE]

[TABLE]

and

[TABLE]

with $\tilde{z}=|q|z$.

The exact charge profile and field absolute value $|\bm{E}|$ are plotted in the (c) and (d) panels of Fig. 3, and compare satisfactorily with the approximate results of the upper panels. The rapid decrease of the field away from the edge is the expected feature of any near field. Although the potential remains finite right at the edge,Volkov and Mikhailov (1988) the field component in the $yz$ plane diverges as $r^{-1/2}$ when $r\to 0$, where $r$ is the distance to the edge.

We are particularly interested in the spin-momentum locking properties of the near-field that enables non-reciprocity, i.e., directional selectivity in the excitation of plasmons by means of local sources (dipoles), with appropriately chosen circular (elliptical) polarization.

In the (c) and (d) panels of Fig. 4, we plot the orientation on the near-field in the $yz$ plane, which together with the $x$ axis (out of plane in the figure) determines the plane where the field rotates. This defines the optimal plane for directional excitation. The ellipticity or degree of circular polarization, $\psi=E_{yz}/iE_{x}$ is plotted on the right panel of the same figure. In general, the ellipticity tends to one (circular polarization) with increasing distance from the edge. On the contrary, the $yz$ component dominates very close to the edge, where it diverges. Again, the comparison between exact and approximate treatments is satisfactory.

Both features of the exact solution are summarized in the (c) and (d) panels of Fig. 5, and again compared with the approximate treatment plotted in the (a) and (b) panels of the same figure. In Fig. 5 (c), we show the rotation of the polarisation plane from $0^{0}$ for $y\to+\infty$ to $\pm 125^{0}$ for $y\to-\infty$ , i.e., the characteristic circling of the polarisation plane around the edge, but with the asymptotic angle $\theta=\arctan(-\tilde{\kappa}_{z}/\tilde{\kappa}_{y})\approx 125^{0}$ inside a graphene plane, not far from $135^{0}$, the result of the Fetter approximation of Fig.5 (a). On the right, the ellipticity as function of $y$ is shown for various distances $z$ with respect to the graphene plane. Far away from the edge, the ellipticity tends to one, a value that coincides with that of an infinite line charge within the instantaneous approximation.

IV Simulations

For the near-field of the edge-plasmon and, in appropriate units, the field along the edge ($x$ axis) is always proportional to the phase $i$, whereas the component in the $yz$ plane is real, a quadrature that leads to elliptical polarization. Therefore, an exciting dipole with an elliptical polarization chosen to match the local polarization of an edge-mode with a given $q$, will preferentially excite that mode and, ideally, not the opposite $-q$ mode, in spite of having the same frequency, leading to directional selectionLodahl et al. (2017).

Here, we will show the results of exact electromagnetic simulations using COMSOL Multiphysics RF Module. We consider a half disk of graphene of radius 800 nm (see Fig. 6). The graphene’s relaxation time and chemical potential are 100 fs and 0.25 eV, respectively. The wavenumber, $q$, of the edge mode, calculated using the dispersion equation Eq. (20), is 64.49 $\mu$m*-1* at the frequency $\nu=40$ THz. We assume that the edge modes are excited by an electric dipole $\mathbf{p}=-i\mathbf{e}_{x}+p_{0}(\mathbf{e}_{y}\cos\phi-\mathbf{e}_{z}\sin\phi)$ A$\cdot$m placed at the point $x=0$ nm and $z=5$ nm ($zq=0.32$). The $y$ coordinate of the dipole is varied as indicated in Fig. 6 (c, d, e, f).

It can be seen from Figs. 6 (a,b) that we can excite either left or right propagating edge modes by choosing the dipole polarization angle, $\phi$. The angle that allows for the most efficient suppression of the right (left) propagating edge mode increases (decreases) with the increase of the dipole coordinate $y$ (see Fig. 6 (c)). This angle does in general not coincide with the angle which provides the most efficient transfer of the dipole energy to the edge modes (compare Fig. 6 (c) and Fig. 6 (d). Particularly, for the dipole located at $y=0$, the edge modes carry the highest energy for the dipole polarized at angle $\phi=90^{\circ}$. On the other hand suppression of the right propagating edge mode requires the dipole polarization to be $\phi=52^{\circ}$.

Finally, we want to extract ellipticity and polarization plane of the edge mode from the numerical simulations. In order to do this, we take into account that for the right propagating edge plasmon the mode field structure is $\mathbf{e}_{m}=-i\mathbf{e}_{x}+E_{0}(-\mathbf{e}_{y}\cos\theta+\mathbf{e}_{z}\sin\theta)$, where $E_{0}$ is an ellipticity and $\theta$ is the polarization plane angle. The condition for the mode suppression is $\mathbf{p}\cdot\mathbf{e}_{m}^{*}=0$, which leads to $p_{0}E_{0}(\cos\phi\cos\theta+\sin\phi\sin\theta)=1$. By calculating numerically the dipole angle $\phi$ that leads to suppression of the right propagating mode for two different dipole momenta, $p_{0}=p_{1,2}$, we can obtain the polarization plane angle and the ellipticity as $\tan\theta=(p_{2}\cos\phi_{2}-p_{1}\cos\phi_{1})/(p_{1}\sin\phi_{1}-p_{2}\sin\phi_{2})$, and $E_{0}=1/(p_{1}(\cos\phi_{1}\cos\theta+\sin\phi_{1}\sin\theta))$. The numerical results are presented in Figs. 6 (e, f) for $p_{1}=1$ and $p_{2}=1.4$. It can be seen that the Wiener-Hopf exact solution provides a very good description of the edge plasmon mode structure in the 2D electron gas as expected, since retardation effects can be neglected.

V Discussion and Conclusion

In this paper, we have investigated the near field of edge modes hosted by general two-dimensional materials together with its spin-momentum characteristics. By an analytical and numerical approach, we were able to identify sweet spots to excite edge-modes uniquely in one or two directions.

We have also extended previous approximate treatments to cover anisotropic and hyperbolic systems, studying the conditions for the existence of edge modes. Furthermore, a detailed assessment of the approximate approach has been carried out, showing it to provide a very good description of the near field of the edge state when compared with exact results.

Our analysis is also backed up by exact electromagnetic simulations. By this, we demonstrate the importance of choosing the optimal polarisation plane and position of exciting dipoles. We hope that our results will help to design more efficient plasmonic circuitries. Future studies can further optimise the excitation of unidirectional edge vs. bulk plasmons.

Acknowledgements

This work has been supported by Spain’s MINECO under Grant No. FIS2017-82260-P and FIS2015-64886-C5-5-P. TL acknowledges support by the National Science Foundation NSF/EFRI grant (#EFRI-1741660).

Appendix A: Derivation of the generalized Fetter solution

The general incoming vector potential inside the half-planes $z<0\;(m=1)$ and $z>0\;(m=2)$ of the dielectric with $\mu_{m}$ and $\epsilon_{m}$ can be written in terms of its Fourier components in $x$-direction as

[TABLE]

where ${\bm{e}}_{n}$ are unit vectors. Maxwell’s equations yield the following wave equation:

[TABLE]

We will now introduce the Fourier transform in $y$-direction

[TABLE]

With this, we can solve Eq. (43), i.e., $a_{m}^{n}(q_{x},q_{y},z)=a_{m}^{n}(q_{x},q_{y})e^{-q_{m}^{\prime}|z|}$ where we defined the perpendicular momentum $q_{m}^{\prime}=\sqrt{q_{x}^{2}+q_{y}^{2}-\mu_{m}\epsilon_{m}(\omega/c)^{2}}$ and the velocity of light in vacuum $c^{2}=1/(\epsilon_{0}\mu_{0})$. The physical fields are obtained as ${\bf E}_{m}=i\omega{\bf A}_{m}$ and $\mu_{m}\mu_{0}{\bf H}_{m}=\nabla\times{\bf A}_{m}$. The sheet currents are also functions of $y,z$ and shall be expanded in their Fourier series with

[TABLE]

Due to the transverse nature of the vector potential in the two half-planes, $\nabla\cdot{\bf A}_{m}=0$, we have the following relations:

[TABLE]

The boundary condition ${\bf n}\times({\bf E}_{2}-{\bf E}_{1})=0$ (with ${\bf n}=(0,0,1)$ the normal to the interface) guarantees the continuity of the parallel components of the vector potential, i.e., $a_{1}^{n}=a_{2}^{n}$ for $n=x,y$. The second boundary condition ${\bf n}\times({\bf H}_{2}-{\bf H}_{1})={\bf j}$ yields the two equations with $n=x,y$:

[TABLE]

With the above equations and ($\mu_{1}=\mu_{2}=1$), this reduces to

[TABLE]

with $n=x,y$ and summation over double indices is implied.

We can then relate the vector potential to the current densities in the following way:

[TABLE]

with $\Omega^{2}=\frac{\epsilon_{1}q_{2}^{\prime}+\epsilon_{2}q_{1}^{\prime}}{q_{1}^{\prime}+q_{2}^{\prime}}\frac{\omega^{2}}{c^{2}}$. Notice that in order to get from Eq. (49) to Eq. (52), one needs to keep the retardation effects because otherwise the determinant of the matrix would vanish.

Let us consider the generating function

[TABLE]

which is related to the modified Bessel function of second kind of $z=0$ and in the non-retarded limit $c\to\infty$. We then have

[TABLE]

with the effective dielectric constant $\epsilon=(\epsilon_{1}+\epsilon_{2})/2$.

Fourier transforming the current in the above equation, we can thus finally write the relation between the vector potential and the current in real space:

[TABLE]

For a finite current only in the half-space $y<0$, we arrive to Eq. (12) of the main text where we have slightly changed the notation with $q_{x}\to\tau q$.

Up to now, all operations have been exact and approximations are necessary to proceed analytically. First, we will neglect retardation effects and second and more crucially, we will approximate the generated function by

[TABLE]

since the first three moments of the two functions are identical (after neglecting retardation effects). $L_{0}(y)$ is now the Green’s function with respect to the operator

[TABLE]

i.e., we have $\mathcal{O}_{0}(y)L_{0}(y)=\delta(y)$ and thus

[TABLE]

The current is generally related to the vector potential as $j^{n}=-\chi_{n,{n^{\prime}}}a^{{n^{\prime}}}$. From Eq. (74), we see that the $y$-component of the gauge potential must be discontinuous. On the other hand, only the first derivative of the $x$-component of the gauge potential is discontinuous. We can thus make the following ansatz with $\kappa_{\pm}>0$:

[TABLE]

From Eq. (74), we then get

[TABLE]

Both expressions yield the same condition, i.e.,

[TABLE]

with the dimensionless quantities $\tau=q_{x}/|q_{x}|$, $\tilde{\kappa}_{\pm}=\kappa_{\pm}/|q_{x}|$ and $\tilde{\chi}_{n,{n^{\prime}}}=\frac{|q_{x}|}{\epsilon_{0}\epsilon\omega^{2}}\chi_{n,{n^{\prime}}}$.

For $y>0$, the resulting matrix has determinant zero and each component has to be zero. This yields $\kappa_{+}=\sqrt{2}|q_{x}|$. Therefore, we can write Eq. 79 as

[TABLE]

which is Eq. 19 of the main text. For $y<0$, we obtain the following condition:

[TABLE]

which coincides with Eq. 20 of the main text.

Appendix B: Discussion of the generalized Fetter solution

Eqs. (80) and (81) provide the conditions for the existence of edge modes and can be combined as

[TABLE]

The evanescence of a putative edge mode requires $\operatorname{Re}(\tilde{\kappa}_{y})>0$, which for a non-absorbing system with TRS, where all entries are real, imply $1-\tilde{\chi}_{yy}<0$. Therefore, $\chi_{yy}$ should be at least positive, justifying Eq. (22). With TRS where $\tilde{\chi}_{xy}=\tilde{\chi}_{yx}$, one can rewrite Eq. (82) as

[TABLE]

The right hand side (RHS) of Eq. (83) vanishes at the edge of the continuum region (see Fig. 2) and, remembering the definition of tilded quantities: $\tilde{\chi}=\frac{q}{\epsilon_{0}\epsilon\omega^{2}}\chi$, a path in $q$ at constant $\omega$ in the (a) panel of Fig. 2 that goes from the continuum edge to $q\to\infty$ drives the RHS of Eq. (83) from [math] to $+\infty$ for a system with $\det(\chi)>0$. This guaranties a solution of Eq. 83 and, therefore, the existence of an edge mode for any ordinary, that is, non-hyperbolic, system with TRS. In the absence of TRS, where $\chi_{xy}=\chi^{\prime}_{xy}+i\chi_{c}$, and $\chi_{yx}=\chi^{\prime}_{yx}-i\chi_{c}$, with the real part is still symmetric $\chi^{\prime}_{xy}=\chi^{\prime}_{yx}$, the previous solution of Eq. (82) simply splits according to the sign of $\tau$, the propagation directionVolkov and Mikhailov (1988). In fact, whenever an edge mode exists with or without TRS, its dispersion relation satisfies

[TABLE]

a result that follows simply from Eq. (82) for the general Drude matrix in non-absorbing systems.

The previous argument for the existence of edge modes fails for a hyperbolic system, at least under TRS: the RHS of Eq. (83) cannot be positive and, at the same time, fulfill the requirements $1-\tilde{\chi}_{yy}<0$ and $\det(\tilde{\chi})<0$. Therefore, no edge mode detaches from the continuum in hyperbolic systems with TRS.

The analysis for edge modes in the Fetter treatment can be extended to the case where the complementary half-plane, $y>0$, is also a material medium with Drude matrix $\chi^{(2)}$, with $\chi^{(1)}$ representing now the response of the original $y<0$ half-plane. For the case where both half-planes respect TRS, $\chi^{(1,2)}_{xy}=\chi^{(1,2)}_{yx}$, an edge mode is guaranteed to exist in the following two cases:

• i)

medium $(1)$ is elliptic and metallic $(\det(\chi^{(1)})>0,\;\operatorname{tr}(\chi^{(1)})>0)$, and medium $(2)$ is elliptic and capacitive Bisharat and Sievenpiper (2017) $(\det(\chi^{(2)})>0,\;\operatorname{tr}(\chi^{(2)})<0)$.

• ii)

medium $(1)$ is elliptic and metallic, and medium $(2)$ is hyperbolic $(\det(\chi^{(2)})<0)$, provided both media share a common region of evanescence in the $(q,\omega)$ plane.

Whenever such an edge mode arises, its dispersion within TRS is given by the following rather simple generalization of Eq. 84

[TABLE]

with $\Delta(\det)=\det(\chi^{(2)})-\det(\chi^{(1)})$, $\Delta\chi_{xx}=\chi^{(2)}_{xx}-\chi^{(1)}_{xx}$, and $\Delta\chi_{yy}=\chi^{(2)}_{yy}-\chi^{(1)}_{yy}$.