1. Introduction

Coherent interaction between a quantum emitter (QE), such as atom, molecule or quantum dot, and the quantized electromagnetic field lies at the heart of quantum optics [1–8] . An incredibly wide range of observable phenomena, such as ultrafast single-photon switches [9–11], trapping atoms by vacuum forces [12–15], photon blockade [16,17], single-molecule sensing [18], single-atom lasers [19,20], enhanced and inhibited spontaneous emission [21–29], giant Lamb shift [30–32], quantum nonlinear optics [33,34], bound state [35–58], etc, have been predicted and demonstrated. Among all, bound state is one particular example, where a photonic excitation is confined to the vicinity of the QE and a discrete eigenstate is formed inside the environmental photonic band gap. The quantum coherence, which plays a prominent role in the fields of quantum physics and technology [59–67], can be protected from environmental noise.

Recently, a bound state between a QE and the surface plasmon polaritons (SPPs) has been predicted [55–57], where an initially excited QE does not decay completely to its ground state and remains in its excited state with a finite population in the steady state even in the presence of the lossy metal. Different from the previous investigation [35–50,68–70], the energy continuum for the photonic environment extends from $0$ to $+\infty$ and there is no photonic band gap in this range. Thus, a bound state will take a discrete negative eigenvalue which is a new quantum phase [38,40]. As stated in Ref. [38], the formation of the bound state actually corresponds to a quantum phase transition where eigenenergy of the ground state changes from zero to a minus. For QE-SPP bound state, the strong field confinement takes great effect. Inspired by the above works, we are interested in the conditions for the existence of a QE-SPP bound state. We will propose some numerical methods for investigating this problem and demonstrate its effect on the non-Markovian decay dynamics, especially for the system undergoing a quantum phase transition around the critical condition.

To determine the existence or absence of a bound state, one usually calculates the long-time values of the excited-state population by solving the non-Markovian master equation, in which a time-convolution integral should be performed. However, this method is time consuming and unfavourable for understanding the underlying physical mechanism. Alternatively, one can determine the existence or absence of a bound state by the resolvent operator method without evaluating the time evolution of the system [1,2,54,55,71–74]. It should be noted that this method depends on the analytical part of the self-energy for frequency below the lower threshold. Very recently, we have proposed an efficient numerical method for calculating the self-energy of a QE in an arbitrary nanostructure for positive frequency [74]. However, it can not be applied for a negative frequency. In this work, we will generalize the above method to the negative frequency domain, in which the eigenfrequency for a bound state is located. We will show that this method is much more efficient than the usual method [see Eq. (13) in this work]. In addition, we will propose a general criterion for identifying the existence of a bound state for a QE located in a plasmonic nanostructure with arbitrary shape. We will show that the photon Green’s function (GF) at zero frequency, which can be calculated by numerous methods [75–86], is sufficient. This helps us to efficiently determine whether and how to form a bound state. It also implies that the bound state condition is independent of the actual kinds of nobel metal-material, since the perfect electric conductor boundary can be assumed for the metal surface in the electrostatic case. We will prove this and also demonstrate that the photon GF at zero frequency is nearly independent of the dimension of the plasmonic nanostructure. Besides, we provide a formalism to easily determine the excited-state population in the steady state without performing a temporal evolution and without calculating the eigenvalue for the bound state, in which one only needs to perform a simple integral of the evolution spectrum.

As demonstrated in Ref. [74] without a bound state, the exact decay 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 in an almost fully numerical way. In this work, we would like to generalize the above method to study the decay dynamics of a QE coupled to surface plasmons when a bound state is formed. We will show that the long-time value of the excited-state population shows a transition from zero to a limited value by slightly altering some parameters (the transition dipole moment in this work) of the system.

The remainder of this paper is organized as follows. In Sec. 2, we first present the theory and derive a general method for calculating the self-energy in the negative frequency domain. Besides, we then provide a general criterion for identifying the existence of a bound state and a new formalism for calculating the excited-state population in the steady state. In Sec. 3, we apply our method to a particular example where a QE is located around a single gold nanosphere with the photon GF obtained analytically. The performance of the above methods for the self-energy in negative frequency domain, the condition for the existence of a bound state, the excited-state population in the steady state and the decay dynamics are shown. Section 4 is devoted to the study of the above methods for a QE located in a general plasmonic nanocavity, where the permittivity of metal is beyond Drude model and the photon GF should be calculated numerically. Finally, we conclude in Sec. 5.

2. 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 [87]

For the above Hamiltonian, the total exciton number operator $\hat {N}=|e\rangle \langle e| + \int _{0}^{+\infty }d\omega \hat {\mathbf {f}}^{\dagger }(\mathbf {r},\omega ) \cdot \hat {\mathbf {f}}(\mathbf {r},\omega )$ is conserved [55]. Thus, the system will remain in the one exciton subspace if the system is initially prepared in a state with only the QE excited. The states of relevance in one single exciton subspace are denoted by $|I\rangle =|e\rangle \otimes |0\rangle$ and $|F_{r,\omega } \rangle =|g\rangle \otimes |1_{r,\omega }\rangle$ with $|1_{r,\omega }\rangle \equiv \hat {\mathbf {f}}_{j}^{\dagger }\left ( \mathbf {r},\omega \right ) |0\rangle$. Here, $|0\rangle$ is the zero photon state.

For a system with an initial state $|I\rangle =|e\rangle \otimes |0\rangle$, its dynamics reads

The probability amplitude for the QE remaining in its excited state $c_{1}(t)$ can be written by the matrix element of the evolution operator

