1. Introduction

It is well-known that there should be trade-offs between the desired confinement effect and the unavoidable dissipation effect when metallic nature is utilized for realizing plasmonic devices [1,2]. For metallic media, the classical Drude model with real relaxation rates needs to be amended in some situations. For instance, thin films containing subwavelength metallic nanoparticles are well modeled by taking the relaxation rate (CRR) γ to be complex-valued [3]. The concept of CRRs has relationships to systems in both classical mechanics [4] and quantum mechanics [5–8]. Recently, superconducting materials themselves are explored for the incorporation into metamaterials [9], thanks to their low losses. The concept of CRRs has been successfully employed in dealing with the vortex dynamics for type-II superconductors through the time-dependent Ginzburg-Landau (TDGL) equation [7], where the longitudinal conductivity and the (perpendicular) Hall conductivity are affected by parameters like ${\gamma }_r$ and ${\gamma }_i$, respectively.

Here, several consequences of incorporating the CRRs in the Drude model will be examined in terms of wave propagations [10–12]. Besides, the resulting symmetry breaking will be examined in terms of the electromagnetic energy flows. Let us consider the conventional Drude model for the dielectric constant of metal as follows, with ${\omega }_p$ as the bulk plasma frequency.

2. Electric conductivity

In view of the circuit theory for nanophotonics [13,14], let us examine the notion of conventional complex electrical conductivity ${\sigma }_0$ defined through $\epsilon \left(\omega \right)\equiv 1+i\left({\sigma }_0/\omega \right)$. In consequence, we are led to the relationship ${\sigma }_0={\left(\gamma -i\omega \right)}^{-1}{\omega }_p^2$ [3]. In order to account for the fact ${\gamma }_i\ne 0$ of the CRR, let us consider the following modified conductivity $\sigma \equiv {\sigma }_r+i{\sigma }_i$ (without any subscript).

In this aspect, we learn that only positive imaginary parts of the complex relaxation time or $\mathrm{Im}\left({\gamma }^{-1}\right)>0$ or ${\gamma }_i<0$ work properly in the case with the TDGL model [6,7]. We ascribe this restriction to the nonlinearity imbedded in the TDGL model, which helps to select only one of the two possible signs. However, both signs of ${\gamma }_i$ work fine in the present model, since the problem under current investigation is essentially linear. For matter of convenience, we take however that ${\gamma }_i<0$ in the ensuing numerical computations.

Meanwhile, we make an interesting comparison to the axial complex conductance examined in [5]. They examined a quantum wire by solving the Schroedinger equation under several assumptions, for instance, prescribed electronic screening functions and Gaussian radial profiles of electric fields. Their model can be represented as a parallel circuit of one $L_1-R_1-C_1$ subcircuit and another $L_0-R_0$ subcircuit. In line with the above-mentioned circuit model [13], they defined another complex electrical conductivity ${\sigma }_c$ (the subscript “c” implying Cuniberti in [5]).

In Fig. 1(a) , we display ${\sigma }_r$ and ${\sigma }_i$ for a set of data: ${\omega }_p=1$, ${\gamma }_r=0.2$, and ${\gamma }_i=0.5$. Returning to our model in Eq. (3), we find that symmetry holds true with respect to ω such that ${\sigma }_r\left(-\omega ,{\gamma }_i,-{\gamma }_i\right)={\sigma }_r\left(\omega ,{\gamma }_i,{\gamma }_i\right)$ and ${\sigma }_i\left(-\omega ,{\gamma }_i,-{\gamma }_i\right)=-{\sigma }_i\left(\omega ,{\gamma }_i,{\gamma }_i\right)$. Therefore, the corresponding behaviors for ${\gamma }_i=-0.5$ can be easily inferred. It is surprising that Fig. 1(a) turns out to be almost identical to Fig. 4 in [5]. Upon closer examinations, the dc limits as $\omega \to 0$ are different. In other words, we have a linear dependence ${\lim }_{\omega \to 0}{\sigma }_r=a_0+a_2\omega$ with $\left(a_0,a_1\right)$ defined in Eq. (4), whereas there is a quadratic dependence ${\lim }_{\omega \to 0}\mathrm{Re}\left({\sigma }_c\right)=b_0+b_2{\omega }^2$ with $\left(b_0,b_1\right)$ as another set of constants defined in [5].

 

