Lasing condition for trapped modes in subwavelength-wired PT-symmetric resonators

Mauro Cuevas; Mauro Cuevas; Mojtaba Karimi Habil; Carlos J. Zapata-Rodríguez

doi:10.1364/OE.414534

1. Introduction

It is known that, under certain conditions, a plasmonic system can support trapped electromagnetics modes, which are electromagnetic non-radiating oscillations that stay localized inside the structure [1]. These plasmon modes are characterized by an ultra-high quality factor provided that the structural symmetry remains unbroken. In fact, the narrow resonances characterizing these mode excitations rely on the breaking degree of the symmetry which allows the trapped mode to couples with free photons [2]. In this sense, plasmon trapped modes can be considered as a kind of symmetry protected states. Due to Ohmic loss in the plasmonic material, the eigenfrequencies of the trapped states are complex valued. However, by introducing gain material elements into the structure it is possible to compensate this material loss and, as a consequence, to reach the lasing condition for which the eigenfrequency associated to an eigenmode is real valued [3].

Because of their fundamental properties as well as their potential applications, the study of hybrid systems composed of gain media and plasmonic materials, which are lossy, is a topic of continuous increasing interest [4–8]. In areas such as condensed matter and surface optics, the amplification of eigenmodes by stimulated emission of radiation has played a key role in the interpretation of a wide variety of experiments, the understanding of various fundamental properties of solids and the engineering of nanolaser devices [9–11]. In particular, new phenomena associated with parity-time (PT) symmetry have been observed in optical systems with balanced loss and gain [12,13]. These optical systems belong to a large family of non-Hermitian systems, which can have a real spectrum provided that the system be invariant under combined operations of parity (P) and time-reversal (T) symmetry [14,15].

Possibilities have been widened towards PT symmetric structures that incorporate graphene as plasmonic material. Doped graphene allows the propagation of surface plasmons with low Ohmic losses, i.e., with high quality factor, from terahertz to near-infrared range [16]. Moreover, the plasmon resonance spectrum, and consequently the optical responses of the system, can be tuned by varying the doped level on graphene. In this context, several works have focused on the influence that long-living and tunable surface plasmons have on graphene based strucures [17,18] and, in particular, on those that present PT-symmetry [19–22].

The present work contains all the ingredients above entered and focusing on the study of a sub-wavelength graphene plasmonic system with balanced volume losses. In particular, our system is composed of two parallel dielectric cylinders, one is with loss and the other with the same level of gain, both wrapped with graphene sheets. Here, the sole effect of graphene is to provide LSP resonances in sub-wavelength wires. Although other more complex structures such as an array of many or infinite cylinders can be designed to present PT-symmetry, our motivation is based in the simplicity of the dimer structure that we have chosen, since it allows us to understand the underlying physics behind the PT-symmetry.

Other dielectric and plasmonic structures (without graphene) with balanced gain and loss have been proposed and studied in the literature. For instance, exceptional points of resonant eigenmodes (whispering gallery eigenmodes) were examined in a system consisting of a finite number of parallel dielectric cylinders in [23]. By exploiting the dynamic characteristics of whispering gallery eigenmodes near the exceptional point, a route to single particle sensing was analyzed for two coupled resonators system [24]. Furthermore, full loss compensation with relative low gain level can also be reached by using the eigenmodes supported by two parallel wires made of plasmonic materials in the visible spectrum [25]. Recently, optical forces have been explored for two interacting parallel plane waveguides, which results from the coupling between the PT-eigenmodes propagating along the symmetry axis of the system [26]. Finally, an overview of recent advancements in the field of non-Hermitian as well as PT-symmetry optical systems can be found in [27], together with a systematic discussion of the key results that are useful to understand exceptional point dynamics.

We use a rigorous formalism based on Mie theory to calculate the eigenmodes dependence with the gain–loss parameter for both trapped and radiating eigenmodes. To do this, we solve the boundary value problem with the corresponding boundary conditions without an incident wave (homogeneous problem). The procedure requires the analytic continuation of the eigenvalue, the frequency in our case, in the complex plane. Analytic continuation is inevitable, even in the case of media with no intrinsic losses, since the open nature of our resonator generates non-null imaginary parts due to radiation losses.

In addition, as it is well established nowadays, the modal lasing analysis in an open plasmonic resonator, as our PT-symmetric structure, can be carried out by applying the lasing eigenvalue problem (see [28] and references therein), in which the modal eigenfrequencies are assumed to be real valued at the lasing condition (or lasing threshold). Following this concept, we solve the homogeneous problem to find the real valued eigenfrequencies, gain–loss parameters and chemical potentials for each modal lasing condition. Moreover, we study the eigenmode influence on the optical response of the system when it is excited by a plane wave near the lasing condition.

By using the quasistatic approximation valid in the long wavelength limit, we show that despite the graphene ohmic losses, a fact that gives rise to a complex spectrum, the structure exhibits a set of properties in common with a PT-symmetric system. For instance, two eigenmode branches coalesce at an exceptional point and, for the gain–loss parameter above certain threshold, these both branches are repelled in the direction of the imaginary part of frequency, making that one of these branches achieves the lasing condition. In addition, we demonstrate that branches eigenmode, which are trapped modes in case of null gain–loss parameter, transform into radiating modes that can be excited by plane wave incidence when the gain–loss parameter is increased. In previous works [29,30] it was reported mechanisms for controlling transitions from trapped to radiating eigenmodes based on producing a difference between the modulus of individual dipole moments on each particle forming the dimer, i.e., by producing a weakly asymmetry with respect to the center of the dimer modifying the dipole moment amplitudes. In addition, in [31] it was experimentally demonstrated that the modulus and phase of the coupling coefficient between the magnetic mode and the electric dipolar mode on a composite metamaterial can be efficiently manipulated for practical applications on radiation control. In particular, we show that the increment of the gain–loss parameter leads to a change in the phase difference between these individual dipole moments, maintaining their modulus values constant.

This paper is organized as follows. In section 2 we present a brief description of the rigorous method used in this work to calculate the scattering of a dimmer composed of two graphene wires. From this method, using the quasistatic approximation, we deduce analytical expressions for eigenfrequencies and eigenvectors as a function of geometrical and constitutive parameters that explain the main features calculated with the rigorous method. In section 3 we present results of two parallel dielectric cylinders tightly coated with a graphene layer, one of them with small inner losses and the other one with the same level of gain. Concluding remarks are provided in section 4. The Gaussian system of units is used and an $\mbox {exp}(-i\, \omega \, t)$ time-dependence is implicit throughout the paper, where $\omega$ is the angular frequency, $t$ is the time coordinate, and $i=\sqrt {-1}$. The symbols Re and Im are used for denoting the real and imaginary parts of a complex quantity, respectively.

2. Theory

2.1 Rigorous description of the fields and scattering efficiencies

We consider a cluster consisting of two parallel and non-overlapping cylindrical dielectric wires, one with gain, $\varepsilon _a=\varepsilon _1-i\varepsilon _i$ ($\varepsilon _i>0$), and the other with equal loss, $\varepsilon _b=\varepsilon _1+i\varepsilon _i$, as shown in Fig. 1. Both wires have the same radius $R_a=R_b=R$ and are wrapped with a graphene sheet. The system is embedded in a lossless and non-magnetic dielectric denoted as medium $v$ with permittivity $\varepsilon _v$. In this case, a PT symmetry around the central axis, denoted by $O$, is fulfilled. We assume that the radius $R$ is sufficiently large to describe the optical properties of the wires as characterized by the same local surface conductivity as planar graphene (see appendix A). Even though nonlocal effects can be observed in our optical system for frequencies lower than characteristic resonance frequencies, this is not interesting for our purposes (see appendix B). We denote by $r_j(\mathbf {r}),\,\phi _j(\mathbf {r})$ ($j=a,\,b$) the polar coordinates of a point at position $\mathbf {r}$ with respect to the local origin $O_j$. A plane wave radiation impinges on the wires with an angle of incidence $\phi _{inc}$ with respect to the $y$ axis. Although some enhanced optical effects related to invisibility modes are observed for $s$ polarization (electric field along the $z$ axis) [32], this work focus on $p$ polarization (magnetic fields along the $z$ axis) for which the electric field in the graphene coating induces electric currents directed along the azimuthal direction $\phi _j(\mathbf {r})$ and LSPs exist in the graphene circular cylinder. In this way, the incident magnetic field (along the $z$ axis) can be written in a system linked to the $j$-cylinder as [30],

