A ug 2 01 9 Bound state and non-Markovian dynamics of a quantum emitter around a surface plasmonic nanostructure

A bound state between a quantum emitter (QE) and surface plasmon polaritons (SPPs) can be formed, where the QE is partially stabilized in its excited state. We put forward a general approach for calculating the energy level shift at a negative frequency ω, which is just the negative of the nonresonant part for the energy level shift at positive frequency −ω. We also propose an efficient formalism for obtaining the long-time value of the excited-state population without calculating the eigenfrequency of the bound state or performing a time evolution of the system, in which the probability amplitude for the excited state in the steady limit is equal to one minus the integral of the evolution spectrum over the positive frequency range. With the above two quantities obtained, we show that the non-Markovian decay dynamics in the presence of a bound state can be obtained by the method based on the Green’s function expression for the evolution operator. A general criterion for identifying the existence of a bound state is presented. These are numerically demonstrated for a QE located around a nanosphere and in a gap plasmonic nanocavity. These findings are instructive in the fields of coherent light-matter interactions.

Recently, a bound state between a QE and surface plasmon polaritons (SPPs) is predicted [61], where the QE does not decay completely to its ground state and part of its excited-state population is trapped in the steady state even in the presence of the lossy metal. Different from the previous investigation [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48], it is an open quantum system and there is no photonic band gap in the positive frequency domain. The formation of a QE-SPP bound state has been attributed to the strong field confinement which can greatly enhance their interaction. Inspired by the above works, we are interested in the conditions for * Electronic address: huang122012@163.com † Electronic address: wxyyun@163.com the existence of a bound state for a QE located around a plasmonic nanostructre.
To determine the existence or absence of a bound state, one usually calculates the long-time values of the excitedstate population by solving the non-Markovian master equation, in which a time-convolution integral should be performed. As discussed in Ref. [62], this method is time consuming. Alternatively, one can determine the existence or absence of a bound state by the resolvent operator method [1,2,52,[62][63][64][65], once the energy level shift (Lamb shift) is known. Very recently, we have proposed a general numerical method for calculating the energy level shift of a QE in an arbitrary nanostructure for positive frequency [62]. In this work, we will generalize the above method to the negative frequency domain, in which the eigenfrequency for a bound state is located. In addition, we will propose a general criterion for identifying the existence of a bound state for a QE located in an arbitrary nonostructure. We will show that the photon Green's function (GF) at zero frequency, which can be calculated by numerous methods [66][67][68][69][70][71][72][73][74][75][76][77], is suffice. This can greatly simplify the problem. Besides, we provide two different formalisms to easily determine the excitedstate population in the steady state without performing a temporal evolution.
As demonstrated in Ref. [62], the exact non-Markovian dynamics can be investigated by the method based on the Green's function expression for the evolution operator, in which there is no need of time convolution. For arbitrary nanostructure, one just needs to calculate the photon GF in the real frequency domain, in which no particular assumption about the permittivity of the material, such as the Drude model, should be made. It removes most limitations, such as the usual tight-binding or a quadratic dispersion assumption for the case of a cavity array or photonic crystal, encountered in previous analytical approaches and allows us to solve the exact decay dynamics for an open system in an almost fully numerical way. In this work, we would like to generalize the above method to study the decay dynamics of a quantum emitter coupled to surface plasmons when a bound state is formed.
The remainder of this paper is organized as follows. In Sec. II, we first present the theory and derive a general method for calculating the energy level shift in the negative frequency domain. Besides, we then provide a general criterion for identifying the existence of a bound state and two new formalisms for calculating the excitedstate population in the steady state. In Sec. III, we apply our method to a particular example where a QE is located around a single gold nanosphere. The performance of the above methods for energy level shift in negative frequency domain, condition for the existence of bound state, excited-state population in the steady state and non-Markovian dynamics are shown. Section IV is devoted to the study of the bound state and non-Markovian dynamics of a QE located in a plasmonic nanocavity. Finally, we conclude in Sec. V.

II. THEORY AND METHOD
Under the dipole and rotating-wave approximations, the total Hamiltonian for a two-level QE coupled to a common electromagnetic reservoir is [78] H = H 0 + H I , Here,f (r, ω) andf † (r, ω) are the annihilation and creation operators for the elementary excitation of the electromagnetic reservoir, respectively. d = g|d|e = dn is the matrix element for the transition dipole moment, with the unit vectorn and its strength d. The electric field vector operatorÊ (r, ω) is given bŷ Here ε(r, ω) = ε R (r, ω) + iε I (r, ω) is the frequency-dependent complex relative dielectric function in space and ε I (r, ω) is its imaginary part. I is the unit dyad and c refers to the speed of light in vacuum.
We assume initially the field is in the vacuum state and the QE is excited. In this case, the states of interest are |I = |e ⊗ |0 and |F r,ω = |g ⊗ |1 r,ω with |e (|g ) the excited (ground) state of the QE and |1 r,ω ≡f † j (r, ω) |0 . |0 is the zero photon state. The state vector of the system evolves as with the initial condition c 1 (0) = 1.
As demonstrated in Ref. [62], the dynamics can be efficiently obtained by the resolvent operator technique, in which the probability amplitude for the excited state is Here, the retarded and advanced Green's functions are G ± (ω) = lim η→0 . By using the rela- in which the spontaneous emission rate Γ(z) and the energy level shift ∆(z) are Here, θ(z) is the step function and the coupling strength is Thus, the probability amplitude for the excited state becomes with the evolution spectrum To evaluate Eq. (4), one should calculate the evolution spectrum S(ω) as well as the energy level shift ∆(ω) in the whole frequency range. For ω ≥ 0, we have demonstrated in Ref. [62] that the energy level shift can be calculated by the subtractive KK method efficiently. For the sake of completeness, we include the derivation of the method here. By using the relations −π Re G(r 0 , r 0 , s) = P +∞ −∞ ds Im G(r 0 , r 0 , s)/(ω − s) and Im G(r 0 , r 0 , −s) = − Im G(r 0 , r 0 , s) , the energy level shift [Eq. (2)] for ω ≥ 0 can be written as [79] ∆(ω) = −π Re g rr (ω) + P +∞ 0 ds Im g rr (s) ω + s .
Subtracting ∆(0) from ∆(ω), we have for ω ≥ 0 ∆(ω) = −π Re g rr (ω) + ∆c(ω), ∆c(ω) = π 2 Re g rr (0) − ω +∞ 0 ds Im g rr (s) (ω + s)s . (8) Here, we have used the relation ∆(0) = −0.5π Re g rr (0), since the second term on the right hand side in Eq. (6) is −∆(0) [see Eq. (2)]. The first term −π Re g rr (ω) is the resonant part arising from the residua at the poles and ∆c(ω) is the correction part which represents the nonresonant part of the dispersion potential. As discussed in Ref. [62], this method is useful and simplifies the numerical integrals in calculating the energy level shift, sicne there is no need to worry about the principal value and the calculation of the GF with imaginary frequency argument. In addition, it converges much more quickly than the method shown in Eq. (2). To generalize the above method to the case for a negative frequency ω < 0, we alternatively use the formalism of Eq. (2) and subtract ∆(0) from ∆(ω) in a similar way. We have for ω < 0 (9) Since the integrand in Eq. (9) decays faster than that in Eq. (2) for large frequency s, it will converge much more quickly than the method shown in Eq. (2). This will be numerically demonstrated in the next section.
Equations (7) and (9) are the main results of our methods for calculating the energy level shift of a quantum emitter for ω ≥ 0 and ω < 0, respectively. It is general and does not imply any specific configuration or system. It should be noted that the energy level shift ∆(ω) in the negetive frequecny range (ω < 0) is just the negative of the nonresonant part of the energy level shift for positive frequency −ω, i.e. −∆c(−ω). Thus, ∆(ω) in the negetive frequecny range (ω < 0) can be obtained directly once the nonresonant part ∆c(ω) is obtained for a positive frequency ω ≥ 0. In addition, it is a monotonically decreasing function and approaching zero as ω → −∞, which can be clearly seen from Eq. (2).
With the energy level shift ∆(ω) calculated, the evolution spectrum S(ω) [Eq. (5)] can be evaluated. For nonzero Γ(ω), S(ω) is of a generalized Lorentzian form But for frequency inside the photonic band gap where Γ(ω) = 0, the evolution spectrum becomes S(ω) = 1 π lim η→0+ η [ω−ω0−∆(ω)] 2 +η 2 . In this case, the evolution spectrum S(ω) is either zero or a delta function depending on the zero of the function If the solution ω b is not inside the photonic band gap, S(ω) = 0. But for ω b inside [33-36, 39, 43, 49-51, 58, 80-82], the evolution spectrum becomes a delta function for frequency around ω b , i.e. S(ω) = 1 π lim η→0+ η [ω−ω0−∆(ω)] 2 +η 2 = Zδ(ω − ω b ), in which Z is a real quantity and can be written as Different from the previous investigations where the energy level shift ∆(ω) and accordingly Z are usually analytically evaluated under a tight-binding or a quadratic dispersion assumption for the case of a cavity array or a photonic crystal, we numerically calculate both by the above methods for an open system in this work. Since there is no photonic band gap in the positive frequency range, the existence of a bound state requires that the eigenfrequency ω b less than or equal to zero, i.e. ω b ≤ 0. Furthermore, there is at most one root ω b for equation (2)]. Thus, when a bound state presents, the probability amplitude in Eq. (4) can be written as in which the second term tends to zero due to out-ofphase interference. The first term indicates a bound state, since lim t→∞ c 1 (t) = Ze −iω b t . This indicates that the quantity Z is the amplitude for the excited state in the steady limit.
If one is only interested in the long-time value of the excited-state population P a (∞) = |c 1 (∞)| 2 = Z 2 , Eq. (12) at t = 0 becomes Equation (11) and (13) are the main results of our methods for calculating the long-time values of the excited-state population P a (∞) without performing a time evolution. One just needs to calculate the photon GF in the real frequency range, which can be used to obtain the energy level shift ∆(ω) and the spontaneous emission rate Γ(ω). For the method shown in Eq. (11), we first solve the transcendental equation (10)] to find the lowest root ω b and then perform a general integral [Eq. (11)] to obtain Z. But for the method shown in Eq. (13), one just needs to perform a general integral about the evolution spectrum S(ω) [Eq. (5)] in the positive frequency range without calculating the root ω b . In the next section, we will demonstrate the performances of both methods [Eq. (11) and (13)].
As stated above, the condition for the appearance of a unique bound state is ω b ≤ 0. This requires f (0) ≥ 0, since the monotonically increasing function f (ω) in the negative frequency range approaches −∞ as ω → −∞. Explicitly, this can be written as where we have used the relation ∆(0) = −0.5π Re g rr (0). Equation (14) is a general criterion for identifying the existence of a bound state when there is no photonic band gap for the electromagnetic reservoir in the positive frequency range. It is only at zero frequency that one can judge whether a bound state exists or not from the real part for the coupling strength Re g rr (0) = Re[d * · G(r 0 , r 0 , 0) · d]/ πε 0 . Here, the photon GF G(r 0 , r 0 , 0) can be obtained by numerous methods [66][67][68][69][70][71][72][73][74][75][76][77]83]. For example, direct electrostatics methods based on solving the Poisson equation, such as the method of images, the method based on separation of variables, the finite element method (FEM), the boundary element method and so on, can be applied. Besides, other methods by extrapolating the solution of the full wave equation near zero frequency can be used [see Ref. [62,66] and references therein]. In this work, we follow the method presented in Ref. [62], in which the analytical method for a nanosphere and the FEM for an arbitrary-shaped nanostructure are adopted.
Before proceeding further, let us give a brief summary about the theory and methods. Based on the photon GF formalism, one can expresses the medium-assisted quantized electromagnetic field by the fundamental bosonic vector field and study the exact quantum decay dynamics of a QE in an open quantum system by the resolvent operator method. For any nanostructure, the photon GF can be obtained by the method shown in Ref. [62,66]. Then, the coupling strength g rr (ω) and accordingly the spontaneous emission rate Γ(ω) can be obtained directly by Eq. (3) and (1), respectively. The energy level shift ∆(ω) can be calculated by Eq. (7) for positive frequency (ω ≥ 0), in which the correction part ∆c(ω) can be calculated by Eq. (8). For ω ≤ 0, ∆(ω) can be obtained directly from Eq. (2) or Eq. (9) with ∆c(−ω) calculated by Eq. (8) for a positive frequency −ω. The non-Markovian decay dynamics can be calculated by Eq. (12). In this equation, we set Z = 0 if there is no bound state. Otherwise, Z can be obtained by Eq. (11) or Eq. (13). Both can be used to determine the long-time value of the excited-state population P e (∞) = |c 1 (∞)| 2 = Z 2 . It should be noted that there is no need to evaluate the bound state eigenfrequency ω b through Eq. (13) as compared with the method by Eq. (11). The bound state eigenfrequency ω b can be obtained by solving the tran- We have provided a general and easy way to judge whether a bound state exists or not [Eq. (14)].
In the following sections, we will first demonstrate the performances of the above methods. For this purpose, we choose a particular example where a two-level QE is located above a metal nanosphere. As shown in Fig. 1(a), a nanosphere with radius a is located at the origin. A QE at a distance h from the surface of the sphere lies on the z-axis. The metal is Gold and characterized by a complex Drude dielectric function [84] ε 2 (ω) = 1 − ω 2 p /ω(ω + iγ p ) with ω p = 1.26 × 10 16 rad/s and γ p = 1.41 × 10 14 rad/s. The background is vacuum with ε 1 = 1. For simplicity, the matrix element for the transition dipole moment is assumed to be polarized along the radial direction of the sphere d =dr and its strength is d = 24 D in this work unless otherwise specified. In section IV , we will apply the method introduced above for a QE located at the center of a plasmonic nanocavity [see Fig. 1(b)], where the permittivity ε 2 for silver is from experimental data [85]. The plasmonic nanocavity is composed of two silver nanorods with a gap distance W . Their radius and height are R = 4 nm and H = 20 nm, respectively. The transition dipole moment is also assumed to be polarized along the axis direction. For simplicity, the transition dipole moment for the QE is thought to be polarized along the z direction. ε1 and ε2 are the permittivities for air and metal, respectively.