Fig. 1 (a) The real and imaginary parts of complex conductance for the Drude model according to Eq. (3). (b) A single planar interface between the upper dielectric and the lower metal-like medium. The in-plane wave number is denoted by k directed in the x-direction. Non-zero ${\gamma }_r$ induces cross-interface energy flows, whereas non-zero ${\gamma }_i$ entails an imbalance between the forward (“F”) and backward (“B”) waves along the interface.

Download Full Size | PDF

 

Fig. 4 The degree of asymmetry e obtained by varying one parameter while fixing the remaining two: $e\left(k\right)$ based on the data of Fig. 2(a) with ${\gamma }_r$ as an abscissa; $e\left({\gamma }_r\right)$ based on the data of Fig. 3(a) with k as an abscissa; and $e\left({\gamma }_i\right)$ vs. ${\gamma }_i$.

Download Full Size | PDF

3. Dissipation-induced stability

For the Drude model in Eq. (1) with a CRR, let us consider the propagations of electromagnetic waves. The problem of surface plasmon resonances (SPRs) occurring across a planar metal-dielectric (M-D) single interface is simplest; yet, it reveals fundamental aspects of SPRs, since diverse geometrical configurations involve a single M-D interface in one way or another [1,15]. In addition, the lifetime of SPRs is of a major concern in applications such as biological sensing [16], in particular, as regards two counter-propagating waves [17].

Let us consider the schematic Fig. 1(b) for air (with $j=a$) and metal-like medium (with $j=m$) separated by a single M-D interface at $z=0$. Here, the “metal-like” zone is more general in implying not only pure metals but also metallic composites. The relative dielectric constant is taken such that ${\epsilon }_a=1$ without loss of generality and ${\epsilon }_m\equiv \epsilon$ for simplicity. With $c_0$ as the light speed in air, we define for nondimensionalization the reference time and length to be $1/{\omega }_p$ and $c_0/{\omega }_p$, respectively. Furthermore, let x and z be the longitudinal (in-plane) and depth (surface-normal) coordinates, respectively. Besides, t is time. Then, the phase factor $\exp \left[i\left(kx-\omega t\right)\right]$ represents propagating electromagnetic waves, where k is the longitudinal wave number. For the purpose of investigating the lifetime of SPRs, it is appropriate henceforth to take a complex-frequency approach [18,19], where $\omega \equiv {\omega }_r+i{\omega }_i$. Instead, k is real-valued.

Owing to stability for $t>0$, it is required that ${\omega }_i<0$ as seen from the expanded phase factor $\exp \left[i\left(kx-{\omega }_{r}t\right)\right]\exp \left({\omega }_{i}t\right)$. As a consequence, a long-life surface plasmon (LLSP) corresponds to the solution with a smaller temporal attenuation $\left|{\omega }_i\right|$ [18–21]. On the other hand, we allow ${\omega }_r$ to take either positive or negative sign. In order to understand the meaning of different signs of ${\omega }_r$, let us take $k>0$ for convenience. For ${\omega }_r>0$, waves are propagating rightward, because $dx/dt={\omega }_r/k>0$ for constant phase of $kx-{\omega }_{r}t$. In comparison, for ${\omega }_r<0$, waves are propagating leftward, because $dx/dt={\omega }_r/k<0$ for constant $kx-{\omega }_{r}t$. Besides, the real-valued k corresponds either to an incident angle for $\left|{\omega }_r\right|>k$ or to the degree of an evanescent coupling for $\left|{\omega }_r\right|<k$ [10,15,17].

Transverse-magnetic (TM) fields are described by a set of non-zero field components $\left(E_x,H_y,E_z\right)$. For such TM waves, the dispersion relation is given by the well-known formula $\epsilon {\left(k^2-{\omega }^2\right)}^{1/2}=-{\left(k^2-\epsilon {\omega }^2\right)}^{1/2}$ [1]. We then square the latter to obtain the following fourth-order (quartic) polynomial in ω.

