Level shift and decay dynamics of a quantum emitter around plasmonic nanostructure

TL;DR

This paper introduces a new, efficient method for calculating quantum energy level shifts of emitters near nanostructures, avoiding complex integrals and requiring less frequency data, with demonstrated applications to nanospheres and nanocavities.

Contribution

A novel approach for computing quantum energy level shifts that simplifies calculations by eliminating principal value integrals and reducing frequency range requirements.

Findings

The new method avoids principal value integrals and imaginary frequency calculations.

It requires a narrower frequency range for accurate results.

Reversible quantum dynamics are observed in the studied systems.

Abstract

We put forward a general approach for calculating the quantum energy level shift for emitter in arbitrary nanostructures, in which the energy level shift is expressed by the sum of the real part of the scattering photon Green function (GF) and a simple integral about the imaginary part of the photon GF in the real frequency range without principle value. Compared with the method of direct principal value integral over the positive frequency axis and the method by transferring into the imaginary axis, this method avoids the principle value integral and the calculation of the scattering GF with imaginary frequency. In addition, a much narrower frequency range about the scattering photon GF in enough to get a convergent result. It is numerically demonstrated in the case for a quantum emitter (QE) located around a nanosphere and in a gap plasmonic nanocavity. Quantum dynamics of the emitter is calculated by the time domain method through solving Schr\"{o}dinger equation in the form of Volterra integral of the second kind and by the frequency domain method based on the Green's function expression for the evolution operator. It is found that the frequency domain method needs information of the scattering GF over a much narrower frequency range. In addition, reversible dynamics is observed. These findings are instructive in the fields of coherent light-matter interactions.

Equations

$$H$$

$$H_0$$

$$H_I$$

$$\hat{\mathbf{E}}(\mathbf{r},\omega )=i\sqrt{\frac{\hslash }{\pi {\epsilon }_0}}\int d\mathbf{s}\sqrt{{\epsilon }_I(\mathbf{s},\omega )}\mathbf{G}(\mathbf{r},\mathbf{s},\omega )\cdot \hat{\mathbf{f}}(\mathbf{s},\omega )$$

$$|\Psi (t)⟩$$

$$+\int dr{\int }_0^{+\infty }d\omega C(r,\omega ,t)e^{-i\omega t}|F_{r,\omega }⟩.$$

$$\stackrel{\cdot }{c_1}(t)=-{\int }_0^{t}K(t-\tau )c_1(\tau )d\tau ,$$

$$g_{rr}(\omega )\equiv \frac{{\mathbf{d}}^{\ast }\cdot \mathbf{G}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\cdot \mathbf{d}}{\hslash \pi {\epsilon }_0}.$$

$$c_1(t)=c_1(0)-{\int }_0^{t}B(t-t^{{}^{\prime }})c_1(t^{{}^{\prime }})d{t}^{{}^{\prime }}$$

$$R_{ii}(z)=\frac{1}{\pi {\epsilon }_0}[{\int }_0^{\infty }d\omega \frac{{\mathbf{d}}^{\ast }\cdot \mathrm{Im}\mathbf{G}({\mathbf{r}}_A,{\mathbf{r}}_A,\omega )\cdot \mathbf{d}}{z-\omega }$$

$$R_{ii}^{\pm }(z)=\underset{\eta \to 0^{+}}{\lim }{R}_{ii}(z\pm i\eta )=\Delta (z)\mp i\frac{\Gamma (z)}{2},$$

$$\Gamma (z)=2\frac{{\mathbf{d}}^{\ast }\cdot \mathrm{Im}\mathbf{G}({\mathbf{r}}_0,{\mathbf{r}}_0,z)\cdot \mathbf{d}}{\hslash {\epsilon }_0}=2\pi \mathrm{Im}{g}_{rr}(z)\theta (z).$$

$$\Delta (z)=\frac{1}{2\pi }\wp {\int }_0^{+\infty }ds\frac{\Gamma (s)}{z-s}.$$

$$c_1(t)=\underset{\eta \to 0_{+}}{\lim }[{\int }_{-\infty }^{+\infty }S(\omega )e^{-i(\omega -{\omega }_0)t}d\omega ,$$

$$S(\omega )=\frac{1}{\pi }\underset{\eta \to 0_{+}}{\lim }\frac{\Gamma (\omega )/2+\eta }{[\omega -{\omega }_0-\Delta (\omega ){]}^2+(\Gamma (\omega )/2+\eta {)}^2]}.$$

$$\Delta (\omega )$$

$$\frac{\omega }{\hslash \pi {\epsilon }_0}\mathbf{d}\cdot {\int }_0^{+\infty }d\xi \frac{\mathbf{G}({\mathbf{r}}_0,{\mathbf{r}}_0,i\xi )}{{\omega }^2+{\xi }^2}]\cdot \mathbf{d}$$

$$=-\pi \mathrm{Re}{g}_{rr}(\omega )+\omega {\int }_0^{+\infty }d\xi \frac{g_{rr}(i\xi )}{{\omega }^2+{\xi }^2}$$

$$\Delta (\omega )$$

