THE JOURNAL OF CHEMICAL PHYSICS 142, 224702 (2015)

Theory of molecule metal nano-particle interaction: Quantum description of plasmonic lasing Yuan Zhanga) and Volkhard Mayb) Institute für Physik, Humboldt-Universität zu Berlin, Netwonstraße 15, D-12489 Berlin, Germany

(Received 20 April 2015; accepted 12 May 2015; published online 8 June 2015) The recent quantum description of a few molecules interacting with plasmon excitations of a spherical metal nano-particle (MNP) as presented in the work of Zhang and May [Phys. Rev. B 89, 245441 (2014)] is extended to systems with up to 100 molecules. We demonstrate the possibility of multiple plasmon excitation and describe their conversion into far-field photons. The calculation of the steady-state photon emission spectrum results in an emission line-narrowing with an increasing number of molecules coupled to the MNP. This is considered as an essential criterion for the action of the molecule-MNP system as a nano-laser. To have exact results for systems with up to 20 molecules, we proceed as recently described by Richter et al. [Phys. Rev. B 91, 035306 (2015)] and study a highly symmetric system. It assumes an equatorial and regular position of identical molecules in such a way that their coupling is dominated by that to a single MNP dipole-plasmon excitation. Changing from the exact computation of the system’s complete density matrix to an approximate theory based on the reduced plasmon density matrix, systems with more than 100 molecules can be described. Finally, nonlinear rate equations are proposed which reproduce the mean number of excited plasmons in their dependence of the number of molecules and of the used pump rate. The second order intensity correlation function of emitted photons is related to the respective plasmon correlation function which approaches unity when the system starts lasing. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921724]

I. INTRODUCTION

The investigation of metal nano-particle (MNP) surface plasmons is of sustained interest (see, for example, the recent reviews1–3). Ongoing intensive research focuses not only on the properties of isolated MNPs, of clusters of interacting MNPs, but also on the behavior of hybrid systems produced when MNPs are combined with molecules, semiconductor nano-crystals or quantum dots, or if they are placed at different types of surfaces. One promising example in this respect is the concept of the plasmonic nano-laser (see Refs. 4 and 5 for a recent review). Based on strong plasmon excitations in the MNP and subsequent photon emission by the MNP, it is expected to overcome the optical diffraction limit of standard laser design. Lasing of such a nano-laser depends on the capability to achieve strong plasmon excitation. Since the plasmons are characterized by an extremely short lifetime of some tens of femtoseconds, strong excitation is necessary. This can only be achieved by a concerted coupling of, for example, many molecules6,7 or many semiconductor quantum dots8–11 to the MNP. The necessary excitation of the different types of quantum emitters themselves is realized either by optical pumping6,7,9,10 or by electrical injection.8,11 Recently, we have theoretically proven that molecular pumping is also possible if the molecules are part of a nano-junction.12 Charging and

a)Electronic address: [email protected] b)Electronic address: [email protected]

0021-9606/2015/142(22)/224702/13/$30.00

discharging of the molecules appearing at an applied voltage may drive them into an excited electronic states. According to Ref. 13, lasing by a plasmonic nano-laser shall be characterized by a narrowing of the emission line shape with increasing pump intensity. As shown by us in Ref. 12, a similar effect can be also realized with a molecular junction where the molecules are not only electrically pumped but also coupled to a MNP. Instead of using a formula of the emission linewidth,13 however, we directly computed the complete far-field photon emission spectrum which includes the necessary plasmon-to-photon conversion. Our calculations show that the emission linewidth is reduced with an increasing number of molecules (coupled to the MNP) in the junction. So far, however, these considerations have been restricted to five molecules in the junction (the direct application of density matrix theory results in an exponential increase of the density matrix elements with an increase of the number of molecules). Therefore, the behavior studied so far does not correspond to that of a nano-laser since only the second and third plasmon excited states are occupied. If the lasing effect is realized, a Poisson-like distribution of the photon state population has to be expected.14 Similarly, the lasing in the plasmonic nano-laser can only be reached if such a distribution across the plasmon states was formed.15 It is the aim of the subsequent considerations to demonstrate lasing behavior of a MNP coupled to more than five molecules. The approach is ready to consider optical pumping of the molecules as well as electrical pumping. Since in the latter case, the molecules are part of a molecular junction;

142, 224702-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-2

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

the whole description becomes more involved and, therefore, is shifted to a forthcoming paper. To achieve a description of up to 100 molecules coupled to the MNP, we proceed in different steps. First, we consider a highly symmetric molecule-MNP model. Here, every molecule couples in the same way to the MNP and the number of density matrix elements to be calculated reduces drastically. This has also been discussed recently in Ref. 15. For us, it serves as a reference computation providing data to justify subsequent approximate considerations. We also note that this treatment is known from cavity-QED theory with many atoms.16 There, the density matrix has been expanded with respect to socalled Dicke states (see also Ref. 17). In calculating the farfield photon emission spectrum, we also consider in detail the plasmon to photon conversion (this also continues the studies of Ref. 15). Although our computations are restricted to MNP dipole plasmons only, we also comment on the possible influence of higher multipole plasmons also excited in the considered spherical MNP. While the former exact description is restricted to 20 molecules, the subsequently presented approximate treatment allows to handle 100 and more molecules, and to deviate from the highly symmetric molecule-MNP system. This theory is based on approximate equations obeyed by the reduced plasmon density matrix. Finally, nonlinear rate equations are introduced which are valid for molecular level populations and the mean plasmon number. They are rather simple but do not allow to compute statistical properties of the excited plasmons as well as the photon emission spectrum. However, in difference to what has been discussed in Ref. 13, rates of spontaneous energy transfer from the molecules to the MNP are also involved. The paper is organized as follows. In Sec. II, we introduce the used model of the molecule-MNP system and the basic density matrix equations. The exact quantum dynamics of the symmetric molecule-MNP system are described in Sec. III. Afterwards in Sec. IV, approximate equations of motion for the plasmon density matrix are introduced. Related computational details can be found in Appendix. There, we also derive the rate equations which are finally used to discuss the properties of the molecule-MNP system in Sec. V. All approximate treatments are related to the exact data of the symmetric system with up 20 molecules. The paper ends with some concluding remarks in Sec. VI. II. MOLECULE-MNP MODEL AND KINETIC EQUATIONS

To highlight the specifications of the system we will introduce in the following, let us start with a rather general model of a MNP coated by molecules (see Fig. 1). The respective Hamiltonian is written as H = Hmol + Hpl + Hmol−pl.

(1)

It covers the molecular part Hmol, the MNP plasmon part Hpl, and the coupling Hmol−pl between both subsystems. The molecules determining Hmol shall be considered as electronic multi-level systems with adiabatic electronic energies Ena . The index n counts the molecules and a is the electronic quantum

FIG. 1. Scheme of some molecules placed in the proximity of a spherical MNP. The molecules are put in an equatorial position to the MNP (x–y plane) with dipole moment in z-direction. An identical coupling of all molecules to the MNP’s dipole moment pointing in z-direction is assumed. The chosen molecule MNP-surface distance of ∆x mol−MNP = 2.5 nm and the coincidence of the molecular excitation energy with the MNP dipole plasmons result in a minor influence of higher multipole plasmons.