Let us examine some properties of the complex-valued ω’s as solutions to Eq. (5) for a complex-valued γ, but with a real-valued k. First of all, we notice a mirror symmetry in Eq. (5) for a given $k^2$, because $\left(-\omega ,-\gamma \right)$ is a solution once $\left(\omega ,\gamma \right)$ is a solution [22]. Furthermore, let us denote the complex conjugate of ω by $\overline{\omega }\equiv {\omega }_r-i{\omega }_i$ and hence ${\left|\omega \right|}^2\equiv \omega \overline{\omega }\equiv {\omega }_r^2+{\omega }_i^2$. When eliminating ${\gamma }_i$ from Eq. (5), we obtain a ${\gamma }_i$-free dispersion equation ${\left|D\right|}^2{\gamma }_r=\mathrm{Im}\left(\overline{D}N\right)$. When further processed, the last relation reduces to

As regards the four complex roots of ω as solutions to Eq. (5) for a given set of real data $\left(k,{\gamma }_r,{\gamma }_i\right)$, it is remarked that a quadratic equation can be analytically solved by the Ferrari-Cardano method. Its geometric interpretation is lucidly illustrated in [23], and it calls for the standard Cardano method in solving an intermediate equation which is cubic in ω [20,21]. Once four complex roots are analytically found, we have checked the validity of the roots so evaluated by ensuring that Eq. (5) is satisfied as well. Besides, we remark that it is non-trivial to sort a set of four roots into proper branches as a parameter is varied. As regards ${\gamma }_r$, we consider both $0<{\gamma }_r<<1$ for low-loss noble metals and $0<<{\gamma }_r=O\left(1\right)$ for high-loss gas plasmas [21]. Here, we call the attention that the value of ${\gamma }_r$ is referred to the plasma frequency ${\omega }_p$.

In order to reveal non-negligible effects of ${\gamma }_i$, we will work with a relatively large value of ${\gamma }_i=-0.2$. On atomic scales, non-zero values of ${\gamma }_i$ may be caused by the structural asymmetry [15]. On larger scales, let us suppose that fine metal nanoparticles, whose sizes are much less than the wavelength of interest, are uniformly dispersed throughout a thin dielectric film [3,13]. The relaxation rate can then be considered as an average value for such a composite film with moderately higher filling ratios. In this case, the way the nanoparticles are dispersed may be responsible for the apparent asymmetry. It is the reason why the lower zone in Fig. 1(b) is denoted as “metal-like” rather than just as “purely metallic”.

4. Lifetimes and Wave Characteristics

Henceforth, let us present typical numerical results. We define a state by an eigenvalue pair on the $\left({\omega }_r,{\omega }_i\right)$-plane. Furthermore, we define a “branch” by a group of such states possessing similar characters as shown in Figs. 2(a) and 3(a) . Each branch is composed of 1000 states in equal intervals of the parametric variations for the current numerical evaluations. According to the naming convention in Eq. (6), we call four branches “L+”, “L-“, “H+”, and “H-”, respectively. The first letter “L” or “H” is decided by the largeness of $\left|{\omega }_r\right|$, whereas the second symbol “+” or “-” is determined by the sign of ${\omega }_r$.

 

Fig. 2 (a) Migrations of four roots on the complex $\left({\omega }_r,{\omega }_i\right)$-plane. (b) The quality factor plotted against ${\omega }_r$. The curves are generated by varying kover the range $0\le k\le 2$, where the arrows mean the direction of increasing k. The filled circles in pink color indicate the starting states at $k=0$, whereas the diamonds indicate the states on the half-way at $k=1$. In the inset of panel (a), we show the symmetric trajectories for the specified data ${\gamma }_r=0.2$ and ${\gamma }_i=0$ (namely, a real-valued γ).

Download Full Size | PDF

 

Fig. 3 (a) Migrations of four roots on the complex $\left({\omega }_r,{\omega }_i\right)$-plane. The arrows imply the direction of increasing ${\gamma }_r$. The filled circles in pink color indicate the starting states at ${\gamma }_r=0$ (where ${\omega }_i=0$, thus being neutrally stable), whereas the diamonds indicate the states on the half-way at ${\gamma }_r=1.5$. In the inset of panel (a), we show the symmetric trajectories for the prescribed data ${\gamma }_i=0$ (a real-valued γ) and $k=0.5$. (b) Migrations of four roots on the complex $\left(v_{gr},v_{gi}\right)$-plane. The thick horizontal arrow in shaded colors indicates the region of the subluminal group velocity, namely, $\left|v_{gr}\right|<k=0.5$.