$$=-\pi \mathrm{Re}{g}_{rr}(\omega )+\wp {\int }_0^{+\infty }ds\frac{\mathrm{Im}{g}_{rr}(s)}{\omega +s}.$$

$$\Delta (0)=-\pi \mathrm{Re}{g}_{rr}(0)+\wp {\int }_0^{+\infty }ds\frac{\mathrm{Im}{g}_{rr}(s)}{s}.$$

$$\Delta (\omega )=-\pi \mathrm{Re}{g}_{rr}(\omega )+\frac{\pi }{2}\mathrm{Re}{g}_{rr}(0)-\omega {\int }_0^{+\infty }ds\frac{\mathrm{Im}{g}_{rr}(s)}{(\omega +s)s}.$$

Full text

one two three

††preprint: HEP/123-qed

year number number identifier Date text]date

LABEL:FirstPage101 LABEL:LastPage#1102

I Introduction

Recently, considerable attention is devoted to the fields of light-matter interaction for both fundamental and applicative purposes Cohen-Tannoudji et al. (1989); Agarwal (1974); Milonni (1994); Buhmann (2012a, b); Oulton et al. (2009); Chikkaraddy et al. (2016); Muskens et al. (2007); Kinkhabwala et al. (2009); Curto et al. (2010); Ringler et al. (2008); Ming et al. (2009); Vecchi et al. (2009); Xue-Hua et al. (2003); Birnbaum et al. (2005); Lukin and Imamoğlu (2000); Khitrova et al. (2006); Yves et al. (2007); Noda et al. (2007); Altug et al. (2006); Schwartz et al. (2011); Ren et al. (2017); Rousseaux et al. (2018); Peng et al. (2017); Zhang et al. (2018); Haran and Chuntonov (2018); Wen et al. (2017); Kneipp et al. (1997). Many novel phenomena have been predicted and related devices have been developed, for example, enhanced and inhibited spontaneous emission, photon blockade, reversible spontaneous emission, unidirectional emission, trapping atoms by vacuum forces, enhanced Raman scattering, LEDs, one-atom maser, low-threshold lasers, etc. Among all, the coherent interaction between a single QE and its electromagnetic environment is of paramount importance, since it constitutes one of the most fundamental aspects in QED and is also a practical need in the fields of quantum information procession as well as single-molecule sensing. Coherent decay dynamics of an excited QE and the energy level shift are two elementary terms in this field.

Plasmonic nanostructure with an ultrasmall optical mode volume is one of the most promising platforms for the above two problems, since strong light-matter interactions at the single-exciton level has been achieved Liu et al. (2017). In addition, precise spatial control of the QE position with respect to the plasmonic nanocavity has been realized Zhang et al. (2017). Further, high-quality metallic nanoparticles with tunable size and controllable shapes can be produced and rationally assembled into well-defined architectures Tan et al. (2011); Shen et al. (2018).

Theoretically, the medium-assisted quantized electromagnetic field can be expressed as the fundamental bosonic vector fields via photon GF Dung et al. (2003). Accordingly, all the information needed to investigate the energy level shift and decay dynamics is contained in the photon GF. For example, the energy level shift is expressed by a principal value integration of the electromagnetic photon GF over the whole frequency range Buhmann (2012a, b). For nanostructure with high symmetry, such as sphere, the photon GF can be obtained by semi-analytical methods Tai (1993). But for nanostructure with arbitrary shape, it is not an easy task. We usually have to resort to numerical methods Gallinet et al. (2015); Zhao et al. (2018a, b); Tian et al. (2019); Chen et al. (2010); Bai et al. (2013); Huang et al. (2012); Vlack and Hughes (2012), such as finite difference time domain method or finite element method, which are time consuming or large memory space requirement. Since the direct principal value integration method needs information for the photon GF over a wide frequency range, one treatment is to transform the principle value into an ordinary integration over the imaginary axis, where the Kramers-Kronig (KK) relations of the GF is exploited Dzsotjan et al. (2011); Buhmann (2012b). In this method, the photon GF as well as the response function for the material with imaginary frequency are needed and hard to be obtained Woods et al. (2016); Rodriguez et al. (2007); Tkatchenko et al. (2013). Differently, we propose a new method by utilizing the subtractive KK relations without transforming into the imaginary frequency and without worrying about the principal value integral. We will show that this will greatly simplifies the calculation and there is no need of the knowledge for the photon GF over a wide imaginary frequency range.

For the decay dynamics, specially in case of strong coupling regime, coherent decay dynamics exhibit significant reversibility and non-Markovian methods should be taken into account Rivas et al. (2014). Usually, this can be numerically resolved by solving a quantum master equation Breuer et al. (2016); Ge and Hughes (2015) or the Schr$\ddot{o}$dinger equation Dung et al. (2003); González-Tudela et al. (2014); Hakami and Zubairy (2016) which leads to the well-known Volterra integral equations of the second kind. Besides the above time domain methods, frequency domain methods based on retarded and advanced Greens function expression for the evolution operator can be used Xue-Hua et al. (2003); Cohen-Tannoudji et al. (1989). As explained in Ref. Huang et al. (2012); Liu et al. (2017); Yang and An (2017); Rousseaux et al. (2018); Gaveau and Schulman (1995); Lambropoulos et al. (2000); Garmon et al. (2019); Cohen-Tannoudji et al. (1989), thit method supplies more obvious physical significance underlying some phenomena and provides distinctive criterion to discriminate between strong coupling and weak coupling. In addition, compared to the time domain method, there is no need of time convolution.