For $\omega \geq 0$, we have demonstrated in Ref. [74] that $\Delta (\omega )$ can be calculated by the subtractive Kramers-Kronig (KK) method efficiently. For the sake of completeness, we include the derivation of the method here. By using the relations $-\pi \operatorname {Re}\mathbf {G}({\mathbf {r}_{0}} ,\mathbf {r}_{0},\omega )=\mathbb {P}\int _{-\infty }^{+\infty }ds\operatorname {Im} \mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0},s)/(\omega -s)$ and $\operatorname {Im}\mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0} ,-\omega )=-\operatorname {Im}\mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0},\omega )$ , $\Delta (\omega )$ [Eq. (13)] for $\omega \geq 0$ can be written as [90]

However, the above method can not be applied for $\omega <0$, since the denominator in Eq. (19) may be zero. To generalize the above method to the case for $\omega <0$, we alternatively use the formalism of Eq. (13) and subtract $\Delta (0)$ from $\Delta (\omega )$. We have for $\omega <0$

Equations (18) and (20) are the main results of our methods for calculating the analytic part of the self-energy $\Delta (\omega )$ for $\omega \geq 0$ and $\omega <0$, respectively. It is general and does not imply any specific configuration or system. It should be noted that $\Delta (\omega )$ for $\omega <0$ is just the negative of the nonresonant part of $\Delta (s)$ for $s=-\omega$, i.e. $-\Delta _c(-\omega )$ [see Eq. (20)]. Thus, $\Delta (\omega )$ for $\omega <0$ can be obtained directly once the nonresonant part $\Delta _c(\omega )$ is obtained for $\omega \geq 0$, and it is a monotonically decreasing function and approaching zero as $\omega \rightarrow -\infty$ [see Eq. (13)].

With $\Delta (\omega )$ calculated in the whole frequency range, the evolution spectrum $S(\omega )$ [Eq. (16)] can be evaluated. For nonzero $\Gamma (\omega )$, $S(\omega )$ is of a generalized Lorentzian form $S(\omega )=\frac {1}{\pi }\frac {\Gamma (\omega )/2}{[\omega -\omega _{0} -\Delta (\omega )]^{2}+[\Gamma (\omega )/2]^{2}}$. But for frequency inside the photonic band gap where $\Gamma (\omega )=0$, the evolution spectrum becomes $S(\omega )=\frac {1}{\pi }\lim _{\eta \rightarrow 0_{+}}\frac {\eta }{[\omega -\omega _{0}-\Delta (\omega )]^{2} +\eta ^{2}}$, and it is either zero or a delta function depending on the zero of the function

Different from the previous system with a limited photonic bandwidth or with a photonic bandgap, there is no photonic bandgap for the plasmonic nanostructure in the positive frequency range $\omega >0$. Thus, the existence of a bound state requires that the eigenfrequency $\omega _{b}$ should be less than or equal to zero, i.e. $\omega _{b}\leq 0$. Furthermore, there is at most one root $\omega _{b}$ for equation $f(\omega )=0$ in the range $\omega \leq 0$, since $f(\omega )=\omega -\omega _{0}-\Delta (\omega )$ is a monotonically increasing function [see Eq. (13)]. Thus, when a bound state is present, the probability amplitude in Eq. (15) can be written as

If one is only interested in the long-time value of the excited-state population $P_{a}(\infty )=\left \vert c_{1}(\infty )\right \vert ^{2}=Z^{2}$, there is another way to determine it besides by Eq. (22). Actually, Eq. (23) at $t=0$ becomes

Compared to the usual method shown in Eq. (22), Eq. (24) provides an efficient numerical method for calculating the long-time value of the excited-state population $P_{a}(\infty )$ without solving the transcendental equation $f(\omega )=\omega -\omega _{0}-\Delta (\omega )=0$ [Eq. (21)] to obtain the lowest root $\omega _{b}$. One just needs to perform a general integral about the evolution spectrum $S(\omega )$ [Eq. (16)] over the positive frequency range. In this work, we will demonstrate the performances of both methods [Eqs. (22) and (24)].

As stated above, the condition for the appearance of a unique bound state is $\omega _{b}\leq 0$. This requires $f(0)\geq 0$, since the monotonically increasing function $f(\omega )$ for $\omega \leq 0$ approaches $-\infty$ as $\omega \rightarrow -\infty$. Explicitly, this can be written as

Equation (25) is a general criterion for identifying the existence of a bound state. 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 $\operatorname {Re} g(0)=\operatorname {Re}[\mathbf {d}^{\ast }\cdot \mathbf {G} ({\mathbf {r}_{0}},\mathbf {r}_{0},0)\cdot \mathbf {d}]/\hbar \pi \varepsilon _{0}$. It should be noted that the photon GF $\mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0},0)$ is nearly independent on the exactly kind of metal, since noble metal such as gold and silver at zero frequency is usually viewed as perfect conductor.

The photon GF $\mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0},0)$ can be obtained by numerous methods [75–86,96]. For example, direct 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 [74,75]. In this work, we follow the latter.

Before proceeding further, let us give a brief summary about the theory and methods. Based on the photon GF formalism, one can express the medium-assisted quantized electromagnetic field by the fundamental bosonic vector field and study the exact quantum decay dynamics of an initially excited QE by the resolvent operator method. For any nanostructure, the photon GF can be obtained by the method shown in Ref. [74,75]. Then, the coupling strength $g(\omega )$ and accordingly $\Gamma (\omega )$ can be easily evaluated through Eqs. (14) and (12), respectively. For $\omega \ge 0$, $\Delta (\omega )$ can be evaluated by Eq. (18), in which the correction part $\Delta _c(\omega )$ can be obtained by Eq. (19). For $\omega \leq 0$, $\Delta (\omega )$ can be obtained by a more efficient method Eq. (20) than by the direct method Eq. (13). The non-Markovian decay dynamics can be obtained by Eq. (23). Similar to Ref. [74] for a system without a bound state, one can easily set $Z=0$. Otherwise, $Z$ can be evaluated by Eq. (22) or Eq. (24). Although both can determine the long-time value of the excited-state population $P_e(\infty ) = |c_1(\infty )|^2= Z^2$, the new form of Eq. (24) avoids the need to explicitly calculate the bound state eigenfrequency $\omega _{b}$ via solving the transcendental equation $f(\omega )=\omega -\omega _{0}-\Delta (\omega )=0$ [Eq. (21)]. We have provided a general and easy way to judge whether a bound state exists or not [Eq. (25)], in which only the photon GF at zero frequency $\mathbf {G}({\mathbf {r}_{0}},\mathbf {r}_{0},0)$ is sufficient.