Download Full Size | PDF

Figure 2(a) displays such four branches in the case of $\gamma =0.2-0.2i$ with a variation over $0\le k\le 2$. The two middle branches L- and L + starting out right from the origin $\left({\omega }_r,{\omega }_i\right)=\left(0,0\right)$ are acoustic by definition, whereas the two side branches H- and H + starting off the origin are optic. Branches L- and H- with ${\omega }_r<0$ exhibit heavier attenuations (i.e., larger $\left|{\omega }_i\right|$) than the respective branches L + and H + with ${\omega }_r>0$. Therefore, the rightward- and leftward-propagating waves have unequal $\left|{\omega }_r\right|$ for a given k. Consequently, there cannot be standing waves formed as a result from two counter-propagating waves of equal amplitude colliding each other [10–12].

The fact that branches with ${\omega }_r<0$ exhibit larger attenuation in Fig. 2(a) is corroborated by the two electrostatic roots marked by “static” with filled symbols, which are located just below the respective branches L- and L + . This electrostatic resonance is described by the condition that $1+\epsilon =0$, which is equivalent to the small-wavelength-limit $k\to \infty$ as easily seen from $k\left(\omega \right)=\omega \sqrt{{\left(1+\epsilon \right)}^{-1}\epsilon }$ [22]. Substituting the Drude model into $1+\epsilon =0$ results in $2\omega \left(\omega +i\gamma \right)-1=0$, which is in turn solved to give rise to two roots $\omega =-\frac{1}{2}i\gamma \pm \frac{1}{2}{\left(2-{\gamma }^2\right)}^{1/2}$ [4]. For $\gamma =0.2-0.2i$, it is found that $\omega =-0.807-i0.114$ and $\omega =0.607-i0.0859$ on the electrostatic resonance. The inset of Fig. 2(a) exhibits four branches in this case of ${\gamma }_i=0$, where they are symmetrically placed with respect to the ${\omega }_i$-axis. In comparison, the four branches in the main panel of Fig. 2(a) are asymmetrically placed with respect to ${\omega }_r=0$, although the degree of asymmetry is still not quite large.

According to the circuit paradigm for metal optics [14], even a planar interface stands for a $L-R-C$ circuit, for which we can define a quality factor by $Q_{\omega }\equiv \frac{1}{2}\left|{\omega }_r/{\omega }_i\right|$ [9,11,18]. In Fig. 2(b), we plot $Q_{\omega }$ against ${\omega }_r$ based on the roots in Fig. 2(a). It shows on both branches H + and H- that $Q_{\omega }$ increases with k (or as $\left|{\omega }_r\right|$ increases). In comparison, on both branches L + and L-, the reverse trend is visible. For both electrostatic roots as mentioned just before, we obtain the same quality factor of $Q_{\omega }=3.54$ as marked by the horizontal broken line in Fig. 2(b). The existence of a minimum in $Q_{\omega }$ at a certain $\left|{\omega }_r\right|$ is in rough agreement with the predictions made for periodic arrays of metal cylinders [24].

As another parametric study, we plot in Fig. 3(a) the four branches on the complex ω-plane as ${\gamma }_r$ is varied over $0\le {\gamma }_r\le 3$. But, we fix the two remaining parameters such that ${\gamma }_i=-0.2$ and $k=0.5$, for an incidence angle of ${\sin }^{-1}\left(1/2\right)={30}^o$. As in Fig. 2(a), a thousand of states are computed for each branch. Figure 3(a) shows that all the four branches start on the ${\omega }_r$-axis at ${\omega }_i=0$ for ${\gamma }_r=0$, as has been proved earlier by Eq. (7). As ${\gamma }_r$ is increased (indicated by the arrows), branches H-, L + , and H + migrate downward first, but they bend upwards. Eventually they almost touch the ${\omega }_r$-axis in the large ${\gamma }_r$-limit. Such a minimum value of $\left|{\omega }_i\right|$ obtained at a finite value of ${\gamma }_r$ has been observed for gas-phase plasmas in the case of isolated nanoparticles as well [21].