In this work, we first present the theory and derive a general method for calculating the energy level shift. We will show that the energy level shift can be expressed by the sum of the real part of the photon GF and an integral part with the integrand a well-behaved function of the imaginary part of the photon GF. For the following two sections, we apply our method to a particular example where a QE is located around a gold nanosphere. For the energy level shift, we numerically compare our method with the direct Hilbert transformation method and the imaginary frequency integration method Buhmann (2012b); Dzsotjan et al. (2011). The good performance of our method will be demonstrated. Then, we show the characteristics of the frequency domain method based on GF and the time domain method of the Volterra integral equation form for calculating dynamics. Section V is devoted to the study of the energy level shift and decay dynamics for QE located in a plasmonic nanocavity. We will demonstrate that frequency domain method for dynamics is more efficient than the time domain method. Finally, a summary is given in Sec. VI.

II Model and Method

Let us consider a two-level QE coupled to a common electromagnetic reservoir. By using the dipole and rotating-wave approximations, the total Hamiltonian for the system is Dung et al. (2003)

[TABLE]

Here, $\hat{\mathbf{f}}\left(\mathbf{r},\omega\right)$ and $\hat{\mathbf{f}}^{\dagger}\left(\mathbf{r},\omega\right)$ are the bosonic vector field annihilation and creation operators for the elementary excitation of the electromagnetic reservoir, respectively. $\mathbf{d}=\langle g|\hat{\mathbf{d}}|e\rangle=d\hat{\mathbf{n}}$ is the element of the transition dipole moment, with the unit vector $\hat{\mathbf{n}}$ and its strength $d$. The electric field vector operator $\hat{\mathbf{E}}\left(\mathbf{r},\omega\right)$ is given by

[TABLE]

where $\mathbf{G}(\mathbf{r},\mathbf{s},\omega)$ is the photon GF defined as [$\nabla\times\nabla\times-\varepsilon(\mathbf{r},\omega)\omega^{2}/c^{2}]\mathbf{G}(\mathbf{r},\mathbf{s};\omega)=\omega^{2}/c^{2}\mathbf{I}\delta(\mathbf{r}-\mathbf{s})$. Here $\varepsilon(\mathbf{r},\omega)=\varepsilon_{R}(\mathbf{r},\omega)+i\varepsilon_{I}(\mathbf{r},\omega)$ is the spatially and frequency-dependent complex relative dielectric function and $\varepsilon_{I}(\mathbf{r},\omega)$ is its imaginary part. $\mathbf{I}$ is the unit dyad and $c$ refers to the speed of light in the 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\rangle=|e\rangle\otimes|0\rangle$ and $|F_{r,\omega}\rangle=|g\rangle\otimes|1_{r,\omega}\rangle$ with $|e\rangle$ ($|g\rangle$) the excited (ground) state of the QE and $|1_{r,\omega}\rangle\equiv\hat{\mathbf{f}}_{j}^{\dagger}\left(\mathbf{r},\omega\right)|0\rangle$. $|0\rangle$ is the zero photon state. The time evaluation for this initial state is given by

[TABLE]

The time-dependent Schrödinger equation leads to the following equation of motion for the probability amplitudes:

[TABLE]

where the kernel function is $K(t-\tau)=\int_{0}^{+\infty}d\omega J(\omega)e^{i(\omega_{0}-\omega)(t-\tau)}$ with the spectral density defined $J(\omega)=\operatorname{Im}g_{rr}(\omega)$. Here, the coupling strength $g_{rr}(\omega)$ is defined by

[TABLE]

One method for dynamics is to take the time integral on both sides of Eq. (1). This leads to the well-known Volterra integral equations of the second kind. Explicitly, the probability amplitude for the QE in the excited states is

[TABLE]

with $B(t-t^{{}^{\prime}})=\int_{0}^{+\infty}d\omega J(\omega)\int_{0}^{t-t^{{}^{\prime}}}e^{-i(\omega-\omega_{0})u}du$, where time-convolution is needed. Note that the kernel $B(t-t^{{}^{\prime}})$ is an integral over the whole positive frequency spectrum.

Another method for dynamics is by the resolvent operator technique Lambropoulos et al. (2000); Xue-Hua et al. (2003); Cohen-Tannoudji et al. (1989). From the GF expression of evolution operator, one can show that $U(t)=\int_{-\infty}^{+\infty}d\omega[G^{-}(\omega)-G^{+}(\omega)]exp(-i\omega t)/(2\pi i)$ where $G^{\pm}(\omega)=\lim_{\eta\rightarrow 0^{+}}G(E\pm i\eta)$ with the resolvent operator $G(z)=(z-H/\hbar)^{-1}$. The probability amplitude $c_{1}(t)$ can be expressed the matrix element of the evolution operator $c_{1}(t)=\langle I|U(t)|I\rangle$. From the operator identity $(z-H_{0}/\hbar)G(z)=H_{I}G(z)/\hbar$, we obtain $G_{ii}(\omega)\equiv\langle I|U(t)|I\rangle=(\omega-\omega_{0}-R_{ii}(\omega))^{-1}$ in which the matrix element of the level-shift operator $R_{ii}(\omega)$ reads