number (ground-state: a = g, first excited state: a = e, etc.). The related electronic wave functions are written as ϕ na . In a first step, we have in mind off-resonant optical pumping of the molecular system, i.e., a higher lying level a = f is resonantly addressed. But, finally, via internal conversion, the first excited molecular level may become strongly populated. This population is of interest since the transition from the first excited state to the ground-state couples to the MNP. Since details of the population of the first excited state will be of no interest, we describe the pumping as an incoherent process moving directly from the molecular ground to the first excited state. The respective rate referring to molecule n is k n,g → e . It enters the density matrix description introduced below. To be complete, we also introduce the excited state decay rate k n,e→ g of the molecules. Moreover, we ignore any higher excited state of the molecule and denote the molecular Hamiltonian as the one of a collection of two-level molecules,   Hmol = En Bn+ Bn + Jnn′ Bn+ Bn′. (2) n, n ′

n

This standard notation (for a recent description, see Ref. 18) uses molecular excitation energies En = Ene − Eng and transition operators Bn+ = |ϕ ne ⟩⟨ϕ ng | and Bn = |ϕ ng ⟩⟨ϕ ne | (note that  the expression assumes Eg = n Eng = 0). Excitation energy exchange among different molecules is accounted for by Jnn′ (standard excitonic coupling18,19). Be also aware of the fact that the used two-level model does not allow to describe excitonexciton annihilation, what represents an important quenching process if many molecules are excited (see, for example, Ref. 18). The collective electronic excitations of the used MNP of spherical shape are multipole plasmons which are described by the Hamiltonian  (3) Hpl = El Cl+mCl m . l, m

The excitation energies are El = (3l/(2l + 1))1/2 × El=1. The dipole plasmon energy El=1 can be fixed by the experiment, but

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-3

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

Because the extension of the molecule-MNP system is in the 20–30 nm range, the internal coupling can be accounted for by the (non-retarded) Coulomb interaction of all involved charges. Moreover, if negative (electronic) and positive (nuclear) charge distributions in the molecules are well balanced locally (absence of permanent dipole moments), an electrostatic coupling of the molecules to the MNP via mirror charge formation is of minor importance. It remains a type of Coulomb-coupling which is responsible for excitation energy exchange between the molecules and the MNP. According to the above given notation, the molecule-MNP coupling reads  ∗ + Hmol−pl = (Vn,l m Bn+Cl m + Vn,l (4) m BnCl m ). n,l, m

The coupling matrix elements take the form20 FIG. 2. Energy transfer rate k n→ l , Eq. (6) between a single molecule and the MNP multipole plasmons l = 1, . . ., 10 drawn versus the molecular excitation energy detuning E mol − E l =1 (parameters according to Table I).

the higher multipole plasmon energies shall follow the wellknown Mie-formula20 (l = 1, 2, . . . and m = −l, . . . , 0, . . . , l are multipolar indices). The Cl+m , Cl m are related harmonic oscillator operators (cf. Ref. 12). Such a harmonic description of quantized plasmon oscillations was demonstrated in Refs. 21–23 as an acceptable approximation for dipole plasmons. We also take the harmonic description for higher multipole plasmons since it is the most plausible approximation at present. Note also that the neglect of any coupling among the different types of plasmons is also not obvious but has been used for simplicity. Fig. 2 gives a graphical representation of the multipole plasmon levels using the parameters collected in Table I (and motivated in Sec. II B). We have drawn the energy transfer rate between the molecule and the various multipole excitations (see below) by varying the molecular excitation energy. The extrema of all curves point at the position of the respective multipole plasmon energy El . All plasmon levels are characterized by the decay rate γpl including radiative and non-radiative contributions (the dephasing rate is γpl/2). In principle, one expects a dependence of the decay rate on the plasmon type. Since nothing is known on this dependence, we introduced a common rate. For the used parameters, we get the following detuning among the lowest multipole levels: El=2 − El=1 = 248 meV and El=3 − El=2 = 100 meV which is larger than the lifetime broadening.

 Vn,l m =

4π Yl m (θ n , φ n ) . eQ l m [dn ∇X] 2l + 1 X nl+1

(5)

It covers the ground-state to the excited state transition dipole moments dn of the molecules n. The multipole moments related to the MNP-plasmon transition are given  2l+1 by eQ l m = r MNP El /2 (r MNP denotes the radius of the MNP). X is a vector connecting the center of the MNP with that of molecule n. The related spherical coordinates are (X n , θ n , φ n ) entering also the spherical harmonics Yl m . Of course, the strength of excitation energy exchange depends on the magnitude of the matrix elements Vn,l m of the molecule-MNP coupling. As indicated in Fig. 2, however, the actual detuning among the molecular energies En and the various multipole plasmon energies is also of importance. This combined effect determining energy exchange is clearly demonstrated via the Golden Rule type energy transfer rate, k n→ l

 2γ lm=−l |Vn,l m |2 = . (En − El )2 + (~γ)2

(6)

The dephasing rate γ is γpl + k n,g → e + k n,e→ g . Fig. 2 displays k n→ l versus En − El=1 and for different multipole plasmons, indicating that they are energetically well separated. This suggests that we can focus on dipole plasmons if the molecular excitation is in resonance to El=1. A. Density matrix theory and emission spectrum

TABLE I. Used parameters (for explanation, see text). E mol d mol ~k ~p ~ω pl ~γ pl d pl Npl ∆x mol−MNP V

2.6 eV 16 D 3 meV 50 meV 2.6 eV 57 meV 2925 D 25 2.5 nm 15.95 meV

The strategy of the following treatment is to handle the system represented by the Hamiltonian Eq. (1) or (10) as the active system in the sense of dissipative quantum dynamics and to introduce molecular pumping as well as molecular and plasmon decay as dissipative processes. This standard treatment has been used recently by us in Ref. 12 and was also applied in Ref. 15. It is based on the following equation for the (reduced) density operator: ∂ ρˆ i = − [H, ρ] ˆ − − D ρ. ˆ ∂t ~

(7)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-4

Y. Zhang and V. May

Here, H is the active system Hamiltonian. The dissipative part is written in the so-called Lindblad-form as  ku
[ Lˆ +u Lˆ u , ρ] ˆ + − 2 Lˆ u ρˆ Lˆ +u . (8) D ρˆ = 2 u Considering molecular excited state decay, the Lˆ +u -operators have to be identified by Bn+ and the k u by k n,e→ g = k n . Offresonant optical pumping is described with k u = k n,g → e = pn and Lˆ +u = Bn . The MNP plasmon damping is described by identifying Lˆ +u with Cl+m and k u with the common decay rate γpl. Since dissipation covers plasmon and molecular decay as well as molecular pumping, the density operator approaches the steady-state expression ρˆ ss at large times which, in particular, involves a non-equilibrium excited plasmon state distribution. A measurable quantity which is of particular interest here is the frequency-resolved emission spectrum. It is written as F(ω) = 4ω3 I (ω) /3πc3~ (cf. Ref. 12) and determines the probability of emitting photons at frequency ω and under steadystate conditions. The related emission line shape function reads  ∞  I(ω) = Re dt e−iωt [d Ad∗B]tr{ Xˆ A+ σ(B; ˆ t)}. (9) 0

A, B

Here, A and B label either the different molecules or the multipole plasmons, and Xˆ A+ and Xˆ B denote the respective transition operators. If A refers to molecule n, we have + (for the inclusion Xˆ A+ = Bn+. Otherwise, it follows Xˆ A+ = Cl=1m of higher multipole plasmons cf. Ref. 24). d A and d∗B abbreviate related transition dipole moments (note that in the case of a spherical MNP, a representation using spherical coordinates becomes necessary). Moreover, we introduced in Eq. (9) the quantity σ(B; ˆ t) = U(t)( Xˆ B ρˆ ss). Time propagation is realized by the time-evolution superoperator U(t), which represents a formal solution of density operator equation (7). The initial value of σ ˆ combines the transition operator Xˆ B with the steady-state (reduced) density operator ρˆ ss. B. Highly symmetric molecule-MNP model