III. PERFROMANCES OF THE ABOVE METHODS FOR BOUND STATE AND NON-MARKOVIAN DYNAMICS
In subsection A, we will first demonstrate the performance of the methods for calculating the energy level shift in the positive frequency range [Eq. (7)] and in the negative frequency range [(9)]. Then, the existence conditions of a bound state [Eq. (14)] will be shown. In subsection B, the methods for the long-time value of the excited-state population [Eq. (11) and (13)] and the non-Markovian decay dynamics [Eq. (12)] are investigated. We adopt the model shown in Fig. 1(a), where the photon GF for the nanosphere can be analytically obtained [62,66,67,83].

A. ENERGY LEVEL SHIFT AND EXISTENCE CONDITIONS OF BOUND STATE
With the photon GF calculated by the analytical method presented in Ref. [62,66,67,83], the spontaneous emission rate Γ(ω) [Eq. (1)] for h = 1nm is shown in Fig. 2(a). The same as the result shown in Ref. [62], great enhancement can be observed for ω in the range of (4 − 7eV ), which can be attributed to the localized surface plasmon resonance of the metal nanosphere.
To evaluate the correction part ∆c(ω) for ω ≥ 0 [Eq. (8)], one should calculate the real part of the coupling strength at zero frequency Re g rr (0). Although this term can be directly obtained from the analytical GF for the nanosphere, we adopt a general method by extrapolating the photon GF near zero frequency to the static case. Different from the previous linear extrapolating method [62], we adopt a linear-quadratic model, in which the the real part of the coupling strength for ω near zero is written in the form of Re g rr (ω) = α + βω + γω 2 . The results for the original data (red circle) and the extrapolating function (black solid line) are shown in Fig. 2(b). We found they agree well with Re g rr (ω) =0.0279−1.1662×10 −6 ω+0.0008ω 2 . Thus, we obtain Re g rr (0) = α = 0.0279 meV , which is consistent with the analytical result at extremely low frequency, for example, Re g rr (ω) = 0.0279 meV with ω = 10 −8 eV .
For the integral part in Eq. (8), the upper limit +∞ in the integral should be replaced by some cut-off frequency ωc to perform a numerical integral. The results are shown in Fig. 2(c) with ωc = 10 eV (the black solid line) and ωc = 200 eV (the red circle). There is no observable difference and their maximum difference is less than 0.006 meV , which means that a small cut-off frequency (ωc = 10 eV ) is enough to obtain a convergent result by Eq. (8).
For a negative frequency ω < 0, the energy level shift ∆(ω) is just the negative of the correction part, i.e. −∆c(−ω) according to Eq. (9) or can be calculated by the method shown in Eq. (2) without performing a principal value integral. The results are shown in Fig. 2(d) by Eq. (9) (the black solid line) and by Eq. (2) (red circle) with ωc = 200 eV . We find that both methods lead to almost the same results. In addition, we observe that ∆(ω) is a monotonically decreasing function in the negative frequency range which is consistent with the theory in the previous section.
To further demonstrate the accuracy, we use the results by Eq. (9) with a large cut-off frequency ωc = 200 eV as a reference and show the absolute difference (∆ error ) for the results by Eq. (2) with a small cut-off frequency ωc = 10 eV (black solid line) or a large cut-off frequency ωc = 200 eV (red dashed line) or by Eq. (9) with ωc = 10 eV (blue dash-dot line) in Fig. 2(e). We find that the method by Eq. (9) converges much more quickly than the method by Eq. (2). Figure 2(f) demonstrates the energy level shift ∆(ω) for our nanosphere model.
With Re grr(0) obtained by the above extrapolating method, one can judge whether a bound state exists or not easily according to Eq. (14). To form a bound state, the strength for the transition dipole moment d should be larger than a critical value dc = [2 ε 0 ω 0 /r · Re G(r 0 , r 0 , 0) ·r] 1/2 [derived by Eq. (14) with Eq. (3)]. Figure 3(a) shows the critical transition dipole moment dc as a function of the distance h between the QE and the surface of the nanosphere. The red circles are for the numerical results and the black solid line is for a cubic fit with dc = 11.177 + 74.563h + 66.083h 2 − 3.060h 3 .  It is found that the critical dipole strength dc increases sharply with the dipole-sphere distance h. We also vary the radius of the nanosphere a with a constant dipolesurface distance h = 1 nm. The results are shown in Fig.  3(b), where the red circles are for the numerical results and the black solid line is for a fit with dc = 138.868 + 16.943 exp(−a/4.14851)+2.472 exp(−a/31.472). Clearly, this shows that the critical dipole strength dc decreases slowly. Thus, we can conclude that dc depends heavily on the dipole-surface distance h but less on the sphere radius a.