[TABLE]

Clearly, by using the relation $1/(z-\omega-i\eta)=\wp(1/(z-\omega))+i\pi\delta(z-\omega)$, this can be written

[TABLE]

with

[TABLE]

Here, $\theta(z)$ is the step function. $\Delta(z)$ is the Hilbert transform of $\Gamma(z)$, which is

[TABLE]

Thus, the probability amplitude reads

[TABLE]

with evolution spectrum

[TABLE]

Although the dynamics for our system can be obtained by the above two methods (Eq. (3) and Eq. (7)), one has to evaluate the kernel $B(t-t^{{}^{\prime}})$ or $\Delta(\omega)$. For both, an integral over the whole frequency range should be performed. However, the Hamiltonian is nonrelativistic and can not be applied in the case of relativistic high frequency. In addition, both $\Delta(\omega)$ and the kernel $B(t-t^{{}^{\prime}})$ are divergent. In the homogeneous case such as vacuum, the divergence can be overcome by the procedure of mass renormalization. In artificial nanostructure, there is no response for the medium in the high-frequency range. Thus, one treatment is to separate the homogeneous and scattering contributions, where the photon GF is decomposed into its bulk and scattered part $\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},\omega)=\mathbf{G}_{0}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},\omega)+\mathbf{G}_{s}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},\omega)$. The scattering photon GF can be used to take place of the total photon GF where the homogeneous-medium contribution is attributed to the definition of the transition frequency of the QE. In this work, we will numerically confirm that this nature renormalization procedure can overcome the divergence.

In the weak coupling limit where both $\Gamma(\omega)$ and $\Delta(\omega)$ are small comparied with $\omega_{0}$ and vary slowly with frequency $\omega$ around $\omega_{0}$, $\Gamma(\omega)$ and $\Delta(\omega)$ in Eq. (8) can be safely replaced by $\Gamma(\omega)=\Gamma$ and $\Delta(\omega)=\Delta$. Then, the spectrum $S(\omega)=\frac{1}{\pi}\frac{\Gamma/2+\eta}{[\omega-\omega_{0}-\Delta]^{2}+(\Gamma/2+\eta)^{2}]}$ is of a Lorentzian form and the result for Eq. (7) is $c_{1}(t)=\exp[-\left(i\Delta+\Gamma/2\right)t]$. This clearly demonstrates that $\Gamma$ and $\Delta$ are the spontaneous emission rate and the energy level shift, respectively. In the following of this paper, $\Gamma(\omega)$ and $\Delta(\omega)$ are termed spontaneous emission rate and energy level shift, even though both of them are maybe highly peaked and vary rapidly with frequency.

The spontaneous emission rate $\Gamma(\omega)$ (Eq. (5)) can be obtained once we known the photon GF (see Eq. (2)). But for the energy level shift $\Delta(\omega)$ (Eq. (6)), there needs a principle integration which is difficult from a numerical view. Alternatively, there is one method based on the contour-integral techniques, where real frequency integral is transformed into ones along the positive imaginary axis plus contributions from the poles Buhmann (2012a, b); Dzsotjan et al. (2011). By utilizing the KK relation of the photon GF, one has Dzsotjan et al. (2011); Dung et al. (2002)

[TABLE]

Besides, we propose a third method by using the subtractive KK relation. By using the relations $-\pi\operatorname{Re}\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},s)=\wp\int_{-\infty}^{+\infty}ds\operatorname{Im}\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},s)/(\omega-s)$ and $\operatorname{Im}\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},-s)=-\operatorname{Im}\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},s)$ , the energy level shift for $\omega\geq 0$ can be written Dzsotjan et al. (2011)

[TABLE]

For $\omega=0$ , the above equation becomes

[TABLE]

Subtracting Eq. (11) from Eq. (10), we have

[TABLE]

Here, we have used the relation $\Delta(0)=-0.5\pi\operatorname{Re}g_{rr}(0)$, since the second term on the right hand side in Eq. (11) is $-\Delta(0)$ ( see Eq. (6) ).

Equation (12) is the central result of our method for evaluating the energy level shift of a QE around arbitrary nanostructure. In this form, there is no need to worry about the principal value. Furthermore, for frequency away from the practical resonance where material becomes transparent, the scattering is weak. Thus, $\operatorname{Im}g_{rr}(s)$ in the integrand is small, which is very useful in calculating the energy level shift by numerical means.