By introducing various specifications, a certain simplification of the general molecule-MNP model explained so far becomes possible. First, we assume only a partial decoration of the MNP with one type of dye molecules. Consequently, the inter-molecular distance shall amount to some nanometer, and the excitation energy transfer coupling Jnn′ may stay small enough to be ignored (exciton-exciton annihilation becomes also of minor importance). We also assume an additional separation between the molecules and the MNP surface (achieved, for example, with a coating of optically inactive molecules). Then, the molecules remain unaffected by so-called MNP spillout electrons. And, we can adjust a distance of the optically active molecule to the MNP surface to let dominate the coupling to dipole plasmons. Finally, the molecules shall be arranged around the MNP’s equator with the transition dipole moment tangential to the MNP surface. This is shown in Fig. 1 and underlines that all the molecules exclusively couple to the dipole plasmon with m = 0, i.e. we set Vn,l m → V = Vl=1, m=0.

J. Chem. Phys. 142, 224702 (2015)

The Hamiltonian, Eq. (1), reduces to   H= Emol Bn+ Bn + EplC +C + (V Bn+C + V ∗ BnC +). n

(10)

n

It describes a one-dimensional quantum oscillator driven by  the composite force n V ∗ Bn . In all the simulations, we take a spherical Au particle with 20 nm diameter. Each dipole plasmon mode has an excitation energy of El=1 = ~ωpl = 2.6 eV, a transition dipole moment of d pl = 2925 D, and the decay rate amounts to ~γpl = 57 meV. The molecular excitation energy shall be in complete resonance to the dipole plasmon, i.e., Emol = ~ωmol = 2.6 eV. The molecular transition dipole moment is chosen as 16 D, which is in line with our previous study.12 We have chosen a uniform molecule MNP-surface distance of ∆x mol−MNP = 2.5 nm. For the configuration shown in Fig. 1, the molecules couple only with one dipole plasmon mode. The molecular excited-state decay rate is taken as ~k = 3 meV (it mimics the influence of the other MNP multipole plasmon excitations, see the discussion in Ref. 12). If not stated, otherwise, the molecular pump rate is chosen as ~p = 50 meV (p = 0.076/fs), what is an optimal value to achieve a multiple plasmon excitation. C. Molecule-MNP model including electrical pumping

We briefly indicate that the used model can be easily modified to change from the case of molecular excitation via optical pumping to the case of electrical pumping. As described in Ref. 12, the molecules shall be a part of a molecular junction. There, they are contacted by two metal nanoelectrodes (left and right X = L, R), one of which may be ready to generate plasmons resonant to the molecular excitation. Due to an applied voltage between the electrodes, the molecules may become charged and discharged. In a simple variant, one assumes that the molecules only undergo a singly negatively charging (acceptance of a single excess electron). If the junction and the applied voltage are properly chosen, the molecules may enter their excited state via discharge. To describe the processes mentioned beforehand, the molecular Hamiltonian, Eq. (2), is extended by  (−) Hmol = En f |ϕ n f ⟩⟨ϕ n f |, (11) n

where En f is the energy of the nth molecule charged state and ϕ n f is the respective wave function. Sequential charging and discharge of the molecules what will be considered as dissipative processes replace the optical pumping in D, Eq. (8). We now identify Lˆ +u by |ϕ na ⟩⟨ϕ n f |, where a = e, g, and k u by   (X ) k n,a→ f = k a→ = Γa(Xf ) f X

×



X

dEel f F(Eel − µ X )δ(En f − Ena − Eel).

(12)

While we labeled the general rate k n,a→ f with the molecular index n, we ignored this for the other parts of the formula. Γa(Xf ) is the molecule-lead coupling function. It refers to lead X and describes charging from the ground or excited state of the neutral molecule. It has been taken to be independent of the lead’s electron energy Eel. Charging becomes possible if

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-5

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

the latter energy combined with Eng or Ene coincides with En f . The Fermi-distribution f F regulates if the electron states of energy Eel are occupied or not. The actual lead chemical potential is µ X =L, R = EF ± |e|V/2 (for the case of a symmetrically applied voltage). In the case of molecular discharge, we take Lˆ +u = |ϕ n f ⟩⟨ϕ na | and k u = k n, f → a with the latter rate obtained from Eq. (12) by replacing f F by 1 − f F (note that the used expression neglects any intermolecular correlation due to charge transmission, cf. the discussion in Ref. 12). This short explanation indicates the possible inclusion of the effect of electrical pumping of the molecules (for more details see our forthcoming paper).

III. EXACT QUANTUM DYNAMICS OF THE SYMMETRIC MOLECULE-MNP SYSTEM

To generate exact data for a molecule-MNP system which is larger than that treated by us in Ref. 12, we will study the steady-state behavior of the symmetric system introduced in Sec. II B, i.e., we consider all molecules as identically coupled to one type of dipole plasmons (m = 0, polarized in z-direction). The inter-molecular distance should be large enough to let become the excitonic coupling as well as excitonexciton annihilation of minor importance. Such a system has been already investigated by us in Ref. 12 for up to 5 molecules (but for the case of electrical pumping of the molecules as introduced in Sec. II C). Recently, it was demonstrated that the description can be extended to 30 molecules if symmetry properties of the related density matrix equations are utilized which are due to the symmetries in the model.15 The main motivation to repeat these computations here is their extension to the consideration of the photon emission spectra. Moreover, we will use these data to confirm a subsequent approximate computation which is ready to consider systems with more than 100 molecules. An exact description of the driven dynamics of the symmetric molecule-MNP system is achieved by introducing the full density matrix. This is possible in using a complete basis which corresponds to the system Hamiltonian, Eq. (10). The most obvious variant of such a basis would be  |α µ⟩ = |α⟩ × | µ⟩ = |ϕ na ⟩ × | µ⟩. (13) n

The molecular product state is characterized by α what comprises the quantum numbers a = g or a = e of the individual molecules n = 1, . . . , Nmol (Nmol is the number of molecules). µ labels the actual excitation of the dipole-plasmon oscillator (pointing in z-direction). Then, the density matrix to be

determined is ρα µ, βν (t) = tr{ ρ(t)| ˆ βν⟩⟨α µ|} ≡ ⟨| βν⟩⟨α µ|⟩.

(14)

Respective equations of motion can be deduced from Eq. (7) (cf. Ref. 12). They resemble those used in the cavity-QED theory formulated for many atoms.16 The only difference to the present consideration is the use of so-called Dicke states as an expansion basis. To solve the equations of motion for ρα µ, βν , we need to introduce a matrix with 4 Nmol × (Npl + 1)2 elements (Npl is the highest plasmon state considered). Obviously, the size of the matrix increases exponentially with Nmol. Consequently, the direct computation of ρα µ, βν is only possible for systems with a few molecules (see Ref. 12). However, when focusing on the system where all molecules are identical and couple in an identical way to the MNP, many density matrix elements coincide. To distinguish between those which are different, we can replace α and β by the set of four numbers (c, d = g, e), n = (nc d ) = (ng g , ng e , neg , nee ),

(15)

where nc d =

N mol 

δ c,a l δ bl, d .

(16)

l=1