In the following sections, we will first demonstrate the performances of the above methods. For this purpose, we first choose a particular example where a two-level QE is located above a metal nanosphere. The photon GF can be obtained analytically. Although the QE-nanosphere system has been widely investigated previously, we are in a new regime where a bound state may exist. 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. Without specification, we set $a=20\,\textrm {nm}$ and $h=1\,\textrm {nm}$. The metal is Gold and characterized by a complex Drude dielectric function [97] $\varepsilon _{2} (\omega )=1-\omega _{p}^{2}/\omega (\omega +i\gamma _{p})$ with $\omega _{p}=1.26\times 10^{16} \, \textrm {\textrm {rad/s}}$ and $\gamma _{p}=1.41\times 10^{14}\,\textrm {rad/s}$. The background is vacuum with $\varepsilon _{1}=1$. For simplicity, the matrix element for the transition dipole moment is assumed to be polarized along the radial direction of the sphere $\mathbf {d=}d \,\mathbf {\hat {r}}$ and its strength is $d=24\,\textrm {D}$ in this work unless otherwise specified. In section 4, we will apply the methods introduced above for a QE located at the center of a plasmonic nanocavity [see Fig. 1(b)], which is a more realistic and general case. The permittivity $\varepsilon _{2}$ for silver is from experimental data [98], which is beyond the Drude model. In addition, the photon GF for such an arbitrary-shaped nanostructure should be obtained by numerical method. The plasmonic nanocavity is composed of two silver nanorods with a gap distance $W=2\,\textrm {nm}$. Without specification, their radius and height are $R=4\,\textrm {nm}$ and $H=20\,\textrm {nm}$, respectively. The transition dipole moment is also assumed to be polarized along the axis direction.

 

Fig. 1. Schematic diagrams. (a) A QE is located around a gold nanosphere with radius $a$. The distance between the emitter and the surface of sphere is $h$. (b) A QE is located at the center of two nanorods with radius $R$ and height $H$. For simplicity, the transition dipole moment for the QE is thought to be polarized along the $z$ direction. $\varepsilon _{1}$ and $\varepsilon _{2}$ are the permittivities for air and metal, respectively.

Download Full Size | PDF

3. Performances of the above methods for bound state and non-Markovian dynamics

In subsection 3.1, we will first demonstrate the performance of the methods for calculating the analytical part for the self-energy $\Delta (\omega )$ for $\omega \geq 0$ [Eq. (18)] and for $\omega <0$ [Eq. (20)]. Then, the existence condition of a bound state [Eq. (25)] will be shown. In subsection 3.2, the methods for the long-time value of the excited-state population [Eqs. (22) and (24)] and the non-Markovian decay dynamics [Eq. (23)] are investigated. We adopt the model shown in Fig. 1(a), where the photon GF for the nanosphere can be analytically obtained [74–76,96].

3.1 Self-energy and existence condition of bound state

With the photon GF calculated by the analytical method presented in Ref. [74–76,96], the nonanalytical part for the self-energy $\Gamma (\omega )$ [Eq. (12)] with $h=1 \textrm {nm}$ is shown in Fig. 2(a). The same as the result shown in Ref. [74,75], great enhancement can be observed for $\omega$ in the range of $(4-7 \textrm {eV})$, which can be attributed to the localized surface plasmon resonance of the metal nanosphere.

 

Fig. 2. $\Gamma (\omega )$ and $\Delta (\omega )$ for a QE around a nanosphere. (a) $\Gamma (\omega )$. (b) $\operatorname {Re}g_{rr} (\omega )$ near zero frequency (red circle) and the fitted results (black solid line). (c) The nonresonant part $\Delta _c(\omega )$ for $\omega \geq 0$ calculated by Eq. (19) with $\omega _c=10 \, \textrm {eV}$ (black solid line) and $\omega _c=200 \, \textrm {eV}$ (red circle). (d) $\Delta (\omega )$ for $\omega <0$ calculated by Eq. (13) (black solid line) and by Eq. (20) (red circle) with $\omega _c=200 \, \textrm {eV}$. (e)The absolute difference between $\Delta (\omega )$ obtained by Eq. (13) (black solid line for $\omega _c=10 \, \textrm {eV}$ and red dash line for $\omega _c=200 \, \textrm {eV}$) or by Eq. (20) (blue dash-dot line for $\omega _c=10 \, \textrm {eV}$) and the results by Eq. (20) with $\omega _c=200 \, \textrm {eV}$. (f) $\Delta (\omega )$ in the range $-10\,\textrm {eV} \leq \omega \leq 10\,\textrm {eV}$.

Download Full Size | PDF

