Signatures of quantized coupling between quantum emitters and localized surface plasmons

TL;DR

This paper investigates the quantized interaction between quantum emitters and localized surface plasmons on a metal nanoparticle, revealing signatures like decay, trapping, and oscillation, which are crucial for quantum device design despite metal losses.

Contribution

It introduces a mechanism for persistent quantum correlations mediated by localized surface plasmons in lossy media, advancing understanding of light-matter interactions at the nanoscale.

Findings

Quantized QE-LSP coupling signatures depend on bound states.

Steady states include decay, trapping, and oscillation.

Results inform quantum device design with plasmonic nanostructures.

Abstract

Confining light to scales beyond the diffraction limit, quantum plasmonics supplies an ideal platform to explore strong light-matter couplings. The light-induced localized surface plasmons (LSPs) on the metal-dielectric interface acting as a quantum bus have wide potential in quantum information processing; however, the loss nature of light in the metal hinders their application. Here we propose a mechanism to make the reversible energy exchange and the multipartite quantum correlation of a collective of quantum emitters (QEs) mediated by the LSPs persistent. Via investigating the quantized interaction between the QEs and the LSPs supported by a spherical metal nanoparticle, we find that the diverse signatures of the quantized QE-LSP coupling in the steady state, including the complete decay, population trapping, and persistent oscillation, are essentially determined by the different number of bound states formed in the energy spectrum of the QE-LSP system. Enriching our understanding on the light-matter interactions in a lossy medium, our result is instructive in the design of quantum devices using plasmonic nanostructures.

Figures

Equations

$$\hat{\mathbf{E}}(\mathbf{r},\omega )=\frac{i{c}^{-2}{\omega }^2}{\sqrt{\pi {\epsilon }_0/\hslash }}\int d^3{\mathbf{r}}^{\prime }\sqrt{\text{Im}[{\epsilon }_{\text{m}}(\omega )]}\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )\cdot \hat{\mathbf{f}}({\mathbf{r}}^{\prime },\omega ),$$

$$\hat{H}$$

$$\dot{\mathbf{c}}(t)+i{\omega }_0\mathbf{c}(t)+{\int }_0^{t}d\tau {\int }_0^{\infty }d\omega e^{-i\omega (t-\tau )}\mathbf{J}(\omega )\mathbf{c}(\tau )=0,$$

$$y_l(\varpi )\equiv {\omega }_0-{\int }_0^{\infty }\frac{D_l(\omega )}{\omega -\varpi }d\omega =\varpi ,\varpi =is.$$

$${\bar{c}}_l(t)=Z_l{e}^{-i{\varpi }_l^{b}t}+{\int }_{iϵ+0}^{iϵ+\infty }\frac{d\varpi }{2\pi }{\bar{c}}_l(-i\varpi )e^{-i\varpi t},$$