Before proceeding further, let us give some discussions about the above three methods to obtain the energy level shift $\Delta(\omega)$. Hereafter, Eq. (9), Eq. (6) and Eq. (12) refer to the imaginary frequency method, the direct Hilbert method and the subtractive KK method respectively. For the imaginary frequency method, real frequency integral is transformed into the positive imaginary axis, where the photon GF $\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},i\xi)$ decays exponentially with $\xi$ and there is no need to worry about the principle value. However, there is extra burden to compute the Green tensor for imaginary frequency $\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},i\xi)$, compared to the other two methods ( Eq. (6) and Eq. (12) ). For nanostructures with high symmetry, $\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},i\xi)$ can be obtained semi-analytically. But for arbitrary nanostructure which is the usual case, it is a difficult task Woods et al. (2016); Rodriguez et al. (2007); Tkatchenko et al. (2013). Comparing the two integral part in Eq. (6) and Eq. (12), we find that the integrand in our subtractive KK method ( Eq. (12) ) decays faster than that in the direct Hilbert method ( Eq. (6) ) with high frequency $s$. This can be clearly seen by comparing the integrands. Thus, our subtractive KK method should converge more rapidly than the direct Hilbert method. Further, there is no need to worry about the principle value as in the imaginary frequency method ( Eq. (9) ). In the next section, we will demonstrate that our subtractive KK method ( Eq. (12)) is powerful and can be used to obtain the exact energy level shift efficiently.

As a demonstration, we will apply the methods introduced above to a particular example where a two-level QE is located above a metal nanosphere (see Fig. 1). The nanosphere with radius $a$ is located at the origin. A QE at a distance $h$ from the surface of the sphere lies on the x-axis of the coordinate system. The metal is chosen to be Gold and characterized by a complex Drude dielectric function Ge and Hughes (2015) $\varepsilon(\omega)=1-\omega_{p}^{2}/\omega(\omega+i\gamma_{p})$ with $\omega_{p}=1.26\times 10^{16}$ $rad/s$ and $\gamma_{p}=1.41\times 10^{14}$ $rad/s$. The background is vacuum with $\varepsilon_{B}=1$. For simplicity, the dipole is polarized along the radial direction of the sphere $\mathbf{d=}d\mathbf{r}$ and its strength is set to $d=24D$. In this case, the photo GF is obtained semi-analytically and and and (1994); Tai (1993)

In the following sections, we first demonstrate the validity and advantages of our method comparing with the other two methods in obtaing the energy level shift. Then, the properties of the the two methods (Eq. (3) and Eq. (7)) for evaluating dynamics are numerically investigated with different truncation conditions. In both cases, we choose the nanosphere system shown in Fig. 1 as an example, which is widely investigated and can be treated analytically. We will show that dynamics by the Fourier-Laplace transformation method with level shift obtained by our subtractive Kramers-Kronig relation can be numerically evaluated with the lease computation resource. At the end, we will apply the method introduced above to investigate the dynamics of a quantum dot in plasmonic nanocavity. Non-Weisskopf-Wigner decay phenomenon will be shown.

III NUMERICAL COMPARISION OF THE ABOVE THEREE METHODS FOR ENERGY LEVEL SHIFT

In this section, the above three methods ( Eq. ( 6 ), Eq. ( 9) and Eq. ( 12 ) ) for obtaining the energy level shift are numerically investigated. For nanosphere, the scattering photon GF can be obtained semi-analytically Zhao et al. (2018a). The properties of the coupling strength defined by $g_{rr}(\omega)\equiv\mathbf{d\cdot G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},\omega)\cdot\mathbf{d/}\hbar\pi\varepsilon_{0}$ ( Eq.(2) ) are shown in Fig. 2 for $a=20nm$ and $h=1nm$. From the inset in Fig. 2(a), we clearly see that there are some peaks for $\operatorname{Im}g_{rr}(\omega)$ in the frequency range between $4eV$ and $6eV$, which stem from the localized surface plasmon (LSP) resonance. Over this frequency range, LSP takes great effect and both the real part and the imaginary part of $g_{rr}(\omega)$ may be large (see Fig. 2(a) and 2(b)). But for frequency away, both $\operatorname{Re}g_{rr}(\omega)$ (real part) and $\operatorname{Im}g_{rr}(\omega)$ (imaginary part) are small. This leads to the usual assumption that higher frequency parts contribute little to the dynamics or energy level shift.

To evaluate the energy level shift, the first method (Eq. (9)) needs information about the coupling strength $g_{rr}(i\xi)$ with imaginary frequency $i\xi$. Although it is not an easy task to obtain the precise photon GF with imaginary frequency in arbitrary nanostructure, its exact value can be obtained by replacing the frequency $\omega$ in the semi-analytical expression of $\mathbf{G}(\mathbf{\mathbf{r}_{0}},\mathbf{r}_{0},\omega)$ with $i\xi$. Figure 2(c) shows that $g_{rr}(i\xi)$ is real and decays with $\xi$, which agree well with the description of Ref. Rodriguez et al. (2007). It is well-behaved and the second part in Eq. (9) can be obtained precisely without difficulty for this particular example. Results from this method are used as a reference for the other two methods (see Eq. (6) and Eq. (12)). Note that $g_{rr}(i\xi)$ spread over a wide frequency range. The inset shows the integral part in Eq. (9) from $\xi$ to $200eV$ $Diff(\xi)=\omega\int_{\xi}^{200eV}dsg_{rr}(is)/(\omega^{2}+s^{2})$ for $\omega=5eV$. It is about $1.1meV$ for $\xi=10eV$.