To evaluate the correction part $\Delta _c(\omega )$ for $\omega \geq 0$ [Eq. (19)], one should calculate the real part of the coupling strength at zero frequency $\operatorname {Re}g(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 [74], we adopt a linear-quadratic model, in which the the real part of the coupling strength for $\omega$ near zero is written in the form of $\operatorname {Re}g(\omega )=\alpha +\beta \omega +\gamma \omega ^2$. The results for the original data (red circle) and the extrapolating function (black solid line) are shown in Fig. 2(b). We find that they agree well with a fit function $\operatorname {Re}g(\omega )=$0.0279$-$1.1662$\times 10^{-6}\omega +$0.0008$\omega ^2$. Thus, we obtain $\operatorname {Re}g(0)=\alpha =0.0279\,\textrm {eV}$, which is consistent with the analytical result at extremely low frequency, for example, $\operatorname {Re}g(\omega )=0.0279\,\textrm {eV}$ with $\omega =10^{-8}\,\textrm {eV}$. This is also consistent with the static result $0.0278885\,\textrm {eV}$ by setting the perfect electric conductor boundary for the nanosphere, which means that $\operatorname {Re}g(0)$ is actually independent on the exactly kind of noble metal.

For the integral part in Eq. (19), the upper limit $+\infty$ should be replaced by some cutoff frequency $\omega _c$. The results are shown in Fig. 2(c) with $\omega _c=10\,\textrm {eV}$ (the black solid line) and $\omega _c=200\,\textrm {eV}$ (the red circle). There is no observable difference and their maximum difference is less than $6\times 10^{-6}\,\textrm {eV}$, which means that a small cutoff frequency ($\omega _c=10\,\textrm {eV}$) is enough to obtain a convergent result by Eq. (19).

According to our new formalism Eq. (20), $\Delta (\omega )$ for $\omega < 0$ is just the negative of the correction part, i.e. $-\Delta _c(-\omega )$. To demonstrate its accuracy, we compare the results with those by the usual method through Eq. (13). The results are shown in Fig. 2(d) by Eq. (20) (red circle) and by Eq. (13) (black solid line) with $\omega _c=200 \, \textrm {eV}$. We find that both methods lead to almost the same results. In addition, we observe that $\Delta (\omega )$ is a monotonically decreasing function for $\omega \leq 0$ which is consistent with the theory in the previous section.

To further demonstrate the efficiency of our new formalism [Eq. (20)], we use the results by Eq. (20) with a large cutoff frequency $\omega _c=200 \, \textrm {eV}$ as a reference and show the absolute difference ($\Delta _{error}$) for the results by Eq. (13) with a small cutoff frequency $\omega _c=10 \, \textrm {eV}$ (black solid line) or a large cutoff frequency $\omega _c=200 \, \textrm {eV}$ (red dash line) or by Eq. (20) with $\omega _c=10 \, \textrm {eV}$ (blue dash-dot line) in Fig. 2(e). We find that the method by Eq. (20) converges much more quickly than the method by Eq. (13). The fast convergence property of the method by Eq. (20) will be more evident when we take into account the response of metal at high frequency. In Sec. 4, we will show that the method by Eq. (13) does not converge even for $\omega _c=1000 \, \textrm {eV}$ while the new method by Eq. (20) converges for $\omega _c=70 \, \textrm {eV}$. Thus, in the following calculation, we adopt the method by Eq. (20). Figure 2(f) demonstrates the final results for $\Delta (\omega )$ in our nanosphere model for $\omega$ in the range of $(-10-10 \textrm {eV})$.

With $\operatorname {Re}g(0)$ obtained, one can judge whether a bound state exists or not easily according to Eq. (25). To form a bound state, the strength for the transition dipole moment $d$ should be larger than a critical value $d_c=[2\hbar \varepsilon _{0}\omega _{0}/\mathbf {\hat {r}}\cdot \operatorname {Re}\mathbf {G}({\mathbf {r}_{0} },\mathbf {r}_{0},0)\cdot \mathbf {\hat {r}}]^{1/2}$ [derived by Eq. (25) together with Eq. (14)]. Figure 3(a) shows the critical transition dipole moment $d_c$ 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 $d_c=11.177+74.563h+66.083 h^2-3.060 h^3$. It is found that the critical dipole strength $d_c$ increases sharply with the dipole-sphere distance $h$. Since there is no natural or artificial QE with the dipole moment as large as several hundred Debye, the dipole-surface distance which is small enough helps to form a bound state. We also vary the radius of the nanosphere $a$ with a constant dipole-surface distance $h=1 \,\textrm {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 $d_c=138.868+16.943 e^{-a/4.14851}+2.472 e^{-a/31.472}$. Clearly, this shows that the critical dipole strength $d_c$ decreases slowly. Thus, we can conclude that $d_c$ depends heavily on the dipole-surface distance $h$ but less on the sphere radius $a$.

 

Fig. 3. The critical transition dipole moment $d_c$ as a function of QE-surface distance $h$ and sphere radius $a$. (a) $d_c$ with $a=20\,\textrm {nm}$. (b) $d_c$ with $h=1\,\textrm {nm}$. The red circles are for the numerical results and the black lines are for fits.

Download Full Size | PDF

3.2 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. (22) and Eq. (24)] and the method for the non-Markovian decay dynamics [Eq. (23)] when a bound state exists.

To evaluate $Z$ by Eq. (22), we first solve the transcendental equation $f(\omega )=\omega -\omega _{0}-\Delta (\omega )=0$ [zero of Eq. (21)] to find the lowest eigenfrequency $\omega _b$. As demonstrated in Fig. 3, there is a critical transition dipole moment $d_c=140.3\,\textrm {D}$ when $h=1\,\textrm {nm}$ and $a=20\,\textrm {nm}$. We find that the root $\omega _b$ is negative when the dipole strength $d$ is beyond the critical value $d_c$. For example, $\omega _b=-0.00352\,\textrm {eV}$ when $d=140.5\,\textrm {D}$. But for dipole moment below the critical one, it is positive. For example, we have $\omega _b=0.00346\,\textrm {eV}$ with $d=140.1\,\textrm {D}$. Figure 4(a) demonstrates the lowest root $\omega _b$ as a function of the dipole strength $d$. It is found that the larger the dipole strength $d$ is, the lower the eigenfrequency $\omega _b$ is.

 

Fig. 4. The lowest eigenfrequency $\omega _b$ and $Z$ as a function of the dipole strength $d$ with $a=20\,\textrm {nm}$. (a) $\omega _b$. (b) $Z$. The inset is for the relative difference of $Z$ obtained by Eqs. (22) and (24).

Download Full Size | PDF

Then, we turn to the calculation of $Z$. Equation (22) can be evaluated once the negative eigenfrequency $\omega _b$ is obtained. The results are shown in Fig. 4(b) (black line). The other method to obtain $Z$ is by Eq. (24), which avoids the calculating of the negative eigenfrequency $\omega _b$. We find that the results by Eq. (24) (red circle) agree well with those by Eq. (22) (black line). The inset is for the relative errors, which is the ratio of their absolute difference to 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 $\omega _b$ for our new method [Eq. (24)].

The non-Markovian dynamics [Eq. (23)] are demonstrated in Fig. 5. Figure 5(a) shows the results when the dipole moments are just below (black line for $d=140.1\,\textrm {D}$) and above (red circle for $d=140.5\,\textrm {D}$) the critical value $d_c=140.3\,\textrm {D}$. For these two dipole moments, the corresponding eigenfrequencies are positive ($\omega _b=0.00346\,\textrm {eV}$) and negative ($\omega _b=-0.00352\,\textrm {eV}$), respectively. We see that they look almost the same at short times. But for longer times (see the inset therein), a slightly smaller dipole strength $d$ (black solid line) leads to a complete decay, which shows the absence of a bound state. 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. This clearly demonstrates that the decay dynamics shows a transition from complete decay to suppressed dissipation by fine tuning the transition dipole moment $d$ [38,40]. In addition, we find that the steady-state population $P_a(\infty )$ obtained here by the dynamics method [Eq. (23)] matches well with the results $Z^2$ obtained by Eqs. (22) or (24).

 

Fig. 5. The excited state population $P_{a}(t)=\left \vert c_{1}(t)\right \vert ^{2}$ and the evolution spectrum $S(\omega )$. (a) $P_{a}(t)$ for $d$ around the critical value $d_c=140.3\,\textrm {D}$. The red circles are for a little larger transition dipole moment $d=140.5\,\textrm {D}$ and the black solid line is for a little smaller dipole moment $d=140.1\,\textrm {D}$. (b) The same as (a) but for a much larger ($d=178\,\textrm {D}$ red dash line) and a much smaller transition dipole moment ($d=100\,\textrm {D}$ black solid line). (c) $S(\omega )$ for $d=140.1\,\textrm {D}$. (d) $S(\omega )$ for $d=100\,\textrm {D}$. The insets in (c) and (d) show the Lorentz fits with widths $w1=0.0453\,\textrm {meV}$ and $w2=4.78\,\textrm {meV}$ for frequency around their first peaks, respectively.

Download Full Size | PDF

To further demonstrate this, we choose another two different values of dipole strength $d$. One is much larger than the critical value $d_c$ and the other is much smaller. The results are shown in Fig. 5(b). The red dashed line ($d=178\,\textrm {D}$ leading to a negative root $\omega _b=-0.65967\,\textrm {eV}$) and the black solid line ($d=100\,\textrm {D}$ leading to a positive root $\omega _b=0.65725\,\textrm {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 for the QE to decay completely to its ground state is much shorter for $d=100\,\textrm {D}$ than that for $d=140.1\,\textrm {D}$ [see the black solid lines in the insets of Figs. 5(a) and 5(b)]. This can be understood by checking the evolution spectrum $S(\omega )$ [see Figs. 5(c) and 5(d) for $d=140.1\,\textrm {D}$ and $d=100\,\textrm {D}$, respectively]. We find that there are two peaks for both, and for frequency away from each peak, $S(\omega )$ is very small. But for frequency around each peak, $S(\omega )$ is large and can be modeled by a perfect Lorentz function. For example, the insets in Figs. 5(c) and 5(d) are for the magnifications around the first peaks, where the red circles are for numerical results and the black solid lines are for their Lorentz fits. Thus, we can divide the evolution spectrum into the sum of two Lorentz functions and a small correction part, i.e. $S(\omega )=S_1(\omega )+S_2(\omega )+S_{diff}(\omega )$ with $S_i(\omega )=\frac {A_i}{\pi }\frac {\Gamma _i/2}{(\omega -\omega _i)^2+(\Gamma _i/2)^2}$ ($i=1,2$).

For the spectrum shown in Fig. 5(c) ($d=140.1\,\textrm {D}$), we have $A_1=0.79$, $\omega _1=0.00346\,\textrm {eV}$, $\Gamma _1=0.0453\,\textrm {meV}$ and $A_2=0.205$, $\omega _2=7.289\,\textrm {eV}$, $\Gamma _2=0.073\,\textrm {eV}$. For both $S_1(\omega )$ and $S_2(\omega )$, the $R_i$ value, which is the ratio of $\omega _i$ to the full width at half maximum $\Gamma _i$, i.e. $R_i=\omega _i/\Gamma _i$, are $R_1=76.38$ and $R_2=100.26$, respectively. Both are much larger than the critical value $1/2$, which means that their Fourier transformations deviate less from the exponentially damped oscillation [99,100], i.e. $\int _{0}^{+\infty }S_i(\omega )e^{-i\omega t}d\omega \simeq A_ie^{-i\omega _it}e^{-\Gamma _it/2}$. Thus, $\Gamma _1$ and $\Gamma _2$ represent the damped rates, and accordingly we can define two time scales $T_1=1/\Gamma _1=14.5\,\textrm {ps}$ and $T_2=1/\Gamma _2=0.009\,\textrm {ps}$ where $T_1\gg T_2$ . The correction part $S_{diff}(\omega )$ is small and contributes much less to the dynamics than those from $S_1(\omega )$ and $S_2(\omega )$ when $t<10T_1$. Thus, the survival probability can be approximated by

We have numerically proved that this expression works well for $P_a(t)>10^{-5}$. For $t<T_2$, the two contributions from $S_1(\omega )$ and $S_2(\omega )$ interference which leads to the oscillation in $P_a(t)$ [see the last term in Eq. (26)]. But for $T_2\ll t<10T_1$, the contribution from $S_2(\omega )$ is small and can be safely ignored. Only the first term in Eq. (26) survivals and it shows a single exponential decay. For times much larger than $10T_1$, all the three terms in Eq. (26) disappear and $P_a(t)$ is dominated mainly by $S_{diff}(\omega )$. However, in this range, $P_a(t)$ is less than $10^{-6}$ and decays to zero due to out-of-phase interference, for $t\rightarrow +\infty$, which is so small and out of the interest of this work. Thus, the decay time is mainly determined by $T_1=1/\Gamma _1$.

The above analysis can also be applied to the case for $d=100\,\textrm {D}$ [the black solid line in Fig. 5(b)]. The spectrum is shown in Fig. 5(d). Similarly, we have $A_1=0.858$, $\omega _1=0.657\,\textrm {eV}$, $\Gamma _1=4.78\,\textrm {meV}$ and $A_2=0.14$, $\omega _2=6.634\,\textrm {eV}$, $\Gamma _2=0.079\,\textrm {eV}$. In this case, the two time scales are $T_1=1/\Gamma _1=0.138\,\textrm {ps}$ and $T_2=1/\Gamma _2=0.008\,\textrm {ps}$ with $T_1\gg T_2$ , respectively. The contribution from the correction part is also very small before $10T_1$. Thus, in this case, the decay time is mainly determined by $T_1=1/\Gamma _1$. The dramatic differences for the two widths shown in the insets of Figs. 5(c) ( $0.0453\,\textrm {meV}$) and 5(d) ($4.78\,\textrm {meV}$) lead to the approximately two orders of magnitude difference for the lifetime of the excited state.

Equation (26) is also helpful to understand the results shown in Fig. 5(a) where the decay dynamics look similar at the initial times for the two dipole moments. Since there is little difference for the two matrix elements of the dipole moment investigated ($d=140.5\,\textrm {D}$ and $d=140.1\,\textrm {D}$), the evolution spectrums $S(\omega )$ for $\omega >0$ look similar except for frequency around the first peak, i.e. $S_1(\omega )$. For a slightly smaller dipole strength $d=140.1\,\textrm {D}$, $S_1(\omega )=\frac {0.79}{\pi }\frac {0.00002265}{(\omega -0.00346)^2+(0.00002265)^2}$. However, when the dipole strength $d=140.5\,\textrm {D}$ is immediately above the critical value, there is a quantum phase transition. $S_1(\omega )$ becomes a delta form, i.e. $S_1(\omega )=Z\delta (\omega -\omega _{b})=0.784\delta (\omega +0.00352)$. We have $Z=0.784$ which is nearly equal to $A_1=0.79$ . Thus, the parameters in Eq. (26) are similar for both dipole moments except for slight differences for $\omega _1$ and $\Gamma _1$, which change from $0.00346\,\textrm {eV}$ to $-0.00352\,\textrm {eV}$ and from $0.0453\,\textrm {meV}$ to $0$, respectively. Since the oscillation period is mainly determined by $\omega _2-\omega _1\simeq \omega _2$, the initial decay dynamics look similar for the two dipole moments.

As demonstrated in Ref. [74], the decay dynamics without a bound state 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. [55]. Although this method is time consuming, we compare it with the resolvent operator method [Eq. (23)] when a bound state appears. Figures 6(a) and 6(b) are the results for $d=142\,\textrm {D}$ and $d=240\,\textrm {D}$, respectively. The lowest eigenfrequencies are $\omega _b=-0.0278\,\textrm {eV}$ and $\omega _b=-1.8209\,\textrm {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. (23) can be used to efficiently and accurately evaluate the non-Markovian decay dynamics when a bound state is present.

 

Fig. 6. Performance of the method [Eq. (23)] for calculating the excited state population $P_{a}(t)=\left \vert c_{1}(t)\right \vert ^{2}$ when a bound state exists with $d$ around and away from the critical value $d_c=140.3\,\textrm {D}$. (a) $d=142\,\textrm {D}$. (b) $d=240\,\textrm {D}$. The black solid lines are the results by the method based on the Green’s function expression for the evolution operator [Eq. (23)]. The red circles are for the results by solving the non-Markovian Schrödinger equation (Volterra integral equation of the second kind) in the time domain (Eq. (3) in Ref. [74]). The results agree well, which means that Eq. (23) is able to obtain the exact non-Markovian dynamic when a bound state is present.

Download Full Size | PDF

The main results for the above two subsections can be summarized as follows. We have demonstrated that $Reg(0)$ can be obtained by the extrapolating method. For $\Delta (\omega )$ with $\omega \leq 0$, we have proven that the method by Eq. (20) is more efficient than the method by Eq. (13). Besides, one can quickly judge whether a bound state exists or not by Eq. (25). We have found that the critical dipole moment is nearly independent of the radius of the nanosphere and the specific kind of noble metals, but heavily dependent on the QE-sphere distance. We have proven that the non-Markovian decay dynamics can be exactly obtained by Eq. (23), in which the parameter $Z$ can be calculated by either Eqs. (22) or (24). In spite of this, one need not to evaluate the lowest eigenfrequency $\omega _b$ [negative zero for the transcendental Eq. (21)] for the method by Eq. (24) compared to the method by Eq. (22). 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.

4. Existence condition 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 $\varepsilon _{2}$ for silver is from experimental data [98]. The photon GF is obtained by the finite element method presented in Ref. [74,75,83]. 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 $d_c$ for the system to form a bound state according to Eq. (25). Then, we show the performances of the methods [Eqs. (22) and (24)] for obtaining the long-time value of the excited-state population. At the end, the non-Markovian dynamics by solving Eq. (23) are demonstrated. For simplicity, the transition frequency for the QE is also $1.5\,\textrm {eV}$.

The calculated $\Gamma (\omega )$ and $\Delta (\omega )$ are shown in Fig. 7 for two different nanocavities with $H=20\,\textrm {nm}$ (the left column) and $H=50\,\textrm {nm}$ (the right column) respectively. Different from the previous section where $\Gamma (\omega )$ takes great value only in a much narrow frequency range, they are relatively large over a wide frequency range [see Figs. 7(a) and 7(b)]. 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 Figs. 7(a) and 7(b), we observe that the lowest plasmonic mode contributes much to $\Gamma (\omega )$, in which the peak position and peak value depend heavily on the rod height $H$.

 

Fig. 7. $\Gamma (\omega )$ and $\Delta (\omega )$. The left (right) column is for $H=20\,\textrm {nm}$ ($H=50\,\textrm {nm}$). (a) and (b) are for $\Gamma (\omega )$. (c) and (d) demonstrate the efficiency of Eq. (20) and Eq. (13) for calculating $\Delta (\omega )$ by choosing different cutoff frequency $\omega _c$. Here, $\omega =-0.1\,\textrm {eV}$. It is found that our new formalism Eq. (20) converges much more quickly than the usual method by Eq. (13). (e) and (f) are for $\Delta (\omega )$. The insets are for frequency around the lowest plasmonic mode.

Download Full Size | PDF

To compare the efficiency of the two methods [Eq. (20) and Eq. (13)] for calculating $\Delta (\omega )$ with $\omega <0$ in this example, we set $\omega =-0.1\,\textrm {eV}$ and vary the cutoff frequency $\omega _c$. The results are shown in Figs. 7(c) and 7(d), in which $\omega _c=70\,\textrm {eV}$ is enough to obtain a convergent result for the method by Eq. (20) while the results by the method of Eq. (13) show a large fluctuation near the average value even for $\omega _c=1000\,\textrm {eV}$. We have performed calculations for other different frequencies $\omega$ and have found similar behavior to the one shown here. Thus, we can conclude that the method by Eq. (20) converges much more quickly than the usual method by Eq. (13).

Figures 7(e) and 7(f) show the calculated $\Delta (\omega )$ by Eq. (18) for $\omega >0$ and by Eq. (20) for $\omega \leq 0$. We have found similar behavior that $\Delta (\omega )$ takes great value over a wide frequency range, and is especially large for $\omega$ around the resonance frequency of the plasmonic mode (see the insets therein).

Although $\Gamma (\omega )$ and $\Delta (\omega )$ look very different for the above two nanocavities, the values of the critical dipole moment $d_c$ are about the same. From the second column of Table 1, we see that they are nearly independent of the radius $R$ or the height $H$ of the nanorod. All are very close to $89.47\,\textrm {D}$ , which is the limiting result for $R\rightarrow \infty$ and $H\rightarrow \infty$ (calculated by using $\mathbf{\hat{r}}\cdot\operatorname{Re}\mathbf{G} ({\mathbf{r}_{0}},\mathbf{r}_{0},0)\cdot\mathbf{\hat{r}}= \sum_{i=1}^{+\infty}1/[8\pi\varepsilon_1(iW/2)^{2}]$ [75]). Physically, the QE is much close to the metal surface and the induced charges are mainly distributed at the undersides of the nanorods around the QE. This mechanism is also helpful to understant the phenomenon that the critical dipole moment for nanocavity composed of larger nanorods is more closer to its limiting value (see Table 1).

Table 1. The results of $d_c$, $\omega _b$ and $Z$ with different heights $H$ and radius $R$ of the nanorod. $d_c$ is the critical strength for the dipole moment obtained by Eq. (25). $\omega _b$ is the negative eigenfrequency for $f(\omega )=\omega -\omega _{0}-\Delta (\omega )=0$ [zero of Eq. (21)]. $Z_1$ and $Z_2$ are the amplitude for the excited state in the steady limit obtained by Eqs. (22) and (24), respectively. Their differences are very small. We observe that $d_c$ and $Z$ are nearly independent of the height $H$ and radius $R$ of the nanorod. But $\omega _b$ depends heavily on them, especially on the transition dipole moment $d$.

View Table

Then we show the performances of the methods [Eqs. (22) and (24)] for obtaining the long-time value of the excited-state population. $Z_1$ and $Z_2$ in Table 1 represent the $Z$ value obtained by Eqs. (22) and (24), respectively. Their differences are very small, which imply that both Eqs. (22) and (24) are able to calculate $Z$. In addition, we observe that $Z$ is also less dependent on the dimensions of the nanorod and the transition dipole moment $d$. This can be understood from $Z=1/(1+x)$ with $x=\int _{0}^{+\infty }\operatorname {Im}g(s)/(\omega _{b}-s)^{2}ds$ [Eq. (22)]. Although both $\omega _{b}$ and $\operatorname {Im}g(s)$ are dependent on the dimensions of the nanocavity, Z varies slowly with $x$ and accordingly is less dependent.

For the lowest eigenfrequency $\omega _b$ , we find that it depends heavily on the dipole strength $d$. This can be understood as follows. According to Eq. (21), ${\omega _b}$ is the solution for $(\omega { - }{\omega _0}{)/}{d^2}{ = }\Delta (\omega )/{d^2}$. Similar the Ref. [93], $\omega _b$ can be obtained through the intersection of the left hand $y_{L} = (\omega { - }{\omega _0})/{d^2}$ and the right hand ${y}_{R} = \Delta (\omega )/{d^2}$. The right hand is independent on the dipole strength $d$ and varies slowly with $\omega$ when $\omega \le 0$ [see Figs. 7(e) and 7(f) or see Eqs. (20) and (19)]. However, the linear function on the left hand depends heavily on the dipole strength $d$. Larger dipole moment shifts the line function more, which intersects the right hand at much smaller frequency $\omega _b$ .

The decay dynamics are shown in Fig. 8. Here, the height and the radius for the nanorod are $H=20\,\textrm {nm}$ and $R=4\,\textrm {nm}$, respectively. The critical dipole strength to form a bound state is $d_c=92.174\,\textrm {D}$ for a transition frequency $\omega _0=1.5\,\textrm {eV}$. Figure 8(a) demonstrates the excited state population $P_{a}(t)$ for the dipole strength $d$ just below and above its critical value $d_c$. 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\,\textrm {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\,\textrm {D}$, the excited state population $P_{a}(t)$ is partially preserved and tends to $0.64$ as $t\rightarrow \infty$.

 

Fig. 8. The excited state population $P_{a}(t)=\left \vert c_{1}(t)\right \vert ^{2}$. (a) $d$ is a little larger ($d=94\,\textrm {D}$ red circles) and smaller ($d=90\,\textrm {D}$ black solid line) than the critical value $d_c$. (b) $d$ is much larger ($d=130\,\textrm {D}$ red dash line) and much smaller ($d=50\,\textrm {D}$ black solid line) than the critical value $d_c$. Here, $d_c=92.1740\,\textrm {D}$ for $H=20\,\textrm {nm}$. (c) and (d) are similar to (a) and (b) except for a much higher rod $H=80\,\textrm {nm}$. In this case, the critical transition dipole moment is $d_c=90.9872\,\textrm {D}$. Here, the radius of the nanorod is $R=4\,\textrm {nm}$.

Download Full Size | PDF

We also consider another two different dipole moments, which are much smaller or larger than the critical value $d_c$. In these cases, the excited state populations $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 also be observed depending on the dipole moment [see the insets in Fig. 8(b)]. For $d=130\,\textrm {D}$, $P_{a}(\infty )=0.593$, which is around $Z^2=0.604$ for $d=120\,\textrm {D}$ (see the last two columns in Table 1). 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. Figures 8(c) and 8(d) demonstrate the results for a higher nanorod $H=80\,\textrm {nm}$ with the same radius $R=4\,\textrm {nm}$. We find that the results shown in Fig. 8(c) look similar to those in Fig. 8(a). Differently, the time for the system to reach its steady state is much smaller in Fig. 8(d) than that in Fig. 8(b), which can also be understood by Eq. (26).

5. Conclusion

In summary, we have presented some numerical methods for systematically investigating the formation of a QE-SPP bound state and its effect on the non-Markovian decay dynamics. Firstly, we have proposed an efficient numerical method [Eq. (25)] to determine the existence or absence of a bound state, in which only the real part of the photon GF $\operatorname {Re}[\mathbf {G} ({\mathbf {r}_{0}},\mathbf {r}_{0},0)]$ at zero frequency is sufficient. This avoids finding the zero of Eq. (10) or solving the decay dynamics in the long-time limit. In addition, it also helps us to understand the decay dynamics around the critical condition, where the long-time value of the excited-state population shows a transition from zero to a limited value by slightly modifying some parameters of the system. Secondly, we have found that the critical dipole moment is much dependent on the distance between the QE and the metal surface but less dependent on the radius of the nanosphere or the size of the nanorod composing the nanocavity. We have numerically proved that the critical condition is independent of the actual kind of nobel metal material, since the perfect electric conductor boundary can usually be assumed for the nobel metal surface in the electrostatic case. Thirdly, we have extended the method for evaluating $\Delta (\omega )$ for $\omega$ in the photonic bandwidth [Eq. (18) together with Eq. (19)] to the case for $\omega$ below the lower threshold [$\omega \leq 0$ Eq. (20) in this work]. We have numerically proved that our method Eq. (20) is much more efficient than the usual method Eq. (13). Fourthly, for a bound state presented, we have proposed a new formalism Eq. (24) to determine the long-time value for the excited state population $P_a(\infty )=Z^2$, in which one only needs to perform a simple integral of the evolution spectrum $S(\omega )$ over the positive frequency range. Compared to the usual method Eq. (22), this method avoids the need to explicitly calculate the energy eigenvalue $\omega _b$ [zero for Eq. (21)] of the bound state, and accordingly the analytical part of the self-energy $\Delta (\omega )$ with $\omega <0$. In addition, we have found that $Z$ is less affected by the dipole strength $d$ with $d\geq d_c$, but the eigenfrequency of the bound state $\omega _b$ is much dependent on it. Lastly, 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 once the energy eigenvalue $\omega _b$ and the parameter $Z$ are obtained. We have clearly shown that $P_a(\infty )$ changes from zero to about 0.6 by slightly modifying the transition dipole moment. In addition, we have found that $P_a(\infty )$ is nearly independent of the sizes of the nanorods and the transition dipole strength $d$ when $d>d_c$, which is consistent with the prediction by Eqs. (22) and (24).

It is envisaged that the decay dynamics with a bound state predicted by the above theory may be experimentally observed by employing J-aggregates as the QE, in which the transition dipole moment is about $d=93.4\,\textrm {D}$ (estimated from $\tau _0=1/\gamma _0=70\,\textrm {ps}$ in Ref. [102]) at low temperature (below $50\,K$) and about $d=34\,\textrm {D}$ at room temperature [5]. It satisfies the requirement of realizing a bound state according to our prediction for the nanocavity system. In addition, the transition dipole moment can be above or below the critical value by tuning the temperature in principle [102]. With currently available nanotechnologies, individual metallic nanoparticle and plasmonic nanocavity have been controllably fabricated [5,7,8]. In addition, sub-nanometre spatial control over the coherent coupling between a single molecule and a plasmonic nanocavity in close proximity is demonstrated [7]. The QE can be excited similarly to the method in Ref. [8], in which the QE is firstly excited to a much higher energy state and then relaxes quickly to the excited state. Thus experimentally demonstrating the predicted phenomenon is possible in principal.