Quantum dynamics of two quantum dots coupled through localized plasmons: An intuitive and accurate quantum optics approach using quasinormal modes
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov Quantum dynamics of two quantum dots coupled through localized plasmons:An intuitive and accurate quantum optics approach using quasinormal modes
Rong-Chun Ge ∗ and Stephen Hughes † Department of Physics, Engineering Physics and Astronomy,Queens University, Kingston, Ontario, Canada K7L 3N6
We study the quantum dynamics of two quantum dots (QDs) or artificial atoms coupled throughthe fundamental localized plasmon of a gold nanorod resonator. We derive an intuitive and efficienttime-local master equation, in which the effect of the metal nanorod is taken into consideration self-consistently using a quasinormal mode (QNM) expansion technique of the photon Green function.Our efficient QNM technique offers an alternative and more powerful approach over the standardJaynes-Cummings model, where the radiative decay, nonradiative decay, and spectral reshaping ef-fect of the electromagnetic environment is rigorously included in a clear and transparent way. Wealso show how one can use our approach to compliment the approximate Jaynes-Cummings modelin certain spatial regimes where it is deemed to be valid. We then present a study of the quantumdynamics and photoluminescence spectra of the two plasmon-coupled QDs. We first explore thenon-Markovian regime, which is found to be important only on the ultrashort time scale of the plas-mon mode which is about 40 fs. For the field free evolution case of excited QDs near the nanorod, wedemonstrate how spatially separated QDs can be effectively coupled through the plasmon resonanceand we show how frequencies away from the plasmon resonance can be more effective for coherentlycoupling the QDs. Despite the strong inherent dissipation of gold nanoresonators, we show thatqubit entanglements as large as 0.7 can be achieved from an initially separate state, which has beenlimited to less than 0.5 in previous work for weakly coupled reservoirs. We also study the superra-diance and subradiance decay dynamics of the QD pair. Finally, we investigate the rich quantumdynamics of QDs that are incoherently pumped, and study the polarization dependent behaviourof the emitted photoluminescence spectrum where a double-resonance structure is observed due tothe strong photon exchange interactions. Our general quantum plasmonics formalism can easilybe extended to include multiple QDs interacting through the QNMs of metallic resonator struc-tures, fully accounting for radiative and non-radiative coupling, as well as nonlinear light-matterinteraction processes.
PACS numbers: 42.50.Nn, 73.20.Mf, 78.67.Bf
I. INTRODUCTION
Quantum nanophotonics is an active field of research,which is driven in part by fundamental ideas in light-plasmon-matter interactions, applications in nanopho-tonics, and by continued advances in nanofabricationtechnologies. In particular, various types of nanoscalestructures have been designed and fabricated to manipu-late the interaction between quantum emitters and localelectric fields, which can be enhanced by tailoring thelocal density of optical states (LDOS) typically throughsome discrete cavity resonance . For coupling to electricfields below the diffraction limit, metallic nanoparticles(MNPs) have been shown to yield an unprecedented con-finement of light at the nm scale due to the spatial struc-ture of the localized surface plasmon (LSP) resonances.For quantum dot (QD) emitters or artificial atoms placedsufficiently close to the MNP, the strong coupling regimehas also been shown to be experimentally accessible .As a consequence of the extreme spatial confinement ofthe LSP, the corresponding effective mode volume of theelectric field is much smaller than more traditional dielec-tric cavity structures, which leads to a strong enhance-ment of the spontaneous emission (SE) rate in the weak-to-intermediate coupling regime . Moreover, fine spa-tial control of the QD dynamics at the single quantum excitation level and processing of the light signal at thenanoscale is possible , resulting in a broad range ofapplications in fields such as high precision quantum in-formation processing and quantum computation, efficientsolar cells , and high precision chemical or biological de-tection .Although long distance and large scale transmissionof information using metallic structures is typically notpractical because of the strong inherent Ohmic losses ofmetals at optical frequencies, improvements can be madeby using hybrid semiconductor-metallic structures ,in which the transmission is carried out by semiconductoroptical technology, while a nanoscale metallic resonatorcan be used as an effective modulator and/or transis-tor . Thus it is of significant fundamental and appliedinterest to study the interaction between dipole emitterssuch as QDs and individual MNPs. Recently, there havebeen a few works studying how quantum emitters cou-ple to LSPs, e.g., to describe enhanced SE (Purcell ef-fect) , entanglement dynamics , and the fluorescencespectrum . In the classical or semiclassical regime, withthe exception of a particular type of simple geometry suchas a spherical cavity , for which classical analytical re-sults are available, most of the nanoplasmonic studiesare carried out by numerical analysis which is numeri-cally cumbersome and not physically intuitive . Inthe quantum optics regime of cavity-QED (cQED), ithas been common to exploit a standard cavitylike masterequation with phenomenological decay rates that are im-plicitly Lorentzian in their decay dynamics ; such anapproach is useful and easy to understand, but it is ul-timately limited since the general non-Lorentzian natureof the LSP is neglected, and it it not clear how to ob-tain the various coupling parameters, e.g., as a functionof QD distance from the MNP resonator.As a result of the fluctuation-dissipation theorem, ingeneral a continuous mode theory instead of the simplesingle mode theory needs to be employed for a quantumoptics description of an inhomogeneous lossy structure,and a quantum noise term can be included phenomeno-logically or by including a continuous reservoir atthe level of a microscopic theory ; both of these ap-proaches result in a powerful framework with the contin-uous response of the medium embedded in the mediumGreen function, which is obtained from an electric dipolesource in Maxwell’s equations. On the other hand, itis highly desirable to be able to describe the physics ofLSPs in terms of one or a few discrete modes, whichhas been the standard approach in dielectric cQED sys-tems. Recently, it is shown that the LSP can be effec-tively described as the quasinormal modes (QNMs) ofthe MNP , which are defined as the eigenfunctions ofwave equation in the frequency domain with open bound-ary conditions . A generalized mode expansion tech-nique of the classical photon Green function based onQNMs has been shown to work extremely well for vari-ous shaped MNPs, and the SE enhancement of an electricdipole located both inside and outside of MNP shows ex-cellent agreement with full numerical calculations overa broad range of frequencies around the the LSP reso-nance . The combination of an insightful QNM ap-proach and a rigorous Green function approach to quan-tum optics is thus highly desired, as MNPs facilitate acoupling regime, in general with a non-Lorentzian spec-tral density, i.e., beyond a dissipative Jaynes-Cummings(JC) model. In certain limits, it can also be used toaid a JC model and justify when such a simpler modelcan work, with clearly identified coupling rates that canbe obtained from QNM theory. Indeed, the MNP yieldsa rich mode coupling regime as a function of positionand polarization, and allows one to explore a complex in-terplay of radiative and nonradiative dynamics that areunique to the metal environment.In this paper we present a quantum optics frameworkto model the quantum dynamics between two QDs cou-pled to the LSP of a MNP system. The extension tomodel more than two QDs is straightforward and also de-scribed. While there have been several papers studyingthe dynamics of two QDs coupled by a MNP , theseapproaches, similar to the methods mentioned above,start from assuming the system could be described bythe standard cQED master equation by adding in phe-nomenological decay parameters by hand; by doing so,they neglect the possible non-Lorentzian features of the LDOS which is important in the case for QDs that aresufficiently close to the LSP resonator , and they donot incorporate the full electromagnetic response of theMNP environment, including both radiative and nonra-diative coupling effects. Instead of assuming a standardLorentzian decay rate of the LSP, we start from a micro-scopic model and derive a master equation that takes intoconsideration the electromagnetic response of the MNPin detail by exploiting a QNM expansion technique forthe photon Green function . As an example applica-tion of this theory, we consider two QDs in the vicinity ofa gold nanorod, as shown in Fig. 1(a). While other MNPshapes can also be used in our theory, including metaldimers , the single nanorod is partly motivated by thefollowing reasons: ( i ) the LSP resonance is around 1.4 eV,which is close to the wavelengths used in optical commu-nication and for many QD emitters, ( ii ) the nanorod is anontrivial geometry for which analytic methods are notreadily available, and ( iii ) , it is dominated by a singlecavity mode, polarized along the axis of the rod. Whilethe technique we exemplify below is a single mode the-ory, it can easily be generalized to include multiple LSPmodes if there are several QNMs in the frequency regimeof interest, and it properly includes the QNM dissipation.Recently, Yang et al. have carried out a somewhat simi-lar effort, to study the simple linear optical properties ofa single dipole coupled to a metal resonator .The layout of our paper is as follows. In Sec. II, wepresent our main theoretical technique and derive a quan-tum master equation based on a rigorous quantum opticsapproach for the medium in terms of the photon Greenfunction, which is obtained from the QNM of the LSP.In Sec. III, we compare the QNM technique to the JCmodel, and present the improvements over the standardJC model with the help of the QNM technique. We alsodiscuss how our approach could be used in conjunctionwith the driven JC model in certain regimes, providinga rigorous definition for the various coupling parameters.In Sec. IV A, our first example studies the simple SEdynamics from a single QD located around the metalnanorod, and shows that the non-Markovian dynamicsis important on a time scale of around the lifetime ofthe LSP. In Sec. IV B, we study the free-field dynamicsof two QDs in a homogeneous background, coupled bythe nanorod LSP and show that two qubit (QD) entan-glement can be established within a few picoseconds forseparate states with a peak value larger than 0.7, de-spite the strong Ohmic losses; we also study the affectof QD pure dephasing on the peak value of the entan-glement evolution, and investigate the concurrence (asa measure of entanglement) for different QD distancesfrom both sides of the nanorod. In Sec. IV C, we studythe incoherent spectrum for the excited two-QD system;in particular, we show explicitly how the real and imag-inary part of the Green function contributes to the cou-pling between spatially separated QDs, and find a richpolarization-dependent behavior of the spectra, includ-ing a double-resonance feature which is mediated by the FIG. 1: (Color online) (a) Schematic of the QD nanorod sys-tem with background refractive index, n B = 1 .
5; QD po-sitions (the dark brown/light grey ellipsoids) are indicatednear both ends of the gold nanorod; the nanorod has the di-mensions L = 100 nm and r c = 15 nm. (b) Norm of theQNM, | ˜ f c ( x, y, z = 0) | , with complex eigenfrequency ˜ ω c / π =( ω c + i γ c ) / π = 324 . − i16 .
58 THz. Yellow (bright) colorindicates the highest intensity regions. (c) Enhancementof the x -projected LDOS/SE (see text), F x [Eq. (11)], foran x -polarized QD, at r a = (60 , ,
0) nm—as is shown byQD a /white arrow in (a)/(b); the orange (solid)/blue (dashed)lines are given by Eq. (8) and full-dipole numerical calcula-tions, respectively. strong photon exchange effects. We present our conclu-sions in Sec. V. II. THEORY: QUANTIZATION SCHEME FOR AGENERAL MEDIUM, MASTER EQUATION,AND GREEN FUNCTION EXPANSION INTERMS OF QUASINORMAL MODES
For our MNP we consider a 3D gold nanorod as shownin Fig. 1(a) with two QDs (dark brown/light grey el-lipsoids), QD a and b , located around both ends ofthe nanorod. We use parameters for a metal nanorod with length L = 100 nm and radius r c = 15 nm, withthe Drude model for the dielectric constant, ε ( ω ) =1 − ω / ( ω + i ωγ ), where ω p = 1 . × rad/s (bulkplasmon frequency) and γ = 1 . × rad/s (collisionrate), similar to the parameters for gold. The backgroundrefractive index is n B = 1 .
5. We will also allow for thepossibility of an incoherent pump field on the QDs.For the medium quantization scheme, we start from theformalism developed by Scheel/Dung et al. , whichhas been widely used to describe the quantum electrody-namics of a quantum emitter around a spherical metallicnanoresonator . This approach, however, is com-pletely general and can be employed for any lossy inho-mogeneous structure as long as the corresponding Greenfunction G ( r , r ′ ; ω ) can be calculated. The photonicGreen function is defined through ∇ × ∇ × G ( r , r ′ ; ω ) − k ε ( r , ω ) G ( r , r ′ ; ω ) = k δ ( r − r ′ ) , for the positiondependent complex permittivity ε ( r , ω ) = ε R ( r , ω ) +i ε I ( r , ω ), where ε ( r , ω ) = ε ( ω ) inside the nanorod and ε ( r , ω ) = ε B = n B elsewhere; here k = ω/c and isthe unit dyadic. The imaginary part of the Green func-tion with the same position arguments, G ( r , r ; ω ) givesthe projected LDOS ∝ Im[ G ii ( r , r ; ω )] /ω , while theGreen function with different position arguments givesthe propagator of the electric field. For a homogeneousdielectric, the imaginary part of the homogeneous Greenfunction is given by Im[ G B ( r , r ; ω )] = ω n B πc .By treating the QDs as two-level systems, and usingthe dipole and rotating-wave approximations, the totalQD-MNP system is described by the Hamiltonian, H = ~ Z d r Z ∞ dω ω f † ( r , ω ) · f ( r , ω ) + X n = a,b ~ ω n σ + n σ − n − (cid:20) X n = a,b σ + n Z ∞ dω d n · E ( r n , ω ) + H . c . (cid:21) , (1)where σ + n /σ − n (with σ n = σ + n σ − n ) are the Pauli matri-ces of the two QDs excitons (electron-hole pairs), and ω n and d n are the resonance frequency and dipole mo-ment of the n -th QD, respectively; f / f † are the bosonfield operators , where the electric field operator is givenby E ( r , ω ) = i R d r ′ G ( r , r ′ ; ω ) · q ~ ε I ( r ′ ,ω ) ε π f ( r ′ , ω ), with ε I ( r , ω ) the imaginary part of ε ( r , ω ).In a rotating frame at the frequency of the QD a , ω a ,the total Hamiltonian becomes H = H S + H R + H I , wherethe system, reservoir, and the interaction terms are re-spectively defined by H S = − ~ ∆ ab σ + b σ − b , (2a) H I = − X n = a,b (cid:18) σ + n e iω a t Z ∞ dω d n · E ( r n , ω ) + H . c . (cid:19) , (2b) H R = ~ Z dr Z ∞ dω ω f † ( r , ω ) · f ( r , ω ) , (2c)where ∆ nm = ω n − ω m . Transforming into the inter-action picture, and using the second-order Born-Markovapproximation, the master equation for the reduced op-erator for QD pair is obtained from ∂∂t ˜ ρ ( t ) = − ~ Z t dτ Tr R n [ ˜ H I ( t ) , [ ˜ H I ( t − τ ) , ˜ ρ ( t ) ρ R ]] o , (3)where ρ R = ρ R (0) is the state of the reservoir;here we have assumed a second-order Born approx-imation, which is valid in the weak-to-intermediatecoupling regime. We assume the temperature ofthe reservoir is 0 K, which is a good approxima-tion at optical frequencies. The field operators satisfythe following relations: Tr R [ f † i ( r , ω ) , f j ( r ′ , ω ′ ) ρ R ] = 0,Tr R [ f i ( r , ω ) , f † j ( r ′ , ω ′ ) ρ R ] = δ ij δ ( r − r ′ ) δ ( ω − ω ′ ). Aftercalculating the integrand on the right hand side of Eq. (3)explicitly, we transform back to the Schr¨odinger picture,and the generalized master equation for the reduced sys-tem is obtained as ∂ρ∂t = 1 i ~ [ H S , ρ ] + Z t dτ (cid:16) X n,m [ σ − n ( − τ ) ρσ + m − σ + n σ − m ( − τ ) ρ ] × J nm ph ( τ ) + H . c . (cid:17) + X n (cid:16) γ ′ n L [ σ n ] + P n L [ σ + n ] (cid:17) . (4)Here J nm ph ( τ ) = R ∞ dωJ nm ph ( ω ) e iτ ( ω a − ω ) , with the photonreservoir function defined through J nm ph ( ω ) = d n · Im[ G nm ( ω )] · d m π ~ ǫ , (5)where, for ease of notation, we have introduced G nm ( ω ) ≡ G ( r n , r m ; ω ) with r n/m the positions of n -th /m -th QD; in addition, d n = d n n n with n n the unitvector of n -th dipole moment, and we have included apure dephasing term, L [ σ n ], on the right hand side ofEq. (4), with a dephasing rate γ ′ n , where L [ O ] = ( OρO † − O † Oρ ) + H . c . is the standard Lindblad superoperator; fi-nally, the last term L [ σ + n ] allows for the possibility of anincoherent pump term on each QD n with pump rate P n .In the following calculations, we assume γ ′ = γ ′ a = γ ′ b and | d n | = 30 D ≈ σ ± n ( − τ ) = e − i H S τ/ ~ σ ± n e i H S τ/ ~ = σ ± n e ∓ i∆ na τ , and for resonant QDs, ∆ ab = 0 , so we have σ ± n ( − τ ) = σ ± n . In order to derive Eq. (4), we haveused the identity R d r ′ G ( r , r ′ ; ω ) · ε I ( r ′ , ω ) G ∗ ( r ′ , r ′′ ; ω ) =Im[ G ( r , r ′′ ; ω )] . In a single QNM picture, note that J ph ( τ ) can be divergent; however, in our master equa-tion, we calculate R t J ph ( τ ) dτ , which is convergent. Ina practical QNM calculation, we compute the frequencyintegral over a finite bandwidth that covers the QNMresonance, with limits at approximately ± second Markov approximation (i.e., t → ∞ )for the reservoir bath sampling, then we obtain the fol-lowing Markovian master equation ∂ρ∂t =i[∆ ab σ b , ρ ] + X n (cid:18) γ ′ L [ σ n ] + γ n L [ σ − n ] + P n × L [ σ + n ] − i∆ ω n [ σ n , ρ ] (cid:19) + L coup [ ρ ] , (6)where we have introduced the QD coupling term L coup [ ρ ] = i P n = m (cid:2) ( σ + n σ − m ρ − σ − m ρσ + n ) g nm − ( ρσ + n σ − m − σ − m ρσ + n ) g ∗ mn (cid:3) , a LSP-induced SE rate γ n = 2 d n · Im[ G nn ( ω n )] · d n ~ ε , a photonic Lamb shift∆ ω n = − d n · Re[ G nn ( ω n )] · d n ~ ε , and a LSP couplingstrength between the QDs, g nm = d n · G nm ( ω m ) · d m ~ ε .When the QDs are resonant with each other,the coupling term can be simplified to L coup [ ρ ] = − i δ mn [ σ + m σ − n , ρ ] + γ nm (2 σ − n ρσ + m − σ + m σ − n ρ − ρσ + m σ − n ),where δ mn = − d m · Re[ G mn ( ω a )] · d n ~ ε describesthe coherent coupling between the QDs, and γ nm = 2 d n · Im[ G nm ( ω a )] · d m ~ ε describes the incoherentcoupling. The master equation (Eq. (6)) then simplifiesto ∂ρ∂t = X n = m γ nm σ − n ρσ + m − σ + m σ − n ρ − ρσ + m σ − n ) − i ~ [ H eff , ρ ]+ X n (cid:18) γ ′ L [ σ n ] + γ n L [ σ − n ] + P n L [ σ + n ] (cid:19) , (7)where the effective Hamiltonian term is defined as H eff = ~ ∆ ba σ b + ~ P n ∆ ω n σ n + ~ P n = m δ mn σ + m σ − n .From the above master equations, it is clear that thedynamics of the coupled QDs will show a strong po-sitional dependence through the Green function terms,and this is fully captured in the theory. Unfortunately,the calculation of the Green functions (apart from verysimple geometries) is generally a very difficult and atime consuming process, even with computations carriedout on clustered computers. In some previous studies,the coupling to the LSP was treated phenomenologicallywithout taking the full detail of the geometry and elec-tromagnetic response into consideration . However,for the MNP, recently it has been shown that the Greenfunction can be accurately obtained in terms of an expan-sion of the QNM, and for the gold nanorod (and indeedmany MNP geometries), a single QNM expansion repre-sents an accurate description of the Green function overbroadband frequencies and spatial positions. For any twospatial points near the MNP, but outside the regime ofOhmic heating, the dyadic Green function is accuratelydescribed through G c ( r , r ; ω ) = ω ω c (˜ ω c − ω )˜ f c ( r )˜ f c ( r ) , (8)where ˜ f c ( r ) and ˜ ω c are the QNM of interest and thecorrespondent complex eigenfrequency, respectively. TheQNMs are normalized through hh ˜ f c | ˜ f c ii = lim V →∞ Z V (cid:18) ω ∂ ( ε ( r , ω ) ω ) ∂ω (cid:19) ω =˜ ω c ˜ f c ( r ) · ˜ f c ( r ) d r + ic ω c Z ∂V p ε ( r )˜ f c ( r ) · ˜ f c ( r ) d r , (9)where in practise we use a computational volume ofabout 1.5 micron cubed. Alternative QNM normaliza-tion schemes are presented in Refs. [32,48], which havebeen shown to be equivalent to the one above. Forour MNP resonator, we have assumed, and verified, thatthere is only one mode in the regime of interest (nearthe LSP), and for the gold nanorod, the resonance ofthe LSP is calculated to be ˜ ω c / (2 π ) = ω c / π + i γ c / π =324 . − i16 .
584 THz (1.344 - i0.0684 eV) with qualityfactor Q = ω c / γ c ≈ .
8; in order to obtain the QNMnumerically, a 6-fs length (Gaussian shape in time do-main) spatial plane wave near 325 THz with polarizationalong the axis of the nanorod is injected, and a run-timeFourier transform with a time window 60 fs is employed; anonuniform conformal mesh scheme is used, with a meshsize of 1 nm cubed is employed around the nanorod. Thespatial dependence of the mode profile, | ˜ f ( x, y, z = 0) | , isshown around the nanorod in Fig. 1(b).In the calculation of the propagator, we use the regu-larized mode ˜ F c ( r , ω ), since it allows one to model spatialregimes from the near to far field regimes , G → G Fc ( r , r ; ω ) = ω ω c (˜ ω c − ω ) ˜ F c ( r , ω ) ˜ F c ( r , ω ) , (10)with the regularized field given by ˜ F c ( r , ω ) ≡ R V G B ( r , r ′ ; ω )∆ ε ( r ′ , ω )˜ f c ( r ′ ) d r ′ , where the volume of theintegral is now confined to the nanorod volume and∆ ε ( r ′ , ω ) = ε ( r ′ , ω ) − ε B . The regularized mode hasa simple physical interpretation: it is the solution to ascattering problem when the nanorod is excited by theQNM, which ensures the correct output characteristicsin the far field. Usually Eq. (8), which uses the divergentQNM field, gives an excellent approximation to the fullGreen function as long as the distance between the twopositions is no more than a few hundred nm away fromthe surface of the nanorod; but as the distance becomessufficiently large, Eq. (10) should be employed to calcu-late both the propagator and enhancement of LDOS (seeRef. [37] for more details).The enhancement of the projected LDOS, in direction n a , is defined as F n a ( ω ) = n a · Im[ G ( r a , r a ; ω )] · n a n a · Im[ G B ( r a , r a ; ω )] · n a , (11)and in terms of the QNM contribution, one simply sim-ply replaces G by G c [i.e., Eq. (8)]. Figure 1(c) shows thecomparison between the enhancement of the x -projected (axis of nanorod) LDOS, F x , at 10 nm [ r a = (60 , ,
0) nm]away from the nanorod [as in shown in Fig. 1(a) by QD a ],calculated via Eq. (8) (orange solid line) and with a fullnumerical dipole calculation using finite-difference timedomain method (FDTD) (blue dashed line). Clearlythe mode expansion technique gives an excellent agree-ment with the full-dipole FDTD calculation, and thusincludes the LSP reservoir function accurately for use inthe presented quantum master equation. The total SErate induced by the QNM, including radiative and non-radiative coupling, is given by γ qnm a ( r a ) = 2 d a · Im[ G c ( r a , r a ; ω a )] · d a ~ ε , (12)where G c is obtained from Eq. (8).From the analysis above, it is clear that two spatiallyseparated QDs could be coupled to each other by the cou-pling term L coup [ ρ ] as shown in Eq. (6). In the following,we will also give the emitted spectrum that can be mea-sured at the detector position r D . For a system contain-ing N QDs, the spectrum at the position of the detec-tor, r D , is given by S ( r D , ω ) = h ( E + S ( r D , ω )) † E + S ( r D , ω ) i ,with E + S ( r D , ω ) = ε P n G ( r D , r d n ; ω ) · d n σ − n (in the ro-tating wave approximation). For continuous wave exci-tation (e.g., from the incoherent pump field), in the pres-ence of just one QD (e.g., QD n ), the spectrum is givenby S p ( r D , ω ) = 1 ε | d n · G ( r n , r D ; ω ) | × lim t →∞ Z ∞ dτ e − i( ω − ω a ) τ h σ + n ( t + τ ) σ − n ( t ) i . (13)However for more than one QD, the the power spectrumis derived to be S p ( r D , ω ) = X n,q S n ( ω ) R qn ( r D , ω )++ X n In the standard JC model both the cavity field andthe quantum emitters (e.g., two level atoms) are treatedas system operators, which makes the model suitable forstudying the physics of strong coupling between the cav-ity mode and the quantum emitters. To include dissipa-tion into the cavity, the JC model assumes an electromag-netic environment with a Lorentzian spectral density, andthis works well for many dielectric cavities. However, thereservoir function from metals resonators can be highlynon-Lorentzian; moreover, it is well known that the plas-monic resonance/spectral density of metallic nanores-onators can change as a function of position around theresonator, which can be probed experimentally by mea-suring the near field electromagnetic response at differ-ent positions . As is shown clearly through Eqs. (5), (8)and (10), the spectral function of the LSP, which is givenby the imaginary part of the Green function at the samespatial point, also depends on the phase of the QNM.Our general model can actually be used to assess whenthe JC may work, with a rigorous definition of the cou-pling parameters, and it can go beyond the Lorentzianlineshape model as well. The JC model, if in a regime ofvalidity, can explore effects beyond the 2nd-order Bornapproximation, e.g., in the strong coupling regime.The Green function in Eq. (8) can be used to obtainthe photon reservoir function [Eq. (5)], e.g., for somedipole position r a (assumed to be near the resonator), G c ( r a , r a ; ω ) = ω ω c (˜ ω c − ω ) ˜ f c ( r a )˜ f c ( r a ). In a rotating waveapproximation, the imaginary part of this function canbe written asIm G c ( r a , r a ; ω ) = ω γ c (˜ f Rc ( r a )) ( ω c − ω ) + γ (cid:20) N ( r a )( ω c − ω ) ω c (cid:21) , (16)where ˜ f Rc and ˜ f Ic are the real and imaginary parts of theQNM function at the dipole position, and we have intro-duced a non-Lorentzian coupling factor defined through N ( r a ) = (˜ f Ic ( r a )) ω c γ c (˜ f Rc ( r a )) . (17)To better quantify the QNM lineshape, consider adipole position r a = (60 , , 0) nm, 10 nm away from theright side of the metal nanorod, as is shown in Fig. 1(b)by the white arrow; the non-Lorentzian coupling factoris N ( r a ) ≈ . 2, which mainly leads to a small frequencyshift of the resonance frequency (which is easily capturedin a Lorentzian function by just moving the resonancefrequency). Thus for this example, the x -projected spec-tral density J xx ph obtained from Eq. (5) (orange solid) is ω (eV) F x (d) ω (eV) F x (c) J xx ph ( µ e V ) (b) J xx ph ( µ e V ) (a) FIG. 2: (Color online) Spectral function and the enhancementof the LDOS/SE for a single QD position. (a) J xx ph at position r a , 10 nm away from the right side of the gold nanorod asshown is Fig. 1(b) by the arrow: with the orange (solid) linegiven by the QNM calculation and the blue (dashed) line is thebest Lorentzian fit. (b) QNM calculation of J xxph at position r b ′ = (0 , , 25) nm (orange solid); the blue (dashed) line showsthe same Lorentzian fit as in (a), but rescaled in amplitude.(c) Enhancement of the LDOS/SE with the QNM calculation(orange solid) and the Lorentzian function (blue dashed) in(a) for an x -polarized QD at r a . (d) QNM calculation of F x for and x -polarized dipole at position r b ′ = (0 , , 25) nm(orange solid); and the blue dashed line uses the Lorentzianspectral function from (a). well described by a Lorentzian line shape (blue dashed)as shown in Fig. 2(a). In general, however, a position de-pendent non-Lorentzian spectral density will be obtainedaround the nanorod; this effect is shown in Fig. 2(b) for J xx ph at the example position r b ′ = (0 , , 25) nm, by thesolid (orange) line, and a slight blueshift of the resonancepeak is also seen which is consistent with the observationin Ref. [51]; the dashed (blue) line is the same Lorentzianfit used in the previous case, which clearly shows that theline shape changes as a function of dipole position. Thecorresponding enhancement of the LDOS/SE is shownin Fig. 2(c). At position r b ′ = (0 , , 25) nm, the non-Lorentzian shape of the spectral density influences theenhanced LDOS/SE as is shown in Fig. 2(d); the non-Lorentzian coupling factor is now N ( r a ) ≈ − . 4, whichhas a more dramatic effect on the spectral line shape. Westress that all the information of the resonance shift andnon-Lorentzian spectral function is included in the spa-tial dependence of the phase factor of the QNM; and thisinformation naturally comes into the calculations belowthrough the QNM normalization condition the analyti-cal Green function. Although the Lorentzian fit of thespectral function is valid at certain spatial locations, thenon-Lorentzian spectrum becomes important when N islarge enough, and one then requires the imaginary partof the QNM as well as the real part [32].As discussed above, and shown in Fig. 2(a), at somepositions the spectral density could be well describedby a Lorentzian line shape (for certain MNPs), so forQDs at these positions the quantum optical interac-tions could be approximately described by a dissipa-tive JC model with the following QD-cavity coherentinteraction terms (in a rotating wave approximation): d n · ˜ f c ( r a ) aσ + + d n · ˜ f ∗ c ( r a ) a † σ − ; however, we see that theseparameters (and the parameters needed to describe QD-QD interactions) actually require the QNM techniquein order to have a rigorous definition of these couplingparameters and the normalized field. The single QD-cavity coupling rate will be given by the usual rate g ,where g ∝ d n | ˜ f c ( r d ) | , while dissipation from the cavitymode is then usually added through a Lindblad operatorthat describes only Lorentzian decay. For a dissipativeLorentzian model to work, we find that the QD positionsmust be located around high symmetry points within thevicinity of the field antinode points, but far enough awayfrom the metal surface. Even when the JC model ap-proximately works, then the decay rates still have to beobtained as a function of position in general. This is pre-cisely what the QNM can provide, if the QD position isin a valid Lorentzian decay regime.We also caution that the JC model still neglects someessential dissipative coupling processes from the metalenvironment. For example, the standard JC model doesnot provide an effective description of the nonradiativeand radiative decay processes; such a description of thenonradiative/radiative decay will be necessary in order tocompute important properties such as the quantum yield(or beta factor), e.g., of a single photon source. Belowwe demonstrate how one can use the QNM technique toachieve the separation of the total decay rate into radia-tive and nonradiative decay channels. Moreover, we willalso show how one can add in Ohmic losses in a clear andsimple way, which is needed for dipole positions very nearthe resonator (e.g., a few nm from the surface) . Impor-tantly, in our approach, all of these physical rates can becomputed analytically using the QNM theory, as a func-tion of space and frequency. We describe and exemplifythese scattering rates below.Without the metal nanorod, the background decay rateis simply γ = d a · Im[ G B ( r a , r a ; ω a )] · d a ~ ε = d a ω a n B ~ ε πc . Whilethe nonradiative decay rate from the QNM is obtainedfrom γ nrq a ( r a ) = 2 ~ ω a ε Z V MNP Re[ j ( r ) · G ∗ c ( r , r a ; ω a ) · d a ] d r , (18)where j ( r ) = ω c ε I ( r , ω a ) G c ( r , r a ; ω a ) · d a is the inducedcurrent density in the nanorod (MNP) at position r a .Thus the radiative decay rate from the QNM is given by γ rq a ( r a ) = γ qnm a ( r a ) − γ nrq a ( r a ) . (19)In addition, the quasistatic decay rate can be obtained h (nm) γ r q a , γ n r q a , γ s t a t a ( γ ) 20 4000.51 h (nm) η r a d , η n r FIG. 3: (Color online) Decay rates (in units of the homoge-neous space radiative decay rate γ ) of an x -polarized QD a induced by quasi-static interaction γ stat a (magenta dashed),radiative contribution of the QNM γ rq a (red solid), and thenon-radiative contribution of the QNM γ nrq a (blue chain),as a function of h [ r a = (50 nm + h, , η rad (lowerred solid), and non-radiative coupling factor, η nr (upper bluesolid) of the decay rate. from γ stat a ( r a ) = 2 d a · Im[ G qs ( r a , r a ; ω a )] · d a ~ ε , (20)with G qs ( r a , r a ; ω a ) = ∓ G B ( r ′ a , − r ′ a ; ω a ) ε ( ω a ) − ε B ε ( ω a )+ ε B ) 37,53 ( ∓ is for s / p -polarized dipoles, respectively). Conse-quently, the total nonradiative decay is given by γ nr = γ nrq a + γ stat a . As is shown above, all of the decay rates arehighly position dependent, but once the QNM is calcu-lated, the decay rates at different positions can be com-puted immediately.Figure 3 shows the decay rates as a function of dis-tance, h , away from the surface of the metal nanorodalong the x -axis; where we see that, at extremely smalldistances h , the system is in quasi-static regime wherethe Ohmic heating effect due to γ stat a (magenta dashed)is strong; as h becomes larger, all of the decay rates de-crease, but the quasi-static decay rate decreases muchfaster than the others; the inset of Fig. 3 shows the radia-tive coupling factor, η rad = γ rq a γ qnm a + γ stat a (lower red solid),and nonradiative coupling factor, η nr = γ stat a + γ nrq a γ qnm a + γ stat a (up-per blue solid), as function of h in the near field regime.We see that the radiative output coupling efficiency isbelow 50%, though this can be increased to about 60%or greater using a dimer configuration . While it is notclear how to include such processes in a JC model, whichwould be further complicated by having different param-eters at different QD positions, they could certainly helpimprove and guide such simpler models. More details onsuch an approach will be reported in a future publicationwhen we will also explore effects beyond a second-orderBorn approximation.With regards to computing the spectrum in a JCmodel, the spatially integrated far-field spectrum fromthe cavity operator is usually given by (assuming arotating frame as the cavity frequency) S cav ( ω ) ∝ κ lim t →∞ R ∞ dτ e − i( ω − ω c ) τ h a † ( t + τ ) a ( t ) i ; however thisassumes that the output coupling rate via the cav-ity (in this case the LSP) is purely radiative. Fora metal resonator, one must include quenching ef-fects by solving the input/output scattering problem,which is exactly what the Green function solution hasdone. In this way one can compliment the JC modelby computing the spatially dependent output spec-trum from the QD system operator dynamics and themedium electric field operators, so that S cav ( r , ω ) ∝ | d · G c ( r , r d ; ω ) | lim t →∞ R ∞ dτ e − i( ω − ω c ) τ h σ + ( t + τ ) σ − ( t ) i .Furthermore, one could obtain the spatially averagedoutput spectrum (e.g., in the far-field) from S radcav ( ω ) ∝ η c κ lim t →∞ R ∞ dτ e − i( ω − ω c ) τ h a † ( t + τ ) a ( t ) i , where κ = γ c and the radiative output coupling factor associated withthe QNM is obtained from η c = γ rq a γ qnm a .To summarize this section, we have discussed howour model can go well beyond the standard JC modelwhile facilitating the simpler JC models in certain spatialregimes. To the extent that the approximate JC couldbe valid, one still has to obtain the coupling parametersfrom a model such as ours, and then carefully includequenching effects into any calculation of emitted fieldsfar away from the system resonator. Thus our modelcan be used to guide and help the simpler JC models incertain regimes as well. IV. RESULTS AND APPLICATIONSA. Localized plasmon induced SE from a singleexcited QD near the nanorod Metal nanoparticles enhance the SE rate of excited sin-gle QDs due to the coupling with the LSP (QNM). TheLDOS at positions around the nanorod changes rapidly inspace compared to the homogeneous dielectric structure,and thus the SE rate of a QD around the nanorod can besignificantly changed , as is shown in Fig. 1(c).In this section, we present an analysis of the SE dynam-ics of a single QD on resonance with the LSP ( ω a = ω c );without loss of generality, we take the case of a QD po-larized along the x -axis at position r a , 10 nm away fromthe nanorod as shown schematically in Fig. 1(a) (QD a ).For the single-QD nanorod system, without an incoher-ent pump field (i.e., P a = 0), the non-Markovian master(Eq. (4)) becomes ∂ρ∂t = Z t dτ (cid:0) [ σ − ρσ + − σ + σ − ρ ] J ph ( τ ) + H . c . (cid:1) + γ ′ L [ σ ] , (21) t (ps) N a t (fs) γ a ( t )( G H z ) FIG. 4: (Color online) Population dynamics of an excitedsingle QD, N a , for an x -polarized dipole with | d | = 30 Dlocated at, r a = (60 , , 0) nm, 10 nm away from the nanorodas shown in Fig. 1(a) by QD a . Exponential decay with γ a = d a · Im[ G ( r a , r a ; ω a = ω c )] · d a ~ ε shown by the blue dashed line, andthe green solid line is the full non-Markovian dynamics givenby Eq. (21). The inset shows the effective exponential decayrate, γ a ( t ), calculated with the full decay dynamics (greensolid), and γ a (blue dashed); at the crossover region, around40 fs ( ≈ π/γ c ), γ a ( t ) agrees with γ a within 2% as shown bythe light grey area). where we have explicitly used the result σ ± ( − τ ) = σ ± since H S = 0 and the kernel function is given by J ph ( τ ) = J aa ph ( τ ). Note the population decay is not affected by puredephasing here so we can neglect pure dephasing for thissingle QD radiative decay study.We assume here that the QD is initially excited.The QD population decay, N a = ρ ee = h e | ρ | e i , us-ing Eq. (4), is shown in Fig. 4 by the light (green)solid line; the dashed (magenta) line shows the resultof a Markovian exponential decay with the rate, γ a = d a · Im[ G ( r a , r a ; ω a )] · d a ~ ε given by Fermi’s golden rule withthe on-resonant projected LDOS. The inset to Fig. 4shows that the SE dynamics is recovered by Fermi’sgolden rule after a characteristic timescale of about 40 fs(shown in the light gray region); this time scale agreesvery well with the corresponding lifetime of the LSP, τ c ≈ π/γ c . B. Localized plasmon induced coupling betweentwo spatially separated QDs in homogeneousbackground For two spatially separated QDs located around thenanorod (as is shown in Fig. 1(a) by QD a and QD b ),these can be effectively coupled to each other by exchang-ing photons via the LSP; in the absence of a pump field(i.e., P n = 0), the non-Markovian master equation be-comes ∂ρ∂t = Z t (cid:18) X n,m [ σ − n e i∆ na τ ρσ + m − σ + n σ − m e i∆ ma τ ρ ] J nl ph ( τ )+ H . c . (cid:19) dτ + X n γ ′ n L [ σ n ] + i[∆ ab σ b , ρ ] . (22)Unless stated otherwise, we will assume the two QDsare resonant with each other (∆ ab = 0), but may beoff resonant with the LSP, where ω a = ω b = ω c + ∆;however, later we also study the case with differentQD resonance frequencies [e.g., in Fig. 11(d)]. Herethe intercoupling between the QDs depends on the pro-jected cross density of optical states (CDOS), ̺ ab ≡ ̺ ( r a , r b ; ω ) ≡ n a · Im[ G ( r a , r b ; ω )] · n b via J ab ph , which givesone part of the characteristic coupling strength betweenthe QDs mediated by the electromagnetic environmentof the nanorod. Since the Green function in use is the re-tarded Green function, the real and imaginary parts arerelated to each other through the Kramers-Kronig rela-tion. The real part of the Green function between the twoQDs yields the coherent coupling, δ ab ( σ + a σ − b + σ + b σ − a ),while the imaginary part gives the incoherent coupling, P n = m γ ab (2 σ − n ρσ + m − σ + m σ − n ρ − ρσ + m σ − n ); the relevant cou-pling strengths are shown in Fig. 5 with the coherentcoupling ( δ ab ) and the incoherent coupling strength ( γ ab )given by the chain (orange) and dashed (blue) lines, re-spectively, and the solid (cyan) line shows γ a . It canbe seen from Fig. 5(a), for x -polarized QDs at positions r a/b = ( ± , , 0) nm, that δ ab may dominate over γ ab and γ a when ∆ = ∆ off , where ∆ off is some offset fre-quency from the real part of the LSP resonance; how-ever, when QD b is z ( y )-polarized, there is almost no cou-pling between the QDs; but due to the complex position-dependent polarization characteristics of the LSP, the off-diagonal element of both the projected CDOS, ̺ ab ( ω ),and the real part of the Green function are non-zero;consequently, QDs with different polarization can beeffectively coupled to each other for certain QD posi-tions. Figure 5(b) shows that an x -polarized QD a at r a = (60 , , 0) nm could be effectively coupled to a z -polarized QD b at r b ′ = ( − , , 23) nm. Below, we willalso look at the effect of the coherent exchange interac-tions in the presence of QD pure dephasing.We first assume that the two QDs are initially in aseparable state | eg i (with the first argument for QD a and the second one for QD b ); as a result of the coherentcoupling, non-classical correlations will be induced be-tween the QDs, and the quantum correlations approacha maximum value at some characteristic interaction time,which eventually decays to zero due to the decoherencecaused by the strong dissipation and the pure dephas-ing of the system. As a measure of the nonlocal quan-tum correlations between the separated QDs, we use theconcurrence C , which is obtained from the eigenvaluesof the flipped density matrix , and its maximum valuein the evolution is denoted as C max . The exciton pop-ulation of QD n is defined as N n = h e n | tr m ρ | e n i with −300 0 300−50050100150 ω − ω c (meV) γ a , γ a b , δ a b ( µ e V ) (a) ∆ off −300 0 300−150−100−50050100150 ω − ω c (meV) (b) ∆ off FIG. 5: (Color online) The various scattering rates as thefunction of detuning from the resonance of the LSP ( ω c ) fortwo resonant QDs with QD a at r a = (60 , , 0) nm (i.e., 10 nmaway from the nanorod surface). (a) γ aa (cyan solid), γ ab (blue dashed), and δ ab (orange dash-dot) with QD b at r b =( − , , 0) nm; the left and right short (red) vertical linesare for later reference when we choose QD detunings of ∆ = ω a − ω c = − 224 meV ≡ ∆ off , and ∆ = 0 meV, respectively.(b) same as (a) but with r b ′ = ( − , , 23) nm [see Fig. 1(a)]. N a , N b , C t (ps)( a ) C m a x γ ′ ( µ eV) t (ps)( b ) C m a x γ ′ ( µ eV) FIG. 6: (Color online) Dynamics of two resonant QDs with x -polarized excited QD a and unexcited QD b , | eg i , at r a =(60 , , 0) nm, r b = ( − , , 0) nm, respectively. (a) Timeevolution of the exciton population of QD a /QD b , N a/b (bluedashed/orange solid), and entanglement C (dark green solid)with pure dephasing rate γ ′ = 10 µ eV for ∆ = 0. The insetshows C max as a function of γ ′ and the circle shows the positionat which the dynamics is studied. (b) same as (a) except for∆ = ∆ off . n = m . Figure 6(a) shows the dynamics of C for x -polarized QDs on resonance with the LSP (∆ = 0) atpositions r a/b = ( ± , , 0) nm [shown in Fig. 1(a) bythe dark brown ellipsoids] by the green (dark) solid linefor γ ′ = 10 µ eV; the exciton populations N a/b are shownby the blue dashed line and orange (light) solid line, re-spectively; the maximum C max as a function of the puredephasing rate γ ′ is shown in the inset and, for this case, C max is always less than 0.5 in agreement with previouswork for entangled atoms in weakly coupled reservoirs .As shown in Fig. 5(a), when the QDs are on resonancewith the LSP, the incoherent coupling rates are muchlarger than the coherent coupling rate | γ ab | , | γ a | ≫ | δ ab | ,but as they are detuned away from the LSP resonance,0 (a) | + i ∆ = 0 N a , N b , C (b) | + i ∆ = ∆ off (c) |−i ∆ = 0 t (ps) N a , N b , C (d) |−i ∆ = ∆ off t (ps) FIG. 7: (Color online) Population decay dynamics of the ini-tial Bell initial states |±i for both x -polarized QDs with thesame resonance frequency, at positions r a/b = ( ± , , 0) nm,respectively; the pure dephasing rate is γ ′ = 10 µ eV. (a) Forthe initial state | + i with ∆ = 0 the dark green solid line is C ,and the blue dashed/orange solid lines show the exciton pop-ulation of QD a /QD b , respectively; (b) same as (a) but with∆ = ∆ off (see Fig. 5); (c) same as (a) but with the initialstate |−i ; (d) same as (c) but with ∆ = ∆ off . the coherent coupling strength begins to dominate overthe incoherent coupling. Figure 6(b) shows the same cal-culation as 6(a) but with ∆ = ∆ off , and now we see that C max could be much larger than the previously limit of0.5 ; indeed our calculations show that it could be evenlarger than 0.7 if the pure dephasing rate were smaller;in addition, we see that the concurrence exhibits an oscil-lating behaviour, which is similar to that of the coherentsystem indicating that the two-QD are effectively coupledthrough the LSP-induced photon exchange.For the detuning value of ∆ = ∆ off , the incoherentrates are around γ a/b ≈ µ eV, the incoherent couplingrate γ ab ≈ µ eV, while the coherent coupling strengthis δ ab ≈ . µ eV; in contrast, for the on resonance case(i.e., ∆ = 0), we have γ a/b ≈ µ eV, γ ab ≈ µ eV, and δ ab ≈ µ eV. It can be seen that the relative coherentcoupling strength with a finite detuning (∆ = ∆ off ) ismuch larger than it is at ω c (neglecting γ ′ ). However,as γ ′ increases, the effective coherent coupling strengthfor ∆ = ∆ off decreases much faster than for ∆ = 0.Thus when γ ′ = 0, C max for ∆ = ∆ off is larger thanthat for ∆ = 0, but it decreases faster as well since γ ′ increases—as shown in the insets of Fig. 6.It is also demonstrated in Fig. 6 that, due to thepresence of QD b , the decay of QD a in the long timelimit slows down. This effect can be explained throughthe effective Hamiltonian, H eff , which in the absenceof dissipation results in four eigenstates | ee i , | gg i , |±i = √ ( | eg i ± | ge i ). The initial state, | eg i , lies in the sub-space composed of |±i which gives the superradiant and subradiant emission depending on the relationship amongthe enhanced SE rate, γ a , and the incoherent coupling, γ ab . Thus the decay of the excited QD may be enhancedat the beginning ( t → 0) due to the faster decay of thecomponent of superradiant state in the initial state; whileat long times, the dynamics is dominated by the slowerdecay of the component of subradiant state in the ini-tial state, which gives a suppressed emission if there is aconsiderable amount of the subradiant component in theinitial state.Figure 7 shows the dynamics of the resonant QDs lo-cated symmetrically at r a/b = ( ± , , 0) nm with theinitial states |±i , with a pure dephasing rate γ ′ = 10 µ eV.It is shown in Fig. 7(a), that when the QDs are onresonant with the LSP, | + i is the superradiant state( γ a/b ≈ γ ab ≫ γ ′ ), and N a/b (blue dashed/orange solid)decay twice as fast than QD a alone. The dynamics withthe initial state |−i is shown is Fig. 7(c), which is nowthe subradiant state. With a detuning of ∆ off , we have γ a/b ≈ γ ab , which are much less than γ a at ∆ = 0,so there is not much difference between the superradi-ant and subradiant states as is shown in Figs. 7(b)-(d).To establish if there are any non rotating-wave effectsnot captured by our master equation approach, we havealso checked that an exact wavefunction method basedon the schr¨odinger equation (with no rotating waveapproximation, but restricted to weak excitation with nopure dephasing) gives the same solution as above with nonoticeable difference.We stress that with our QNM formulation, one doesnot need to calculate additional Green function simula-tions for different QD positions, which makes the ap-proach convenient for exploring the position-dependentbehaviour of QDs (as we have demonstrated earlier forthe position dependent decay rates). For the initialstate | eg i , numerical calculations (with ∆ ab = ∆ = 0, γ ′ = 10 µ eV) show that the maximum achievable en-tanglement, C max , is not a monotonic function of dis-tance h from the QDs to the both sides of the nanorod, r a/b = ( ± ± h, , 0) nm. It is found that, at first C max increases as h becomes larger, and reaches its max-imum around h = 10 nm; then, it decreases as h in-creases further; for example, C max ( h = 2 nm) ≈ . C max ( h = 10 nm) ≈ . 45, and C max ( h = 18 nm) ≈ . h , the system is in the quasi-static cou-pling regime where Ohmic losses due to γ stat a (magentadashed) are strong; as h becomes larger, the Ohmic lossesbecomes smaller and the effective coupling between theQDs becomes larger and thus C max increases; however, asthe spatial distance increases further, the effective cou-pling strength becomes weaker and weaker with respectto the pure dephasing rate, γ ′ , which causes C max todecrease again.1 −200 0 20000.51 S p ( a r b . un i t s ) (b) P = 0 . µ eV −200 0 20000.51 ω − ω L ( µ eV) S x / y p ( a r b . un i t s ) (c) P = 0 . µ eV −200 0 20000.51 S ( a r b . un i t s ) (a) P = 0 . µ eV −200 0 20000.51 ω − ω L ( µ eV) S p ( a r b . un i t s ) (d) P = 10 µ eV FIG. 8: (Color online) Incoherent spectra for QDs at r a/b =( ± , , 0) nm, respectively ( ω a/b = ω c ) with γ ′ = 1 µ eV.(a) S with P = 0 . µ eV: the black solid/green dashed areresults with/without the presence of QDs, respectively; (b)Incoherent spectra S p at position r D = (0 , , . / . / µ m(magenta/blue dashed/green) with P a = 0 . µ eV. (c) Polar-ization dependent spectra S xp (orange chain) and S zp (cyan)at r D = (0 , , . µ m with P a = 0 . µ eV. (d) S P at position r D = (0 , , . / . µ m (magenta/blue dashed); S xp (orangechain) at r D = (0 , , . µ m with P a = 10 µ eV. −200 0 20000.51 S p ( a r b . un i t s ) ω − ω a ( µ eV) (a) P a/b = 0 . µ eV −200 0 20000.51 ω − ω a ( µ eV) S p ( a r b . un i t s ) (b) P a/b = 10 µ eV FIG. 9: (Color online) Incoherent spectra, S P at position r D = (0 , , . / . µ m (magenta/blue dashed) for QDs at r a/b respectively ( ω a/b = ω c ) with γ ′ = 1 µ eV. (a) P a/b =0 . µ eV. (d) P a/b = 10 µ eV. C. Emitted spectrum from an incoherent pump As is analysed in Sec.II and shown explicitly inSec. IV B, two spatially separated QDs can be effec-tively coupled to each other due to the characteristicsof the CDOS (incoherent coupling) and the real part ofthe Green function (coherent coupling) of the nanorod.In the following, we will concentrate on the spectrumthat can be measured using excitation from an incoher-ent pump field.For our first incoherent pump investigation, we as-sume both the two x -polarized QDs are resonant withthe LSP of the nanorod (∆ ab = ∆ = 0), and they aresymmetrically located at 10 nm ( r a/b ) away from the both sides of the nanorod ( γ a = γ b , ∆ ω a = ∆ ω b ) asis shown in Fig. 1(a). We first assume only QD a is in-coherently pumped ( P b = 0), and we will compare thisresult with the spectrum emitted when both QDs are in-coherently excited. The emitted spectra are shown inFig. 8, with γ ′ = 1 µ eV and P b = 0. Without the pres-ence of QD b , the bare spectrum of QD a , S a , is shownin Fig. 8(a) by the black solid line, and its full-width athalf maximum (FWHM) is much larger than for a ho-mogeneous medium due to the LSP coupling; the green(dark) dashed line is S a from QD a including the pres-ence of QD b . The linear spectrum, S p , at difference posi-tions are shown in Fig. 8(b) for r D = (0 , , . / . / µ m(magenta solid, blue dashed, green solid). It is in-teresting that the spectrum in the near field regime, r D = (0 , , . µ m, shows a sharp spectral peak; butas the detector position is located further away from thenanorod, at r D = (0 , , . µ m, a broadened peak witha sharp peak located at the same spectral position as inthe near field is observed; as the spectrum propagates tothe far field regime, at r D = (0 , , µ m, then the sharppeak develops into a dip, which indicates that there isinterference between the sharp peak and the broad reso-nance which is a Fano resonance effect. From the analysisin Sec. IV B, |±i are eigenstates of the H eff , which arethe superradiant and subradiant states, respectively. Thesharp peak is the result of decay from |−i to | gg i , andthe broad peak is the decay from | + i to | gg i , while S p is the total contributions from the two including inter-ference effects; the separation between the peaks is givenby 2 ~ | δ ab | in the linear regime, and the asymmetry ofthe position with respect to ω a depends on the inducedLamb-shift ∆ ω a/b . In Figs. 8(a) and (b), the bare spec-trum is almost the same as the sharp peak, which meansthe population of |−i is much larger than | + i in the lin-ear regime where P a ≪ γ a/b , γ ′ , γ ab . In fact, under thissituation, the rate equations of ρ ++ , and ρ −− are simplygiven by dρ ++ dt = γ ′ ρ −− − ρ ++ ) + γ ( ρ ee − ρ ++ )+ P a ρ gg − ρ ++ ) + γ ab ( ρ ee − ρ ++ ) , (23a) dρ −− dt = γ ′ ρ ++ − ρ −− ) + γ ( ρ ee − ρ −− )+ P a ρ gg + ρ −− ) + γ ab ( ρ −− − ρ ee ) . (23b)Since at the steady state ρ gg ≈ ρ ee ≈ P a , then the ratio of steady state population ρ −− to ρ ++ is given by ρ −− ρ ++ ≈ γ ′ + γ + γ γ ′ + γ − γ ≫ S zp (cyan solid), and S xp (orange chain); it is seen thatthe sharp peak is mainly z -polarized, while the broadpeak is primarily x -polarized. So, at the far field, the S p displays mainly the broad peak as a result of the dipoleradiation that is observed in Fig. 8(b). In the presence ofa nonlinear pump field, with P a = 10 µ eV, the computed2 −100 −50 0 5000.20.40.60.81 S p ( a r b . un i t s ) ω − ω a ( µ eV) (a) P = 0 . µ eV −100 −50 0 5000.20.40.60.81 ω − ω a ( µ eV) S p ( a r b . un i t s ) (b) P = 10 µ eV FIG. 10: (Color online) Incoherent spectra for reso-nant QDs ( ω a/b = ω c + ∆ off ) at positions r a/b =( ± , , 0) nm respectively with γ ′ = 1 µ eV. (a) S p at po-sition r D = (0 , , . / . / µ m (thick magenta solid/thickblue dashed/thin grey solid) with P a = 0 . µ eV. (b) S p (ma-genta) and S xp (orange dashed) at position r D = (0 , , . µ mwith P a = 10 µ eV. spectrum S p , is shown in Fig. 8(d) at r D = (0 , , . µ mby the magenta solid line, and the orange chain is the x -polarized spectrum S xp .We next consider both QDs incoherently excited bythe same pump field, with P a = P b . In the linear regime(weak pump limit), we get basically the same result asthe case with one QD incoherently pumped, as shown inFig. 9(a). As long as the pump rate is small, the inco-herent spectrum is similar with one or two pump fields.However, as the pump field is increased, the power broad-ening with two QDs excited is notably larger, as depictedin Fig. 9(b). While there are many pumping scenariosthat we could study, in what follows, we will concentrateon the case that only QD a is incoherently excited.As we detune the QDs from the LSP resonance, using∆ = ∆ off , a double-peak is observed for the total spec-trum, S p , at r D = (0 , , . µ m as is shown in Fig. 10(a)by the magenta solid line; as before, the right (higherfrequency) peak is suppressed as we evolve to the farfield regime. The corresponding high pump spectrum S p is shown in Fig. 10(b) by the magenta solid line at r D = (0 , , . µ m, and the two peaks are now less ac-cessible than it in the low pump (linear) regime; for apump rate around P a = 40 µ eV, the steady state pop-ulations are around N a ≈ . , N b ≈ . 5, and the doublepeaks merge into a single resonance peak; in contrast,the populations of both QDs are negligible for the pumprate as low as P a = 0 . µ eV. It is interesting to notethat similar physics occurs in an incoherently pumpedquantum-dot–cavity system , though in that case therewas no incoherent coupling contribution ( γ = 0), sothe two peaks (vacuum Rabi splitting peaks) were at thesame height. In the present case, the doublet feature isentirely due to photon exchange effects, mimicking thewell known vacuum Rabi doublet.As is discussed above, by using the QMN technique wecan efficiently conduct a detailed study of the positionaldependence on the dynamics of the system. By way of anexample, we next study the incoherent spectra of orthog- −200 −100 0 10000.51 S p ( a r b . un i t s ) (a) −100 −50 0 5000.51 ω − ω a ( µ eV) S p ( a r b . un i t s ) (c) −100 −50 0 5000.51 (b) −50 0 5001 S z p ω − ω a ( µ eV) −100 −50 0 5000.51 ω − ω a ( µ eV) (d) FIG. 11: (Color online) Incoherent spectra observed at r D = (0 , . , µ m for x -polarized QD a at position r a =(60 , , 0) nm and z -polarized QDb at r b ′ = ( − , , 23) nm(a) S p for resonant QDs (∆ = 0) with P a = 0 . µ eV (greendashed) and P a = 10 µ eV (magenta solid); γ ′ = 1 µ eV; (b)same as (a) but with ∆ = ∆ off ; insert shows S zp with thesame color scheme. (c) S p (magenta solid) and S zp (orangedashed) for resonant QDs (∆ = ∆ off ) with P a = 10 µ eV and γ ′ = 5 µ eV. (d) S p (magenta solid) and S zp (orange dashed)for off resonant QDs (∆ ab = − µ eV, ω a = ω c + ∆ off ) with P a = 10 µ eV and γ ′ = 1 µ eV. onal QDs at the detector position r D = (0 , . , µ m,with an x -polarized QD a at r a , and using a z -polarizedQD b at r b ′ = ( − , , 23) nm; for this configuration,note that QD b (polarized in the z direction) obtains aneven larger SE enhancement induced by the QNM cou-pling, γ qnm b ≈ 500 ( γ ), but a slightly smaller nonra-diative contribution of γ nrq b ≈ 212 ( γ ) when ∆ ab = 0;the spectra, S p , with P a = 0 . / µ eV, γ ′ = 1 µ eVand ∆ = 0 are shown in Fig. 11(a) by the (dark) greendashed/magenta solid lines, respectively. While the spec-tra, S p , at ∆ = ∆ off are shown in Fig. 11(b) by thegreen dashed line ( P a = 0 . µ eV) and magenta solidline ( P a = 10 µ eV); the inset shows S zp with the samecolor scheme; we see that it is now much easier to ac-cess the double-peak structure with the polarization de-pendent spectrum. For a larger pure dephasing rate of γ ′ = 5 µ eV, S p is shown in Fig. 11(c) by the magentasolid at ∆ = ∆ off , and it shows the double-peak struc-ture is less visible as the pure dephasing rate increase;the orange dashed line displays S zp . Finally, we have alsostudied the effect of detuning between the QDs on thespectra. Figure 11(d) shows the spectra with QD-QDdetuning ∆ ab = − µ eV, P a = 10 µ eV and γ ′ = 1 µ eV:the detuning changes both the separation and positionof the peaks, and the double peak structure is seen tobe robust against detuning as long as it is not too large(with respect to δ ab ) as is shown by the magenta solidline ( S p ). But as the detuning becomes larger and largerthe double-peak inevitably begins to disappear. How-3ever, this robustness is in general much larger than fornarrowband dielectric cavity systems.In general the splitting of the incoherent spectrum S p could be effectively controlled by the coherent couplingstrength between the QDs ( δ ab ∝ Re[ G ˆ n a ˆ n b ab ( ω a )]), whichcan be achieved by changing both the location and po-larization of the QDs (or moving the nanorod); the rel-ative height of the double peak will be changed at thesame time since the position dependent behaviour of theplasmon-induced decay rate ( γ n ∝ Im[ G ˆn n ˆn n kk ]( ω n )) andthe cross decay rate ( γ ∝ Im[ G ˆ n a ˆ n b ab ( ω a )]) will alsochange. V. CONCLUSION In summary, we have presented an efficient masterequation formalism to include the effect of coupling ar-tificial atoms (QDs) to the dissipative electromagneticresponse of a gold nanorod or general shaped metal res-onator, which was aided through a QNM expansion of themedium Green function. Using the derived master equa-tion we studied the dynamics of two QDs, and showedthat due to the complicated position-dependent polariza-tion characteristic of the LSP, QDs could be effectivelycoupled together even with orthogonal polarization. Ourresults first show that the non-Markovian regime can beimportant for time scales of the order of the LSP lifetime,after which the dynamics of the SE decay is well describedby the exponential decay law with decay rate given by theimaginary part of the Green function at the frequency ofthe QD. We have also discussed how our model differs andcompares with a simpler JC approach, and the potentiallimitations of the JC model are highlighted. In certain regimes where a JC model could work, then the requiredparameters can also be obtained directly from the QNMtheory. Using our more general theory, we then presenteda selection of examples to study the quantum dynamicsof two QDs coupled to a gold nanorod, and discussed thevarious radiative and nonradiative coupling rates as afunction of QD position. For separate initial states withone of QDs excited, maximum entanglements of greaterthan 0.7 could be achieved within a few ps, which alsoshows a non-monotonic behavior as a function of distancefrom the QDs to the nanorod; we also showed that in or-der to get the QDs more effectively coupled, the QDsshould be detuned away from the resonance of the LSP.With an incoherent pump, Fano resonance features arepredicted in emitted spectrum, with a rich polarizationdependent behaviour and a double-peak structure thatsignals strong photon exchange effects between the LSPcoupled QDs. Importantly, our theory can quickly treatthe coupling dynamics between multiple QDs at variousspatial locations over a wide range of frequencies and al-lows an intuitive understand of the underlying physicsof LSP coupling, including a proper decoupling of ra-diative and nonradiative decay channels. As shown byKewes et al. , the ability to separate such processes isimportant to accurately model emerging nanoplasmonicdevices such as SPASERS. Acknowledgements This work was supported by the Natural Sciences andEngineering Research Council of Canada and Queen’sUniversity. ∗ Electronic address: [email protected] † Electronic address: [email protected] K. J. Vahala, Nature , 839 (2003). M. S. Tame, K. R. McEnery, S. K. ¨Ozdemir, J. Lee, S. A.Maier, and M. S. Kim, Nat. Phys. , 329 (2013). C. Van Vlack, P. T. Kristensen, and S. Hughes, Phys. Rev.B , 075303 (2012). A. Tr¨ugler and U. Hohenester, Phys. Rev. B , 115403(2008). F. Alpeggiani, S. D’Agostino, and L. C. Andreani, Phys.Rev. B , 035421 (2012). A. Delga, J. Feist, J. Bravo-Abad, and F. J. Garcia-Vidal,Phys. Rev. Lett. , 253601 (2014). P. Anger, P. Bharadwaj, and L. Novotny, Phys. Rev. Lett. , 113002 (2006). M. Frimmer and A. F. Koenderink, Phys. Rev. Lett. ,217405 (2013). C. Belacel, B. Habert, F. Bigourdan, F. Marquier, J.-P.Hugonin, S. M. Vasconcellos, X. Lafosse, L. Coolen, C.Schwob, C. Javaus, B. Dubertret, J.-J. Greffet, P. Senel-lart, and A. Maitre, Nano Lett. , 1516 (2013). A. V. Krasavin, and A. V. Zayats, Phys. Rev. Lett. , 053901 (2012). H. W. Lee, G. Papadakis, S. P. Burgos, K. Chander, A.Kriesch, R. Pala, U. Peschel, and H. A. Atwater, NanoLett. , 6463 (2014). P. Bharadwaj, B. Deutsch, and L. Novotny, Adv. Opt.Photon. , 438 (2009). K. A. Willets, and R. P. Van Duyne, Ann. Rev. Phys.Chem. , 267 (2007). K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan,R. R. Dasari, and M. S. Feld, Phys. Rev. Lett. , 1667(1997). R. Zhang, Y. Zhang, Z.-C. Dong, S. Jiang, C. Zhang, L. G.Chen, L. Zhang, Y. Liao, J. Aizpurua, Y. Luo, J.-L. Yang,and J.-G. Hou, Nature , 82 (2013). R. F. Oulton, V. J. Sorger, D. A. Genov, D. F. P. Pile, andX. Zhang, Nat. Photon. 2, 496 (2008). V. J. Sorger, Z. Ye, R. F. Oulton, Y. Wang, G. Bartal, X.Yin, and X. Zhang, Nat. Commun. 331 (2011). D. E. Chang, A. S. Sørensen, E. A. Demler, and M. D.Lukin, Nat. Phys. , 807 (2007). A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L.Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Phys. Rev. Lett. , 020501 (2011). R.-C. Ge, C. Van Vlack, P. Yao, Jeff. F. Young, and S.Hughes, Phys. Rev. B , 205425 (2013). G. Sun and J. B. Khurgin, Appl. Phys. Lett. , 113116(2011). D. M. Sullivan, Electromagnetic simulation using theFDTD method, IEEE Press, (2000). K. Stannigel, M. K¨onig, J Niegemann, and K. Busch, Op-tics Express , 14934 (2009). J. M. McMahon, A.-I. Henry, K. L. Wustholz, M. J. Natan,R. G. Freeman, R. P. Van Duyne, and G. C. Schatz, Anal.Bioanal. Chem. , 1819 (2009). E. Waks, and D. Sridharan, Phys. Rev. A , 043845(2010). A. Ridolfo, O. Di Stefano, N. Fina, R. Saija, and S.Savasta, Phys. Rev. Lett. , 263601 (2010). R. A. Shah, N. F. Scherer, M. Pelton, and S. K. Gray,Phys. Rev. B , 075411 (2013). T. Gruner and D.-G. Welsch, Phys. Rev. A , 3246(1995). S. Scheel, L. Kn¨oll, and D.-G. Welsch, Phys. Rev. A ,700 (1998). B. Huttner, S. M. Barnett, Phys. Rev. A , 4306 (1992). L. G. Suttorp, and M. Wubs, Phys. Rev. A , 013816(2004). C. Sauvan, J. P. Hugonin, I. S. Maksymov, and P. Lalanne,Phys. Rev. Lett. , 237401 (2013). P. T. Kristensen, and S. Hughes, ACS Photonics , 2(2014). P. T. Leung, S. Y. Liu, and K. Young, Phys. Rev. A ,3982 (1994). K. M. Lee, P. T. Leung, and K. M. Pang, J. Opt. Soc. Am.B , 1409 (1999). R.-C. Ge, and S. Hughes, Opt. Lett. , 4235 (2014). R.-C. Ge, P. T. Kristensen, Jeff. F. Young, S. Hughes, NewJ. Phys. , 113048 (2014). R.-C. Ge, Jeff Young, and S. Hughes, Optica , 246 (2015). R. D. Artuso and G. W. Bryant, Phys. Rev. B , 125423(2013). K. S lowik, R. Filter, J. Straubel, F. Lederer, and C. Rock-stuhl, Phys. Rev. B , 195414 (2013). Jianji Yang, Mathias Perrin, and Philippe Lalanne, Phys.Rev. X , 021008 (2015). H. T. Dung, L. Kn¨oll, and D.-G. Welsch, Phys. Rev. A ,3931 (1998). S. Scheel, L. Kn¨oll, and D.-G. Welsch, Phys. Rev. A X.-W. Chen, V. Sandoghdar and M. Agio, Phys. Rev. Lett.110, 153605 (2013) L. Novotny and B. Hecht, Principles of Nano Optics, Cam-bridge University Press, (2006). E. N. Economou, Green’s Functions in Quantum Physics,Springer, (Piscataway, NJ 2006). Note that the real part of the Green function diverges asthe two space points approaching each other, which con-tributes to the homogeneous space Lamb shift of the QDs;this contribution has been included in the definition of theresonant frequencies ω n . E. A. Muljarov, M. D. Doost, and W. Langbein, e-print:arXiv:1409.6877 (2014). P. T. Kristensen, R. Ge, and S. Hughes, e-print:arXiv:1501.05938 (2015), Phys. Rev. A., in press. G. W. Bryant, F. J. Garc´ıa de Abajo, and J. Aizpurua, Nano Lett. , 631 (2008). P. Anger, P. Bharadwaj, and L. Novotny, Phys. Rev. Lett. , 113002 (2006) P. Gay-Balmaz and O. J. F. Martin, Appl. Opt. , 4562(2001). W. K. Wootters, Phys. Rev. Lett. , 2245 (1998). H. T. Dung, S. Scheel, and D.-G. Welsch, and L. Kn¨oll, J.Opt. B , S169 (2002). S. Hughes, Phys. Rev. Lett. , 227402 (2005). H. T. Dung, L. Kn¨oll, and D.-G. Welsch, Phys. Rev. A ,063810 (2002). P. Yao, P. K. Pathak, E. Illes, S. Hughes, S. M¨unch, S.Reitzenstein, P. Franeck, A. L¨offler, T. Heindel, S. H¨ofling,L. Worschech, and A. Forchel, Phys. Rev. B , 033309(2010).59