Figure 2(d) shows the integrand $\omega_{0}\operatorname{Im}g_{rr}(\omega)/(\omega_{0}+\omega)\omega$ as a function of frequency $\omega$ in our subtractive KK method ( see Eq. (12)) with $\omega=4.41eV$, which is around the LSP dipole mode. Also shown, in the inset of Fig. 2(d), is a plot over the frequency range where LSP takes great effect. We find that this integrand is also well-behaved and can be handled easily and precisely. This is true for other transition frequency $\omega$. Note that there is no need about the knowledge of the photon GF on the imaginary axis, which is different from the first method. In addition, the integrand in our method (Eq. (12)) is peaked over a narrower frequency range around the resonance frequency of the nano-sphere, which is different from the integrand in Eq. (9) (see Fig. 2(c) where it remains over a much wider frequency range).

The energy levle shift obtained by the above three methods are shown in Fig. 3. From a computational perspective, the upper limit of integration $+\infty$ should be replaced by some cut-off frequency $\omega_{max}$. Figure 3(a) shows the energy level shift $\Delta(\omega)$ by the first method (Eq. (9)), where $\omega_{\max}=200eV$. The results are convergent, as the difference between $\omega_{\max}=100eV$ and $\omega_{\max}=200eV$ is less than $10^{-7}eV$. The inset shows the integral term for the imaginary frequency method ( Eq. (9)) ). Although, it is in the range of a few $meV$ for parameters considered here, it can not be overlooked where its appearance and importance are discussed in Ref. Dzsotjan et al. (2011).

As a demonstration of our subtractive method, Fig. 3(b) shows the results from Eq. (12) with the integral part ( the last term ) shown in the inset. The black solid line and the red dots are results for integral on the range between $(0eV,200eV)$ and $(3eV,8eV)$, respectively. It is found that they agree well, which means that our method needs knowledge over a narrow frequency range which is different from the imaginary frequency method (see Fig. 2(c) and inset therein). In addition, this term is also in the range of a few $meV$. Different from that shown in the inset in Fig. 3(a), it grows with increasing transition frequency. In the visible range, this integral part is lower than that in the inset of Fig. 3(a). For example, at $\omega=2eV$, it is below $1.5meV$, while it is above $3.5meV$ for the imaginary frequency (see the inset in Fig. 3(a)).

To test the accuracy, result from the imaginary frequency method ( Eq. (9)) is thought to be precise and served as a reference. We define the relative errors $\Delta_{i\_error}(\omega)=\Delta_{i}(\omega)-$ $\Delta(\omega)$ $(i=2,3)$, where $\Delta(\omega)$ is the energy level shift obtained by the imaginary frequency method (results in Fig. 3(a) by Eq. (9)) and $\Delta_{2}(\omega)$ ($\Delta_{3}(\omega)$) represents the results by the direct Hilbert method ( our subtractive KK method ) through Eq. (6) ( Eq. (12) ). Numerically, the upper limit of both integrations should be truncated to some value $\omega_{\max}$. Figure 3(c) are the results for $\Delta_{2\_error}(\omega)$ with $\omega_{\max}=10eV$ (black dash dot), $20eV$ ( purple dash ), $50eV$ ( blue dot ) and $200eV$ ( red solid line ). We find that $\left|\Delta_{2\_error}(\omega)\right|$ in the frequency range is a few $meV$, when the cut-off frequency is less than $\omega_{\max}=20eV$. A much higher value for $\omega_{\max}$, for example, $\omega_{\max}>20eV$ is needed in order to meet $\left|\Delta_{2\_error}(\omega)\right|<1meV$ . But for our subtractive method, $\left|\Delta_{3\_error}(\omega)\right|$ shown in Fig. 3(d) are less than $0.008meV$ for all the above four different cut-off frequency $\omega_{\max}$. This clearly confirms that knowledge of the photon GF over a narrow frequency range is enough to get the precise energy level shift by our subtractive KK method ( Eq. (12) ).

It should be stressed that numerical evaluation of the photon GF is not an easy task for arbitrary nanostructure. Recently, we have proposed a finite element method to exactly calculate the scattering photon GF. For one frequency point, two different runs are needed. In addition, for a single run, the typical computational time is about half an hour on our workstation with processor ”Intel(R) Xeon(R) E5-2697 v3” and the memory usage is about $30GB$, where the simulation domain has been even reduced to one quarter for symmetry consideration of the nanosphere system. Since the frequency range by our method is much narrower than by the direct Hilbert method, we can conclude that our method is much better. In addition, once the dynamics is calculated by the frequency domain method through Eq. (7), one has to resolve the photon GF on the imaginary axis. Compared to our subtractive KK method where knowledge of photon GF with real frequency is enough, the imaginary frequency method needs extra simulation.

To further demonstrate the good performance of our subtractive KK method, we also investigate the case for $a=10nm,h=2nm$ and $a=20nm,h=1nm$. Figure 4(a) and 4(b) show similar phenomena as those in Fig. 3(c) and 3(d), where the relative error is a few $\mu eV$ for our method and is a few $meV$ for the direct Hilbert method with $\omega_{max}\leq 20eV$. The relative error for the direct Hilbert method increases with the emitter approaching the surface of nanosphere, which can be clearly seen by comparing Fig. 4(c) with Fig. 4(a) or Fig. 3(c). However, the relative error for our method remains at a extremely low level, for example, less than $10\mu eV$ for $\omega_{max}\geq 10eV$.