B. LONG-TIME VALUE OF THE EXCITED-STATE POPULATION AND NON-MARKOVIAN DYNAMICS
In this subsection, we demonstrate the performance of the two methods for calculating the amplitude of the excited state in the long time limit Z [Eq. (11) and Eq. (13)] and the method for the non-Markovian decay dynamics [Eq. (12)].
To evaluate Z by Eq. (11), we first solve the transcendental equation f (ω) = ω − ω 0 − ∆(ω) = 0 [zero of Eq. (10)] to find the lowest eigenfrequency ω b . As demonstrated in Fig. 3(b), there is a critical transition dipole moment dc = 140.3 D when h = 1 nm. We find that the root ω b is negative when the dipole strength d is beyond the critical value dc. Then, we turn to the calculation of Z. Equation (11) can be evaluated once the negative eigenfrequency ω b is obtained. The results are shown in Fig. 4(b) (red circle). The other method to obtain Z is by Eq. (13). We find that the results [the black solid line in Fig. 4(b)] agree well with those by Eq. (11) (red circle). The inset is for the relative errors, which is the ratio of their absolute difference and their average value. It is less than 0.25% for the considered transition dipole moment. It should be noted that there is no need to calculate the eigenfrequency ω b for the method by Eq. (13). The non-Markovian dynamics [Eq. (12)] are demonstrated in Fig. 5. Figure 5 are for slightly below and above the critical value dc = 140.3 D, respectively. In these cases, the corresponding eigenfrequencies are positive (ω b = 0.00346 eV )and negative (ω b = −0.00352 eV ), respectively. We see that they look almost the same at short times. But for longer times (see the inset therein), a slight smaller dipole strength d (black solid line) leads to an almost perfect decay when the bound state is absent. Importantly, for the case of a little larger dipole strength d (red circle), the QE experiences a partial decay and becomes dissipationless after a long time. In addition, the steady-state population P a (∞) matches well with the results Z 2 obtained by Eq. (11) or (13). A suppressed dissipation dynamics is observed. To further demonstrate this, we choose another two different values of the dipole strength d. One is much larger than the critical value dc and the other is much smaller. The results are shown in Fig. 5(b). The red dashed line (d = 178 D leading to a negative root ω b = −0.65967 eV ) and the black solid line (d = 100 D leading to a positive root ω b = 0.65725 eV ) demonstrate a totally different decay characteristics. Different from the case shown in Fig. 5(a), the short time behaviors of the decay dynamics are different for these two values of dipole strength d. In addition, the time needed for the complete decay is much shorter for d = 100 D than that for d = 140.1 D [see the black solid lines in the insets of Fig. 5(a) and 5(b)]. This can be understood by checking the evolution spectrum S(ω) in the positive frequency domain (see Fig. 5(c) and 5(d) for d = 140.1 D and d = 100 D, respectively). We find that S(ω) around the peak value can be modeled by a perfect Lorentz function in both cases [see the insets in Fig. 5(c) and (d), where the red circles are for numerical results and the black solid lines are for their Lorentz fits]. The half widths, which represent the spontaneous emission rates, are w1 = 0.0453 meV and w2 = 4.78 meV for the above two different dipole moments, which is approximately two orders of magnitude difference. As demonstrated in Ref. [62], the decay dynamics can also be investigated by solving the non-Markovian Schrödinger equation in the time domain, which is equivalent to the method adopted in Ref. [61]. Although this method is time consuming, we compare it with the resolvent operator method [Eq. (12)] when a bound state appears. Figure 6(a) and 6(b) are the results for d = 142 D and d = 240 D, respectively. The lowest eigenfrequencies are ω b = −0.0278 eV and ω b = −1.8209 eV respectively, which are around and far away from zero. We find that the results for the excited state population P a (t) by both methods agree well in both cases. The insets are for the long-time results. It should be noted that the resolvent operator method requires much less computation time than the method by solving the Schrödinger equation in the time domain. This clearly demonstrates that Eq. (12) can be used to efficiently and accurately evaluate the non-Markovian decay dynamics when a bound state presents.
The main results of the above two sections can be summarized as follows. We have demonstrated that the real part of the coupling strength at zero frequency Reg rr (0) can be obtained by the extrapolating method. For the energy level shift, we have proven that the method by Eq. (9) is more efficient than the method by Eq. (2). Thus, the exact energy level shift can be obtained efficiently by Eq. (7) for positive frequency and by Eq. (9) for negative frequency. Besides, one can quickly judge whether a bound state exists or not by Eq. (14) with Reg rr (0) obtained. We have proven that results for the non-Markovian decay dynamics by Eq. (12) or by solving the non-Markovian Schrödinger equation (see Ref. [61,62]) are the same. Thus, the non-Markovian decay dynamics can be exactly obtained by Eq. (12), in which the parameter Z can be either calculated by Eq. (11) or (13). In spite of this, one need not to evaluate the lowest eigenfrequency ω b [negative zero for the transcendental Eq. (10)] for the method by Eq. (13) compared to the method by Eq. (11). We will use these methods to investigate the existence condition of a bound state and the non-Markovian decay dynamics when a QE is located in a gap plasmonic nanocavity.