(1)$$H_{inc}(\mathbf{r}) = e^{i k_v r^j \sin(\phi_{inc}-\phi^j)} \sum_{m={-}\infty}^{+\infty} ({-}1)^m J_m(k_v r_j(\mathbf{r})) e^{i m \phi_l(\mathbf{r})} e^{{-}i m \phi_{inc}},$$

where $r^j,\,\phi ^j$ are the polar coordinates of the $j$-cylinder, $k_v=\sqrt {\varepsilon _v} \frac {\omega }{c}$, is the modulus of the photon wave vector in medium $v$, $\omega$ is the angular frequency, $c$ is the vacuum speed of light and $J_m(x)$ is the $n$th Bessel function. The scattered magnetic field in medium $v$ ($r_j(\mathbf {r})>R$) can be written as a superposition of the field scattered by each of the cylinders,

(2)$$H^{(v)}_s(\mathbf{r}) = \sum_{j=a,b} \sum_{m={-}\infty}^{+\infty} b_{j\,m} H_m(k_v r_j(\mathbf{r})) e^{i m \phi_j(\mathbf{r})},$$

where $H_m(x)$ is the $n$th Hankel functions of the first kind. Note that the $j$th term of the summation corresponds to the field scattered by the $j$th cylinder linked to the local system with origin $O_j$. In the region inside the cylinders, $r_j(\mathbf {r})<R$, the transmitted field is written as

(3)$$H^{(j)}(\mathbf{r}) = \sum_{m={-}\infty}^{+\infty} a_{j\,m} J_m(k_j r_j(\mathbf{r})) e^{i m \phi_j(\mathbf{r})},$$

where $j=a,\,b$. To find the unknown complex amplitudes of the reflected $b_{j m}$ and transmitted $a_{j m}$ fields (2) and (3), we use the usual boundary conditions and the addition theorem for Bessel and Hankel functions [33]. This theorem allow us to write one of the terms in (2), associated to the scattered field of one of the cylinders (for example, the $j=b$ cylinder) in the other local coordinates (the $j=a$ cylinder). In this way, the scattered field (2) will be represented in the form of expansions in Hankel functions written in the local coordinate $j=a$. By replacing this expression and Eq. (3) with $j=a$ into the boundary conditions along the surface $r_a=R$ of the $j=a$ cylinder, one obtain a set of $2$ equations for the $2 \times 2$ unknown amplitudes. Similarly, we can write the scattered field (2) in the local coordinate $j=b$ and use the boundary conditions on the surface of the $j=b$ cylinder to obtain other set of 2 equations for the $2 \times 2$ unknown amplitudes. However, we closely follow a variant of the method, developed in [34], that allows to reduce to half the dimension of the system of equations. The detailed of this implementation has been given in [30], and leads to the following system of equations for the amplitudes $b_{j m}$

(4)$$\left[ \begin{array}{llll} \,\,\,\,\,\,\,\,\,\,\,\,\overline{\overline{\mathbf{I}}} & -\overline{\overline{\mathbf{S}}}_a\cdot\overline{\overline{\mathbf{T}}}_{a\,b} \\ -\overline{\overline{\mathbf{S}}}_b\cdot \overline{\overline{\mathbf{T}}}_{b\,a} & \,\,\,\,\,\,\,\,\,\,\,\,\overline{\overline{\mathbf{I}}} \\ \end{array} \right] \left( \begin{array}{lll} \overline{\mathbf{b}}_a\\ \overline{\mathbf{b}}_b\\ \end{array} \right)= \left( \begin{array}{lll} \overline{\overline{\mathbf{S}}}_a\cdot \overline{Q}_a\\ \overline{\overline{\mathbf{S}}}_b\cdot \overline{Q}_b\\ \end{array} \right),$$

where $\overline {\mathbf {b}}_l$ and $\overline {\mathbf {Q}}_l$ are vectors whose coordinates are the elements $b_{l m}$ and

(5)$$Q_{j m}= e^{i k_v r^j \sin(\phi_{inc}-\phi^j)} ({-}1)^m e^{{-}i m \phi_{inc}},$$

respectively, $\overline {\overline {\mathbf {T}}}_{lj}$ is the matrix with elements

(6)$$T_{l j n m}= H_{n-m}(k_v r_l^j) e^{i (n-m) \phi_l^j},$$

and $\mathbf {S_j}$ is the matrix with elements $S_{j m n}=s_{j m} \,\delta _{m n}$, where $s_{j m}$ are the elements of the scattering matrix associated to the $j$th cylinder [35]

(7)$$s_{j\,m} = \frac{[k_j \varepsilon_v J_m(y) J^{'}_m(x)-k_v \varepsilon_j J_m^{'}(y) J_m(x)-\frac{4\pi \sigma}{c k_0} i k_j k_v J^{'}_m(x) J^{'}_m(y)]}{k_v\varepsilon_j J_m(x) H^{'}_m(y)-k_j \varepsilon_v J^{'}(x)H_m(y)+\frac{4\pi \sigma}{c k_0} i k_j k_v J^{'}_m(x) H^{'}_m(y)},$$

where $x=k_j R$, $y=k_v R$ and the prime denotes the derivative with respect to the argument.

Fig. 1. Schematic illustration of the system composed by $2$ dielectric cylinders wrapped with graphene sheets.

Download Full Size | PDF

Knowing the total electromagnetic field allows us to calculate the scattering cross sections. The time-averaged scattered power is calculated from the integral of the radial component of the complex Poynting vector flux through an imaginary cylinder of length $L$ of radius $\rho _0$ which envelops the graphene wire system (see Fig. 1),

(8)$$P_s=\rho_0 L \int_0^{2 \pi} \left\langle \mathbf{S}\right\rangle \cdot \hat{r}\, d\phi,$$

(9)$$\left\langle \mathbf{S}(\rho_0,\phi)\right\rangle \cdot \hat{r}=\frac{c}{8 \pi} \textrm{Re} \left\{ E_{s, \phi} H_s^{*}\right\}= \frac{c}{8 \pi} \textrm{Re} \left\{\frac{1}{i k_0 \varepsilon_v} \frac{\partial \, H_s}{\partial r} H_s^{*}\right\}.$$

It is convenient to express each of the terms in Eq. (2) in a same coordinate $O$ [30]. Substituting the obtained expression into Eq. (9) we obtain