From the above results, we can conclude that Eq. (12) is an efficient method to obtain the energy level shift of a QE in arbitrary nanostructure, where the argument of the frequency in photon GF is real. The integration can be made by common techniques and converges much quickly.

IV CHARACTERISTICS OF THE TWO METHODS FOR DYNAMICS

In this section, we numerically demonstrate the characteristics for the above two methods shown in Eq. (3) and Eq. (7) used in calculating dynamics. For simplicity, the parameters about the system are the same as those in Fig. 4(c) and 4(d). The transition frequency is set to $\omega_{0}=5eV$. Figure 5(a) are the results by the time domain method through solving the Volterra integral equation ( Eq. (3) ) with the cut-off frequency $\omega_{max}=10eV$ ( blue solid line ) and $\omega_{max}=20eV$ ( red circles ). We see that results are the same for both cases, which means that a small cut-off frequency, for example, $\omega_{max}=10eV$ is enough to obtain a convergent result for this particular example.

For the frequency domain method based on Green’s function expression for the evolution operator ( Eq. (7) ), we need the knowledge of the energy level shift $\Delta(\omega)$. In the above section, we have shown that our method converges more quickly than the direct Hilbert transform method ( Eq. 6 ). In addition, there are some minor errors for relative low cut-off frequency. Figure 5(b) shows how these errors influence the dynamics of the emitter. From the inset, we find that results when $\Delta$ is obtained from the direct Hilbert method ( Eq. 6 ) deviate from those when $\Delta(\omega)$ is obtained by the imaginary frequency domain method ( Eq. (9) ) or our subtractive KK method ( Eq. (12) ). In addition, the larger the cut-off frequency $\omega_{max}$ is, the less the error is.

Thus, the dynamics by the frequency domain method based on the Green’s function expression of the evolution operator is sensitive to the energy level shift $\Delta(\omega)$, which implies that $\Delta(\omega)$ should be precise to obtain a convergent results. In addition, Fig. 5(c) shows that both methods produce the same results as long as $\Delta(\omega)$ is convergent. This phenomenon remains for different parameters, such as the transition frequency $\omega_{0}$, the location of the emitter and the radius of the nanosphere, which are not shown here.

Although the dynamics of an excited dipole can be addressed by both methods in this example, the frequency domain method (Eq. (7)) provides direct information about the spontaneous emission spectrum, which is helpful to understand the time evolution of the system in various regimes (Detail discussion can be found in Chapter III of Ref. Cohen-Tannoudji et al. (1989), where the spontaneous emission spectrum changes progressively from a Lorentizan form to a set of two delta functions are demonstrated.). Recently, bound state with decoherence dynamics in photonic crystal John and Wang (1990); John and Quang (1994); González-Tudela and Galve (2019); Liu and Houck (2016), cavity arrays Lombardo et al. (2014) and plasmonic nanostructure have been discussed. To form a bound state, there should be a discrete eigenstate with eigenenergy in the photonic band gap or in the negative axis $\omega_{b}<0$, where $\omega_{b}$ is one root for equation $\omega-\omega_{0}-\Delta(\omega)=0$. With this in hand, it is helpful to explain the above phenomena and it is instructive to classify the coupling in various regimes. In addition, we will demonstrate in the next section that a much narrower frequency range for the photon GF is needed by the frequency domain method, when more realistic permittivity beyond the Drude model assumption for material is taken into account.

The main results of this section can be summarized as follows. Decay dynamics can be obtained by the time domain method by solving the Volterra integral equation ( Eq. (3) ) or by the frequency domain method based on the Green’s function expression of the evolution operator ( Eq. (7) ) with the energy level shift $\Delta(\omega)$ calculated by our subtractive KK method ( Eq. (12) ), in which information over a narrow real frequency range for the photon GF is enough. If the energy level shift $\Delta(\omega)$ is calculated by the imaginary frequency method ( Eq. (9)) or by the direct Hilbert method ( Eq. (6)), knowledge of the photon GF either over a wide imaginary frequency range or a wide real frequency range is demanded to ensure a convergent result.

V LEVEL SHIFT AND DYNAMICS OF A QUANTUM EMITTER IN A PLASMONIC

NANO-CAVITY

In this section, we apply Eq. (12), Eq. (7) and Eq. (3) to investigate the level shift and dynamics of a QE in a plasmonic nanocavity (see Fig. 6). Here, the nanocavity is composed of a silver nanorod above a silver substrate with a gap distance $b=3nm$. The diameter and height of the nanorod are $a=10nm$ and $h=30nm$ respectively. For simplicity, a QE with transition dipole moment $d=72D$ polarized along the z-axis is at the center of the gap. Permittivity for the air is $\varepsilon_{1}=1$. Different from the previous example where Drude model is used for the permittivity of metal over the whole frequency range, permittivity for silver $\varepsilon_{2}$ is from experiment data Palik (1997) which is beyond the Drude model.