In contrast, branch L- appears to be stretched downward on this display window. Indeed, we found numerically that branch L- asymptotically approaches the ${\omega }_r$-axis in the large ${\gamma }_r$-limit. For this purpose, let us examine Fig. 3 in more detail for the asymptotic behaviors of the four branches as ${\gamma }_r\to \infty$. In this limit, Eq. (5) lends itself to the following four component solutions, being labeled with the corresponding branch L-, H-, H + , and L + as shown in Fig. 4.

Therefore, the relative magnitude of $\left|{\omega }_r\right|$ is exchanged between branches L- and H- as ${\gamma }_r$ increases. As in Fig. 2(a), Fig. 3(a) shows an asymmetry in $\left|{\omega }_r\right|$. All the more clearly, Fig. 3(a) depicts that the accompanying attenuation $\left|{\omega }_i\right|$ is widely different among the four branches. As a consequence, branch H- appears to live longest among the four, thanks to its smallest $\left|{\omega }_i\right|$.

Figure 3(a) displays two additional curves denoting the electrostatic states as found for Fig. 2(a). In particular, Fig. 3(a) shows that branch L- and L + run asymptotically parallel to their corresponding electrostatic branches. There are always propagating waves with ${\omega }_r\ne 0$ as long as ${\gamma }_i\ne 0$, irrespectively of the largeness of ${\gamma }_r$. In comparison, it is found that there are no propagating waves for large ${\gamma }_r$ in the case of ${\gamma }_i=0$ because ${\omega }_r=0$, which is seen within the inset of Fig. 3(a) on the central portion of the branch (indicated by the upward arrow within the inset) [20,21].

Although the meaning of group velocity gets less clear for complex-valued γ, we can still formally define it by

In order to summarize the results in Fig. 3(a), we define the degree of asymmetry e as follows.

There are three conceivable levels of asymmetry as regards the SPRs for the M-D configuration shown in Fig. 1(b). First, the zero-level asymmetry is concerned with the direction normal to the M-D interface. For $\gamma =0$, the time-averaged Poynting vectors are directed along the M-D interface in both media, thus disallowing cross-interface energy flow. In this case, the energy flows are directed along the M-D interface, although the flow in air is opposite to that in the metal-like medium. On the other hand, a non-zero positive material loss $\gamma >0$ causes non-zero cross-interface energy flow [22]. This first-level asymmetry is therefore associated with the real part of the relaxation rate ${\gamma }_r$, and the resulting effect takes the form of energy flow in the direction normal to the M-D interface. The second-level asymmetry is incurred if ${\gamma }_i\ne 0$, for which there will be no standing waves in the direction along the M-D interface, because the forward- and backward-propagating waves viewed along the M-D interface do not balance among them (${\omega }_r$’s being asymmetric with respect to ${\omega }_r=0$).

Other symmetry-breaking phenomena occur with the nonlinear second-harmonic generation due to applied magnetic fields [15], which is closely related to the asymmetry with magnetic fields as occurs with the TDGL for superconductors [6,7]. The occurrence of complex eigenvalues can also be found in the Hatano-Nelson model of superconductors due to the asymmetric effect of hopping electrons [8].

5. Complex frequency versus complex wave number

We have already taken frequency to be real-valued in Sec. 2, whereas frequency is complex-valued in Sec. 4. In general, whether frequency is real- or complex-valued depends on how to properly to describe physical phenomena under investigation [20–22].

For a real-valued k, the complex frequency ω describes both the direction of wave propagations by ${\omega }_r$ and the wave attenuation by $-\left|{\omega }_i\right|$ as shown by the numerical simulations ensuing after pulsed excitations in [10,25]. Besides, the complex frequencies are sought as solutions to the generalized eigenvalue problems for light propagations through photonic crystals where metallic wires or filaments are periodically embedded [19,24,26,27]. In particular, positive values of ${\omega }_i$ are suitable for describing optical gains due to external pumping [9,18,26]. In general, system responses to short-duration signals such as encountered in biological or chemical sensing [16,17] are better to be described in terms of complex frequencies rather than complex wave numbers (the latter being discussed shortly). Transient responses after the application of pulsed lasers can also be adequately described via complex frequencies [15,17], although the laser pulses themselves are better described by real-valued frequencies. In addition, complex frequency is associated with the aforementioned “time-dependent” Ginzburg-Landau (TDGL) equation [7], instead of its steady-state version.