(10)$$\left\langle \mathbf{S}(\rho_0,\phi)\right\rangle \cdot \hat{r}= \frac{c}{8 \pi} \textrm{Re} \left\{\frac{1}{i \sqrt{\varepsilon_v}} \sum_{m,n} B_m B_n^{*} H'_m(k_v \rho_0) H^{*}_m(k_v \rho_0) e^{i(m-n)\phi}\right\},$$

where

(11)$$B_m=\sum_{j=a,b} \sum_{q={-}\infty}^{+\infty} b_{j m} J_{m-q}(k_v r^j) e^{i(q-m) \phi^j}.$$

Inserting Eq. (10) into (8) and taking into account the wronskian $W\left \{J_m(x),Y_m(x)\right \}=2/(\pi x)$ ($x=k_v \rho _0$), after some algebraic manipulations, we obtain the power scattered by the cylinders,

(12)$$P_s=\frac{c^2 L}{2\pi \omega \varepsilon_v} \sum_{m={-}\infty}^{+\infty} |B_m|^2.$$

The scattering cross section is defined as the ratio between the total power scattered by the cylinders, given by Eq. (12), and the incident power $P_{inc}$ intersected by the area of all the cylinders.

It is known that the scattering cross section and the near to the cylinders field are strongly affected by complex singularities in the field amplitudes $b_{j m}$. Singularities occur at complex locations and they represent the frequency of the eigenmodes supported by the cylinders system. Complex frequencies of these modes are obtained by solving the homogeneous problem, i.e., by imposing the vectors $\overline {\mathbf {Q}}_a$ and $\overline {\mathbf {Q}}_b$ in Eq. (4) be zero. Then, a condition to determine eigenfrequencies is to require the determinant of the matrix in Eq. (4) to be zero,

(13)$$\mbox{det}\left[ \begin{array}{llll} \,\,\,\,\,\,\,\,\,\,\,\,\overline{\overline{\mathbf{I}}} & -\overline{\overline{\mathbf{S}}}_a\cdot\overline{\overline{\mathbf{T}}}_{a\,b}\\ -\overline{\overline{\mathbf{S}}}_b\cdot\overline{\overline{\mathbf{T}}}_{b\,a} & \,\,\,\,\,\,\,\,\,\,\,\,\overline{\overline{\mathbf{I}}} \\ \end{array} \right]=0.$$

This condition corresponds to the full retarded dispersion equation (FR) for eigenmodes and it determines the complex frequencies $\omega =\omega _R+i\omega _I$ ($\omega _I<0$) in terms of all the geometrical and constitutive parameters of the system.

2.2 Quasistatic approximation: a simple model based on two coupled electric dipoles

Although the rigorous treatment represented by Eq. (13) gives us all the kinematics and dynamics characteristics of the PT-symmetric eigenmodes, the method lacks of analytical expressions that explain the main dependencies with both geometrical and constitutive parameters. For the purpose of showing this dependencies, by applying the quasistatic method, here, we reduce the full treatment provided by the homogeneous part of Eq. (4) to a simple $2\times 2$ matrix description as follows.

Assuming that the radius $R$ of cylinders is much smaller than the wavelength $\lambda =2\pi c/\omega$, the problem can be treated using the dipole approximation where the dimer eigenfunctions calculated with (4) are, as a good approximation, a superposition of single plasmons with angular momentum $m=\pm 1$ linked to each cylinder. In this way, only four coordinates, $b_{j\pm 1}$ ($j=a,\,b$), define the amplitudes vector.

We first consider the case in which the induced dipole moments are along the $\pm x$ directions, i.e., the local magnetic field associated to each of the cylinders has a dependence $\approx \sin \phi _j$ ($j=a,\,b$). As a consequence, the amplitudes $b_{j 1}=b_{j -1}$ ($j=a,\,b$). Here, the subscript $j -1$ stand for cylinder $j$ and angular momentum $m=-1$. In this way, the matrix equation for the modal amplitudes reduces to a 2$\times$2 matrix system for amplitudes $b_{a 1}$ and $b_{b 1}$,

(14)$$\left[ \begin{array}{llll} 1 & - 2\,s_a\,H_1(z)/z \\ - 2\,s_b\,H_1(z)/z & 1 \\ \end{array} \right] \left( \begin{array}{lll} b_{a 1}\\ b_{b 1}\\ \end{array} \right)= \left( \begin{array}{lll} 0\\ 0\\ \end{array} \right),$$

where we have used $2 H_1(z)/z=H_0(z)+H_2(z)$, $z=k_v R_{ab}$ ($R_{ab}$ is the distance between cylinder centers). Using the small argument asymptotic expansions for Bessel and Hankel functions [33], it follows that $s_{j}$ can be written as,

(15)$$s_j=\frac{\pi\, k_v^2}{4} i \alpha_j,$$

$j=a,\,b$, where

(16)$$\alpha_j= R^2 \frac{\varepsilon_1\pm i\varepsilon_i-\varepsilon_v+g(\omega)}{\varepsilon_1\pm i\varepsilon_i+\varepsilon_v+g(\omega)},$$

are the dipolar polarizability of the cylinders $j=a$ or $j=b$, respectively, $g(\omega )=-\frac {\omega _g^2}{\omega ^2+i \omega \gamma _g}$, $\omega _g^2=\frac {4 e^2 \mu }{ \hbar ^2 R}$ is the effective plasma frequency for the dipolar mode [36]. Note that the value of $\alpha _a$ differs from $\alpha _b$ in the sign of the gain–loss parameter $\varepsilon _i$, signs $+$ and $-$ correspond to $j=a$ and $j=b$ respectively. Near the resonance frequency, the polarizability $\alpha _j$ ($j=a,\,b$) can be written as

(17)$$\alpha_j= R^2 \frac{\varepsilon_1\pm i\varepsilon_i-\varepsilon_v+g(\omega_j)}{g'(\omega_j) (\omega-\omega_j)}=\frac{A_j}{\omega-\omega_j},$$

where $\omega _j=\omega _{jR}+i\omega _{jI}$ ($\omega _{jI}<0$) is the complex pole of $\alpha _j$, i.e., the eigenfrequency of the single graphene cylinder $j$, $g'(\omega )$ is the derivative of $g(\omega )$ and

(18)$$A_j=R^2\frac{\varepsilon_1\pm i\varepsilon_i-\varepsilon_v+g(\omega_a)}{g'(\omega_a)}.$$

The single cylinder eigenfrequencies for cylinders $a$ and $b$ are written as [36]

(19)$$\omega_{j}= \frac{\omega_g}{\sqrt{\varepsilon_j+\varepsilon_v}}-i\frac{\gamma_g}{2}\approx \omega_0+i\omega_{jI},$$

where

(20)$$\omega_{0}= \frac{\omega_g}{\sqrt{\varepsilon_1+\varepsilon_v}},$$

(21)$$\omega_{jI}={-}(\Gamma_{sup}\pm\Gamma_{core}),$$

the signs $+$ and $-$ correspond to $j=a$ and $j=b$, respectively, and

(22)$$\begin{aligned} \Gamma_{sup}=\frac{\gamma_g}{2},\\ \Gamma_{core}=\frac{1}{2}\varepsilon_i \frac{\omega_g}{(\varepsilon_1+\varepsilon_v)^{3/2}}, \end{aligned}$$

are the damping rates corresponding to the ohmic loss in graphene covers and dielectric cores, respectively. We have used the fact that $\varepsilon _i<<\varepsilon _1+\varepsilon _v$ in the last equality in Eq. (19). From Eq. (20) we see that the real part of the eigenfrequencies do not depend on $\varepsilon _i$, a fact that also results by solving the fully retarded dispersion equation for a single graphene cylinder [37].

By taking into account the small argument $z=k_v R_{ab}$ in the Hankel functions, the off diagonal elements in the matrix in Eq. (14) are written as

(23)$$-2\frac{\,H(k_v R_{ab})}{k_v R_{ab}} s_j={-}2\frac{-2 i}{\pi (k_v R_{ab})^2} \frac{\pi k_v^2 i}{4} \alpha_j={-}\frac{\alpha_j}{R_{ab}^2}.$$

Therefore, the matrix Eq. (14) takes the form

(24)$$\left[ \begin{array}{llll} \omega-\omega_a & -A_a/R_{ab}^2 \\ - A_b/R_{ab}^2 & \omega-\omega_b \\ \end{array} \right] \left( \begin{array}{lll} b_{a 1}\\ b_{b 1}\\ \end{array} \right)= \left( \begin{array}{lll} 0\\ 0\\ \end{array} \right),$$

where $A_j$ $j=a,\,b$ are given by Eq. (18). Equation (24) gives us a simple description for the dimer dynamic. For large separations, $A_j/R^2_{ab}<<\omega _a$ (or $R/R_{ab}<<1$), the matrix (24) is diagonal and thus the eigenfrequencies correspond to that of each individual graphene wires composing the dimer, i.e., the plasmonic wires do not interact between them. Taking into account the PT parameters, the real parts of both frequencies are the same whereas their imaginary parts differ in $2\Gamma _{core}$. For small enough values of $R_{ab}$, extra diagonal terms take appreciable values and a splitting between the real parts of eigenfrequencies occurs. From Eq. (C8) we can see that this splitting is proportional to $R^2/R^2_{ab}\omega _0$. For system (24) to have a non-trivial solution, its determinant must be equal to zero, a condition which can be written as

(25)$$(\omega-\omega_a)(\omega-\omega_b)=\frac{A_a A_b}{R_{ab}^4}.$$

A detailed developed of the right hand side of this equation can be seen in appendix C. By replacing expression (C8) into Eq. (25) and solving for $\omega$ eigenfrequencies,

(26)$$\begin{aligned} \omega_{{\pm}}={-}i\frac{\gamma_g}{2}+\omega_0 \pm \sqrt{-\frac{\gamma^2}{4}+\omega_{aI}\omega_{bI}+\frac{A_a A_b}{R_{ab}^4}}=\\ \omega_0-i\frac{\gamma_g}{2} \pm \frac{\omega_0}{(\varepsilon_1+\varepsilon_v)}\sqrt{\frac{-\varepsilon_i^2}{4}+\frac{R^4}{R_{ab}^4} \varepsilon_v^2 }, \end{aligned}$$

where in the last equality we have used Eq. (21). This equation shows the dependence of the eigenfrequencies with the gain–loss optical parameter $\varepsilon _i$. For $\varepsilon _i=0$, the real parts of the eigenfrequencies split by $\Delta \,\omega =\omega _+ -\omega _- = 2 \frac {R^2}{R_{ab}^2} \frac {\omega _0}{(\varepsilon _1+\varepsilon _v)} \varepsilon _v$ while their imaginary parts remain degenerated at the value $-\gamma _g/2$. From Eq. (24), we see that the eigenvector associated to the upper branch verify $b_{a 1}=-b_{b 1}$ pointing out that the induced dipole moments on each of the cylinders move out of phase, i.e., the state corresponds to a trapped mode, while the eigenvector associated to the lower branch verify $b_{a 1}=b_{a 1}$ pointing out that this mode corresponds to a radiating mode that can be excited, for example, by an incident plane wave [30]. We use the notation $+x-x$ and $+x+x$ to refer to the upper and lower branches, respectively. Since the terms inside the root have opposed signs between them, the splitting between upper ($+x-x$) and lower ($+x+x$) branches monotonically decreases as $\varepsilon _i$ is increased until reaching the value

(27)$$\varepsilon_{ep}= 2 \frac{R^2}{R_{ab}^2} \varepsilon_v,$$

for which the eigenfrequencies coalesce at the exceptional point $\omega =-i\frac {\gamma _g}{2}+\frac {\omega _g}{\sqrt {\varepsilon _1+\varepsilon _v}}$. At the same time, the imaginary parts of both eigenfrequencies remain equal. Beyond the exceptional point and for the gain–loss parameter $\varepsilon _i>\varepsilon _{ep}$, the imaginary parts of the eigenfrequencies bifurcate, while their real parts remain degenerated. In Fig. 2 we illustrate in the complex plane the above described trajectories for the upper and lower branches (26) as parametric functions of the gain–loss parameter $\varepsilon _i$. We have taken the straight segment $z=0,\,\mbox {with}\,\textrm {Re}\,z>0$ as the cut line for the square root function in Eq. (26) so that $\sqrt {-w^2}=i w$ for $w$ real and positive.

Fig. 2. Trajectory in the complex plane of the eigenfrequencies (26) as a parametric function of the gain–loss parameter $\varepsilon _i$. Two branches, (a) $+x-x$ and $+x+x$ and (b) $+y-y$ and $+y+y$, are observed. As $\varepsilon _i$ increases, the two eigenfrequencies approach each other until they coalesce in an exceptional point. After passing the exceptional point, they are repelled in their imaginary parts. The real parts of the permittivity of the cylinders are $\textrm {Re} \, \varepsilon _1=\textrm {Re} \, \varepsilon _2=2.13$, $\varepsilon _v=1$, the radius $R_a=R_b=0.03\mu$m and the gap $\Delta =0.04\mu$m. The graphene parameters are $T=300$K, $\gamma _a=\gamma _b=0.1$meV and $\mu _a=\mu _b=0.5$eV.

Download Full Size | PDF

We now consider the case of polarization along the $y$ axis, in which each of the induced dipole moments are along $\pm y$ direction, i.e., the local magnetic field associated to each of the cylinders has a dependence $\approx \cos \phi _j$ ($j=a,\,b$). As a consequence, the amplitudes $b_{j 1}=-b_{j -1}$ ($j=a,\,b$). In this case, the matrix equation for these amplitudes is

(28)$$\left[ \begin{array}{llll} 1 & - 2\,s_a\,H'_1(z) \\ - 2\,s_b\,H'_1(z) & 1 \\ \end{array} \right] \left( \begin{array}{lll} b_{a 1}\\ b_{b 1}\\ \end{array} \right)= \left( \begin{array}{lll} 0\\ 0\\ \end{array} \right),$$

where we have used $2 H'_1(z)/z=H_0(z)-H_2(z)$, and $H'_1(z)$ is the derivative of $H_1(z)$. Following the same steps as in the case of $x$ polarization, the matrix Eq. (28) takes the form

(29)$$\left[ \begin{array}{llll} \omega-\omega_a & A_a/R_{ab}^2 \\ A_b/R_{ab}^2 & \omega-\omega_b \\ \end{array} \right] \left( \begin{array}{lll} b_{a 1}\\ b_{b 1}\\ \end{array} \right)= \left( \begin{array}{lll} 0\\ 0\\ \end{array} \right),$$

where $A_j$ are given by Eq. (18). Note that the matrix in Eq. (29) differs from that in Eq. (24) only in the signs of the non diagonal terms. Therefore, the eigenfrequencies corresponding to induced dipole moments along $y$ direction are formally given by Eq. (26). Unlike the $x$ polarization case, for which the upper frequency branch for $\varepsilon _i=0$ corresponds to a trapped mode, From Eq. (29) and taking $\varepsilon _i=0$, we see that the eigenvector associated to the upper branch verify $b_{a1}=b_{b1}$, i.e., it corresponds to a radiating mode. Conversely, it is straightforward to verify that the lower frequency branch corresponds to a trapped mode for which $b_{a1}=-b_{b1}$. It is worth noting that in this case, $\pm y$ oscillations, it is convenient to take the cut of the complex square root function $\sqrt {z}$ as the straight line $z=i w$ ($w$ real and positive) so that $\sqrt {-w^2}=-i w$. In this way, beyond the exceptional point the upper branch $+y+y$ moves away from the real axis whereas the lower branch $+y-y$ reaches the real axis, as shown in Fig. 2(b).

3. Results

We consider a system of two dielectric cylinders with permittivities $\varepsilon _a=2.13+i \varepsilon _i$ ($\varepsilon _i>0$), $\varepsilon _b=2.13-i\varepsilon _i$ for lossy and gain cylinders, respectively. The radii $R_a=R_b=R=0.03\mu$m and both cylinders are coated with a graphene monolayer and immersed in vacuum ($\varepsilon _v=1$). The graphene parameters are: temperature $T=300$K, chemical potentials $\mu _1=\mu _2=0.5$eV and the carriers scattering rates $\gamma _1=\gamma _2=0.1$meV. The positions of the cylinders are $\mathbf {r}^a=-0.05\mu$m$\hat {x}$, $\mathbf {r}^b=0.05\mu$m$\hat {x}$ (center to center distance $R_{ab}=0.1\mu$m) and the gap between them is $\Delta =0.04\mu$m. Fig. 3(a) shows the trajectory of the eigenfrequencies in the complex plane as a parametric function of the gain–loss parameter $\varepsilon _i$ calculated by solving the full retarded dispersion equation (FR). To solve this equation, we use a Newton–Raphson method adapted to treat complex variable. Four branches are observed, two of them ($+x+x$ and $+x-x$) corresponds to both cylinders polarized along the $x$ axis (dipole moments along the $x$ axis) while the other two branches ($+y-y$ and $+y+y$) corresponds to the case in which both cylinders are polarized in the $y$ axis (dipole moments along $y$ axis). We are using the notation of the asymptotic case when $\varepsilon _i=0$ for naming the dimer surface plasmons branches, so that $+y+y$ ($+x+x$) branch corresponds to the curve starting at the point for which both dipole moments move in phase on the $y$ axis ($x$ axis), point O” (point O”’) in Fig. 3(a), and the $+y-y$ ($+x-x$) branch corresponds to the curve starting at the point for which the dipole moments oscillating in opposite phase on the $y$ axis ($x$ axis), point O’ (point O) in Fig. 3(a). For instance, the branch $+x+x$ starts at frequency $\omega /c=0.8589226-i0.5789274\,10^{-3}\mu$m$^{-1}$ for $\varepsilon _i=0$, where both dipole moments are in phase, and it moves to the right side leaving away from the real axis. Moreover, the real part of this trajectory approaches asymptotically to the value $\omega /c=0.8868\mu$m$^{-1}$ corresponding to the real part of the eigenfrequency of a single graphene cylinder [36]. On the contrary, the branch $+x-x$ starts at $\omega /c=0.9107941-i0.253260710\,10^{-3}\mu$m$^{-1}$ for $\varepsilon _i=0$, with both dipole moments in opposite phase, and it approaches to the real axis where the lasing threshold at $\omega _{crit}/c \approx 0.8954615$m$^{-1}$ is reached for $\varepsilon _i=\varepsilon _{crit}=0.166$ (point A in Fig. 3(a)).

Fig. 3. (a) Trajectory of the four branches as a parametric function of the gain–loss parameter $\varepsilon _i$. The lasing threshold is achieved for branches $+x-x$ and $+y-y$ at points A and C, respectively. The gain–loss parameter values are $\varepsilon _i=0.166$, for point A, and $\varepsilon _i=0.1715275$, for point C. Points B and D correspond to states on the $+x+x$ and $+y+y$ branches for $\varepsilon _i=0.166$ and $\varepsilon _i=0.1715275$, respectively. Dashed lines correspond to QA branches plotted in Fig. 2. (b) Scattering cross section for $p$-polarized incident waves for illumination direction along ($\phi _{inc}=90^\circ$) and perpendicular ($\phi _{inc}=0$) to the axis joining the cylinder centers. The permittivity of the cylinders are $\varepsilon _a=2.13+i \varepsilon _i$, $\varepsilon _b=2.13-i \varepsilon _i$, the radius $R_a=R_b=R=0.03\mu$m and the gap $\Delta =0.04\mu$m. The graphene parameters are $T=300$K, $\gamma _a=\gamma _b=0.1$meV and $\mu _a=\mu _b=0.5$eV.

Download Full Size | PDF

On the other hand, Fig. 3(a) also shows the trajectory of the branch corresponding to the polarization along the $y$ axis. The branch $+y-y$ starts at frequency $\omega /c=0.8594583-i0.25349871\,10^{-3}\mu$m$^{-1}$ for $\varepsilon _i=0$, where both dipole moments are in opposed phase, and it moves toward the right side, reaching the real axis at point C where the critical eigenfrequency $\omega _{crit}/c=0.8755407\mu$m$^{-1}$ for a critical value of the gain–loss parameter $\varepsilon _{crit}=0.1715275$. On the other hand, the $+y+y$ branch starts at frequency $\omega /c=0.9102954-i0.5560595\,10^{-3}\mu$m$^{-1}$ for $\varepsilon _i=0$ and it moves to the left side leaving away from the real axis. As in the $x$ polarized case, both branches are repelled in the direction of the imaginary axis allowing the $+y-y$ branch to achieve the lasing threshold at point C.

In order to understand the gain–loss compensation near the critical points at which the eigenfrequencies are almost real, in Fig. 3(b) we plot the scattering cross sections for a plane wave impinging at an angle $\phi =0$ (electric field along the $x$ axis) and $\phi =90^\circ$ (electric field along the $y$ axis). The corresponding gain–loss parameter is $\varepsilon _i=0.166$ for $\phi =0$, and $\varepsilon _i=0.1715$ for $\phi =90^\circ$, i.e., it values are near to the critical values at which lasing conditions are achieved. In the $\phi =0$ case, we observe that the scattering cross section is enhanced at frequency near $\omega /c \approx 0.895$m$^{-1}$ that agree well with the lasing frequency for the $+x-x$ branch calculated by solving the eigenmode problem (point A in Fig. 3(a)). Moreover, we observe another peak (less intense) near $\omega /c=0.873\mu$m$^{-1}$, a value that falls near the real part of the eigenfrequency of a state with $\varepsilon _i=0.166$ but corresponding to the $+x+x$ branch (point B in Fig. 3(a)). A similar behavior presents the scattering curve for $\phi =90^\circ$ and $\varepsilon _i=0.1715$. From Fig. 3(b) we observe a very sharp peak at frequency that coincides with the lasing frequency $\omega _{crit}/c \approx 0.8755407\mu$m$^{-1}$ calculated by solving the eigenmode problem (point C in Fig. 3(a)). Moreover, a second peak is observed at a frequency $\approx 0.894\mu$m$^{-1}$ associated to the excitation of the state on the $+y+y$ branch for $\varepsilon _i=0.1715$ (point D in Fig. 3(a)).

The question that arises from the above results is how eigenmodes on the $+x-x$ and $+y-y$ branches, which correspond to trapped modes for $\varepsilon _i=0$, can be excited with a plane wave by varying the gain–loss parameter $\varepsilon _i$. Furthermore, these branches reach their lasing threshold for a critical value of the gain–loss parameter. To find a response, we have calculated the eigenvectors, containing all the field amplitudes of the eigenmodes. In particular, we have verified that coefficients with $|m|\not =1$ are orders of magnitude less than those corresponding to $m=\pm 1$, suggesting that the dimer plasmons can be considered, as a good approximation, as a superposition of single plasmon with $m=\pm 1$ linked to each cylinder. As a consequence, to gain further insight into the underlying physics of these branches excitations, we applied the QA as follows. Without loss of generality, we consider the case for which the induced dipole moments on cylinders are in $\pm x$ direction. By replacing Eq. (26) into Eq. (24) and using Eq. (19), we find the following relation between the dimer amplitudes

(30)$$b_{a 1}={\mp} \frac{\varepsilon_v (R/R_{ab})^2}{\sqrt{-\frac{\varepsilon_i^2}{4}+\varepsilon_v^2 (R/R_{ab})^4} + i \frac{\varepsilon_i}{2}} \, b_{b 1},$$

where the $-$ and $+$ signs correspond to the $+x-x$ and $+x+x$ branches, respectively. From this equation we see that the modulus of the $b_{a 1}$ and $b_{b 1}$ amplitudes are equal providing that the gain–loss parameter be less than $\varepsilon _{ep}$. Taking into account the fact that $b_{a -1}=b_{a 1}$ and $b_{b -1}=b_{b 1}$ and using Eq. (2), we can write the field scattered by the cylinders as

(31)$$H^{(v)}_s(\mathbf{r}) = 2 i b_{a 1} H_1(k_v r_a) \sin(\phi_a) + 2 i b_{b 1} H_1(k_v r_b) \sin(\phi_b).$$

Comparing this expression with that corresponding to a single dipole moment $p$ along the $x$ axis and centered at the origin ($H \approx p H_1(k_v r) \sin (\phi )$), we deduce that Eq. (31) corresponds to a superposition of two fields, one of them due to a dipole moment of amplitude $p_a \approx b_{a 1}$ centered at the cylinder $a$ and other due to a dipole moment of amplitude $p_b \approx b_{b 1}$ centered at the cylinder $b$. Since $|b_{a 1}|=|b_{b 1}|$, the induced dipole moments $|p_a|=|p_b|$ and, as a consequence, both emitted fields for each of the dipoles have the same intensity. Focusing our attention on the $x-x$ eigenmode, from Eq. (30) we see that the phase difference $\Phi$ between $b_{a 1}$ and $b_{b 1}$ amplitudes, and thus between $p_a$ and $p_b$, is shifted from $\pi$ to $\pi /2$ when $\varepsilon _i$ increases from $0$ up to above the value $\varepsilon _{ep}$. This fact implies that the eigenmodes on the $+x-x$ trajectory pass from trapped to bright by increasing the gain–loss parameter. Our calculation confirm this expectation, as can be seen in Fig. 4 where we show plots of the phase $\Phi$ as a function of the gain–loss parameter $\varepsilon _i$ by applying the FR dispersion rigorous method (continuous line) and using Eq. (30) (dashed line). From this Figure we see that the curve calculated with the QA agree well with that calculated by using the FR dispersion method. In particular, we see that the phase $\Phi =\pi$ for $\varepsilon _i=0$, which means that the eigenmodes are dark states, and monotonically decreases for until reach the lasing modes A (for $+x-x$ brach) and B (for $+y-y$ branch).

Fig. 4. Phase difference between the induced dipoles on each cylinders as a function of the gain–loss parameter $\varepsilon _i$ for the $+x-x$ and $+y-y$ branches. The calculations have been carried out by applying the FR dispersion method (continuous line) and the QA method (dashed line). Constitutive and geometrical parameters are the same as in Fig. 3.

Download Full Size | PDF

In order to visualize the above behavior, in Figs. 5(a) and 5(b) we have plotted the spatial distribution of the magnetic field $H_s^{v}(\mathbf {r})$ (calculated by using the FR method) corresponding to the $+x-x$ mode for $\varepsilon _i=0$ (point O in Fig. 3(a)) and $\varepsilon _i=0.166$ (point A in Fig. 3(a)), respectively. The eigenfrequencies are $\omega /c=(0.9107941-i 0.2532607\,10^{-3})\mu$m$^{-1}$ (Fig. 5(a)) and $\omega /c=(0.8954615-i 0.2850863\,10^{-5})\mu$m$^{-1}$ (Fig. 5(b)). The map of Fig. 5(b) is characterized by four zones, two of them are in red (blue) pointing out that the magnetic field is in $-z$ ($+z$) direction. This field distribution corresponds to that of two identical dipole moments placed at the center of each cylinders and oscillating with opposite phase ($\Phi =\pi$). On the other hand, from Fig. 5(b) we observe two regions, one of them is red ($y>0$) and the other is blue ($y<0$). Moreover, the field intensity near the first cylinder is notably less than that near the second cylinder, indicating that, at this time, the dipole moment on the first cylinder is less than that corresponding to the second cylinder, as shown in Fig. 5(b). This behaviour is in accordance with the fact that the phase difference between both dipole moments $p_a$ and $p_b$ is near $\Phi =\pi /2$ for the lasing mode A and, as a consequence, this mode can be excited by plane wave incidence. Similarly, in Figs. 5(c) and Fig. 5(d), we plotted the magnetic field map for the $+y-y$ mode for $\varepsilon _i=0$ (point O’ in Fig. 3(a)) and $\varepsilon _i=0.1715$ (point C in Fig. 3(a)), respectively. The magnetic field map shown in Fig. 5(c) is similar to that of two dipole moments of equal amplitude oscillating with opposite phase ($\Phi =\pi$). On the contrary, as the gain–loss parameter reach the critical value $\varepsilon _i=0.1715$, the magnetic field distribution looks like that of two dipole moments oscillating out of phase (but not in opposed phase) and, as a consequence, the lasing mode B can be excited by plane wave incidence.

Fig. 5. Spatial distribution of the magnetic field $H_s^{v}(\mathbf {r})$ at a fixed time for plasmon eigenfrequencies corresponding to points O, O’, A and C in Fig. 3(a). The resonance frequency is $\omega /c=(0.9107941-i 0.2532607\,10^{-3})\mu$m$^{-1}$ (a), $\omega /c=(0.8954615-i 0.2850863\,10^{-5})\mu$m$^{-1}$ (b), $\omega /c=(0.8594583-i 0.2534987\,10^{-3})\mu$m$^{-1}$ (c) and $\omega /c=(0.875540-i 0.6965941\,10^{-7})\mu$m$^{-1}$ (d). The arrow plotted at the center of each cylinder denotes the direction and relative amplitude of the induced electric dipole.

Download Full Size | PDF

It is worth noting some similarities and differences between the eigenmode branches calculation by using the FR dispersion equation and those calculated by using the QA. On the one hand, the $+x+x$ and $+x-x$ branches in Fig. 3(a) are repelled in the direction of the imaginary axis, as predicted by the QA for a PT-symmetric system (Fig. 2), a feature that allows the $+x-x$ branch to achieve the lasing threshold at point A. Furthermore, the critical value (27) calculated by using the QA results $\varepsilon _{ep}=0.18$, a value that agree well with the gain–loss parameter for which the lasing threshold is achieved in Fig. 3(a). Moreover, the phase difference between the induced dipole moments on each cylinders calculated by FR and QA methods matches quite well. On the other hand, $+x+x$ and $+x-x$ branches in Fig. 3(a) do not start at points with the same value of their imaginary parts as occur for the branches in Fig. 2 and, as a consequence, these trajectories does not coalesce at an exceptional point as in Fig. 2. This is true because the lack of a radiation losses term in the QA, i.e., the radiation losses would prevent the system from having all the properties of a full loss compensated PT-symmetric structure, shown in Fig. 2, such as the existence of an exceptional point.

In order to study the behavior with the chemical potential on graphene, we set $\mu _a=\mu _b=\mu$ and vary the values of $\mu$. In Fig. 6 we have plotted the lasing frequency $\omega _{crit}$ and the corresponding gain–loss parameter $\varepsilon _{crit}$ as functions of $\mu$, calculated with the rigorous FR method. From Fig. 6(a) we can see that the lasing frequency for both $+x-x$ and $+y-y$ branches are increasing functions of $\mu$. This fact can be understood by taking into account that the real part of the eigenfrequency corresponding to a single graphene cylinder, which falls between the lasing frequencies for $+y-y$ and $+x-x$ branches, is proportional to $\sqrt {\mu }$ (see Eq. (20)). On the other hand, from Fig. 6(b) we see that the critical value of the gain–loss parameter for which the lasing condition is achieved is a decreasing function of the chemical potential $\mu$. This behaviour has a similarity with that presented by a single graphene cylinder for which has been demonstrated that the gain level to achieve the lasing condition decreases with the chemical potential value [37].

Fig. 6. Lasing frequency (a) and critical gain–loss parameter (b) as functions of the chemical potential $\mu$ ($\mu _a=\mu _b=\mu$) for branches $+x-x$ and $+y-y$. The calculations have been carried out by applying the FR dispersion method. Geometrical and other constitutive parameters are the same as in Fig. 3.

Download Full Size | PDF

4. Conclusions

In conclusion, we have analytically studied the scattering and the eigenmode problems for a dimer composed of two graphene coated dielectric cylinders, one of them with loss and the other with the same level of gain. We have demonstrated the existence of two branches, corresponding to trapped modes when $\varepsilon _i=0$, that reach the lasing conditions for suitable values of the gain–loss parameter. While the phase difference between the induced dipole moments on individual cylinders changes from $\pi$ to a value near $\pi /2$ when the gain–loss parameter is incremented, a fact that provides the mode transformation from trapped to radiating modes, the modulus of each individual dipole moments maintains equals in between.

Other mechanisms to transform a trapped mode into a resonant observable which can be excited by a plane wave have been reported in other works. All these methods are based on the introduction of a small asymmetry with respect to the center of the dimer by producing dissimilar dipole moment modulus on each of the cylinders. Interestingly, here, we found that in the transformation from trapped to radiating eigenmodes on both $+x-x$ and $+y-y$ branches, the modulus of the individual dipole moments does not change, while it is changed the phase between them. A distinction between these two kind of mechanism, the phase variation and modulus variation mechanisms, to transform a trapped mode into a resonant observable has not been reported before to our knowledge. We believe that our results will be usefull for a deeper understanding of the PT-symmetric LSP characteristics and that the LSP effects we have demonstrated opens up possibilities for practical applications involving sub-wavelength laser structures.

Appendix A. Graphene conductivity

We consider the graphene layer as an infinitesimally thin, local and isotropic two-sided layer with frequency-dependent surface conductivity $\sigma (\omega )$ given by the Kubo formula [38], which can be read as $\sigma _{loc}= \sigma ^{intra}+\sigma ^{inter}$, with the intraband and interband contributions being

(A1)$$\sigma^{intra}(\omega)= \frac{2i e^2 k_B T}{\pi \hbar^2 (\omega+i\gamma_g)} \mbox{ln}\left[2 \mbox{cosh}(\mu_g/2 k_B T)\right],$$

(A2)$$\begin{aligned} \sigma^{inter}(\omega)&= \frac{e^2}{\hbar} \Bigg\{ \frac{1}{2}+\frac{1}{\pi}\mbox{arctan}\left[(\hbar \omega-2\mu_g)/2k_BT\right]-\\ &\frac{i}{2\pi}\mbox{ln}\left[\frac{(\hbar \omega+2\mu_g)^2}{(\hbar \omega-2\mu_g)^2+(2k_BT)^2}\right] \Bigg\}, \end{aligned}$$

where $\mu _g$ is the chemical potential (controlled with the help of a gate voltage), $\gamma _g$ the carriers scattering rate, $e$ the electron charge, $k_B$ the Boltzmann constant and $\hbar$ the reduced Planck constant.

Appendix B. Nonlocal response in graphene conductivity

Nonlocal, or spatial dispersion effects, in graphene conductivity can present deviations from the results predicted by the local conductivity (A1) and (A2). In general, these effects appear when the graphene-based structure size is of a few nanometers. At this size scale, the surface plasmon propagates with low phase velocity $v_{SP}$ approaching to the electron Fermi velocity $v_{F} \approx 10^6\,m/s$. We can estimate the value of $v_F/v_{SP}$ using the quasistatic approximation as follows. Since the surface plasmon effective momentum $k_{SP}$, which is along the azimutal angle ($\phi$ axis), can be written as [36],

(B1)$$k_{SP}=\frac{m}{R},$$

if we consider $m=1$ (dipolar order), the factor

(B2)$$\frac{v_F}{v_{SP}}=\frac{v_F}{R \, \omega} \approx \frac{10^{12} \mu \mbox{m}/s}{ 0.03\mu \mbox{m} \, 3 \times 10^{14} 1/s } \approx \frac{1}{9},$$

where we have consider a working frequency (frequency range for plasmon resonances in considered examples) $\omega \approx 1 \mu$m$^{-1} \, c \approx 10^{14}\mu \mbox {m}/s$ and $R=0.03\mu$m. With this value of $v_F/v_{SP}$, the graphene surface conductivity slightly deviates from the value corresponding to the local approximation. This fact can be viewed by considering the graphene conductivity including non local response, as in [39], through the hydrodynamic approximation,

(B3)$$\sigma(\omega,k)=\frac{\sigma_{loc}(\omega)}{1-\frac{\beta^2 \, k^2}{\omega^2}},$$

where $\beta =\sqrt {\frac {3}{4}} v_F$. By replacing the value of $k_{SP}$ given by Eq. (B1) with $m=1$, we obtain a value which differ less than $0.8$% from the local conductivity $\sigma _{loc}(\omega )$. As a consequence, we have neglected the non local response on graphene conductivity.

Appendix C. Developing of the right hand side in Eq. (25)

The right side of Eq. (25) can be calculated by using Eq. (18) as follows.

(C1)$$\frac{A_a A_b}{R_{ab}^4} = \frac{a^4}{R_{ab}^4} \frac{[\varepsilon_1+i\varepsilon_i-\varepsilon_v+g(\omega_a)]}{g'(\omega_a)} \frac{[\varepsilon_1-i\varepsilon_i-\varepsilon_v+g(\omega_b)]}{g'(\omega_b)},$$

where

(C2)$$g(\omega_j)={-}\frac{\omega_g^2}{\omega_j^2+i\gamma_g \omega_j},$$

and

(C3)$$g'(\omega_j)=\frac{\omega_g^2 (2 \omega_j+i \gamma_g)}{(\omega_j^2+i\gamma_g \omega_j)^2}.$$

By taking into account that $\gamma _g<<|\omega _j|$, we can write this equations as,

(C4)$$g(\omega_j)={-}\frac{\omega_g^2}{\omega_j^2}\frac{1}{1+i\gamma_g/\omega_j}\approx{-}\frac{\omega_g^2}{\omega_j^2} \left(1-i\frac{\gamma_g}{\omega_j}+\frac{\gamma_g^2}{\omega_j^2}+\cdots\right),$$

and

(C5)$$\begin{aligned} g'(\omega_j) &\approx \frac{ 2 \omega_g^2}{\omega_j^3} \left(1-i\frac{\gamma_g}{\omega_j}+\frac{\gamma_g^2}{\omega_j^2}+\cdots\right)-\frac{\omega_g^2}{\omega_j^2} \left(i\frac{\gamma_g}{\omega_j^2}-2 \frac{\gamma_g^2}{\omega_j^3}+\cdots\right) =\\ &\frac{2 \omega_g^2}{\omega_j^3} \left( 1 - \frac{3}{2}i\frac{\gamma_g}{\omega_j}+2\frac{ \gamma_g^2}{\omega_j^2}+\cdots\right). \end{aligned}$$

Considering the lowest order in $\Gamma _{sup}/\omega _0$ and $\Gamma _{cover}/\omega _0$, functions $g(\omega _j)$ and $g'(\omega _j)$ are written as

(C6)$$g(\omega_j)\approx{-}\frac{\omega_g^2}{\omega_0^2} \left[1 \pm \frac{i}{\omega_0}2\Gamma_{core}\right],$$

and

(C7)$$g'(\omega_j) \approx \frac{2 \omega_g^2}{\omega_0^3} \left[1 \pm i\frac{3}{\omega_0}\Gamma_{core}\right].$$

Signs $+$ and $-$ correspond to $j=a$ and $j=b$, respectively. In this approximation, Eq. (C1) reduces to

(C8)$$\frac{A_a A_b}{R_{ab}^4} = \frac{R^4}{R_{ab}^4} \frac{\omega_g^2}{(\varepsilon_1+\varepsilon_v)} \varepsilon_v^2.$$

Disclosures

The authors declare no conflicts of interest.

References

1. S. Prosvirnin and S. Zouhdi, “Resonances of closed modes in thin arrays of complex particles,” in Advances in Electromagnetics of Complex Media and Metamaterials, (Springer, 2002).

2. V. Fedotov, M. Rose, S. Prosvirnin, N. Papasimakis, and N. Zheludev, “Sharp trapped-mode resonances in planar metamaterials with a broken structural symmetry,” Phys. Rev. Lett. 99(14), 147401 (2007). [CrossRef]

3. N. I. Zheludev, S. Prosvirnin, N. Papasimakis, and V. Fedotov, “Lasing spaser,” Nat. Photonics 2(6), 351–354 (2008). [CrossRef]

4. D. J. Bergman and M. I. Stockman, “Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems,” Phys. Rev. Lett. 90(2), 027402 (2003). [CrossRef]

5. M. I. Stockman, “Spasers explained,” Nat. Photonics 2(6), 327–329 (2008). [CrossRef]

6. A. Krasnok and A. Alù, “Active nanophotonics,” Proc. IEEE 108(5), 628–654 (2020). [CrossRef]

7. D. Wang, W. Wang, M. P. Knudson, G. C. Schatz, and T. W. Odom, “Structural engineering in plasmon nanolasers,” Chem. Rev. 118(6), 2865–2881 (2018). [CrossRef]

8. S. I. Azzam, A. V. Kildishev, R.-M. Ma, C.-Z. Ning, R. Oulton, V. M. Shalaev, M. I. Stockman, J.-L. Xu, and X. Zhang, “Ten years of spasers and plasmonic nanolasers,” Light: Sci. Appl. 9(1), 90 (2020). [CrossRef]

9. R. F. Oulton, V. J. Sorger, T. Zentgraf, R.-M. Ma, C. Gladden, L. Dai, G. Bartal, and X. Zhang, “Plasmon lasers at deep subwavelength scale,” Nature 461(7264), 629–632 (2009). [CrossRef]

10. M. Moccia, G. Castaldi, A. Alù, and V. Galdi, “Harnessing spectral singularities in non-hermitian cylindrical structures,” IEEE Trans. Antennas Propag. 68(3), 1704–1716 (2020). [CrossRef]

11. V. M. Alvarez, J. B. Vargas, M. Berdakin, and L. F. Torres, “Topological states of non-hermitian systems,” Eur. Phys. J. Spec. Top. 227(12), 1295–1308 (2018). [CrossRef]

12. M.-A. Miri, P. LiKamWa, and D. N. Christodoulides, “Large area single-mode parity–time-symmetric laser amplifiers,” Opt. Lett. 37(5), 764–766 (2012). [CrossRef]

13. F. Monticone, C. A. Valagiannopoulos, and A. Alù, “Parity-time symmetric nonlocal metasurfaces: all-angle negative refraction and volumetric imaging,” Phys. Rev. X 6(4), 041018 (2016). [CrossRef]

14. C. M. Bender, R. Tateo, P. E. Dorey, T. C. Dunning, G. Levai, S. Kuzhel, H. F. Jones, A. Fring, and D. W. Hook, PT symmetry: In quantum and classical physics (World Scientific Publishing, 2018).

15. C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, “Observation of parity–time symmetry in optics,” Nat. Phys. 6(3), 192–195 (2010). [CrossRef]

16. R. A. Depine, Graphene Optics: Electromagnetic solution of canonical problems (Morgan & Claypool Publishers, 2017).

17. Y. Fan, N.-H. Shen, T. Koschny, and C. M. Soukoulis, “Tunable terahertz meta-surface with graphene cut-wires,” ACS Photonics 2(1), 151–156 (2015). [CrossRef]

18. Y. Fan, N.-H. Shen, F. Zhang, Q. Zhao, H. Wu, Q. Fu, Z. Wei, H. Li, and C. M. Soukoulis, “Graphene plasmonics: a platform for 2d optics,” Adv. Opt. Mater. 7(3), 1800537 (2019). [CrossRef]

19. P.-Y. Chen and J. Jung, “P t symmetry and singularity-enhanced sensing based on photoexcited graphene metasurfaces,” Phys. Rev. Appl. 5(6), 064018 (2016). [CrossRef]

20. X. Lin, R. Li, F. Gao, E. Li, X. Zhang, B. Zhang, and H. Chen, “Loss induced amplification of graphene plasmons,” Opt. Lett. 41(4), 681–684 (2016). [CrossRef]

21. O. Zhernovnykova, O. Popova, G. Deynychenko, T. Deynichenko, and Y. V. Bludov, “Surface plasmon-polaritons in graphene, embedded into medium with gain and losses,” J. Phys.: Condens. Matter 31(46), 465301 (2019). [CrossRef]

22. W. Zhang, T. Wu, and X. Zhang, “Tailoring eigenmodes at spectral singularities in graphene-based pt systems,” Sci. Rep. 7(1), 11407 (2017). [CrossRef]

23. A. Abdrabou and Y. Y. Lu, “Exceptional points for resonant states on parallel circular dielectric cylinders,” J. Opt. Soc. Am. B 36(6), 1659–1667 (2019). [CrossRef]

24. W. Chen, J. Zhang, B. Peng, Ş. K. Özdemir, X. Fan, and L. Yang, “Parity-time-symmetric whispering-gallery mode nanoparticle sensor,” Photonics Res. 6(5), A23–A30 (2018). [CrossRef]

25. V. Klimov, I. Zabkov, D. Guzatov, and A. Vinogradov, “Loss compensation symmetry in dimers made of gain and lossy nanoparticles,” Laser Phys. Lett. 15(3), 035901 (2018). [CrossRef]

26. M.-A. Miri, M. Cotrufo, and A. Alù, “Anomalous optical forces in pt-symmetric waveguides,” Opt. Lett. 44(14), 3558–3561 (2019). [CrossRef]

27. M. Parto, Y. G. Liu, B. Bahari, M. Khajavikhan, and D. N. Christodoulides, “Non-hermitian and topological photonics: optics at an exceptional point,” Nanophotonics 10(1), 403–423 (2020). [CrossRef]

28. D. M. Natarov, T. M. Benson, and A. I. Nosich, “Electromagnetic analysis of the lasing thresholds of hybrid plasmon modes of a silver tube nanolaser with active core and active shell,” Beilstein J. Nanotechnol. 10, 294–304 (2019). [CrossRef]

29. A. S. Kupriianov, Y. Xu, A. Sayanskiy, V. Dmitriev, Y. S. Kivshar, and V. R. Tuz, “Metasurface engineering through bound states in the continuum,” Phys. Rev. Appl. 12(1), 014024 (2019). [CrossRef]

30. M. Cuevas, S. H. Raad, and C. J. Zapata-Rodríguez, “Coupled plasmonic graphene wires: theoretical study including complex frequencies and field distributions of bright and dark surface plasmons,” J. Opt. Soc. Am. B 37(10), 3084–3093 (2020). [CrossRef]

31. F. Zhang, C. Li, Y. Fan, R. Yang, N.-H. Shen, Q. Fu, W. Zhang, Q. Zhao, J. Zhou, T. Koschny, and C. M. Soukoulis, “Phase-modulated scattering manipulation for exterior cloaking in metal–dielectric hybrid metamaterials,” Adv. Mater. 31(39), 1903206 (2019). [CrossRef]

32. M. Naserpour, C. J. Zapata-Rodríguez, S. M. Vuković, H. Pashaeiadl, and M. R. Belić, “Tunable invisibility cloaking by using isolated graphene-coated nanowires and dimers,” Sci. Rep. 7(1), 12186–14 (2017). [CrossRef]

33. M. Abramowitz and I. A. Stegun, “Handbook of mathematical functions with formulas, graphs, and mathematical table,” in US Department of Commerce, (National Bureau of Standards Applied Mathematics series 55, 1965).

34. D. Felbacq, G. Tayeb, and D. Maystre, “Scattering by a random set of parallel cylinders,” J. Opt. Soc. Am. A 11(9), 2526–2538 (1994). [CrossRef]

35. M. Cuevas, “Enhancement, suppression of the emission and the energy transfer by using a graphene subwavelength wire,” J. Quant. Spectrosc. Radiat. Transfer 214, 8–17 (2018). [CrossRef]

36. M. Cuevas, M. A. Riso, and R. A. Depine, “Complex frequencies and field distributions of localized surface plasmon modes in graphene-coated subwavelength wires,” J. Quant. Spectrosc. Radiat. Transfer 173, 26–33 (2016). [CrossRef]

37. L. Prelat, M. Cuevas, N. Passarelli, R. Bustos Marún, and R. Depine, “Theoretical study of a graphene based–localized surface plasmon sapser,” Nanophotonics and Nanoplasmonics, Frontiers in Optics, To be published (2020).

38. L. A. Falkovsky, “Optical properties of graphene and iv–vi semiconductors,” Phys.-Usp. 51(9), 887–897 (2008). [CrossRef]

39. T. Christensen, A.-P. Jauho, M. Wubs, and N. A. Mortensen, “Localized plasmons in graphene-coated nanospheres,” Phys. Rev. B 91(12), 125414 (2015). [CrossRef]

Lasing condition for trapped modes in subwavelength-wired PT-symmetric resonators

Abstract

1. Introduction

2. Theory

2.1 Rigorous description of the fields and scattering efficiencies

2.2 Quasistatic approximation: a simple model based on two coupled electric dipoles

3. Results

4. Conclusions

Appendix A. Graphene conductivity

Appendix B. Nonlocal response in graphene conductivity

Appendix C. Developing of the right hand side in Eq. (25)

Disclosures

References

Cited By

Figures (6)

Equations (44)

Optics Express