IV. EXISTENCE CONDITIONS OF A BOUND STATE AND NON-MARKOVIAN DYNAMICS OF A QE IN A PLASMONIC NANOCAVITY
In this section, we adopt the model shown in Fig. 1(b). Different from the previous section where the Drude model for the permittivity of metal over the whole frequency range is assumed, the permittivity ε 2 for silver is from experimental data [85]. The photon GF for the nanocavity is numerically obtained by COMSOL MUL-TIPHYSICS software with the method presented in Ref. [62,66,74]. We first vary the geometric parameters of the nanocavity, such as the height H and the radius R of the nanorod, to find the critical dipole moment dc for the system to form a bound state according to Eq. (14). Then, we investigate the long-time value of the excitedstate population by Eq. (11) and (13). At the end, the non-Markovian dynamics by solving Eq. (12) are demon-strated. For simplicity, the transition frequency for the QE is also set to be 1.5 eV . The spontaneous emission rate Γ(ω) and the energy level shift ∆(ω) are shown in Fig. 7 for two different nanocavities with H = 20nm and H = 50nm. Different from the above section where both Γ(ω) and ∆(ω) take great value only in a much narrow frequency range, they are relative large over a wide frequency range. This is because the permittivity for metal is from experimental data, which is beyond the Drude model and takes into account the realistic response of metal at high frequency. From the insets in Fig. 7, we observe that the lowest plasmonic mode contributes much to the spontaneous emission rate Γ(ω) and the energy level shift ∆(ω), in which the peak position and peak value depend heavily on the rod height H.
Although Γ(ω) and ∆(ω) look very different for the above different nanocavities, the critical dipole moment dc for the system to form a bound state is nearly independent of the radius R and height H of the nanorod. This can be clearly seen from the second column of Table I. Z1 and Z2 are the amplitude for the excited state in the steady limit obtained by Eq. (11) and (13), respectively. Their differences are very small (see Table I), which imply that Eq. (11) and (13) are able to calculate Z. It should be noted that there is no need to calculate the eigenfrequency ω b by Eq. (13). In addition, we observe that Z is nearly independent of the sizes of the nanorod and the transition dipole moment d. But for the negative eigenfrequency ω b , it depends heavily on them, especially on the dipole strength d.
The decay dynamics are shown in Fig. 8. Here, the height and the radius for the nanorod are H = 20 nm and R = 4 nm, respectively. The critical dipole strength to form a bound state is dc = 92.174 D for a transition frequency ω 0 = 1.5 eV . Figure 8(a) demonstrates the excited state population P a (t) for the dipole strength d just below and above its critical value dc. At the very beginning, there is little difference, which is similar to the phenomenon shown in Fig. 5(a). The difference increases with time. For a little smaller dipole strength d = 90 D , we observe a complete decay where the excited state population P a (t) tends to zero [see the black solid line in the inset of Fig. 8(a)]. But for a slightly larger dipole strength d = 94 D, the excited state population P a (t) is partially preserved and tends to 0.64 as t → ∞.
We also consider another two different dipole moments, which are much smaller or larger than the critical value dc. In this case, the excited state population P a (t) for the above two different dipole strengths behaves differently at the beginning. This implies that the decay dynamics is much affected by the dipole strength at the beginning. For the long time performance, a complete decay and a partial limited decay can be observed depending on the dipole moment [see the insets in Fig. 8(b)]. For d = 130 D, P a (∞) = 0.593, which is around Z 2 = 0.604 for d = 120 D (see the last two columns in Table I). This is also consistent with the argument that the dipole strength d takes little effect on the excited state population in the long time limit.
The above phenomena can also exist for another plasmonic nanocavity composed of nanorods with different geometric parameter. Figure 8(c) and 8(d) demonstrate the results for a higher nanorod H = 80 nm with the same radius R = 4 nm. We find that the results shown in Fig. 8(c) looks similar to that in Fig. 8(a). Differently, the time for the system to reach its steady state is much smaller for Fig. 8(d) than that in Fig. 8(b).