The coupling strength $g$ is obtained by COMSOL Multiphysics software with the method in Ref. Zhao et al. (2018a), where the scattering GF is expressed by the difference of the electric fields of an oscillating electric point dipole with and without nanostructre. The real part and imaginary part for the coupling strength $g$ are shown in Fig. 7 (a) and 7(b) respectively. Different from the results shown in Fig. 2(a) and Fig. 2(b) where the coupling strength $g$ is near zero for frequency away from the plasmonic resonance, there are some response for relative large frequency (see the inset in Fig. 7(b)). Thus, to evaluate the third part in our subtractive KK method $g(0)$ (Eq. (12)) $\Delta_{cor}(\omega)=\omega\int_{0}^{\omega_{max}}ds\operatorname{Im}g_{zz}(s)/(\omega+s)s$ , we should choose a relative large cut-off frequency $\omega_{max}$. Figure 7(c) shows $\Delta_{cor}(\omega)$ with $\omega_{max}=50eV$ (black solid line) and $\omega_{max}=10eV$ (red dashed line). Their difference is shown in the inset. We find that their difference is about $10meV$ for $\omega\approx 2eV$. For the second part in our subtractive KK method $g(0)$ (Eq. (12)), linearly extrapolating method is used. Results for $Reg(\omega)$ with $\omega$ in the range $[0.125eV,0.2eV]$ is shown in Fig. 7(d). A liner function $Reg_{i}(\omega)=178.685+2.44152\omega$ is obtained. This result agrees well with that when silver is thought to be perfect conductor.

Figure 8 shows the performances for the time domain method and for frequency domain method. The transition frequency is $\omega_{0}=2.5eV$. Figure 8(a) are results by the frequency domain method (Eq. (7)). The black solid line is for the cut-off frequency $\omega_{max}=50eV$ while the red circle is for $\omega_{max}=10eV$. We find that they agree well (see the inset therein). This means that knowledge about the scattering GF over a narrow frequency range is enough. But for the results by the time domain method (Eq. (3)) shown in Fig. 8(b), we find that results with $\omega_{max}=20eV$ (red square) differs much from that with $\omega_{max}=50eV$ (black solid). This means that $\omega_{max}=20eV$ is not enough for this method. To compare the results by the above two methods, we show the results by both methods with a relative high cut-off frequency $\omega_{max}=50eV$ in Fig. 8(c). They agree with each other. Thus, we can conclude that a much narrow frequency range for the photon GF is enough to get a convergent results by the frequency domain method than by the time domain method.

VI SUMMARY

In summary, we have proposed a general numerical method for calculating the energy level shift of a QE in arbitrary nanostructure. By subtracting the expression for the energy level shift $\Delta(0)$ and using the Kramers-Kronig relations for the scattering photon GF, we have shown that the principal value integral in calculating $\Delta(\omega)$ is transforms into an ordinary integration(see Eq. (12)). We have made numerical comparisons with the method of direct Hilbert transformation over the positive frequency axis (see Eq. (6)) and the method by transferring the integration to the imaginary frequency axis (see Eq. (9)) for emitter located around a gold nanosphere and at the center of a gap plasmonic nonacavity. In the gold nanosphere case, permittivity for metal is supposed to be Drude model, although it can not be extended to a wide frequency range in reality. By the method of integration on the imaginary frequency axis (Eq. (9))), $g(i\xi)$ spread over a wide frequency range and contribution from $\xi>10eV$ is in the order of $meV$. For the method of direct Hilbert method (Eq. (6)), errors are also in the orders of $meV$ with a cut-off frequency $\omega_{max}=10eV$. But for our method, numerical errors is several $\mu eV$ with a cut-off frequency $\omega_{max}=10eV$. This clearly demonstrate that a much narrower frequency range about the scattering GF is enough for our subtractive KK method (Eq. (12)), which is very useful when calculating the GF by numerical means. In addition, we have demonstrated that dynamics by the frequency domain method based on the Greens function expression for the evolution operator sensitively depends on the energy level shift.

For the gap plasmonic nonacavity case, permittivity for the metal is from experiment data which is beyond the Drude model. We have found that the coupling strength $Img(\omega)$ is relatively strong and response can not be ignored in the high frequency range. We have found that the integral part over the frequency range $[10eV,50eV]$ of our subtractive method (Eq. (12)) is about $10meV$ for $\omega$ around $2eV$. Nevertheless, dynamics by the frequency domain method is less affected in the high frequency range. There is no visible difference for the dynamics of the excited-state population between $\omega_{max}=10eV$ and $\omega_{max}=50eV$. Differently, the time domain method in the form of Volterra integral of the second kind is strongly is affected. We have found that cut-off frequency with $\omega_{max}=20eV$ is not enough to get a convergent results. In addition, we have observed that results gradually approached to those by the frequency domain method with the cut-off frequency $\omega_{max}$ increasing.

Acknowledgements.

This work was financially supported by the National Natural Science Foundation of China (Grants No. 11464014, 11347215, 11564013, 11402096, 11464013) and Hunan Provincial Innovation Foundation For Postgraduate (Grants No.CX2018B706).