The al describe the electronic state of molecule l and the whole set {al } is identical with α. The quantum number β is represented by the bl . In this manner, nc d indicates how many molecules of product state |α⟩ (cf. Eq. (13)) are in the single molecule state ϕc and, simultaneously, how many molecules of product state | β⟩ are in the single molecule state ϕ d . When establishing equations of motion for ρα µ, βν , the transition frequencies (Eα − E β )/~ appear (Eα , E β denote the energy of the respective molecular product state). For the symmetric system, the energy Eα (E β ) is identical for all states |α⟩ (| β⟩) where the number of molecules in the excited state is identical. Consequently, we may write the transition frequencies alternatively as (neg − ng e )ωmol. The term with neg gives the number of all molecules which are in state |α⟩ in their excited state and which are in state | β⟩ in their ground-state. Those contribute to the transition frequency. The term with ng e gives the other number which also contributes to the transition frequency (molecules which are in state |α⟩ in their ground-state and which are in state | β⟩ in their excited state). In a similar manner, one can translate other terms of the original density matrix equations. Accordingly, ρα µ, βν is replaced by ρ(µ,ν) (the plasmon n quantum numbers remain unchanged) and the equations of motion turn into the following form:

)  γpl ( ∂ (µ,ν) (µ + ν) ρ(µ,ν) ρn = −i(neg − ng e )ωmol ρ(µ,ν) − iωpl(µ − ν)ρ(µ,ν) − − 2 (ν + 1)(µ + 1)ρ(µ+1,ν+1) n n n n ∂t 2 ( ) p( )

(µ,ν) k (µ,ν) − neg + 2nee + ng e ρ(µ,ν) − 2ng g ρ(µ,ν) − 2n + n + n ρ − 2n ρ g g g e eg ee n n (n g g −1, n ee +1) (n ee −1, n g g +1) 2 2  (√ ) √  √ (µ+1,ν) (µ,ν+1) (µ−1,ν) + iv νnag ρ(µ,ν−1) − µ + 1n ρ + ν + 1n ρ − µn ρ . ea ae g a (n ag −1, n ae +1) (n ea −1, n g a +1) (n ae −1, n ag +1) (n g a −1, n ea +1)

(17)

a

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-6

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

Since all molecules are identical, the molecular decay rate k, the molecular pump rate p, and the molecule-MNP coupling V = ~v do not depend on the molecular index n. The also accounts for a change of two given notation of ρ(µ,ν) n components of n. We indicate such change by explicitly denot-
ing the changed components,
for example, ng g − 1, ng e + 1 ≡ ng g − 1, ng e + 1, neg , nee . The vector n satisfies the following conditions. Each component is an integer in the range [0, Nmol] and the sum of the four components equals Nmol. The number of different vectors n is obtained by asking on the number of possibilities to distribute Nmol indistinguishable balls (units) on four different places (the four components of n). The result is given N +3 by the binomial coefficient (3 mol ) ≡ C(Nmol + 3, 3) = (1/6) (Nmol + 3) (Nmol + 2) (Nmol + 1). Noting also the number Npl + 1 of different considered plasmon states, the number of elements forming
ρnµ,ν is ntot = (1/6) (Nmol + 3) (Nmol + 2) 2 (Nmol + 1) Npl + 1 (polynomial increases with increase of Nmol). In order to analyze the non-equilibrium steady state of the system, we introduce the following excited state population: µ) Pm µ (t) = ρ(µ, (t). (N −m,0,0, m) mol

(18)

It is the probability to realize a single excited-state configuration with m = nee excited molecules and µ excited plasmons. The probability to find m molecules in their excited state is  Pm (t) = Pm µ (t). (19) µ

We further note that the total number of excited-state configurations with m excited molecules is given by the binomial coefficient C(Nmol, m). Therefore, the plasmon state population follows as Pµ (t) =

N mol 

C(Nmol, m)Pm µ (t).

(20)

m=0

In a further step, we adopt the formula for the emission line shape, Eq. (9), to the symmetric molecule-MNP system considered here as follows:  ∞ 2 Re I(ω) = d pl dt e−iωt 0  √ µ−1, µ µC(Nmol, ng g )σ(n (t). (21) × −n g g ) g g ,0,0, N µ, n g g

mol

The expression notices the dominance of plasmon dipole moments in Eq. (9). Therefore, the respective molecular contri2 butions have been ignored and it only results the prefactor d pl . The σ-matrix replaces the general type σα µ, βν (t) as in the case of the density matrix explained beforehand. And, σn(µ,ν) is solution of Eq. (17) but with the initial condition  also a (µ+1,ν) µ + 1ρn . The involved density matrix is the steady state solution of Eq. (17). In order to display the emission linewidth, we fitted F(ω) by a Lorentzian. A. Variation of pump rate

To demonstrate the influence of the molecular pump rate, it suffices to consider the case of five molecules. Fig. 3 displays

FIG. 3. steady-state properties of five identical molecules coupled to the MNP. Variation of the molecular pump rate p (other parameters according to Table I). Upper panel: plasmon state population P µ . Lower panel, blue squares: linewidth (FWHM) of the photon emission spectrum and red triangles: emission maximum.

the respective results. Using small pump rates p, only a small population of the first excited plasmon state appears. Higher excited plasmon states become populated when increasing ~p up to about 50 meV. The plasmon state population Pµ decreases if p is further increased. An optimal value of p is also suggested by the lower panel of Fig. 3 displaying the maximum of the emission and its width. The decrease of Pµ with increasing p takes place if the excitation energy transferred back to the molecules dominates and further pumping is hindered (saturated; see also the respective analysis in Sec. IV where an approximate treatment of the system is introduced). Notice that the emission linewidth is much smaller than both the plasmon decay rate ~γpl = 57 meV (it determines the emission linewidth of an isolated MNP) and the molecular pump rate ~p = 50 meV (it determines the emission linewidth of an isolated molecule). The decrease of the linewidth and a later increase with increasing p goes along with an increase and decrease of plasmon excitation, respectively. Here, the increase of plasmon excitation can be understood as an increase of an overall lifetime. As a result, the energy

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-7

Y. Zhang and V. May

uncertainty of the emitted photons and thus the emission line-broadening are decreased.25 Although the optimal value of p somewhat changes with the number of molecules coupled to the MNP (not shown), we use in all following computations ~p = 50 meV. B. Variation of the number of molecules

Having fixed the rate which pumps the molecules into their excited state, we now investigate the behavior of the system with an increasing number Nmol of involved molecules. Such a dependence is demonstrated in Fig. 4. If Nmol is increased, the plasmon state population Pµ gets maximal in a certain range of µ, and this range moves to higher excited plasmon states. For system with 20 molecules involved, even the fourteenth excited state is populated. Pµ resembles a Poisson distribution indicating the strong driving of the MNP plasmons by the action of the molecules. Something like a coherent state of plasmons is formed. This statement is in line with the narrowed photon emission shown in the lower panel of Fig. 4 (cf. also Ref. 14). In particular, the lower panel demonstrates that for a

J. Chem. Phys. 142, 224702 (2015)

system with 20 molecules, the linewidth is almost one sixth smaller and the emission maximum is 255 times larger as in the case of a single molecule. We expect that for systems with a further increased number of molecules, the center of the plasmon state population further shifts to higher excited states. At the same time, the emission peaks get larger and the linewidth is narrowed. Next, we will introduce an approximate theory which is ready to consider much larger numbers of molecules and which confirms this expectation.

IV. APPROXIMATE EQUATION OF MOTION FOR THE PLASMON DENSITY MATRIX

The consideration of the symmetric molecule-MNP system in Sec. III has been carried out to perform exact computations for up to 20 molecules coupling to the MNP. To further increase the number, an approximate description becomes necessary. It is based on the (reduced) plasmon density matrix, ρ µν (t) = tr{ ρ(t)|ν⟩⟨µ|}. ˆ

(22)

Its introduction and the following treatment was inspired by textbook laser theory, where a somewhat similar consideration of the photon density matrix is described.26 Noting general density operator equation (7), one easily verifies the following equation for ρ µν :
∂ ρ µν = −i (µ − ν)ωpl − i(µ + ν)γpl/2 ρ µν ∂t  + γpl (µ + 1)(ν + 1)ρ µ+1,ν+1   √ (n) −i µρe µ−1,gν vn µ + 1ρ(n) gµ+1,eν + n

√

− ν ρ(n) gµ,eν−1 −

√


ν + 1ρ(n) e µ,gν+1 .

(23)

The molecule-MNP coupling matrix elements ~vn = Vn,l=1m=0 may differ from molecule to molecule. We can also allow for varying molecular excitation energies En = ~ω n . The new type of density matrices emerging in Eq. (23) reads (a, b = g, e) ρ(n) = tr{ ρ(t)|ϕ ˆ nb ⟩⟨ϕ na | × |ν⟩⟨µ|}. a µ,bν

(24)

This quantity correlates a plasmon excitation with a molecular transition or with a molecular state population. The subsequent equations for the ρ(n) generate a hierarchy of dena µ,bν sity matrices including contributions of two molecules, three molecules, etc. Since the appearance of such higher density matrices is connected with the coupling matrix elements vn as a prefactor, the treatment of the hierarchy corresponds to a perturbation theory with respect to the molecule-MNP coupling. As demonstrated below, a second-order consideration would be sufficient. Details of the computations are collected in Appendix A. As a final result valid again in a steady-state situation, we obtain the following balance equation obeyed by the plasmon state populations Pµ = ρ µ µ : (k µ + µγpl)Pµ = pµ Pµ−1. FIG. 4. steady-state properties of up to 20 identical molecules coupled to the MNP (parameters according to Table I). Upper panel: plasmon state population P µ . Lower panel, blue squares: linewidth (FWHM) of photon emission spectrum, red triangles: emission maximum (both curves versus number of molecules).

(25)

The decay of plasmon level µ with rate k µ + µγpl is compensated by plasmon excitation with pump rate pµ . These new rates are formed by the molecular pumping rates pn , the decay rates k n , and the molecule MNP energy transfer rates, cf.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-8

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

Eqs. (A22) and (A23). In Appendix A, we also claimed that k µ combines energy transfer from the MNP to the molecules and subsequent molecular decay, while pµ is a combination of molecular pumping and energy transfer into the MNP. Here, we only quote these rates for the case where the molecular excitation is in complete resonance with the dipole plasmon (note δ n = k n/2 + pn/2), kµ =

 n

µvn2 k n δ n (δ n + γpl/(4µ)) + 2µvn2

(26)

µvn2 pn . δ n (δ n + γpl/(4µ)) + 2µvn2

(27)

and pµ =

 n

In a first step, we used Eq. (25) to reproduce Pµ versus the pump rate p as discussed in Sec. III and drawn in Fig. 3. In order to analyze the behavior we specify balance equation (25) to a recursion relation valid for five identical molecules, which also identically couple to the MNP 2

Pµ 20v p . = Pµ−1 20v 2 k + γpl(k + p)2 + 4µv 2γpl

(28)

This population ratio passes through a maximum with increasing pump rate p what explains the behavior of the distribution Pµ as shown in Fig. 3. When calculating the plasmon state populations for different numbers of molecules Nmol, we first concentrate on the cases Nmol = 1, . . . , 20 already treated within the exact computation scheme (upper panel of Fig. 4). Since the resulting data are perfectly reproduced, we do not present them again but change to the cases Nmol = 20, . . . , 100 (see Fig. 5). The population distributions are shifted to higher plasmon excited states. In addition, they become broader and their maximum decreases. For example, if Nmol = 20, the distribution covers about 20 states which are noticeably populated and the maximum is located at µ = 6. For Nmol = 100, the population is distributed across 40 states and has its maximum at µ = 40. The behavior of the plasmon state population can be easily justified by an analysis of Eq. (25) which simplifies in the

FIG. 5. steady-state plasmon state population P µ for up to 100 identical molecules coupled to the MNP (parameters according to Table I). Computations based on the plasmon density matrix introduced in Sec. IV.

present case to (note κ = (2v/(k + p))2), Pµ Nmol p = . Pµ−1 Nmol k + (1/κ + 2µ)γpl

(29)

Let us start with a single molecule coupled to the MNP. According to the parameters given in Table I, the population ratio for µ = 1 is smaller than one and the ratios for larger µ are negligible (cf. upper panel of Fig. 4). If Nmol is increased, the ratios Pµ /Pµ−1 increase but all of them are still smaller than one (see the curve of Fig. 4, upper panel for Nmol = 5). A further increase of Nmol results in the ratio close to one for some critical number µc . The ratios for other µ , µc are much smaller than one. We could estimate this number µc by setting the ratio equal to one. Obviously, this critical number increases with increasing Nmol, which explains why the population distribution shifts to higher plasmon excited states (see Figs. 4 and 5 for Nmol ≥ 10). A. Second order intensity correlation function of emitted photons

The steady-state emission spectra of the molecule-MNP system introduced so far are related to the photon’s first-order intensity correlation function.27 To obtain more information about the light source formed by the molecule-MNP system, we consider the second-order (time-equal and normalized) intensity correlation function g (2)(0) = ⟨a+a+aa⟩/⟨a+a⟩2 (a and a+ denote photon operators of a particular mode). Following Ref. 28, we assume proportionality to the respective plasmon correlation function gpl(2)(0) = ⟨C +C +CC⟩/⟨C +C⟩2   = µ(µ − 1)Pµ /( Pν )2. µ

(30)

ν

This quantity drawn versus the pump rate p and the number of molecule Nmol is shown in Fig. 6. gpl(2)(0) is larger than one

FIG. 6. Equal time second-order plasmon correlation function versus molecular pump rate ~p and the number of molecules identically coupled to the MNP (parameters according to Table I).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-9

Y. Zhang and V. May

for all the chosen values of ~p and Nmol. The function gpl(2)(0) approaches one in the range ~p < 10 meV and Nmol < 20 and also in the range ~p > 10 meV and Nmol > 20. In the former case, the system is not lasing what is also indicated by the plasmon state population (not shown). Therefore, the equality gpl(2)(0) = 1 alone does not indicate lasing operation. For a fixed molecular pump rate, gpl(2)(0) increases first and then decreases with increasing number of molecules. The behavior is similar for a fixed number of molecules and an increase of the pump rate. In the latter case, gpl(2)(0) decreases with increasing pumprate and for a fixed number of molecules. It decreases with increasing number of molecules and with fixed pump-rate.

V. RATE EQUATION APPROACH

An approximation scheme somewhat alternative to that of Sec. IV is based on the use of diagonal density matrices

J. Chem. Phys. 142, 224702 (2015)

describing the molecular level populations Pna , Eq. (B1), and the mean number Npl of excited plasmons, Eq. (B4). All these quantities obey the approximate nonlinear rate equations as derived in Appendix B. A comparison of Npl calculated within different approximation schemes is given in Fig. 7. Additionally, we also present data based on the rate equations derived in Ref. 13. There, the molecular dephasing does not depend on the molecular pump rate. This leads to a linear increase of the plasmon mean number with increasing pump rate. However, in our equations, the behavior is different since the molecular dephasing rate is given by (k + p)/2. Compared to Eqs. (B6) and (B7), the equations of Ref. 13 do not include the spontaneous energy transfer contribution k n Pne . This is equivalent to a neglect of plasmon quantum fluctuations. The upper panel of Fig. 7 compares the change of Npl with the molecular pump rate as obtained within different approximations. Here, five molecules are considered. Since the green squares are always larger than the black dots, the solution of the plasmon state population balance equation overestimates the plasmon excitations somewhat. Although the blue solid line (solution of nonlinear rate equations) does not completely coincide with the black dots (exact result), they follow the same tendency. On the contrary, the red dashed line (nonlinear rate equations ignoring the term k n Pne ) starts to increase from a pump rate with ~p ≈ 5 meV. This demonstrates the so-called threshold behavior. This approximate form of Npl versus p (red dashed line) arrives at its maximum around 38 meV and then starts to decrease afterwards. It becomes zero when the molecular pump rate is larger than 70 meV. Such a behavior is due to the assumption ⟨C +C⟩ = ⟨C +⟩ ⟨C⟩ taken to arrive at rate equations used in Ref. 13. The assumption implies that the MNP is strongly excited, i.e., Npl ≫ 1. For the results shown here, this assumption is not valid since Npl < 1. The lower panel of Fig. 7 displays the change of Npl with the number of molecules Nmol (the molecular pump rate is taken as ~p = 50 meV). Again, different approximations are compared with the exact result. The solution of the plasmonstate population balance equation overestimates a little bit the exact result. The data obtained by a solution of the nonlinear rate equations nearly coincide with the exact result as long as Nmol < 5. For larger Nmol, the exact results are somewhat underestimated. VI. CONCLUSIONS

FIG. 7. Mean plasmon number Npl calculated within different approximations (parameters according to Table I). Upper panel: variation of the molecular pump rate p for the case of five molecules identically coupled to the MNP, lower panel: variation of the number of molecules. Black dots: use of the exact density matrix equations (17), green squares: use of Eq. (25) obtained in calculating the plasmon density matrix, blue solid curve: use of rate equations (B6) and (B7), red dashed curve: as blue solid curve but neglecting the terms k n P ne .

We have generalized our recent quantum description restricted to a few molecules interacting with MNP plasmon excitations of Ref. 12 to systems with up to 100 molecules. While the former work concentrated on the case of electrical pumping of the molecules here, we focused on optical pumping. A symmetric system with up to 20 molecules could be described exactly similar to Ref. 15. In extending this recent study, we also describe the conversion of plasmon excitations into far-field photons. The steady-state photon emission spectrum displays emission line-narrowing with an increasing number of molecules coupled to the MNP. This is considered as an essential criterion for the action of the molecule-MNP system as a nano-laser.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-10

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

Systems with up to 100 molecules were described in changing from the exact computation of the system’s complete density matrix to an approximate theory based on the reduced plasmon density matrix. The quality of the used approximation was confirmed by a reproduction of the results derived in the framework of the exact theory. This approximation has been also used to compute the second-order intensity correlation function of the emitted photons. It was related to the respective plasmon correlation function which approaches one when the system starts lasing. Finally, nonlinear rate equations are proposed which reproduce the mean number of excited plasmons in their dependence of the number of molecules and the used pump rate. To confront the exact description of a symmetric molecule-MNP system with some approximate treatments paves the way to exclude some assumptions we took here. For example, we may consider changes of the molecular excitation energies from molecule to molecule. The excitonic coupling, which is of some importance for the case of a closer coating of the MNP with molecules, can be also accounted for. We may introduce some higher multipole plasmons into the description. And, our theory can be easily extended to the case of

three-level molecules. This becomes necessary when studying plasmonic lasing due to electron injection in a molecular junction.12 ACKNOWLEDGMENTS

Financial support by the Deutsche Forschungsgemeinschaft through Sfb 951, as well as by the China Scholarship Council (Y.Z.) are gratefully acknowledged. We also thank Yaroslav Zelinskyy, Dirk Ziemann, and Thomas Plehn for several illuminating discussions. APPENDIX A: THE REDUCED PLASMON DENSITY MATRIX

We explain in somewhat more detail the computation of the plasmon density matrix ρ µν (t), Eq. (22), also discussed in Sec. IV. The respective equation of motion Eq. (23) generates the new type of density matrix, Eq. (24), combining molecular and plasmon excitations. To proceed further, we need to set up equations of motion for this new type of density matrices. We start with a presentation of the following equation of motion:

(√ )  ∂ (n) √ (n) ρgµ,eν−1 = i(ω˜ ∗n + ω µν−1)ρ(n) µρe µ−1,eν−1 − γpl([(µ + ν − 1)/2]ρ(n) ) + ivn ν ρ(n) − (µ + 1)ν ρ(n) gµ,gν − gµ,eν−1 gµ+1,eν gµ,eν−1 ∂t  √  √ √ (n, n′) n ′) (n, n ′) n ′)
µρgeµ−1,egν−1 − ν ρ(n, (A1) −i vn′ µ + 1ρ(n, geµ,egν . ggµ+1,eeν−1 − ν − 1ρggµ,eeν−2 − n′

The plasmon transition frequencies are defined as ω µν = (µ − ν)ωpl. Note also the introduction of ω˜ n = En/~ − iδ n , where δ n = (k n + pn )/2. There also appear more complex density matrices which are defined as (n ′ , n) ′

n) ρ(n, = trS{ ρ(t)|ϕ ˆ nc ⟩⟨ϕ na | × |ϕ n ′d ⟩⟨ϕ n ′b | × |ν⟩⟨µ|}. ab µ,c dν

(A2)

(n) Be aware of the relation ρ µν = ρ(n) gµ,gν + ρe µ,eν which results from the simple completeness relation of the two considered molecular electronic levels. The equations of motion for the other density matrices entering the equation for ρ µν read

( )  ∂ (n) (n) (n) ρe µ−1,gν = −i(ω˜ n + ω µ−1ν )ρ(n) µ(ν + 1)ρ − γ [(µ − 1 + ν)/2]ρ − pl e µ−1,gν e µ,gν+1 e µ−1,gν ∂t  ( ) √ √ (n) ′ √ (n) n) (n, n ′) − ν + 1ρ + ivn ( ν ρe µ−1,eν−1 − µρgµ,gν ) − i vn′ µ − 1ρ(n, ee µ−2,ggν ee µ−1,ggν+1 n′

−i



vn′

(√

(n, n ′)

µρegµ,geν −

√

(n, n ′)

) ν ρegµ−1,geν−1 ,

(A3)

n′

( )  ∂ (n) (n) (n) ρgµ,gν = −(2δ n + iω µν )ρ(n) gµ,gν − γpl [(µ + ν)/2]ρgµ,gν − (µ + 1)(ν + 1)ρgµ+1,gν+1 ∂t  ( ) √ √ (n, n′) √ (n) n ′) + ivn ( ν ρ(n) µρe µ−1,gν ) + k n ρ µ,ν − i vn′ µ + 1ρ(n, gµ,eν−1 − ggµ+1,geν − ν ρggµ,geν−1 n′

−i



vn′

(√

(n, n ′)

µρgeµ−1,ggν −

√

(n, n ′)

)

ν + 1ρgeµ,ggν+1 ,

(A4)

n′

and ∂ (n) √ (n) ρ = −(2δ n + iω µ−1ν−1)ρ(n) µν ρ(n) e µ,eν ) e µ−1,eν − γpl([(µ + ν − 2)/2]ρe µ−1,eν−1 − ∂t e µ−1,eν−1  ( ) √ √ √ (n) √ (n, n′) (n, n ′) ′ + ivn ( ν ρ(n) − µρ ) + p ρ − i v µρ − ν − 1ρ n µ−1,ν−1 n e µ−1,gν gµ,eν−1 egµ,eeν−1 egµ−1,eeν−2 n′

−i



vn′

(

(n, n ′)

µ − 1ρee µ−2,egν−1 −

√

(n, n ′)

)

ν ρee µ−1,egν .

(A5)

n′

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-11

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

Restricting to the lowest order with respect to the moleculeMNP′ coupling vn , we neglect density matrices of type n) ρ(n, , Eq. (A2). Then, the remaining equations are ab µ,cdν closed and can be solved numerically. However, we would like to go one step further and will demonstrate that an analytical solution for the ρ(n) can be constructed. This a µ,bν is only possible if we remove the density matrices which correspond to the next higher level of plasmon excitations. For example, considering the equation of motion for ρ(n) gµ,eν−1, there also appears ρ(n) . Respective terms describing gµ+1,eν plasmon decay from the next higher states are also present in the remaining equations. To avoid overestimation of plasmon decay by simply neglecting these terms which increase the plasmon population of the actual levels, we use the following approximation in Eq. (A1): ) (  (n) − γpl [(µ + ν − 1)/2]ρ(n) gµ,eν−1 − (µ + 1)ν ρgµ+1,eν ) (  (A6) ≈ −γpl [(µ + ν − 1)/2] − µ(ν − 1) ρ(n) gµ,eν−1. Such an approximation is also applied with respect to the other equations. Together with a neglect of higher order density matrices, it results in ∂ (n) ρ = i(ω˜ ∗n + ω˜ µν−1)ρ(n) gµ,eν−1 ∂t gµ,eν−1 √ (n) √ + ivn ( ν ρgµ,gν − µρ(n) e µ−1,eν−1), (A7) ∂ (n) ρ = −i(ω˜ n + ω˜ µ−1ν )ρ(n) e µ−1,gν ∂t e µ−1,gν √ (n) √ + ivn ( ν ρe µ−1,eν−1 − µρ(n) gµ,gν ), (A8) ∂ (n) ρgµ,gν = −(2δ n + iω˜ µν )ρ(n) gµ,gν + k n ρ µ,ν ∂t √ √ (n) µρe µ−1,gν ), + ivn ( ν ρ(n) gµ,eν−1 −

(A9)

∂ (n) ρ = −(2δ n + iω˜ µ−1ν−1)ρ(n) e µ−1,eν ∂t e µ−1,eν−1 √ (n) √ + ivn ( ν ρe µ−1,gν − µρ(n) gµ,eν−1) + pn ρ µ−1,ν−1.

(A10) √ Be aware of the definition ω˜ µν = ω µν − i((µ + ν)/2 − µν). Solving the equations of motion for the ρ(n) introduces a µ, bν some memory effects into the equation for the original density matrix ρ µν . We will assume that such memory effects are of minor importance and will neglect them. This is equivalent to a solution of the equation of motion for the ρ(n) by ignoring a µ,bν the time-derivative. Thus, it is assumed that the respective decay rates dominate the time evolution (approximation of instantaneous dephasing; for further justification, see also the discussion of the numerical results). Consequently, we arrive at the respective solution of the Eqs. (A7)–(A10). Here, we only quote those results which are necessary to determine ρ µν √ (n)√ a(n) µν µ + c µν ν (n) (n) ρgµ,eν−1 = d µ−1ν−1 pn ρ µ−1,ν−1 λ (n) µν √ (n)√ a(n) µν ν + c µν µ (n) − d µν k n ρ µ,ν (A11) λ (n) µν and

√ (n)√ b(n) µν ν − c µν µ

ρ(n) e µ−1,gν =

λ (n) µν −

d (n) µ−1ν−1 pn ρ µ−1,ν−1

√ (n)√ b(n) µν µ − c µν ν λ (n) µ,ν

d (n) µν k n ρ µ,ν .

(A12)

In the above expressions, we introduced (n) a(n) ˜ n + ω˜ µ−1ν − i(d (n) µν = ω µν ν + d µ−1ν−1 µ)vn ,

(A13)

(n) b(n) ˜ ∗n − ω˜ µν−1 + i(d (n) µν = ω µν ν + d µ−1ν−1 µ)vn , √ (n) (n) c(n) µν = i(d µν + d µ−1ν−1) µνvn ,

(A14)

d (n) ˜ µν + k n + pn ) µν = vn /(i ω

(A16)

(n) (n) (n) 2 λ (n) µν = a µν b µν + (c µν ) .

(A17)

(A15)

and Inserting the above results into Eq. (23), we arrive at ∂ ρ µν = −iω µν ρ µν ∂t  γpl − (µ + ν)ρ µν + γpl (µ + 1)(ν + 1)ρ µ+1ν+1 2  (n) (n) (n) − (h(n) µ+1ν+1 d µν pn + g µν d µν k n )ρ µν −

n 

(n) f µ+1ν+1 d (n) µ+1ν+1 k n ρ µ+1ν+1

n

−



(n) (n) f µν d µ−1ν−1 pn ρ µ−1ν−1,

(A18)

n

including the additional abbreviations (n) (n) √ (b(n) µν − a µν ) µν − c µν (µ + ν) (n) f µν = ivn λ (n) µν

(A19)

and (n) g µν

= ivn

(n) (n)√ a(n) µν ν − b µν µ + 2c µν µν

λ (n) µν

.

(A20)

The resulting equations for the diagonal matrix elements Pµ = ρ µ µ , i.e., the plasmon-state populations, follow as: ∂ Pµ = −(µγpl + k µ )Pµ + ((µ + 1)γpl + k µ+1)Pµ+1 ∂t − pµ+1 Pµ + pµ Pµ−1. (A21) To derive these equations, we have noticed that d (n) µ µ = vn /(k n + pn ). In addition, we have introduced the following rates:   µκ n µ k n kµ = g µ(n)µ d (n) (A22) µ µ kn = p + k n + 2µκ n µ n n n and pµ = −



f µ(n)µ d (n) µ µ pn =

n

 n

µκ n µ pn . pn + k n + 2µκ n µ

(A23)

k µ acts as a plasmon decay rate (beside the term µγpl in Eq. (A21)) and plasmon pumping is realized via pµ . Both types of rates are of composite character, since the molecular rates k n and pn contribute together with the energy transfer rate κ n µ ≡ 2Im

vn2 . ω˜ ∗n − ω˜ µ µ−1

(A24)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-12

Y. Zhang and V. May

J. Chem. Phys. 142, 224702 (2015)

This rate describes energy exchange between the molecule and the plasmon level µ and µ − 1. Accordingly, plasmon pumping via pµ is constituted by the product of the molecular pump rate pn and the energy transfer rate µκ n µ into the µ’th excited plasmon level. In order to get the correct dimension of a rate, the product is divided by the total rate pn + µκ n µ + k n + µκ n µ affecting the molecular states. Since different molecules couple to the MNP, we finally have to sum up with respect to all molecules. The same interpretation is possible with respect to the plasmon decay rate k µ . It describes plasmon decay via energy transfer to the molecules. Considering steady state conditions Eq. (A21) further simplifies and can be easily translated into the balance relation given in Eq. (25).

APPENDIX B: RATE EQUATIONS IN MEAN-FIELD APPROXIMATION

So far, the dynamics of the molecule-MNP system has been studied by means of the exact density matrix approach using Eq. (17) as well as in utilizing the treatment explained in Appendix A. While the latter consideration has been based on the computation of the reduced plasmon density matrix, the subsequent discussion will exclusively use diagonal density matrix elements, i.e., state populations. We first derive an equation of motion for the molecular excited state population, + + Pne (t) = tr{ ρ(t)B ˆ n Bn } ≡ ⟨Bn Bn ⟩.

(B1)

Using basic density operator equation (7), one immediately arrives at ∂ Pne = pn − (k n + pn )Pne − 2vn Im⟨BnC +⟩. ∂t

(B2)

The equation of motion for the mixed molecule-plasmon function takes the form ∂ ∗ ⟨BnC +⟩ = −i(ω˜ n − ω˜ pl )⟨BnC +⟩ ∂t
+ ivn Pne + ⟨(Bn+ Bn − Bn Bn+)C +C⟩  +i vn′⟨Bn+′ Bn ⟩. (B3) n,n ′

We introduced ω˜ n = ω n − i(k n + pn )/2 and ω˜ pl = ωpl − iγpl/2. In a next step, we ignore the last term describing higher intermolecular correlations and solve the equation in the limit of instantaneous dephasing by neglecting the time-derivative (this approximation implies that the memory effect appearing when solving the equation exactly is of less importance; see Appendix A and the recent application18). Further on, we introduce the factorization ⟨(Bn+ Bn − Bn Bn+)C +C⟩ ≈ (Pne − Png )Npl, where Npl = ⟨C +C⟩

(B4)

is the mean plasmon number. It results in ⟨BnC +⟩ =


ivn ∗ Pne + (Pne − Png )Npl . ω˜ n − ω˜ pl

(B5)

If the above expression is inserted in Eq. (B2), we arrive at the following nonlinear rate equation:

∂ Pne = pn − (k n + pn )Pne ∂t
− κ n Pne (1 + Npl) − Png Npl .

(B6)

Note the introduction of the energy transfer rate κ n = 2vn2 γn/((ωpl − ω n )2 + γn2 ), where γn = (γpl + pn + k n )/2. The combination of the this rate with Pne describes the spontaneous and plasmon induced energy transfer from molecule n to the MNP. The subsequent term proportional to Png characterizes the reverse energy transfer from the MNP to the molecule (be aware of Png = 1 − Pne ). Following the same line of approximations which led us to the above given rate equation, we get 
∂ Npl = −γplNpl + κ n Pne (1 + Npl) − Png Npl . (B7) ∂t n We obtained the same structure of the molecule-MNP energy transfer as in the rate equation for Pne . Noting also the alternative form Pne + (Pne − Png )Npl, one realizes that a positive molecular population inversion amplifies the mean plasmon number. Additionally, rate equation (B7) also indicates the possibility to increase the mean plasmon number by an increase of the number of molecules. If all molecules are identical we may simply combine the number of molecules Nmol with the energy transfer rate to get Nmol κ. If this is much larger than the plasmon damping, the molecules can increase the plasmon number. In the reverse case the molecules cannot. 1M. L. Brongersma, N. J. Halas, and P. Nordlander, Nat. Nanotechnol. 10, 25

(2015). Ringe, B. Sharma, A.-I. Henry, L. D. Marks, and R. P. Van Duyne, Phys. Chem. Chem. Phys. 15, 4110 (2013). 3M. S. Tame, K. R. McEnery, S. K. Özdemir, J. Lee, S. A. Maier, and M. S. Kim, Nat. Phys. 9, 329 (2013). 4Y. Yin, T. Qiu, J. Li, and P. K. Chu, Nano Energy 1, 25 (2012). 5P. Berini and I. D. Leon, Nat. Photonics 6, 16 (2012). 6M. A. Noginov, G. Zhu, A. M. Belgrave, R. Bakker, V. M. Shalaev, E. E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, Nature 460, 1110 (2009). 7X. G. Meng, A. V. Kildishev, K. Fujita, K. Tanaka, and V. M. Shalaev, Nano. Lett. 9, 4106 (2013). 8M. T. Hill, M. Marell, E. S. P. Leong, B. Smalbrugee, Y. Zhu, M. Sun, P. J. van Veldhoven, E. J. Geluk, F. Karouta, Y. S. Oei, R. Nötzel, C. Z. Ning, and M. K. Smit, Opt. Express 17, 11107 (2009). 9R. M. Ma, R. F. Oulton, V. J. Sorger, G. Bartal, and X. Zhang, Nat. Mater. 10, 110 (2011). 10Y. J. Lu, J. Kim, H. Y. Chen et al., Science 337, 450 (2012). 11K. Ding, Z. C. Liu, L. J. Yin, M. T. Hill, M. J. H. Marell, P. J. van Veldhoven, R. Noetzel, and C. Z. Ning, Phys. Rev. B 85, 041301(R) (2012). 12Y. Zhang and V. May, Phys. Rev. B 89, 245441 (2014). 13M. I. Stockman, J. Opt. 12, 02404 (2010). 14P. Meystre and M. Sargent, Elements of Quantum Optics (Springer-Verlag, Berlin, Heidelberg, 1990). 15M. Richter, M. Gegg, T. S. Theuerholz, and A. Knorr, Phys. Rev. B 91, 035306 (2015). 16U. Martini and A. Schenzle, Lect. Notes Phys. 561, 238 (2001). 17H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, 2002). 18V. May, J. Chem. Phys. 140, 054103 (2014). 19V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley-VCH, Weinheim, 2000, 2004, 2011). 20Y. Zelinskyy, Y. Zhang, and V. May, J. Phys. Chem. A 116, 11330 (2012). 21L. G. Gerchikov, C. Guet, and A. N. Ipatov, Phys. Rev. A 66, 053202 (2002). 22G. Weick, R. A. Molina, D. Weinmann, and R. A. Jalabert, Phys. Rev. B 72, 115410 (2005). 23G. Weick, G.-L. Ingold, R. A. Jalabert, and D. Weinmann, Phys. Rev. B 74, 165421 (2006). 24Y. Zhang, Y. Zelinskyy, and V. May, J. Nanophotonics 6, 063533 (2012). 2E.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

224702-13 25The

Y. Zhang and V. May

well-known relation ∆E∆t ≥ ~ relates in the present case the energy uncertainty of the emitted photons, i.e., their linewidth, to the timeuncertainty, i.e., the lifetime of the system. Since an increasing plasmon excitation results in an increase of the lifetime and in an increase of the photon emission intensity, the latter quantity behaves reversely to the linewidth. It decreases with increasing photon emission.

J. Chem. Phys. 142, 224702 (2015) 26M.

Sargent II, M. O. Scully, and W. E. Lamb, Laser Physics (AddisonWesley Publishing Company, 1974). 27W. Vogel and D. G. Welsch, Lectures on Quantum Optics (Akademie Verlag, 1994). 28T. S. Theuerholz, A. Carmele, M. Richter, and A. Knorr, Phys. Rev. B 87, 245313 (2013).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 61.16.135.116 On: Mon, 15 Jun 2015 06:39:43

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.