V. CONCLUSIONS
In summary, we have presented an efficient numerical method [Eq. (14)] to determine the existence or absence of a bound state for a QE coupled to surface plasmon polaritons without calculating the time evolution of the system. We have found that the critical dipole strength dc is heavily dependent on the dipole-sphere distance h but much less dependent on the sphere-radius a for the QE located around a nanosphere. For the QE at the center of a plasmonic nanocavity composed of two nanorods, we have found that dc is nearly independent of the radius and height of the rod. For a bound state presented, we have proposed two different methods to determine the long-time value for the excited state population P a (∞) = Z 2 [Eq. (11) and (13)]. To evaluate Eq. (11), we have proposed a new formalism to calculate the energy level shift for a negative frequency [Eq. (9)], which is more efficient than the method by Eq (2). This is helpful to determine the exact eigenfrequency ω b [zero for Eq. (10)] for the bound state, which is an important parameter to evaluate Z by Eq. (11). For a QE located The results of dc, ω b and Z with different heights H and radius R of the nanorod. dc is the critical strength for the dipole moment obtained by Eq. (14). ω b is the negative eigenfrequency for f (ω) = ω − ω0 − ∆(ω) = 0 [zero of Eq. (10)]. Z1 and Z2 are the amplitude for the excited state in the steady limit obtained by Eq. (11) and (13), respectively. Their differences are very small. We observe that dc and Z are nearly independent of the height H and radius R of the nanorod. But ω b depends heavily on them, especially on the transition dipole moment d.  around a nanosphere or in a plasmonic nanocavity, we have shown that both methods [Eq. (11) and (13)] lead to the same results. In addition, we have found that Z is less affected by the dipole strength d with d ≥ dc, but the lowest eigenfrequency ω b is much dependent on it. By comparing with the time-domain method via solving the non-Markovian Schrödinger equation, we have shown that the non-Markovian decay dynamics in the presence of a bound state can be exactly and efficiently obtained by the method based on the Green's function expression for the evolution operator. It is found that the excited state population P a (∞) matches well with Z 2 obtained by Eq. (11) or (13). In addition, we have found that P a (∞) is nearly independent of the sizes of the nanorod and the transition dipole strength d when d is larger than the critical value dc, which is consistent with the prediction by Eq. (11) and (13).