On the other hand, conventional waveguide theory has been discussed in terms of the real-valued frequency, say, with $\omega >0$ such as in the case of external illuminations by continuous-wave (cw) lasers [11–13,22]. In this approach, the complex-valued k is sought as a solution to the dispersion relation, which is obtained as follows for Eq. (5).

Hence, the solution pairs $\left(k_r,k_i\right)$ and $\left(-k_r,-k_i\right)$ refer to the right- and left-propagating waves. Of course, both type of propagating waves undergo attenuations in their respective propagation directions. Therefore, we could establish standing waves by illuminating two continuous waves of equal amplitude from opposite directions.

We plot in Fig. 5 the positive pair $\left(k_r,k_i\right)$ against ω for the complex relaxation rate γ as indicated within each panel. The range $0\le \omega <k_r$ and $k_r<\omega$ refer to the evanescent and guided waves, respectively. In all the cases, $k_r\to \infty$ and $k_i\to 0$ as $\omega \to \infty$ according to the limiting relation $k_r\to 2^{-1/2}\omega$ from Eq. (13). Moreover, $k_r$ exhibits both a local maximum and a local minimum with respect to ω in all the cases. In comparison, $k_i$ exhibits one single local maximum. In addition, the character of acoustic branch is displayed in all the cases, namely, $\omega \to 0$ as $k_r\to 0$.

 

Fig. 5 The real and imaginary parts of $k=\omega \sqrt{{\left(1+\epsilon \right)}^{-1}\epsilon }$ according to Eq. (13), where both of $\left(k_r,k_i\right)$ are positive. Three values of the complex relaxation rate are examined: (a) ${\gamma }_i=0.2-i0.2$, (b) ${\gamma }_i=0.03-i0.2$, and (c) ${\gamma }_i=0.2-i$.

Download Full Size | PDF

From comparison of Fig. 5(a) with 5(b), the local maxima get sharper for both of $\left(k_r,k_i\right)$ with a decrease in ${\gamma }_r$. In particular, we could infer from Fig. 5(b) that there is a mini-frequency-gap (where $k_r\to 0$ as indicated by the large hatched horizontal arrows) approximately near $0.8<\omega <1$ as ${\gamma }_r\to 0$ with ${\gamma }_i=-0.2$ fixed [19,24,26,27]. On the other hand, the increase in $\left|{\gamma }_i\right|$ does not alter the overall trends of $\left(k_r,k_i\right)$ with ${\gamma }_r=0.2$ fixed, as seen from Figs. 5(a) and 5(c).

In order for standing waves to be established in the case of complex relaxation rates, we require that continuous-wave beams of equal frequency and amplitude are carefully illuminated from exactly opposite directions. However, in the case of pulsed illuminations, it is quite natural that counter-propagating waves hardly lead to standing waves. This notion is valid for linear geometry. In comparison, the case with non-linear geometry is not quite simple as shown by the next case with circular geometry.

An interesting measurement setting is that a continuous-wave laser beam introduced from an optical fiber through a subwavelength-sized nanoparticle bead is incident tangentially along a certain circumference of a micro-sphere resonator [11]. It is because a whispering-gallery-mode (WGM) is established by the interaction of a rotating wave (directed into the external laser beam) and a counter-rotating wave (advancing into its opposite direction) aroused by back scattering. As a result from one-sided external excitation by a laser beam, the resulting waves are composed of two components: one running (propagating) and the other standing. Moreover, both waves are attenuated with time due to both resonator and coupling. In fact, this notion of propagating and standing waves are crucial in properly designing whispering-gallery-resonator (WGR) lasers [12,28].

6. Summary

To summarize, by examining the effects of the complex relaxation rates, we are a little closer to understanding the true nature of the lifetime and asymmetry of the surface plasmon resonances in the context of the most fundamental single-interface configuration. Extension to more complicated geometries can be carried out as well. Ways of controlling the relative magnitudes among ${\omega }_p$, ${\gamma }_r$, and ${\gamma }_i$ should depend on the choice of materials and the formulation of composite materials, which is under our current investigation.