Interaction of single molecules with metallic nanoparticles

0.64 MB·English

by  Ulrich Hohenester

Interaction of single molecules with metallic nanoparticles
Ulrich Hohenester∗ and Andreas Tru¨gler
Institut fu¨r Physik, Karl–Franzens–Universita¨t Graz, Universita¨tsplatz 5, 8010 Graz, Austria
(Dated: February 3, 2008)
Wetheoreticallyinvestigatetheinteractionbetweenasinglemoleculeandametallicnanoparticle.
We develop a general quantum mechanical description for the calculation of the enhancement of
radiative and non-radiative decay channels for a molecule situated in the nearfield regime of the
metallic nanoparticle. Using a boundary element method approach, we compute the scattering
rates for several nanoparticle shapes. We also introduce an eigenmode expansion and quantization
scheme for the surface plasmons, which allows us to analyze the scattering processes in simple
8 physicalterms. Anintuitiveexplanationisgivenforthelargequantumyieldofquasione-andtwo-
0 dimensionalnanostructures. Finally,webrieflydiscussresonantFo¨rsterenergytransferinpresence
0 of metallic nanoparticles.
2
PACSnumbers: 42.50.Nn,42.50.Ct,73.20.Mf,33.50.-j
n
a
J
I. INTRODUCTION establish, in linear response, a convenient link between
5
quantum electrodynamics and classical Maxwell theory.
2
Metal nanoparticles can sustain local surface plasmon In this work we start by reviewing the description

] excitations,particleplasmons,whicharehybridmodesof scheme for a molecular dipole placed in the vicinity of
l
l a light field coupled to a coherent electron charge oscil- a metallic nanoparticle. We show how to use the dyadic
a
Green function of Maxwell’s theory to express the elec-
h lation [1, 2]. The properties of these excitations depend
tricfieldfluctuationsintermsofthecurrentfluctuations
- strongly on particle geometry and interparticle coupling,
s in the metal, and how to derive the molecule scatter-
e and give rise to a variety of effects, such as frequency-
ing rate within linear response theory. In the quasistatic
m dependent absorption and scattering or near field en-
regime, where the light wavelength is much larger than
hancement. Particle plasmons enable the concentration
.
at oflightfieldstonanoscalevolumesandplayakeyrolein the nanoparticle size and the molecule-MNP distance,
thedecayratecanbeseparatedintoaradiativeandnon-
m surfaceenhancedspectroscopy[3]. Thiscanbeexploited
radiative contribution, which both show a simple scal-
for a variety of applications, such as surface-enhanced
-
ing behavior. We next use a boundary element method
d Raman scattering [4, 5] or biochemical sensorics [6].
n approach [26] to compute these decay rates for realistic
Improvednanofabricationmethodsnowadaysallowad-
o nanoparticle geometries. The results are shown to be
vancedcontrolofthenanoparticleshapeandthearrange-
c in nice agreement with Mie theory. We additionally in-
[ ment of nanoparticle ensembles [7], and open the possi-
troduce for a Drude description of the metal electrons
bilitytoflexiblytailorspecificmolecule-nanoparticlecou-
1 a quantization scheme for the surface plasmons [27, 28].
plings. Anger et al. [8] recently investigated the fluores-
v Although this approach could be used to study the non-
cence of a single molecule and a single laser-irradiated
0 linear optical response of metallic nanoparticles, in this
spherical gold nanoparticle, and demonstrated by vary-
0 work we use this framework to provide an intuitive pic-
9 ing the molecule-nanoparticle distance the continuous
ture for the radiative and non-radiative molecular de-
3 transitionfromfluorescenceenhancementtofluorescence
cay channels only, and give a simple explanation for the
. quenching [9]. This finding was supported by Ku¨hn and
1 higher quantum yield achievable in quasi one- and two-
0 coworkers [10] using a similar setup, and stimulated a
dimensionalnanoparticles. Finally,webrieflydiscussres-
8 large number of related experiments [11, 12].
onant F¨orster energy transfer between donor and accep-
0 The theoretical description of the modified molecule
tor molecules placed in the nanoparticle nearfield.
:
v dynamics in presence of metallic nanoparticles (MNPs)
i has been addressed by several authors in the past (see,
X
e.g. [4,13,14,15]andreferencestherein). Morerecently,
r thequestionofhowtotreatquantumelectrodynamicsin II. THEORY
a
the presence of absorbing media has seen anewed inter-
est[16,17,18,19],motivatedbyrelatedworkonthermal As for the description of the molecule, in this work
nearfieldradiation[20,21]andthermallyinducedscatter- we follow Refs. [14, 15] who considered a generic two-
ing processes of atoms located in the vicinity of metallic level system. This approach is also best suited for other
or superconducting bodies [22, 23, 24, 25]. This work quantum emitters, such as colloidal quantum dots. We
has identified the Green tensor as the central element to describe the molecule in terms of a generic two level
scheme,with0thegroundstateand1theexcitedmolecu-
larstatewithenergyE . Other,morerefineddescription
1
schemes,considering,e.g.,thedifferentvibronicstatesof
∗Electronicaddress: [email protected] the molecule [29, 30] or details of the molecular orbitals
2
[31], could be treated in a similar fashion. wherewehavemadeuseoftheEinsteinsumconvention.1
We first consider the optical decay from the excited SomedetailsofitsderivationareprovidedinappendixA.
state 1 to the groundstate. Associated to this transition Theexpressioninbracketsdescribesthefieldfluctuations
is a dipole operator in thermal equilibrium, to be determined at the molecu-
lar transition frequency ω and position r . These fluc-
m
dˆ=µ(σˆ +σˆ )≡dˆ++dˆ−, (1) tuations are induced by current noise in the metal. Let
− +
jˆ(r)=jˆ+(r)+jˆ−(r) denote the current operator in the
where µ is the molecular dipole moment, assumed to metal, which we have again split into positive and nega-
be real, σˆ = |1(cid:105)(cid:104)0| is the molecular excitation operator tive frequency components. If j were a classical current,
+
whichbringsthemoleculefromthegroundtotheexcited evolving with frequency e−iωt, the electric field could be
state, and σˆ = |0(cid:105)(cid:104)1| the corresponding de-excitation determined from the wave equation [35]
−
operator. In an interaction representation according to
the free Hamiltonian Hˆ0 = E1|1(cid:105)(cid:104)1| the operators σˆ− ∇×∇×E(r,ω)−k2(cid:15)(r,ω)E(r,ω)= 4πiω j(r,ω), (4)
and σˆ oscillate with positive and negative frequencies c2
+
e−iωt and eiωt, respectively, as indicated by the dipole
de-excitation and excitation operators dˆ± in Eq. (1). It where E(r,ω) and j(r,ω) are the Fourier transform of
the electric field and current, respectively, c is the vac-
isclearthatthelatteroperatorscanbeeasilygeneralized
uum speedof light, k =ω/c the wavenumberin vacuum,
in case of more complicated level schemes.
and (cid:15)(r,ω) the space and frequency dependent dielec-
In a similar fashion, we split the electric field operator
tric function. We set the magnetic permeability µ = 1
Eˆ(r) = Eˆ+(r) + Eˆ−(r) into two contributions evolv-
throughout. A convenient way to compute the electric
ing with positive and negative frequencies [32, 33]. The
field for a given current source is by means of the dyadic
molecule-lightcouplingisdescribedwithinthedipoleap-
Green tensor G(r,r(cid:48),ω) [2], to be determined from the
proximation
wave equation
(cid:16) (cid:17)
Hˆop =−dˆ·Eˆ ∼=− dˆ+·Eˆ−+dˆ−·Eˆ+ , (2) ∇×∇×G(r,r(cid:48)ω)−k2(cid:15)(r,ω)G(r,r(cid:48),ω)=4πδ(r−r(cid:48))11,
(5)
where Eˆ(r ) has to be taken at the position r of the togetherwithappropriateboundaryconditions. Thecal-
molecule. Wmehaveneglectedcontributionswhermebothdˆ culation of G(r,r(cid:48),ω) is a common problem in classical
and Eˆ oscillate with either positive or negative frequen- electrodynamics,andwillbedescribedinmoredetailfur-
therbelow. Forthemomentitsufficestoassumethatthe
cies, which corresponds to the common rotating-wave
Greenfunctionisatourhand. BycomparisonofEqs.(4)
approximation [32, 33]. In quantum optics one usually
and (5) we readily observe
expands the field operators Eˆ± in the eigenmodes of
Maxwell’s equation to obtain photon annihilation and iω (cid:90)
creation operators. The first expression in parentheses E(r,ω)= G(r,r(cid:48)ω)j(r(cid:48),ω)dr(cid:48). (6)
c2
thendescribesamolecularde-excitationandthecreation
of a photon, and the second expression the reversed pro- Thus, for a linear material response, described by the
cess. dielectric function (cid:15)(r,ω), the Green function provides
Things change considerably if the field operators Eˆ thelinkbetweenthecurrentsourceandtheelectricfield.
evolveinacomplexdielectricenvironmentwhereabsorp- Eq. (6) can be also used for the corresponding electric-
tion takes place, as is the case for MNPs. Owing to ab- field and current operators [16], where Eˆ±(r,ω) is in-
sorption the eigenmodes acquire a complex energy, the duced by jˆ±(r(cid:48),ω). Note that in the linear response
imaginarypartassociatedwiththedampingoftheeigen-
regime positive and negative frequency components do
modes,whichspoilstheusualquantizationprocedurefor
not mix. We can now use Eq. (6) to relate the field fluc-
theelectromagneticfields. Aswewillshownext,inlinear
tuations through
response one can make use of the fluctuation-dissipation
theorem [34]toavoidsuchaneigenmodeexpansionofEˆ. (cid:68) (cid:69) ω2 (cid:90)
Eˆ+(r ,ω)Eˆ−(r ,ω) = G (r ,r,ω)
i m j m c4 ik m
(cid:68) (cid:69)
A. Molecule–MNP coupling × ˆjk+(r,ω)ˆjl−(r(cid:48),ω) G∗jl(rm,r(cid:48),ω)drdr(cid:48) (7)
to the current fluctuations in the dielectric. Here, the
1. Green function approach
Green functions describe how the field propagates from
the current source to the position of the molecule. In
OurstartingpointisprovidedbytheFermi-golden-rule
expression for the molecule decay rate
(cid:68) (cid:69)
γ =2πµi Eˆi+(rm,ω)Eˆj−(rm,ω) µj (3) 1 WeuseGaussandatomicunitse=m=(cid:126)=1throughout.
3
appendix B we show that for local and isotropic di- molecular dipole, as computed in the quasistatic limit.
electricsthecurrentcorrelationfunctioncanbeexpressed Wehavetoadditionallyaccountforradiativelosses. This
as δ δ(r − r(cid:48))ω2(cid:15)(cid:48)(cid:48)(r,ω)/(2π), with (cid:15)(cid:48)(cid:48)(ω) being the is done by using the standard expression
kl
imaginary part of the dielectric function. Together with
the integral equation [16] |µ+µ |2
γ =γ MNP (12)
(cid:90) r 0 µ2
G (r,s,ω)k2(cid:15)(cid:48)(cid:48)(s,ω)G∗ (r(cid:48),s,ω)ds
ik jk
=4π(cid:61)m[G (r,r(cid:48),ω)] , (8) for the composite molecule-MNP dipole radiator. Here
ij
γ = 4n µ2k3isthefree-spacedecayrateofthemolecule,
0 3 b
which is obtained by multiplying Eq. (5) with G∗ and with n being the refractive index of the embedding
jk b
subtracting the equation obtained by complex conjuga- medium. The second term on the right-hand side of
tion, we then arrive at the expression for the molecular Eq. (12) is the enhancement of the radiative molecular
decay rate decay in presence of the MNP, where µ is the MNP
MNP
dipole induced by the molecule. Below we will show how
γ =2k2µ· (cid:61)m[G(rm,rm,ω)] ·µ. (9) to compute Eqs(rm) and µMNP within a boundary ele-
ment method approach.
Equation (9) is our final result. It shows that the decay
rate of a molecule in proximity to a MNP is fully deter-
mined by the dyadic Green function of classical Maxwell
theory. This is a huge simplification because all the dy-
3. Scaling
namicsofcurrentfluctuationsinthemetalareembodied
in the dielectric function, which can be obtained from
either experiment or first-principles calculation. In the general case there are two length scales which
Eq. (9) can be brought into an even more transparent define the problem, namely the size of the nanoparticle
formbyusingtherelationj =−iωµbetweenthecurrent andthewavelengthofthelight. Inthequasistaticregime,
and the molecular dipole. Inserting this expression into where the light wavelength is assumed to be much larger
Eq. (6), one immediately obtains for Eq. (9) the form thantheMNPsize,itisasinglelengthscaleenteringthe
problem. Inthisregimeonecanfindasimplescalinglaw
γ
=−(cid:61)m[µ·E(r )] , (10) for γ and γ on the size of the MNP. Suppose that we
2 m nr r
have computed the two decay rates for a given nanopar-
whereE(r )istheelectricfieldinducedbythedipoleµ ticlegeometry. Wecanthenscalethewholegeometryby
m
itself. Equation (10) describes a self-interaction, where ascalingfactorλ, accordingtor →λr, whilekeepingall
the dipole polarizes the MNP, and the total electric field otherquantities,suchastransitionfrequencyωormolec-
of the dipole and of the polarized MNP acts, in turn, ular dipole moment µ, fixed. One can show that such a
back on the dipole. It is the imaginary part of this self scaling transforms the scattering rates according to
interaction which accounts for the decay of the excited
molecular state. γ (λ)=γ /λ3, γ (λ)=γ . (13)
nr nr r r
We will find that Eq. (13) is an extremely useful ex-
2. Quasistatic approximation
pression for estimating the relative importance of non-
radiative and radiative decay rates for realistic MNPs.
In many cases of interest the size of the MNP is much
smaller than the light wavelength. It is then possible
to employ the so-called quasistatic approximation. For
a given source, the electric field is computed from the
B. Mie theory
scalarpotentialΦ(r)accordingtoE(r)=−∇Φ(r). This
approach differs from the truly static case in that one
uses the frequency-dependent dielectric function (cid:15)(r,ω), Theelectromagneticresponseofsphericalmetallicpar-
at the molecular transition frequency ω, rather than the ticles can be solved exactly within Mie theory. This the-
static limit. The imaginary part of (cid:15)(r,ω) is associated ory provides closed expressions in both the retarded and
with Ohmic losses of electromagnetic fields inside the quasistatic case. Indeed, it is a very rare situation that
metal. The self-interaction expression in the quasistatic onecanfindanexactsolutionforanexperimentallyinter-
limit thus accounts only for nonradiative losses esting, nontrivial problem. In this work we will be using
the Mie theory primarily as a test case for our numerical
γ
2nr =−(cid:61)m[µ·Eqs(rm)] , (11) scheme described below. As there exists a vast amount
of literature on Mie theory [36], in appendix C we only
as indicated by the subscript on γ. E (r ) is the elec- briefly sketch the quasistatic approach and provide the
qs m
tric field at the position of the molecule, induced by the explicit expressions used in our calculations.
4
C. Boundary element method (BEM) where the different signs depend on whether the surface
derivative is taken inside or outside the nanoparticle.
In the quasistatic limit the wave equation (4) for the Atthispointitisconvenienttochangefromboundary
electric field transforms to the Poisson equation integrals to boundary elements, suitable for a compu-
tational implementation. Within the boundary element
∇2Φ(r)=−4πρ(r) (14) method (BEM)approach,thesurfaceofthenanoparticle
is approximated by small surface elements of triangular
for the scalar potential. Here ρ(r) is the charge distri-
shape, as shown in Fig. 1 for the example of a spheri-
bution of the molecular dipole. The problem can be sig-
cal particle. In its most simple form, which we will be
nificantly simplified for a MNP described by a homoge-
using in this work, the surface charges are assumed to
neous dielectric function (cid:15) (ω) embedded in a medium
m be located at the center of each triangle s , and the ma-
i
with dielectric constant (cid:15) . As long as thereis no danger
b trices G and F connecting different surface elements
ij ij
of confusion, we will suppress the frequency dependence are obtained from G(s ,s ) and F(s ,s ).3 The integral
i j i j
of the metal dielectric function and assume that (cid:15) is
m equation (17) then reduces to two matrix equations
evaluated at the molecular transition frequency. We can
nowusethestaticGreenfunctionG(r,r(cid:48))=1/(|r−r(cid:48)|), Φ(cid:48) =(F +2π11)σ+Φ(cid:48) , Φ(cid:48) =(F −2π11)σ+Φ(cid:48) ,
m ext b ext
which is the solution of the Poisson equation (18)
where F is the matrix with elements F(s ,s ) and the
∇2G(r,r(cid:48))=−4πδ(r−r(cid:48)) (15) subscript of Φ(cid:48) indicates the side on whicih tjhe surface
derivative is taken. We can now use the boundary con-
for a point source located at position r(cid:48). The solution of
dition (cid:15) Φ(cid:48) −(cid:15) Φ =0 to compute, for a given external
Eq. (14) can then be expressed as [26, 37] m m b b
charge distribution, the surface charges according to
(cid:90)
Φ(r)= G(r,s(cid:48))σ(s(cid:48))ds(cid:48)+Φext(r), (16) σ =−(cid:8)2π((cid:15) +(cid:15) )11+((cid:15) −(cid:15) )F(cid:9)−1((cid:15) −(cid:15) )Φ(cid:48) . (19)
∂Ω m b m b m b ext
where ∂Ω denotes the surface of the metallic nanopar- Equation (19) is the main result for our BEM ap-
ticle, σ is a surface charge that has to be chosen such proach. It shows that the surface charge can be com-
that the appropriate boundary conditions are fulfilled, putedthroughnumericalinversionofamatrixhavingthe
and Φ is the scalar potential for the external charge size of the number of triangular surface elements, which
ext
distributionofthemoleculardipole. Indeed,Eq.(16)ful- is typically of the order of thousands. Once we have
fills the Poisson equation (14) everywhere except at the computedσ foragiveninhomogeneityΦ(cid:48) ,wecancom-
ext
boundaries. Theboundaryconditionsat∂Ωarethecon- pute the potential and electric fields at any given space
tinuityofthetangentialelectricfieldandnˆ·(D −D )= pointfromEq.(16),andtheinduceddipolemomentfrom
b m (cid:80)
0, where nˆ is the unit vector normal to the surface and µMNP = isiσi. In our computational approach we de-
pointing away from the MNP. D is the dielectric dis- scribe the molecular dipole by two oppositely charged
b
placementattheoutersurface,andD thedisplacement particlesseparatedbyasmalldistance,andcomputethe
m
at the inner surface. The continuity of the electric field electricfieldfromEq.(16)bymeansofafinitedifference
is guaranteed for a potential which is continuous at the scheme.
boundary, as is the case in Eq. (16). To implement the
second boundary condition, we have to take the surface
derivative of Φ at either side of the boundary. The limit D. Eigenmode expansion
r → s in Eq. (16) has to be treated with some care.2
With the abbreviations F(s,s(cid:48)) = (nˆ · ∇)G(s,s(cid:48)) and 1. Boundary element method approach
Φ(cid:48) =(nˆ ·∇)Φ, we then get
Things can be formulated differently when the dielec-
(cid:90)
lim Φ(cid:48)(r)= F(s,s(cid:48))σ(s(cid:48))ds(cid:48)±2πσ(s)+Φ(cid:48) (s), tric function is approximately of Drude form
ext
r→s
∂Ω
(17) ω2
(cid:15) (ω)=(cid:15) − p . (20)
m 0 ω(ω+iγ)
2 Letusconsiderlimr→snˆ·∇R G(r,s(cid:48))σ(s(cid:48))ds(cid:48) foracoordinate
system where nˆ = eˆz, r = (0,0,z), and s(cid:48) = ρ(cosφ,sinφ,0) 3 ForthediagonalelementsofGwechangetopolarcoordinatesρ
is given in polar coordinates ρ and φ. We compute the bound-
andφ,andperformtheintegrationwithinatrianglewherer(φ)
ary integral within a small circle with radius R, within which
denotestheupperlimitofρforagivenφ. Wethenobtain
the surface charge σ can be approximated by a constant. The
integralthenbecomes Z 2π Z r(φ) 1 Z 2π
dφ ρdρ = r(φ)dφ,
limnˆ·Z r−s(cid:48) ds(cid:48)→ lim2πzZ Rρdρ(ρ2+z2)−23 =2π. 0 0 ρ 0
z→0 |r−s(cid:48)|3 z→0 0 which can be easily calculated numerically. The diagonal ele-
mentsofF areallzero.
5
Here (cid:15) is a dielectric constant, accounting for screening β is a real number. We can now expand the surface
0 λ
of d-band electrons in transition metals, ω is the bulk charge distribution in terms of these eigenfunctions viz.
p
plasma frequency, and γ the Landau damping of plas- [27, 28, 38]
mons. In gold the above Drude form is valid for energies
bd-ebloawndtheexctihtarteisohnosldsevtaliune. oTfyappicparloxviamluaetselfyor2AeuV, wwhheicrhe σ =(cid:88)(cid:18)ω2nβ0 (cid:19)12 uλ(r)(cid:0)aλe−iωλt+a∗λeiωλt(cid:1) , (24)
λ λ
we will be using in this work, are (cid:15) = 10, ω = 9 eV, λ
0 p
and 1/γ =10 fs. Instead of the Drude form we can also
wherea areexpansioncoefficientsfortheplasmonmode
λ
use a hydrodynamic model where electrons with particle
λ. With Eq. (24) the Hamilton function can be brought
densityn0 movefreelyinamediumwiththebackground to the intriguing form H = 1(cid:80) ω (a∗a +a a∗). The
dielectric constant (cid:15) . This model is completely equiv- 2 λ λ λ λ λ λ
0 quantization may now be performed through the substi-
alent to the framework of the Drude dielectric function. tution a , a∗ → aˆ , aˆ†, where aˆ and aˆ† are the oper-
Aswewillshowbelow,throughthehydrodynamicmodel λ λ λ λ λ λ
ators for annihilation and creation of a surface plasmon
we can establish a microscopic description for the elec-
mode λ, respectively. They obey the usual bosonic com-
tron dynamics, which will prove useful for interpreting
mutation relations. With these quantized surface plas-
our results in physical terms and for performing a quan-
mon modes Eq. (23) can be brought to the form
tization of surface plasmons.
The energy of a classical electron plasma is the sum of Hˆ =(cid:88)ω (cid:16)aˆ†aˆ + 1(cid:17) . (25)
kinetic and electrostatic energy [27, 28] pl λ λ λ 2
λ
(cid:90)
H = 1 (cid:8)n (∇Ψ)2+ρΦ(cid:9)d3r. (21) This is the Hamilton operator for surface plasmons. In
2 0
principle it can be used to describe the non-linear prop-
ertiesofmetallicnanoparticles, althoughinthisworkwe
Hereρ(r)isthechargedensitydisplacementfromequilib-
shall only be interested in linear response. What still re-
rium, Φ(r)istheelectrostaticpotentialinducedbyρ(r),
mains to be done is to show how Landau damping can
andΨ(r)isthevelocitypotential,whosederivativegives
be included, and how to obtain the nonradiative and ra-
the velocity density v =−∇Ψ [28]. For the surface plas-
diative scattering rates starting from Eq. (25).
mons of our present concern we consider surface charge
Let us first consider the coupling between the molec-
distributions σ which are nonzero only at the surface of
ular dipole and the surface plasmons (24). Using the
theMNP.AsdetailedinappendixD,theHamiltonfunc-
expressions given in appendix D, we can compute for a
tion(21)canberewritteninaboundaryelementmethod
given surface mode u the scalar potential Φ and its
approach as λ λ
surface derivative Φ(cid:48) according to
λ
1 (cid:110) (cid:16) (cid:17)−1
H = σ˙T 2π11+F˜ Gσ˙ (cid:104) (cid:105)−1
2n0 Φλ = −4π 2π((cid:15)0+(cid:15)b)11+((cid:15)0−(cid:15)b)F˜ Guλ
(cid:104) (cid:105)−1 (cid:111)
+ ωp2σT 2π((cid:15)0+(cid:15)b)11+((cid:15)0−(cid:15)b)F˜ Gσ . (22) Φ(cid:48) = −G−1(cid:16)2π11−F˜(cid:17)Φ . (26)
λ λ
Here σ is the vector of the surface charges within the
discretizedsurfaceelements, ωp =(4πn0)12 istheplasma The scalar potential Φ˜λ at the position of the molecu-
frequency, and F˜ is the surface derivative of the Green lar dipole is then obtained by means of Green’s second
identity [Eq. (D1)], and the coupling between the molec-
function with respect to the second argument. Equa-
ulardipoleandthesurfaceplasmonu fromthestandard
tion (22) accounts for undamped plasma oscillations. λ
expression g = ρ Φ˜ . In our numerical approach we
Damping as well as coupling to the molecular dipole has λ dip λ
again describe the charge distribution ρ of the molec-
tobeincludedadditionally,aswewilldiscussfurtherbe- dip
ular dipole by means of two oppositely charged particles
low.
separated by a small distance. Within this approach we
Before doing so, we show that Eq. (22) allows us to
then obtain for the coupling between the molecule and
compute eigenmodes for the surface charge oscillations.
the surface plasmons the interaction Hamiltonian
As the matrices appearing in the Hamilton function are
real and symmetric, they can be diagonalized simultane- (cid:88)(cid:16) (cid:17)
ously. Let ωλ2 and uλ denote the eigenvectors and eigen- Hpl−mol = gλσˆ+aˆλ+gλ∗σˆ−aˆ†λ , (27)
values of the generalized eigenvalue problem λ
(cid:16) (cid:17)−1 wherewehavemadeuseoftherotating-waveapproxima-
ωλ2 2π11+F˜ Guλ tion. σˆ± are the molecular excitation and de-excitation
operators introduced in Eq. (1). The first term on the
(cid:104) (cid:105)−1
=ω2 2π((cid:15) +(cid:15) )11+((cid:15) −(cid:15) )F˜ Gu . (23) right-hand side of Eq. (27) describes the creation of a
p 0 b 0 b λ
surface plasmon through de-excitation of the molecule,
The eigenvectors u can be chosen real and are orthog- and the second term the reversed process of plasmon an-
λ
onal in the sense uTλ (2π11+F˜)−1Guλ(cid:48) =βλδλλ(cid:48), where nihilation and molecule excitation.
6
whereldenotestheusualdegreesofsphericalharmonics.
For the coupling constants we get
(cid:115)
ω3
glm = lnl nbal+12 alm. (29)
0
Hereaistheradiusofthesphere,anda aretheexpan-
lm
sioncoefficientsofthemoleculardipolegiveninappendix
C.
3. Surface plasmon dynamics
Withthequantizationofthesurfaceplasmonswehave
opened the quantum optics toolbox, which would allow
us to investigate strong coupling or other nonlinear ef-
fects [33]. In this work, however, we will use the surface
plasmon eigenmodes only to analyze the behavior of the
radiativeand non-radiativedecay ratesγ and γ . Pur-
r nr
suing a similar open-system approach as in appendix A,
but for the molecule-plasmon coupling (27), we can ex-
press the nonradiative decay rate as
γnr =(cid:88)|g |2(cid:90) 0 e−iωt(cid:68)aˆ (0)aˆ†(t)(cid:69) dt. (30)
2 λ λ λ
λ −∞
Heretheexponentialaccountsforthefreepropagationof
themolecule,andthetermsinbracketsdescribetheplas-
moncorrelationfunctiontobeevaluatedatzerotemper-
ature. The latter is of the form exp[i(ω −iγ)t], where
λ 2
γ is the Landau damping of the plasmons. In a more
general approach we could have also introduced damp-
ing through a master equation approach with suitable
Lindblad operators. We can finally solve the scattering
FIG. 1: (a) Comparison of BEM and Mie calculations for
integral of Eq. (30) to arrive at
scattered light intensity of a spherical gold nanoparticle, as
cdfooiremledpciufftrteiercdefnuwtnisctuhtriiofnancteohfedRiqsecufr.aes[t3iisz9ta]a.ttiTiocnhasepwspyirtmohxbi(omcl)satr1ie4op4nrevasenerndtticfreoesrs,ut(lhdtes) γ2nr ∼=(cid:88)λ |gλ|2 (ωλ−ω)γ22+(cid:0)γ2(cid:1)2 . (31)
400vertices,(e)400verticesandafinergridaroundthenorth
Again we have neglected the frequency renormalization
pole, and (f) 900 vertices. (b) Nonradiative decay rate of a
due to plasmon coupling. Eq. (31) shows that, through
molecule in vicinity of the spherical nanoparticle, in units of
aneigenmodeexpansion,thenonradiativedecayratecan
free-space decay rate γ0 and for a molecular dipole oriented
r be decomposed into the different contributions for each
in radial direction. The transition energy of the molecule
is 2.4 eV [see arrow in panel (a)], and the diameter of the plasmonmode. Eachcomponentisgivenbythesquareof
nanosphereis10nm. TheinsetreportstheratiooftheBEM the coupling constant together with a Lorentzian whose
scattering rate with the exact Mie result. broadening is given by the Landau damping. Note that
in the limit γ → 0 the Lorentzian would reduce to
πδ(ω −ω). In addition to Eq. (31), the radiative scat-
λ
teringratecanbeobtainedbycomputing,withinlowest-
2. Mie theory
orderperturbationtheory,theoccupationoftheplasmon
modes,andmultiplyingitwiththedipolemomentofthe
Forsphericalparticlesitbecomespossibletoobtainan-
corresponding mode.
alyticexpressionsforthesurfaceplasmonenergiesω and
λ
coupling constants g by means of Mie theory. We here
λ
provide for comparison the corresponding results which III. RESULTS
werefirstobtainedbyRitchie[27]. Theplasmonenergies
read
A. Spherical particles
(cid:115)
l
Figure 1(a) reports the scattered light intensity from
ω =ω , (28)
l p (l+1)(cid:15)b+l(cid:15)0 a spherical gold nanoparticle, as computed within our
7
spherediscretizationwith144vertices[seepanel(c)]sig-
nificantly deviates already at a distance of 2 nm, while
the results for sphere discretizations with 400 and 900
vertices [panels (d) and (e)] show nice agreement down
to values of about 1 nm. Things improve considerably
for the special sphere discretization shown in panel (e).
Here we have introduced a finer mesh around the north
pole where the molecule approaches the sphere. As evi-
dent from the inset of panel (b), the difference between
the BEM and Mie calculations is at most ten percent,
even for such small distances as 0.2 nm.
B. Other particle shapes
In the following we consider different particle shapes
where the scattering properties cannot be obtained an-
alytically. Figure 2 reports optical spectra for different
particleshapes,whichareshowninthelowerpanel. Inall
cases the light polarization is along the long axes of the
nanoparticles. Theplasmonenergyofthedipolemodein-
creases from the quasi one-dimensional rod (a) and cigar
(b),overthequasitwo-dimensionalellipsoid(c)anddisk
FIG. 2: Scattered light intensity for different nanoparticle (d),tothesphere(e). Thesolidanddashedlinesshowre-
shapes of (a) rod, (b) cigar, (c) ellipsoid, (d) disk, and (e) sults as obtained from the dielectric function of Ref. [39]
sphere. The ratio of height : diameter is 5 : 1 for (a,b) and
and the Drude form (20), respectively. With the excep-
1 : 5 for (c,d). The solid and dashed lines show results ob-
tionofthesphere,bothresultsareinniceagreementand
tainedforthedielectricfunction[39]andtheDrudeform(20),
thusjustifytheDrudedescriptionforquasione-andtwo-
respectively. Lightpolarizationischosenalongthelongaxes.
dimensional nanoparticles. From the comparison of the
The dots above the spectra indicate the positions of the sur-
face plasmon eigenmodes, as obtained from the solutions of results for the rod and cigar, as well as for the ellipsoid
Eq.(23). Thecorrespondingeigenmodesoflowestenergyare and disk, we observe that the detailed shape of the par-
shown at the bottom of the figure. For the considered par- ticlehasnodramaticimpactonthespectra. Finally, the
ticle shapes only the eigenmodes with lowest energy have a heightofthepeakshastheapproximateratio1:4:9for
non-zero dipole moment and can couple to light. the particles of different dimension. Indeed, this is the
behavior one would expect for a dipole radiator where
the oscillator strength is distributed between (a,b) one,
quasistatic approach, and for the dielectric function of (c,d) two, and (e) three dipole modes, where only one is
Ref. [39]. The solid line is the result of Mie theory, optically excited.
which is summarized in appendix C, and the symbols The dots in the upper panel indicate the plasmon
show results of our boundary element method approach eigenenergies ω obtained from the generalized eigen-
λ
describedinSec.IIC.Thebroadpeakataphotonenergy value problem (23). The lowest dipole modes are pre-
of approximately 2.3 eV is associated to the dipole res- cisely at the positions of the corresponding peaks in the
onance of the nanoparticle. We observe from the figure optical spectra. In the lower panel of the figure we show
thatthedifferencebetweenBEMandMieresultsreduces theeigenmodesu forthetwosurfaceplasmonsoflowest
λ
with increasing number of surface elements. energy. For the particle shapes under consideration only
Panel (b) shows the nonradiative decay rate for a the mode of lowest energy has a finite dipole moment
nanosphere,with10nmdiameter,andamoleculeplaced and can thus couple to light. A striking feature of the
at a given distance away from the particle. We assume eigenenergiesisthatwithdecreasingenergyofthedipole
that the dipole is oriented in radial direction, and set mode the energy separation to the next mode increases.
the molecular transition frequency to 2.4 eV. One ob- It is of the order of 0.4 eV for the quasi one dimensional
serves that with decreasing distance the decay enhance- particles (a,b), about 0.2 eV for the quasi two dimen-
ment drastically increases, reaching a value of approxi- sionalparticles(c,d),anddecreasesbyafurtherfactorof
mately107 forthesmallestdistanceof0.2nm. Thesym- twoforthesphere. Forthespherewealsoobserveaquasi
bols are the results of our BEM calculations. From the continuum of surface plasmon modes at energies around
inset, whichreports theratio betweenthe BEMandMie 2.5 eV, in agreement with the result (28) of Mie theory.
decay rates, one observes that the BEM results give reli- As we will show next, these different energy separations
ableresultsonlydowntomolecule-nanoparticledistances betweenplasmonstateshaveimportantconsequencesfor
comparabletothediscretizationlength. Theresultsfora the non-radiative decay rates.
8
FIG. 4: Enhancement of FRET transfer rate for two differ-
ent donor positions. The transition frequency of the donor
moleculeisinresonancewiththedipolesurfaceplasmon. We
perform an average over all dipole orientations of the donor
andacceptormolecules. Notethatinthequasistaticlimitthe
FRET enhancement does not depend on the particle size.
torof,e.g. 2or4,γ woulddecreasebyafactorof8and
nr
64, respectively, and the quantum yield increase by the
same factor. Thus, for a quasi one- or two-dimensional
metallic nanoparticle whose length is around 50 nm, the
quantum yield can easily be of the order of several ten
percent.
We can now use Eq. (31) to decompose the scattering
rate into contributions for different plasmon modes λ.
Thisiswhatisshowninthelowerpanelsofthefigurefor
(c)thecigar,(d)thedisk,and(e)thesphere. Eachpoint
in the density plots corresponds to a given molecule-
particle distance and plasmon energy, and we have in-
troducedabroadeningoftheorderofγ foreachplasmon
mode. The colors indicate the relative importance of the
different plasmon modes. For the cigar-shaped particle,
weobserveinpanel(c)thatuptoverysmalldistancesit
isprimarilythedipolesurfaceplasmonmodethatcouples
tothemolecule. Thisisbecauseofthelargeenergysepa-
ration between the plasmon modes of lowest energy, and
the resulting weak coupling to the off-resonant modes.
FIG. 3: (a) Non-radiative decay rate in units of radiative Onlyforthesmallestdistancesthecouplingconstantsg
free-space decay rate γ0 for different nanoparticle shapes of λ
r dominate over the detunings ω −ω. Similar behavior is
sphere (circles), disk (triangles), and cigar (squares). The λ
observedinpanel(d)forthedisk-shapedparticle. Things
length of the particles along the long axes are 28 nm for the
are different in panel (e) showing results for the sphere.
cigar, 16 nm for the disk, and 10 nm for the sphere. These
Because of the small energy separation between the dif-
dimensions have been chosen such that all particles have the
samevolume. Theorientationofthemoleculardipoleandthe ferentsurfaceplasmonmodes,alreadyatrelativelysmall
directionalongwhichthedipoleapproachesthenanoparticle distances,sayaround0.1inunitsofthespherediameter,
areindicatedintheinset,andthemoleculartransitionenergy themoleculestartstosignificantlyexcitethehigher-lying
is assumed to be in resonance with the dipole modes. (b) plasmonmodes. Asconsequence,thequantumyielddra-
Quantum yield γr/γnr for the particles shown in panel (a), matically drops at these distances.
andfordifferentmolecule-nanoparticledistancesmeasuredin
units of the sphere diameter of 10 nm. (c) Decomposition of
thenon-radiativedecayrateintoeigenmodes,asdescribedin
C. F¨orster resonant enery transfer (FRET)
text.
We conclude this section by briefly discussing F¨orster
Figure 3 shows (a) the non-radiative decay rate, and energy resonance transfer (FRET) in presence of
(b) the quantum yield γ /γ , for the different particle nanoparticles. Here a donor and acceptor molecule ex-
r nr
shapes shown in panel (a). In all cases the molecular change energy through dipole-dipole interaction. This
transition frequency is assumed to be in resonance with processbecomedrasticallyenhancedifthemoleculescan
the dipole modes. Note that the volume of the particles benefit from the nearfield enhancement of the metal-
corresponds to that of a relatively small sphere of 10 nm lic nanoparticle [40, 41]. For given donor and acceptor
diameter. Byincreasingthesizeoftheparticlesbyafac- dipole moments µ and µ , the enhancement of the
don acc
9
FRET process can be computed from APPENDIX A: MASTER EQUATION
|µ ·G(r ,r )·µ |2
η = acc acc don don . (32) Inthisappendixweprovidesomedetailsforthederiva-
FRET |µ ·G (r ,r )·µ |2
acc 0 acc don don tionoftheFermi-golden-ruleresult(3). Aswearedealing
r and r denote the positions of the donor and ac- withanopensystem,i.e.,amoleculeinteractingwiththe
don acc
ceptor molecule, and G and G are the dyadic Green continuumofphotonmodes,wehavetoadoptadensity-
0
functionsinpresenceandabsenceofthenanoparticle. In matrix description [33]. We assume that the molecule is
our approach, they are computed in the quasistatic ap- initiallyinthestate1. Thecorrespondingdensityopera-
proximation and for the molecular transition frequency, tor is of the form ρˆ0 =|1(cid:105)(cid:104)1|. The total density operator
in accordance to the prescription of Sec. IIC. Note that wˆ has to additionally account for the degrees of freedom
the denominator in Eq. (32) is precisely the free-space of the electromagnetic environment. Let us make the
dipole-dipole coupling. usual assumption that at time zero the molecule is de-
Figure 4 shows the FRET enhancement for a disk- coupledfromtheelectromagneticenvironment. Thetotal
shaped nanoparticle and for a molecular transition fre- density operator then factorizes wˆ0 =ρˆ0⊗ρˆR, where ρˆR
quency equal to the dipole surface plasmon energy. The is the density operator of the electromagnetic reservoir
positions of the donor molecules are indicated in panels whichweconsidertobeinthermalequilibrium. Thetime
(a) and (b). For simplicity, in the plot we assume that evolutionofwˆ(t)followsfromtheLiouvillevon-Neumann
theacceptormoleculesareplacedonasheet,andwehave equation wˆ˙(t) = −i[Hˆ (t),wˆ(t)], which we have writ-
op
averaged over all orientation of the donor and acceptor ten in the interaction representation with H the light-
op
molecules. The FRET transfer becomes enhanced by up matter coupling of Eq. (2). As we are only interested
to four orders of magnitude, in particular for transfer in the molecular dynamics, we remove the explicit de-
processes across the nanoparticle. Here, the donor and pendence on the reservoir by tracing out the degrees of
acceptormoleculeexchangetheirenergyviatheresonant freedom of the electromagnetic environment through
surfaceplasmonmode. Notethatinordertoestimatethe
(cid:16) (cid:17)
FRETefficiencyweadditionallyhavetoconsidertheen- ρˆ˙(t)=tr wˆ˙(t)=−itr [Hˆ (t),wˆ(t)] . (A1)
R R op
hancement of the radiative and nonradiative decay pro-
cesses, as discussed in Ref. [41].
Sofarwehavenotachievedtoomuchsincetheexpression
ontheright-handsidestillinvolvesthefulldensityopera-
torwˆ. Toproceed, wedescribethelight-mattercoupling
IV. SUMMARY
in lowest order time-dependent perturbation theory, also
known as the Born–Markov approximation [33],
In summary, we have developed a general framework,
based on a suitable description of quantum electrody- (cid:90) t(cid:104) (cid:104) (cid:105)(cid:105)
namics in the presence of absorbing media and a general ρˆ˙(t)∼=−trR Hˆop(t), Hˆop(t(cid:48)),ρ(t)⊗ρR dt(cid:48).
boundary element method approach, for the description 0
(A2)
ofthemodifiedmoleculardynamicsinpresenceofmetal-
Equation (A2) is correct up to second order in the light-
lic nanoparticles. Within the quasistatic limit the decay
matter coupling Hˆ . There is a subtle point regard-
rate has been separated into radiative and non-radiative op
ing the time argument of ρ(t) on the right-hand side,
decay channels, whichbothobeyasimple scalingbehav-
for which we refer the interested reader to the literature
ior. We have discussed different nanoparticle discretiza-
[42, 43]. We next insert the light-matter coupling (2) in
tions,andhaveshownthatreliableresultscanonlybeob-
rotating-wave approximation into Eq. (A2), to arrive at
tained for properly defined surface discretizations. From
the generalized master equation
a comparison with Mie theory we have estimated an er-
rorbar of less than ten percent for our numerical results. (cid:90) t (cid:104)
We have performed an eigenmode expansion for the sur- ρˆ˙(t)∼=− dt(cid:48)
face plasmons, and have shown that the higher quantum 0
(cid:68) (cid:69)
yield for quasi one- and two-dimensional nanoparticles Eˆ+(t)Eˆ−(t(cid:48)) dˆ−(t)dˆ+(t(cid:48))ρˆ(t)
i j i j
canbeattributedtothelargeenergysplittingofthecor-
(cid:68) (cid:69)
respondingeigenmodes. Finally,wehaveshownthatsig- + Eˆ+(t(cid:48))Eˆ−(t) ρˆ(t)dˆ−(t(cid:48))dˆ+(t)
i j i j
nificant enhancements of the F¨orster transfer rates can
(cid:68) (cid:69)
be expected for donor and acceptor molecules placed in − Eˆ+(t(cid:48))Eˆ−(t) dˆ+(t)ρˆ(t)dˆ−(t(cid:48))
i j j i
the vicinity of the nanoparticle.
(cid:68) (cid:69)
− Eˆ+(t)Eˆ−(t(cid:48)) dˆ+(t(cid:48))ρˆ(t)dˆ−(t)
i j j i
(cid:105)
Acknowledgments + (±)←→(∓) . (A3)
We are grateful to Joachim Krenn and Alfred Leitner Here we have introduced the shorthand notation (cid:104).(cid:105) =
for most helpful discussions. tr (ρˆ .), have exploited the cyclic permutation of the
R R
10
fieldoperatorsunderthetrace,andhaveusedthattheex- the explicit form (1) for the dipole operators and tak-
pectation values (cid:104)Eˆ±Eˆ±(cid:105) vanish in thermal equilibrium. ing the matrix elements ρ = (cid:104)1|ρˆ|1(cid:105) of the operator
i j 11
The four terms explicitly given on the right-hand side equation (A6), we find that the upper-state population
of Eq. (A3) correspond to photon emissions, and those of the molecule decays according to ρ˙ = −γρ , with
11 11
resulting from the exchange of (±) with (∓) to photon the decay rate γ given in Eq. (3).
absorptions. Asinthisworkwewillbedealingwithelec-
tromagneticfieldsinthevisibleregimeandroomtemper-
ature,wecansafelyignorethethermaloccupationofthe
APPENDIX B: CURRENT CORRELATION
radiation modes, and thus neglect absorption processes
FUNCTION
of thermal photons throughout.
Let us assume that the dipole operators oscillate
In this appendix we provide some details of how to re-
with the molecular transition frequency ω according to
dˆ±(t)=e∓iω0tdˆ±. In addition, we introdu0ce the Fourier late the current correlation function (cid:104)ˆjk+(r,ω)ˆjl−(r(cid:48),ω)(cid:105)
to the imaginary part of the dielectric function of the
transform of the field operators
metal. We will only consider local and isotropic me-
Eˆ±(t)=(cid:90) ∞e∓iωtEˆ±(ω)dω . (A4) dia, such that the expression given above reduces to
2π δ δ(r −r(cid:48))(cid:104)ˆj+(ω)ˆj−(ω)(cid:105). The calculation of such cor-
0 kl
relation functions is a common problem in solid state
We will use the same symbol Eˆ for the field operators in
physics [45, 46].
timeandfrequencyspaceaslongasthereisnodangerof
Westartwithageneralresultoflinearresponsetheory
confusion. For sufficiently large times t, where ω t (cid:29) 1
0
is fulfilled, we can replace the lower limit 0 in the time 4π (cid:90) ∞ (cid:68) (cid:69)
integration (A3) by −∞. This replacement is called the (cid:15)(cid:48)(cid:48)(ω)= (cid:61)m eiωt [ˆj+(t),ˆj−(0)] dt, (B1)
ω2
adiabaticapproximation [33,43]andisusuallyconsidered 0
tobevalidunderbroadconditions. Thetimeintegrations
which relates the imaginary part of the dielectric func-
then reduce to integrals over exponentials, which can be
tion to the retarded current correlation function [34].
solved, with the relative time τ =t−t(cid:48), according to
A common link between the ordered and retarded cor-
(cid:90)−0∞e−i(ω0−ω+i0+)τdτ = ω0−ωi+i0+ ∼=πδ(ω−ω0). rρe(lta)t=ion(cid:68)[fˆju(nt)c,tˆjio†(n0)i]s(cid:69).prToovicdoemdpbuytetρh(et)swpeecwtroaulldfunnecetdiona
(A5)
microscopicmodelfortheelectrondynamicsinthemetal,
The infinitesimally small and positive quantity 0+ has
which is certainly beyond the scope of our work. How-
been introduced to provide a damping of the exponen-
ever, inordertoobtain,inthermalequilibrium, therela-
tial at early times. Its appearance has nothing to do
tionbetweenthedifferentcorrelationfunctionsitsuffices
with a real damping process. Rather it allows to per-
to assume that the states |m(cid:105) and energies E of the
m
form the adiabatic limit in the integration prior to the
full many-body hamiltonian Hˆ could be determined in
thermodynamic limit in the expectation value (cid:104).(cid:105) [44].
principle. In fact, we never have to compute these states
Finally, in the last expression on the right hand side of
explicitly. Inserting in (cid:104).(cid:105) = Z−1tr(e−βHˆ .) a complete
Eq. (A5) we have only considered Dirac’s delta function,
set of states, where β =1/(k T) is the inverse tempera-
which accounts for the energy conservation in the scat- B
tureandZ thepartitionfunction,wethenobtainforthe
tering processes, and have neglected the Cauchy princi-
Fourier transform of the spectral function [45, 46]
pal value. The latter contribution is responsible for an
energy shift of the molecular transition, usually known
as the Lamb shift [33], which could be incorporated into ρ(ω) = (cid:0)1−e−βω(cid:1)Z−1(cid:88)e−βEm(cid:12)(cid:12)(cid:104)m|ˆj+|n(cid:105)(cid:12)(cid:12)2
(cid:12) (cid:12)
a renormalized ω . In thermal equilibrium, the expec-
0 m,n
tation value of the field operators (cid:104)Eˆ+(ω)Eˆ−(ω(cid:48))(cid:105) =
i j ×2πδ(ω−[E −E ]). (B2)
n m
2πδ(ω−ω(cid:48))(cid:104)Eˆ+(ω)Eˆ−(ω)(cid:105)isdiagonalinfrequencyspace.
i j
Putting together all the results and changing back from From Eq. (B2) one then obtains (cid:104)ˆj+(ω)ˆj−(ω)(cid:105) =
the interaction to the Schr¨odinger representation, then ρ(ω)/(1−e−βω). Similarly, upon insertion of the many-
brings us to the master equation
body statesinto theretarded correlationfunctionwe get
ρˆ˙ ∼= −i(cid:104)Hˆ ,ρˆ(cid:105)−π(cid:68)Eˆ+(ω )Eˆ−(ω )(cid:69) ρ(ω) = ω2(cid:15)(cid:48)(cid:48)(ω)/(2π). Putting together all results we
0 i 0 j 0 arrive at our final expression
(cid:16) (cid:17)
× dˆ−dˆ+ρˆ+ρˆdˆ−dˆ+−2dˆ+ρˆdˆ− . (A6)
i j i j j i (cid:68) (cid:69) ω2(cid:15)(cid:48)(cid:48)(ω)
ˆj+(ω)ˆj−(ω) =[n¯ (ω)+1] ,, (B3)
The first term on the right hand side accounts for the th 2π
free propagation of the molecule, and the second term
for scattering processes due to the coupling to the elec- withn¯ (ω)=1/(eβω−1)theBose-Einsteindistribution
th
tromagnetic fields. The terms in parentheses are remi- function. For the problem of our present concern we can
niscent of the Lindblad master equation. Finally, using safely neglect the thermal occupation n¯ of photons.
th