$$\underset{t\to \infty }{\lim }|P(t){|}^2=\{\begin{array}{rl}& 0,M=0\\ & Z^2/4,M=1\\ & [Z_0^2+Z_1^2+D(t)]/4,M=2,\end{array}$$

$${\mathbf{G}}^0(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$$

$${\mathbf{G}}^{\text{R}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$$

$${\mathcal{R}}^H=\frac{{\tau }_2\partial {\tau }_1-{\tau }_1\partial {\tau }_2}{{\kappa }_1\partial {\tau }_2-{\tau }_2\partial {\kappa }_1},{\mathcal{R}}^V=\frac{k_1^2{\tau }_1\partial {\tau }_2-k_2^2{\tau }_2\partial {\tau }_1}{k_2^2{\tau }_2\partial {\kappa }_1-k_1^2{\kappa }_1\partial {\tau }_2},$$

$${\mathbf{M}}_{mn}^e(k\mathbf{r})$$

$${\mathbf{M}}_{mn}^o(k\mathbf{r})$$

$${\mathbf{N}}_{mn}^e(k\mathbf{r})$$

$${\mathbf{N}}_{mn}^o(k\mathbf{r})$$

$${\mathbf{G}}_{rr}({\mathbf{r}}_l,{\mathbf{r}}_j,\omega )=-\frac{\delta ({\mathbf{r}}_l-{\mathbf{r}}_j)}{k_1^2}+\frac{i{k}_1}{4\pi }\sum _{n=1}^{\infty }\sum _{m=0}^n{c}_{mn}\cos [\frac{2\pi m(l-j)}{N}]\frac{h_n^{(1)}(k_{1}r)[j_n(k_{1}r)+{\mathcal{R}}^V{h}_n^{(1)}(k_{1}r)]}{[k_{1}r/P_n^m(0){]}^2},$$

$$J_{lj}(\omega )=\frac{3{\gamma }_0{\omega }^3\sqrt{{\epsilon }_{\text{d}}}}{4\pi {\omega }_0^3}\text{Re}\left[\sum _{n=1}^{\infty }\sum _{m=0}^n{c}_{mn}\cos [\frac{2\pi m(l-j)}{N}]\frac{h_n^{(1)}(k_{1}r)[j_n(k_{1}r)+{\mathcal{R}}^V{h}_n^{(1)}(k_{1}r)]}{[k_{1}r/P_n^m(0){]}^2}\right].$$

$$\underset{\rho \to 0}{\lim }{j}_n(\rho )$$

$$\underset{\rho \to 0}{\lim }{\partial }_{\rho }[\rho j_n(\rho )]$$

$$\underset{\rho \to 0}{\lim }{h}_n^{(1)}(\rho )$$

$$\underset{\rho \to 0}{\lim }{\partial }_{\rho }[\rho h_n^{(1)}(\rho )]$$

$${\mathcal{R}}_n^V=\frac{-i(k_{1}R{)}^{2n+1}(n+1)}{(2n+1)!!(2n-1)!!}\frac{{\epsilon }_{\text{d}}-{\epsilon }_{\text{m}}(\omega )}{n{\epsilon }_{\text{m}}(\omega )+(n+1){\epsilon }_{\text{d}}}.$$

$${\mathbf{G}}^{\text{R}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\sum _{n=1}^{\infty }{\mathbf{G}}_n^{\text{R}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega ),$$

$$n{\epsilon }_{\text{m}}({\omega }_n)+(n+1){\epsilon }_{\text{d}}=0,$$

$${\dot{c}}_l(t)=-i{\omega }_l{c}_l(t)-\int d\omega \int d^3\mathbf{r}\frac{c^{-2}{\omega }^2}{\sqrt{\pi {\epsilon }_0\hslash }}\sqrt{\text{Im}[{\epsilon }_m(\omega )]}{\mu }_{l\hat{j}}{G}_{\hat{j}\hat{i}}({\mathbf{r}}_l,\mathbf{r},\omega )d_{\mathbf{r},\omega }(t),$$

$${\dot{d}}_{\mathbf{r},\omega }(t)=-i\omega d_{\mathbf{r},\omega }(t)+\sum _l\frac{c^{-2}{\omega }^2}{\sqrt{\pi {\epsilon }_0\hslash }}\sqrt{\text{Im}[{\epsilon }_m(\omega )]}{\mu }_{l\hat{k}}^{\ast }{G}_{\hat{k}\hat{i}}^{\ast }({\mathbf{r}}_l,\mathbf{r},\omega )c_l(t),$$

$$d_{\mathbf{r},\omega }(t)$$

$${\dot{c}}_l(t)+i{\omega }_l{c}_l(t)+\sum _j{\int }_0^{t}d\tau {\int }_0^{\infty }d\omega e^{-i\omega (t-\tau )}{J}_{lj}(\omega )c_j(\tau )=0,$$

$$E{c}_l$$

$$E{d}_{\mathbf{r},\omega }$$

$$(E-\hslash {\omega }_0)c_l-{\hslash }^2\sum _{j=0}^{N-1}\int d\omega \frac{J_{lj}(\omega )}{E-\hslash \omega }{c}_j=0,$$

$$(E-\hslash {\omega }_0)\mathbf{c}-{\hslash }^2\int \frac{\mathbf{J}(\omega )d\omega }{E-\hslash \omega }\mathbf{c}=0,$$

$$[E-\hslash {\omega }_0-{\hslash }^2\int \frac{\mathbf{D}(\omega )}{E-\hslash \omega }d\omega ]\bar{\mathbf{c}}=0,$$

$$E=\hslash {\omega }_0+{\hslash }^2\int \frac{D_l(\omega )}{E-\hslash \omega }d\omega .$$

$$|\Psi (t)⟩=\sum _{\alpha =1}^M{x}_{\alpha }^b{e}^{-i{\varpi }_{\alpha }^{b}t}|{\Phi }_{\alpha }^b⟩+\sum _{E\in \text{CB}}{x}_E{e}^{-iEt/\hslash }|{\Phi }_E⟩,$$

$$⟨\Psi (0)|\Psi (t)⟩=\sum _{\alpha =1}^M|x_{\alpha }^b{|}^2{e}^{-i{\varpi }_{\alpha }^{b}t}+\sum _{E\in \text{CB}}|x_E{|}^2{e}^{-iEt/\hslash }.$$

$$P(t)$$

Full text

I Introduction

Hybrid systems composed of metal nanoparticles (MNPs) and quantum emitters (QEs) have drawn intense attention in physics, chemistry, and materials and life sciences Giannini et al. (2011); Atwater and Polman (2010); Kabashin et al. (2009); Tame et al. (2013); Lee et al. (2015). By confining light within regions far below the diffraction limit in modes of localized surface plasmons (LSPs), the strong light-matter interaction is realizable in the vicinity of the MNPs Trügler and Hohenester (2008); Wersäll et al. (2017); Baranov et al. (2018); Chikkaraddy et al. (2016); Savasta et al. (2010); Santhosh et al. (2016); Chikkaraddy et al. (2016); Matsuzaki et al. (2017); Kewes et al. (2018); Baranov et al. (2018); Vasa and Lienau (2018). Recently, dramatic progress has been made to reveal the modified radiative properties of QEs by the LSPs in quantum plasmonics. Fascinating effects, including the superradiance of an ensemble of dipoles Pustovit and Shahbazyan (2009), the surface plasmon amplification by stimulated emission of radiation Noginov et al. (2009), the quantum statistics control of photons Ridolfo et al. (2010), and the suppression of quantum fluctuations of light Martín-Cano et al. (2014), have been found. These effects have led to a wide application of the LSPs in quantum information processing and quantum device designing. However, the dissipation of the LSPs induced by the loss nature of light in metal severely restricts their practical applications Tame et al. (2013); Hugall et al. (2018).

It has been found that a QE residing near the metal is quenched by its decay through the nonradiative electromagnetic modes absorbed by the metal Dulkeith et al. (2002); Anger et al. (2006); Kühn et al. (2006); Gurlek et al. (2018). Such quenching hampers the complete quantum control in plasmonic systems, where a persistent quantum coherence is of importance Gurlek et al. (2018); Peng et al. (2017). In the systems of a collective of QEs, the cooperative effect makes the strong coupling between the QEs and the radiative mode dominate the metal absorption Delga et al. (2014a) and suppresses quenching to the QEs Kongsuwan et al. (2018). It endows the multiple-QE system coupled to metal nanostructures with a promising route to suppress the loss of LSPs in metal Sáez-Blázquez et al. (2017, 2018). Going beyond the weak-coupling description of QE-LSP interactions Morshed Behbahani et al. (2016); Thanopulos et al. (2017); Iliopoulos et al. (2018), it has been found that the LSPs can act as a quantum bus to mediate the coherent interactions and generate the entanglement among QEs Nerkararyan and Bozhevolnyi (2015); Hakami and Zubairy (2016); Amooghorban and Aleebrahim (2017). However, such quantum coherence is dynamically transient and tends to vanish in the long-time limit. In terms of practical applications, persistent quantum coherence and entanglement of the QEs are desired. On the other hand, a widely used description of strong QE-LSP coupling is based on the pseudomode method Waks and Sridharan (2010); Delga et al. (2014b); Peng et al. (2017); Delga et al. (2014a); Ge and Hughes (2015); Li et al. (2016), which decomposes the spectrum into a sum of discrete resonant modes with Lorentzian expansion and succeeds in mapping the non-Markovian dynamics into a Markovian one González-Tudela et al. (2014). When the coupling is strong enough, the QEs and LSPs are highly hybridized, and thus the pseudomode method is insufficient and a rigorous continuous-mode theory is needed.

In this paper, going beyond the pseudomode method, we study exactly the dissipative dynamics of a collective of QEs interacting with the LSPs supported by a MNP. A mechanism to overcome the loss effect of the LSPs in the metal is discovered. We find the diverse signatures of the strong QE-LSP couplings, including complete decay, population trapping, and persistent oscillations, in the long-time steady state. Our analyses reveal that they are determined by the formation of different numbers of QE-LSP bound states. We also find that, as a consequence of the suppression of loss effect of the LSPs, a persistent entanglement among the QEs can be mediated by the LSPs. Such bound-state-favored persistent entanglement among the QEs plays a constructive role in applying the LSPs as a quantum bus in quantum information processing.

II System and quantization

The system is composed of a MNP surrounded by $N$ QEs. The QEs labeled by $l$ are positioned at $\mathbf{r}_{l}$ on the equator plane of the MNP (see Fig. 1). Each QE is modeled as a two-level system with frequency $\omega_{l}$ and dipole moment $\boldsymbol{\mu}_{l}$. The MNP has a radius $R$ and a dielectric permittivity denoted by a complex Drude model $\varepsilon_{\text{m}}(\omega)=\varepsilon_{\infty}-\omega_{p}^{2}/[\omega(\omega+i\gamma_{p})]$, where $\omega_{p}$ is the bulk plasma frequency, $\varepsilon_{\infty}$ is the high-frequency limit of $\varepsilon_{\text{m}}(\omega)$, and $\gamma_{p}$ is the Ohmic loss of light in the MNP Johnson and Christy (1972). The whole system is embedded in a homogeneous medium with dielectric constant $\varepsilon_{\text{d}}$. We consider that both of the dielectric and the metal are nonmagnetic and thus their permeability $\mu_{\text{d}}=\mu_{\text{m}}\equiv 1$.

Besides propagating into the dielectric as a radiative mode and being absorbed by the MNP as a nonradiative mode, the optical field emitted by the QE also induces a confined hybrid mode which consists of LSPs localized near the metal-dielectric interface Pitarke et al. (2007). The LSPs enable a confinement of light within the subwavelength areas on the interface, which supplies an ideal platform to explore the strong quantized light-matter coupling González-Tudela et al. (2014); Yang and An (2017). A quantization method of light in the absorbing medium has been proposed based on the dyadic Green’s function, where the absorption of the medium to light is described by a Langevin noise Gruner and Welsch (1996); Dung et al. (1998). Then the electric field reads

[TABLE]

where $\varepsilon_{0}$ is the vacuum permittivity, $c$ is the speed of light, and $\hat{\bf{f}}(\bf{r},\omega)$ satisfying $[\hat{\bf{f}}(\bf{r},\omega),\hat{\bf{f}}^{\dagger}(\bf{r}^{\prime},\omega^{\prime})]=\delta(\bf{r}-\bf{r}^{\prime})\delta({\omega-\omega^{\prime}})$ is the annihilation operator of light. The Green’s function $\bf{G(\bf{r},\bf{r}^{\prime},\omega)}$ satisfying the Helmholtz equation $[{\boldsymbol{\nabla}}\times{\boldsymbol{\nabla}}\times-\omega^{2}c^{-2}\varepsilon_{\text{m}}(\omega)]\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\mathbf{I}\delta(\mathbf{r}-\mathbf{r}^{\prime})$, with $\mathbf{I}$ being the identity matrix, denotes the field in frequency $\omega$ evaluated at $\bf{r}$ due to a point source at $\bf{r}^{\prime}$. The spatial distribution of all of the three modes has been incorporated in $\bf{G(\bf{r},\bf{r}^{\prime},\omega)}$ by solving the Helmholtz equation subject to the boundary condition of the system geometry. It allows for a complete description of the quantized light-matter coupling by calculating $\bf{G(\bf{r},\bf{r}^{\prime},\omega)}$. For a spherical MNP, the Green’s function is analytically solvable. For more details see Appendix A.

The Hamiltonian of the full QE-MNP system under the dipole and rotating-wave approximations reads Dung et al. (2000)

[TABLE]

where $\hat{\sigma}_{l}=|g_{l}\rangle\langle e_{l}|$ is the transition operator from the excited state $|e_{l}\rangle$ to the ground state $|g_{l}\rangle$ of the $l$th QE. The validity of the rotating-wave approximation in a related system was revealed in Iliopoulos et al. (2018). The dipole approximation works when the QE size is sufficiently small Gonzalez-Tudela et al. (2010); Stobbe et al. (2012); Yang and An (2016). Conventionally, the LSPs are viewed as a few discrete pseudomodes with Lorentzian expansion. Then one can use the standard cavity QED method to describe the QE-LSP coupling Waks and Sridharan (2010); Delga et al. (2014b); Peng et al. (2017); Delga et al. (2014a); Ge and Hughes (2015); Li et al. (2016). It neglects the non-Lorentzian features of the spectrum and may be insufficient when the QE is close to the interface Yang and An (2017), where the hybridization of the QEs and the LSPs dominates.

III Exact dynamics

We can see that the total excitation number $\hat{\mathcal{N}}=\sum_{l}\hat{\sigma}^{{\dagger}}_{l}\hat{\sigma}_{l}+\int d^{3}\mathbf{r}\int d\omega\mathbf{\hat{f}}^{{\dagger}}(\mathbf{r},\omega)\cdot\mathbf{\hat{f}}(\mathbf{r},\omega)$ is conserved. In the single-excitation subspace, the time-evolved state can be expanded as $|\Psi(t)\rangle=[\sum_{l}c_{l}(t)\hat{\sigma}_{l}^{\dagger}+\int d^{3}\mathbf{r}\int d\omega d_{\mathbf{r},\omega}(t)\hat{\mathbf{f}}^{\dagger}(\mathbf{r},\omega)]|G;\{0_{\omega}\}\rangle$, where $|G\rangle$ denotes all the QEs in the ground state and $|\{0_{\omega}\}\rangle$ is the vacuum state of the total modes. It can be derived that $c_{l}(t)$ obeys (see Appendix B)

[TABLE]

where $\mathbf{c}(t)=(c_{0}(t),\ldots,c_{N-1}(t))^{T}$ is a column vector with $c_{l}(t)$ being the excited-state probability amplitude of $l$th QE, and $\mathbf{J}(\omega)$ is a matrix, with $J_{lj}(\omega)=\omega^{2}\boldsymbol{\mu}_{l}\cdot\textrm{Im}[\mathbf{G}(\mathbf{r}_{l},\mathbf{r}_{j},\omega)]\cdot\boldsymbol{\mu}_{j}^{\ast}/(\pi\hbar\varepsilon_{0}c^{2})$ the correlated spectral densities between the $l$th and $j$th QEs. Thus all the actions of the metal-dielectric structure on the QEs have been collected in $\mathbf{J}(\omega)$. We have chosen the QEs having identical frequency $\omega_{l}=\omega_{0}$ and used $\int d^{3}\mathbf{s}\frac{\omega^{2}}{c^{2}}$Im$[\varepsilon_{m}(\omega)]\mathbf{G}(\mathbf{r},\mathbf{s},\omega)\mathbf{G}^{\ast}(\mathbf{r}^{\prime},\mathbf{s},\omega)=\textrm{Im}[\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)]$ Dung et al. (2000). The convolution in Eq. (2) renders the QE dynamics non-Markovian. The correlation of different $c_{l}(t)$ indicates that, although direct couplings of QEs in Eq. (1) are absent, their indirect couplings can be effectively induced by exchanging the virtual excitations of the photons.

The solution of Eq. (2) can be analyzed by a Laplace transform, which yields $\tilde{\mathbf{c}}(s)=\mathbf{V}\bar{\mathbf{c}}(s)\mathbf{V}^{-1}\mathbf{c}(0)$, with $\bar{\mathbf{c}}(s)=[s+i\omega_{0}+\int_{0}^{\infty}d\omega\frac{\mathbf{D}(\omega)}{s+i\omega}]^{-1}$. We have used the Jordan decomposition of $\mathbf{J}(\omega)=\mathbf{V}\mathbf{D}(\omega)\mathbf{V}^{-1}$, with $\mathbf{V}$ and $\mathbf{D}(\omega)=\text{diag}[D_{0}(\omega),\ldots,D_{N-1}(\omega)]$ its similarity matrix and Jordan canonical form, respectively. Then $\mathbf{c}(t)$ is obtainable by in inverse Laplace transform to $\bar{\mathbf{c}}(s)$, which can be done by finding its poles from

[TABLE]

It can be proven that the roots $\varpi$ multiplied by $\hbar$ are just the hybrid eigenenergies of the QEs and the LSPs in the single-excitation subspace (see Appendix C). Since $y_{l}(\varpi)$ is a monotonically decreasing function when $\varpi<0$, each one of Eqs. (3) has one discrete root $\varpi^{b}_{l}$ if $y_{l}(0)<0$. We call the discrete eigenstates with eigenenergy $\hbar\varpi_{l}^{b}$ the bound state. In the region $\varpi>0$, it has an infinite number of roots, which form a continuous energy band. Determined by the system parameters, at most $N$ independent bound states could be formed. Using Cauchy’s residue theorem, we readily have $\mathbf{c}(t)=\mathbf{V}\bar{\mathbf{c}}(t)\mathbf{V}^{-1}\mathbf{c}(0)$ with the elements of $\bar{\mathbf{c}}(t)$ as Yang and An (2017)

[TABLE]

where the first term with $Z_{l}=[1+\int_{0}^{\infty}\frac{D_{l}(\omega)}{(\varpi_{l}^{b}-\omega)^{2}}d\omega]^{-1}$ is from the bound state and the second term is from the energy band. Oscillating with time in continuously changing frequencies, the second term behaves as a decay and tends to zero due to out-of-phase interference. Thus, if the bound state is absent, then $\lim_{t\rightarrow\infty}\mathbf{c}(t)=0$ characterizes a complete decay; while if the bound states are formed, then $\lim_{t\rightarrow\infty}\mathbf{c}(t)=\mathbf{V(Z}e^{-i\boldsymbol{\varpi}^{b}t}\mathbf{)V}^{-1}\mathbf{c}(0)$, with $\bf{x}=\text{diag}(x_{0},\ldots,x_{N-1})$ for $\bf{x}=\bf{Z}$ and $\boldsymbol{\varpi}^{b}$, implies decoherence suppression. This indicates that the dynamics of the QEs in the long-time limit is intrinsically determined by the energy-spectrum characters of the whole QE-LSP system. Generally, solving $\mathbf{V}$ and $\mathbf{D}(\omega)$ needs numerical calculations. Here, for concreteness, we choose that all the QEs have identical dipole moments and uniform coordinates $\mathbf{r}_{l}=(r,\pi/2,2\pi l/N)$ such that $\mathbf{J}(\omega)$ is a symmetric circulant matrix with $J_{lj}(\omega)=J_{mn}(\omega)\equiv J_{|l-j|}(\omega)$ for $|l-j|=|m-n|$ (see Appendix A). Because $\mathbf{J}(\omega)$ is a symmetric circulant matrix, we readily have $D_{l}(\omega)=\sum_{j=0}^{N-1}J_{j}(\omega)\lambda_{l}^{N-j}$ and $\mathbf{V}=(\upsilon_{0},\ldots,\upsilon_{N-1})$ with $\upsilon_{l}=\frac{1}{\sqrt{N}}(1,\lambda_{l},\ldots,\lambda_{l}^{N-1})^{T}$ and $\lambda_{l}=\exp(-2\pi il/N)$ Gray (2006).

IV Results and discussion

It was previously found that the reversible energy exchange between the QEs induced by a common surface plasma tends to vanish in the long-time limit under the Born-Markovian approximation Gonzalez-Tudela et al. (2011). Different from that result, we will show that such mediated coherent coupling can induce a persistently reversible energy exchange between the QEs even in the steady state when the approximation is relaxed. We choose silver for the metal with $\hbar\omega_{p}=9.01$ eV, $\varepsilon_{\infty}=3.718$, and $\hbar\gamma_{p}=0.09$ eV in the interested frequency range Scholl et al. (2012) and the QEs with $\hbar\gamma_{0}=0.1$ meV. We focus on the QE dynamics by studying the initial-state fidelity $P(t)=|\langle\Psi(0)|\Psi(t)\rangle|^{2}$.

First, taking $N=2$, we consider that only one of the QEs is excited initially, i.e., $|\Psi(0)\rangle=\hat{\sigma}_{0}^{{\dagger}}|G;\{0_{\omega}\}\rangle$. We can calculate that with time evolution the fidelity reads $P(t)=|c_{0}(t)|^{2}$. Figure 2(a) shows the evolution of $P(t)$ in three characteristic values of $r$. As a result of the near-field enhancement of the LSPs, a significant oscillation appears in the dynamics for all three cases. Absent in the Born-Markovian approximate result, this is entirely the non-Markovian effect, which represents a reversible energy exchange and thus manifests the strong coupling between the QEs mediated by the LSPs González-Tudela et al. (2014). It is interesting to see that the non-Markovian effect manifests its action on the QEs not only in its transient dynamics, but also in its steady state. When $r=9.5$ nm, $P(t)$ tends to zero accompanying the QEs decay completely to the ground state, which is consistent with the previous results He and Zhu (2012); Otten et al. (2015). However, a remarkable difference appears with further decreasing $r$. One can see that $P(t)$ tends to a nonzero value when $r=9.0$ nm, which represents a stable population trapping in the system, while when $r=8.0$ nm, $P(t)$ tends to a lossless oscillation with a constant frequency, which is quite like the Rabi oscillation Dudin et al. (2012) and represents a persistent energy exchange among QEs caused by the QE-LSP interaction. These diverse signatures can be explained by our bound-state analysis. From Eq. (4) we have (see Appendix D)

[TABLE]

where $M$ is the number of formed bound states and $D(t)=2Z_{0}Z_{1}\cos[(\varpi_{1}^{b}-\varpi_{0}^{b})t]$ is the interference between the two bound states. This conclusion can be confirmed by the energy spectrum shown in Fig. 2(b). The two branches of bound states formed in the band gap divide the spectrum into three regions: without bound state when $r\gtrsim 9.0$ nm, one bound state when $8.5\lesssim r\lesssim 9.0$ nm, and two bound states when $r\lesssim 8.5$ nm. The regions match well with the ones where $P(\infty)$ shows different behaviors [see Fig. 2(c)], i.e., complete decay, population trapping, and persistent oscillation, as expected from Eq. (5). Such bound-state-favored behaviors are constructive to generate entanglement between the QEs. Different from the asymptotic vanishing in the Born-Markovian approximation Gonzalez-Tudela et al. (2011) and in the absence of the bound state, the generated entanglement can be preserved as long as the bound states are formed [see Fig. 2(d)]. This is helpful for utilizing plasmonic nanostructures in designing quantum devices. Our results can be generalized to the case of a large number of QEs. With more of the bound states being formed in the large-$N$ case, the persistent oscillations will be complicated, but the mechanism is the same as in the present case. In Appendix E, the dynamics for $N=4$ is provided. Note that the similar bound-state-induced decoherence suppression for the single-QE case has been found in Refs. Thanopulos et al. (2017); Yang and An (2017).

Next we consider that the QEs are initially in a $|\Psi(0)\rangle=\frac{1}{\sqrt{N}}\sum_{l=0}^{N-1}\hat{\sigma}^{{\dagger}}_{l}|G;\{0_{\omega}\}\rangle$, which is a multipartite entangled state widely used in quantum information processing Joo et al. (2003); Gorbachev et al. (2003). The canonical transformation $\mathbf{V}$ can convert Eq. (2) into $\dot{\bar{\mathbf{c}}}(t)+i\omega_{0}\bar{\mathbf{c}}(t)+\int_{0}^{t}d\tau\int d\omega e^{-i\omega(t-\tau)}\mathbf{D}(\omega)\bar{\mathbf{c}}(\tau)=0$ with $\bar{\mathbf{c}}(t)\equiv\mathbf{V}^{-1}\mathbf{c}(t)$. Its initial condition can be calculated as $\bar{\mathbf{c}}(0)=(1,0,\ldots,0)$, under which only the $\bar{c}_{0}(t)$ component of this matrix equation has a nonzero solution. Thus its dynamics has the same equation of motion as the one of a single QE coupled to the LSPs Yang and An (2017). This indicates that the $N$ QEs collectively act as a two-level superatom to interact with the LSPs with the spectral density characterized by $D_{0}(\omega)$. This notion of a superatom is a powerful concept in designing single-photon quantum sources Saffman et al. (2010); Paris-Mandoki et al. (2017); Zeiher et al. (2015); Kumlin et al. (2018). We can calculate the initial-state fidelity $P(t)=|\bar{c}_{0}(t)|^{2}$.

In the same mechanism as the case of $N=2$, the entanglement of the QEs can be preserved in the steady state due to the formation of the bound state. Figure 3(a) shows the evolution of $P(t)$ for a different number $N$ of QEs. It shows that $P(t)$ tends to a finite value for large $N$, where the QEs remain entangled. It can be understood from the bound-state analysis. As discussed above, we readily obtain $\lim_{t\rightarrow\infty}\bar{c}_{0}(t)=Z_{0}e^{-i\varpi_{0}^{b}t}$ when Eq. (3), with $l=0$, has an isolate root in the region $\varpi<0$. Figures 3(b) and 3(c) show that the region where $P(t)$ tends to a stable value matches well with the one where a bound state is formed in the energy spectrum of the whole system. This verifies again our conclusion that it is the formation of a bound state that preserves the entanglement in the steady state. We also plot in Fig. 3(d) the spectral density $D_{0}(\omega)$, which measures the coupling strength of the QEs and the LSPs. We can see that the contribution of the resonant dipole mode $\omega_{1}=3.77$ eV is entirely canceled, while the one of the quadrupole mode $\omega_{2}=3.94$ eV is enhanced by increasing $N$ (see Appendix A). This is due to the destructive interference of the undistinguished coupling channels between different QEs and the LSPs Delga et al. (2014a); Kumlin et al. (2018); Xiang et al. (2016).

We note that, although we consider only the case that the dipole moments of the QEs are polarized along the radial direction, our result can be generalized to other cases. Some quantitative difference might occur, but the constructive role played by the bound states in overcoming the loss effect of the LSPs in the MNP does not change. We emphasize that our finding is realizable in the state-of-art technique of experiments. The parameters used in our calculation are near the ones of silver as the MNP and the $J$ aggregates as the QEs. Their strong coupling has been studied Schlather et al. (2013); Sugawara et al. (2006); Fidder et al. (1990); Aberra Guebrou et al. (2012). The bound state and its distinguished role in the non-Markovian dissipative dynamics have recently been observed in both photonic crystal Liu and Houck (2017) and ultracold-atom systems Krinner et al. (2018). This means that our finding is completely realizable in quantum plasmonics system, where the strong light-matter coupling is more manifest than in other systems.

V Conclusion

We have proposed a mechanism to overcome the loss effect of LSPs in metal by investigating the exact dynamics of $N$ QEs coupled to LSPs supported by a MNP. It has been found that, in sharp contrast to the previous approximate result that the reversible energy exchange and the entanglement of the QEs mediated by the LSPs exclusively tends to vanish due to the loss effect of LSPs in metal, the persistent quantum coherence and entanglement can be established among the QEs by the LSPs. Our analysis indicates that it is the formation of hybrid bound states in the energy spectrum of the QE-LSP system that governs this lossless behavior. Such bound-state-assisted behavior is helpful in the application of LSPs as a quantum bus. The further study of the multipartite $W$-class state demonstrates the collective suppression of the resonant dipole mode and the enhancement of the quadrupole mode in the QE-LSP coupling. Within the present experimental state of the art, our finding supplies a guideline for experiments to design quantum devices using the plasmonic nanostructures.

VI Acknowledgments

The work was supported by the Natural Science Foundation of China (Grants No. 11704103, No. 11875150, and No. 11834005), by the Doctoral Scientific Research Foundation of Henan Normal University (Grant No. 5101029170296), and by the Fundamental Research Funds for the Central Universities of China.

Appendix A Green’s function of the spherical metal nanoparticle

In this appendix, we give the derivations of the Green’s function of the spherical MNP in calculating the exact dynamics of QEs coupled to LSPs.

Given a spherical MNP with permittivity $\varepsilon_{\text{m}}(\omega)$ and radius $R$ embedded in a homogeneous medium with dielectric constant $\varepsilon_{\text{d}}$, the Green’s functions contributed by the free-space radiation sources and by the MNP-QE interaction are given by Hakami and Zubairy (2016); Van Vlack et al. (2012); Hakami et al. (2014)

[TABLE]

where $\mathcal{R}^{H}$ and $\mathcal{R}^{V}$ are the scattering coefficients corresponding to the transverse electric field $\mathbf{M}_{mn_{o}^{e}}$ and the transverse magnetic field $\mathbf{N}_{mn_{o}^{e}}$ with even and odd contributions. According to the boundary conditions at the surface, $\mathcal{R}^{H}$ and $\mathcal{R}^{V}$ are given by

[TABLE]

where $\tau_{i}=j_{n}(k_{i}R)$, $\kappa_{i}=h_{n}^{(1)}(k_{i}R)$, $\partial\tau_{i}=\partial_{\rho}[\rho j_{n}(\rho)]_{\rho=k_{i}R}$, and $\partial\kappa_{i}=\partial_{\rho}[\rho h_{n}^{(1)}(\rho)]_{\rho=k_{i}R}$. Here, $j_{n}(x)$ and $h_{n}^{(1)}(x)$ are the spherical Bessel functions and the Hankel functions of the first kind, respectively, with $k_{1}=\omega\sqrt{\varepsilon_{\text{d}}}/c$ and $k_{2}=\omega\sqrt{\varepsilon_{\text{m}}(\omega)}/c$ the wave vectors in the dielectric and the metal. The vector functions in spherical coordinates are defined as

[TABLE]

where $P_{n}^{m}(x)$ are the associated Legendre polynomials. In Eqs. (8) and (9), the superscript $(1)$ denotes that $j_{n}(x)$ has to be replaced by $h_{n}^{(1)}(x)$.

In the case that the dipole moments of the QEs are polarized along the radial direction, only the $rr$ component of the Green’s function contributes to the interactions. In the structure studied, the QEs labeled by $l$ are located at $\mathbf{r}_{l}=(r,\pi/2,2\pi l/N)$, with $l=0,\ldots,N-1$. We obtain

[TABLE]

where $c_{mn}=(2-\delta_{0m})n(n+1)(2n+1)(n-m)!/(n+m)!$ and contributions from both the free-space field and the scattered field have been incorporated. From the definition, the spectral density characterizing the coupling strength between QEs and LSPs can be calculated as

[TABLE]

Defining $J_{lj}(\omega)\equiv J_{|l-j|}(\omega)$, we can verify that the spectral densities are periodic with $J_{l}(\omega)=J_{N-l}(\omega)$.

When the radius of the MNP is very small compared to the wavelength, i.e., $|k_{2}R|\ll 1,|k_{1}R|\ll 1,$, the Green’s function can be further simplified. Substituting the limits

[TABLE]

into Eq. (10), we can readily have $\mathcal{R}^{H}=0$ and $\mathcal{R}^{V}=\sum_{n=1}^{\infty}\mathcal{R}_{n}^{V}$, with

[TABLE]

Then the scattered Green’s function can be decomposed into

[TABLE]

with the scattering coefficient in Eq. (9) replaced by $\mathcal{R}_{n}^{V}$. The poles of $\mathcal{R}_{n}^{V}$ determines the resonance frequency of the LSPs. In this manner, the LSPs are expressed as a series of resonant modes labeled by $n$ with eigenfrequency determined by

[TABLE]

from which the contributions of the different resonant modes of LSPs to the light-matter interaction can be studied. In the low-frequency condition, the resonant frequencies can be determined by $\text{Re}[\varepsilon_{\text{m}}(\omega_{n})]=-(n+1)\varepsilon_{\text{d}}/n$ due to $\textrm{Re}[\varepsilon_{m}(\omega)]\gg\textrm{Im}[\varepsilon_{m}(\omega)]$ Johnson and Christy (1972). The first resonant mode is called the dipole mode and the second one is the quadrupole mode Delga et al. (2014a). Using the parameters in our system, we can calculated the frequencies of the dipole and quadrupole modes $\omega_{1}=3.77$ eV and $\omega_{2}=3.94$ eV.

Appendix B Derivation of the evolution equations

In this appendix, we give the derivation of Eq. (2). In the single-excitation subspace, the time-evolved state can be expanded as $|\Psi(t)\rangle=[\sum_{l}c_{l}(t)\hat{\sigma}_{l}^{\dagger}+\int d^{3}\mathbf{r}\int d\omega d_{\mathbf{r},\omega}(t)\hat{\mathbf{f}}^{\dagger}(\mathbf{r},\omega)]|G;\{0_{\omega}\}\rangle$. According to the Schrödinger equation $i\hbar|\dot{\Psi}(t)\rangle=\hat{H}|\Psi(t)\rangle$, we have

[TABLE]

with $l,j=0,\ldots,N-1$ and $\hat{i},\hat{j},\hat{k}=x,y,z$. Using $d_{\mathbf{r},\omega}(0)=0$, Eq. (25) can be formally solved as

[TABLE]

Substituting this into Eq. (24), we have

[TABLE]

with $J_{lj}(\omega)=\omega^{2}\boldsymbol{\mu}_{l}\cdot\textrm{Im}[\mathbf{G}(\mathbf{r}_{l},\mathbf{r}_{j},\omega)]\cdot\boldsymbol{\mu}_{j}^{\ast}/(\pi\hbar\varepsilon_{0}c^{2})$ being the spectral density. We have used $\int d^{3}\mathbf{s}\frac{\omega^{2}}{c^{2}}$Im$[\varepsilon_{m}(\omega)]\mathbf{G}(\mathbf{r},\mathbf{s},\omega)\mathbf{G}^{\ast}(\mathbf{r}^{\prime},\mathbf{s},\omega)=\textrm{Im}[\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)]$ Dung et al. (2000).

For simplify, we choose the QEs having identical frequency $\omega_{l}=\omega_{0}$. Introducing a column vector $\mathbf{c}(t)=(c_{0}(t),c_{1}(t),\ldots,c_{N-1}(t))^{T}$ and a spectral density matrix ${\bf J}(\omega)=J_{lj}(\omega)$, we obtain the evolution equation as Eq. (2).

Appendix C Eigenenergies of the system

In this appendix, we give the derivation of the energy spectrum of the whole system in the single-excitation subspace and the proof that they are exactly the same as the poles in the evolution equation under the Laplace transform in the main text.

The eigenstate $|\Phi\rangle$ of the QE-LSP system in the single-excitation subspace can be expanded as $|\Phi\rangle=[\sum_{l=0}^{N-1}c_{l}\hat{\sigma}_{l}^{{\dagger}}+\int d^{3}\mathbf{r}\int d\omega d_{\mathbf{r},\omega}\hat{\mathbf{f}}^{{\dagger}}(\mathbf{r},\omega)]|G;\{0_{\omega}\}\rangle$. According to the stationary Schrödinger equation $\hat{H}|\Phi\rangle=E|\Phi\rangle$, with $E$ the eigenenergy, we have

[TABLE]

with $l,j=0,\ldots,N-1$ and $\hat{i},\hat{j},\hat{k}=x,y,z$. Solving $d_{\mathbf{r},\omega}$ and substituting it into Eq. (28), it is easy to obtain

[TABLE]

or

[TABLE]

expressed in a matrix form. Using the Jordan decomposition of $\mathbf{J}(\omega)=\mathbf{V}\mathbf{D}(\omega)\mathbf{V}^{-1}$, with $\mathbf{V}$ and $\mathbf{D}(\omega)=\text{diag}[D_{0}(\omega),\ldots,D_{N-1}(\omega)]$ its similarity matrix and Jordan canonical form, Eq. (31) can be expressed as

[TABLE]

where $\mathbf{\bar{c}}=\mathbf{V}^{-1}\mathbf{c}$. The equations have nontrivial solutions if and only if the determinant of the coefficient matrix is zero. Therefore, the eigenvalues of the QE-LSP system in the single-excitation subspace are determined by

[TABLE]

Equation (33) takes the same form as the equation to determine the bound state obtained in the main text. This clearly demonstrates that the dynamics of QEs essentially depends on the energy-spectrum character of the whole QE-LSP system.

Appendix D Solution in the steady state

Suppose $M$ bound states form outside the continuous energy band. Then, according to the completeness of the eigenstates, the time evolution of $|\Psi(0)\rangle=\hat{\sigma}_{0}^{\dagger}|G;\{0_{\omega}\}\rangle$ can be expanded as

[TABLE]

where $x^{b}_{\alpha}=\langle\Phi_{\alpha}^{b}|\Psi(0)\rangle$ and $x_{E}=\langle\Phi_{E}|\Psi(0)\rangle$. The first term is contributed by the bound eigenstates and the second one is from the continuous-band eigenstates. Its overlap with the initial state reads

[TABLE]

The initial-state fidelity defined as $P(t)=|\langle\Psi(0)|\Psi(t)\rangle|^{2}$ reads

[TABLE]

Both the second and the third term contain the oscillating frequencies $E/\hbar$, which are continuously summed in the continuous energy band. Such terms in the continuous energy band tend to vanish due to the out-of-phase interference of the different components in the long-time limit. Thus only the isolated bound states survive in the long-time limit. For the $N=2$ case, at most two bound states can be formed. Therefore, we have

[TABLE]

with $F(t)=2|x^{b}_{1}|^{2}|x^{b}_{2}|^{2}\cos(\varpi^{b}_{1}-\varpi^{b}_{0})t$.

To determine $x^{b}_{\alpha}$, we solve the stationary Schrödinger equation $\hat{H}|\Phi\rangle=E|\Phi\rangle$. From Eqs. (28)$-$(30) we have

[TABLE]

Equation (40) leads to

[TABLE]

Equations (38) and (39) have nontrivial solutions if and only if

[TABLE]

If $y_{\pm}(0)<0$, two bound states with the eigenenergies $E^{b}=\hbar\varpi^{b}_{+}$ and $\hbar\varpi^{b}_{-}$ determined by Eq. (42) can be formed in the band-gap area. Focusing on these bound states, we calculate their corresponding excited-state populations $|c_{0,\pm}^{b}|^{2}$ in the first QE. Substituting Eqs. (38), (39), and (41) into the normalization condition $\sum_{l=0}^{1}|c_{l}|^{2}+\int d^{3}\mathbf{r}\int_{0}^{\infty}d\omega|d_{\mathbf{r},\omega}|^{2}=1$ and repeatedly using $E^{b}-\hbar\omega_{0}-\int_{0}^{\infty}d\omega{\hbar^{2}J_{0}(\omega)\over E^{b}-\hbar\omega}=\pm\int_{0}^{\infty}d\omega{\hbar^{2}J_{1}(\omega)\over E^{b}-\hbar\omega}$ obtained from Eq. (42) for the bound states, we obtain

[TABLE]

It can be verified that $x_{1,2}^{b}=c_{0,\pm}^{b\ast}$. According to the forms of $Z_{l}$ obtained by the Laplace transform in the main text, we can readily see that $|x_{1,2}^{b}|^{2}={Z_{0,1}\over 2}$.

The above process gives the analytical proof to Eq. (5) from the bound states. From this proof we can clearly see the distinguished role of the formed bound states in lossless steady-state behaviors. Such suppression to the decay in the lossy medium is guaranteed by the characters of the bound states as stationary states with isolated eigenenergies of the whole system.

Appendix E Exact Dynamics for $N=4$

The spectral density matrix for $N=4$ reads

[TABLE]

where the periodic condition $J_{l}(\omega)=J_{N-l}(\omega)$ has been used. As a symmetric and circulant matrix, $\mathbf{J}(\omega)=\mathbf{V}\mathbf{D}(\omega)\mathbf{V}^{-1}$, where $\mathbf{D}(\omega)={\text{diag}}[J_{0}(\omega)+2J_{1}(\omega)+J_{2}(\omega),J_{0}(\omega)-J_{2}(\omega),J_{0}(\omega)-2J_{1}(\omega)+J_{2}(\omega),J_{0}(\omega)-J_{2}(\omega)]$ and \mathbf{V}=\frac{1}{2}\left[\begin{array}[]{cccc}1&1&1&1\\ 1&i&-1&-i\\ 1&-1&1&-1\\ 1&-i&-1&i\end{array}\right] Gray (2006). Note that the eigenvalues are degenerate with $D_{1}(\omega)=D_{3}(\omega)$. The canonical transform $\mathbf{\bar{c}}(t)=\mathbf{V}^{-1}\mathbf{c}(t)$ can convert the integro-differential equation into

[TABLE]

with

[TABLE]

and $\mathbf{\bar{c}}(0)=\frac{1}{2}\left[1,1,1,0\right]^{T}$ under the initial condition $|\Psi(0)\rangle=\hat{\sigma}_{0}^{{\dagger}}|G;\{0_{\omega}\}\rangle$, where $c_{1}(t)=c_{3}(t)$ has been used. Equation (45) is analytically solvable by the Laplace transform. As shown in the main text, its solution in the long-time limit reads

[TABLE]

This clearly shows that the dynamics of the system in the long-time limit is determined by the formation of a bound of the whole system. It is not easy to find that the initial-state fidelity equals $P(t)=|c_{0}(t)|^{2}$, with $c_{0}(t)=\frac{1}{4}[\bar{c}_{0}(t)+2\bar{c}_{1}(t)+\bar{c}_{2}(t)]$.

Figure 4(a) plots the evolution of $P(t)$ in different $r$. The different behaviors, i.e., complete decay, population trapping, and persistent oscillation, are present depending on the value of $r$. Details on the long-time behaviors are shown in Fig. 4(b). Such phenomena are associated with the formation of the bound state of the QE-LSP system. Figure 4 (c) shows the energy spectrum of the whole system. If no bound state is formed, then $P(t)$ tends to zero, which characterizes the complete decoherence. If one bound state is formed, then $P(t)$ tends to a finite value, which describes the population trapping. If two or more bound states are formed, then $P(t)$ tends to the Rabi-like persistent oscillations in the long-time limit. Such behaviors coincide with our analytical analysis.