Decay dynamics of Localised Surface Plasmons: damping of coherences and populations of the oscillatory plasmon modes

Properties of plasmonic materials are associated with surface plasmons - the electromagnetic excitations coupled to coherent electron charge density oscillations on a metal/dielectric interface. Although decay of such oscillations cannot be avoided, there are prospects for controlling plasmon damping dynamics. In spherical metal nanoparticles (MNPs) the basic properties of Localized Surface Plasmons (LSPs) can be controlled with their radius. The present paper handles the link between the size-dependent description of LSP properties derived from the dispersion relation based on Maxwell's equations and the quantum picture in which MNPs are treated as"quasi-particles". Such picture, based on the reduced density-matrix of quantum open systems ruled by the master equation in the Lindblad form, enables to distinguish between damping processes of populations and coherences of multipolar plasmon oscillatory states and to establish the intrinsic relations between the rates of these processes, independently of the size of MNP. The impact of the radiative and the nonradiative energy dissipation channels is discussed.

The optical features of metal nanoparticles (MNPs) are dominated by their ability to resonate with the EM radiation. Excitation of collective surface charge oscillations in form of the standing waves known as Localized Surface Plasmons (LSPs) gives rise to a variety of effects such as the near-field concentration to the region well below the diffraction limit or the far-field resonant absorption and scattering in the spectral ranges which can be manipulated by MNP's dimensions, shape and the dielectric properties of the environment [27][28][29][30][31].
Size-dependent spectral properties of noble MNPs (e.g. [32][33][34][35]) and plasmon dynamics are of basic importance in many studied plasmonic issues. The example can be harnessing hot electrons and holes, resulting from the decay of LSPs which recently attracted considerable attention because of their promising applications in photodetection, energy-harvesting, hot-carrierinduced chemistry, photocatalysis, etc. [36][37][38][39][40][41][42][43]. The performance of plasmonic devices correlate with numerous parameters that need to be studied to reach the optimum and desired properties. In all of them, LSP dynamics tailored by MNPs size play a crucial role.
The shapes and the spectral widths of the MNPs spectra are usually studied on the ground of the formalism known as Lorentz-Mie scattering theory with the solutions in the form of a sum of an infinite series of the spherical multipole partial waves. Such fully classical electrodynamic description of scattering of a plane wave by a sphere (e.g. [44][45][46][47][48]) have been extended for the case of non-absorbing hosts [11], multi-layered spheres [49,50], and arbitrary incident light beams [51]. LSP damping rates manifested in the broadening of the NPs spectra were experimentally studied in the case of the dominant dipole contribution to the spectra for NPs with available, limited sizes [33,34,[52][53][54][55] (and references therein). The resulting damping times and the frequencies corresponding to the maxima in the far-field intensity signals) are often suggested to directly characterize the LSP damping process and LSP resonance position.
However, the desired parameters (both for fundamental reasons and with a view to applications) such as size-dependent oscillation frequencies of multipolar modes corresponding to LSP's resonances, or the decay rates (times) of the excited oscillations, are not the explicit parameters of Lorentz-Mie scattering theory. In [56][57][58][59] such functions were derived from considering the classical dispersion relation [60,61] for the surface localized EM (SLEM) fields basing on the self-consistent Maxwell divergent-free equations. The modelling [56][57][58][59] provides the explicit size dependence of the oscillation frequencies of SLEM fields and of the damping rates of such oscillations for the dipole and higher order multipolar surface modes and goes beyond the limitations related to the MNP's size including those which result from the often used quasistatic ap-proximation (e.g. [35]) which kills the size dependence (see also [17] in case of cavity QED formalism using the electrostatic electric field potential).
However, the such classical description of the dynamics of LSPs damping imposes constraints to fully understand the LSP dynamics including the damping processes, because it does not distinguish between the dephasing and population damping. In case of the dipole LSP such quantities were suggested to be connected to each other [34,[53][54][55] and composed of radiative and nonradiative decay into electron-hole excitations. However, more detailed analysis of the problem was not given.
The present paper is aimed to present such an analysis by adopting the theory of quantum open systems and formulating the conclusions which apply to the issue of the LSP damping.
Plasmonic phenomena are inherently quantum (see e.g. [62][63][64]). Quantum plasmonics has recently emerged as a new fascinating field of research with a view of observing quantum phenomena in light-matter interactions at the nanoscale ( [22,63,[65][66][67] for reviews). In particular, there have been a large number of theoretical studies of interactions between a (dipole) emitters and confined plasmonic structures in weak or also in strong-coupling regime using different approaches ('macroscopic' QED using Green's functions, quasinormal decomposition, resonant-state-expansion) [17-26, 68, 69]. The confinement of light field to scales far below that met in the case of the conventional optics by metalic nanostructures enables to describe the interactions between atomic systems and plasmonic structures in the formalism of cavity QED and (see e.g. [17]) to explore the potential of such systems [70] for developing future quantum technologies like single-photon transistor or for carrying quantum information (see [66,67,71] for reviews).
In this paper we describe the intrinsic dynamics of LSP which includes the relaxation pathways leading to the dephasing and population damping of plasmon oscillations within a quantum dynamical description of open systems. The master equation in the Lindblad form within the Born-Markov approximations [72] is used to describe the evolution of an open quantum system of N electrons confined in an MNP in absence of the driving EM field. Conclusions allow to generalize the common understanding of the LSP damping as the process which consists not only from the dephasing of LSP oscillations, but also from the population damping at the twice larger rate, irrespective the MNPs size and LSP multipolarity.
The results of such modeling have been related to the previously derived multipolar damping rates versus MNPs radius which describe the LSP decoherence process. Such rates resulted from the dispersion relation for the surface localized electromagnetic (SLEM) fields [56][57][58][59] which we shortly reconsider for completeness of the modelling. Derived in absence of the illuminating light, the dephasing and population damping rates tailor the transient LSP dynamics but manifest also in the amplitudes and widths of the spectra. The size dependence of studied parameters allows predicting the optimum and desired properties of LSPs, being a useful tool in tailoring MNPs plasmonic perfor-mance in experiments and applications also in the size regions still practically unavailable.

Plasmon oscillatory eigenstates (classical description)
Energy levels of atoms and molecules manifest in transitions from one energy level to another when they absorb and emit light. We use a similar picture of energy levels for a plasmonic system of N free-electrons confined in a spherical MNP. The starting point is based on the results of rigorous classical electrody- Fig. 1 The scheme leading to the dispersion relation for SLEM fields (left) and the example of the resulting discrete complex eigenvalueshω l − ihΓ l /2 for the consecutive modes l = 1, 2, 3... (right) for Au MNPs embedded in water [59]. namic description based on the self-consistent divergent free Maxwell equations (with no external sources). This problem was completely solved in a classical paper by Mie (1908) [44]. However, in the present study based on the formulations [21, 56-59, 61, 73-75], unlike in the case of the popular Lorentz-Mie scattering theory (e.g. [44][45][46][47][48]), the radiation illuminating the MNP is absent. A spherical metal/dielectric interface forms a cavity [21] which allows excitation of diverse modes of the SLEM fields in the form of the standing waves adjusted to the MNP's dimensions. Such modes can be excited after MNP is illuminated by the light within the appropriate frequency range. The continuity relation at the spherical interface leads to a set of separate dispersion relations for each multipole mode of the partial wave. The solutions exist for TM polarized modes only, because contrary to TE modes, TM (or p polarized electric waves) posses the non-zero normal to the surface (radial) components of the electric field, which are able to couple with the surface free charges.
The dispersion relations for SLEM fields define the complex, discrete eigenvaluesh(ω l − iΓ l /2) and connect the allowed oscillation frequency ω l and the corresponding damping rate Γ l of the l'th mode oscillations to the (inverse Fig. 2 Ascribing the energy levels (Fig. b)) to the oscillation energieshω l (R) (Fig. a)) which resulted from the dispersion relation for SLEM fields.
of) NS's radius R. The size dependence of ω l (R) and Γ l /2(R) define the dynamics of SLEM fields: exp (iω l − Γ l /2)t) and is unambiguously determined by the material properties of the MNP and its dielectric environment. Found in absence of the illuminating radiation, ω l and Γ l inherently characterize an MNP of the radius R in the same way as the energy levels and the inverse of lifetimes characterize an atom or a molecule. In both cases, these quantities manifest in the spectra, when the systems are illuminated. The example of ω l (R) and Γ l (R) dependence for gold MNPs embedded in water [59] is shown in Figure 1.
Let us note, that the problem of eigenmodes for plasmonic resonators recently attracted new interest resulting in a number of theoretical studies of quasinormal modes for plasmonic resonators and open cavities. In particular, based on the concept of the resonant-state expansion (e.g. [21,26,[76][77][78]), it was confirmed, that an open optical system like a plasmonic confined structure can be characterized by the complex eigenfrequencies with the real and imaginary parts corresponding to, respectively, the spectral positions of the resonances and defining the spectral linewidths of resonances. However, such a classical description does not include populations of the plasmon modes nor their damping. The next step is to replace the oscillation energieshω l (orhω r l ) of the classical modes by the discrete energy levels, which are distinct from the zeroenergy non-oscillatory level by the energieshω l (R) (orhω r l (R)) (see Figure  2a),b). The corresponding states of the plasmonic systems S in the Hilbert space are |ψ l with l = 1, 2, 3... (see Figure 2c)). The only possible transitions are those with the absorption or emission of a photon with the energyhω l . Such transitions occur between the state |ψ l and the non-oscillatory state |ψ 0 .
3 The density matrix for N electrons of the plasmonic system and the quantum master equation The density matrix (the density operator) is an alternate representation of the state of a quantum system which is very convenient for systems in the mixed states and in time-dependent problems. One of the reasons to consider the mixed states is the entanglement of the systems with the environment. The diagonal elements of the density matrix correspond to the probabilities p n = N n /N of occupying a quantum states |ψ n , n =0,1,2... , so they describe the relative populations of these states. The complex off-diagonal elements of the density matrix in the basis {|ψ l , |ψ 0 } contain a time-dependent phase factors that describe the evolution of coherent superposition of the states.
The plasmonic system S we consider consists of N electrons which in general can be distributed over the states |ψ n , n = 0, l = 1, 2... with |ψ 0 for the ground, nonoscillatory states and |ψ l for the excited oscillatory states of the plasmon modes l (see Figure 2). Coherences between the states |ψ 0 and |ψ l can be created after the interaction of the system e.g. with the external EM field. As no physical system is absolutely isolated from its surroundings, the plasmonic system S has to be considered as an open quantum system which is a subsystem of a larger combined quantum system S + E, where E represents the environment to which the open system S is coupled. Following the main assumption of the basic theory of open quantum systems [72], the environment is assumed to be a large system with an infinite number of degrees of freedom, (with a continuous and wide spectrum of characteristic frequencies) in thermal equilibrium. Therefore, the state of the system E is practically unaffected by coupling to the system S. The interaction of the open system S with the environment causes an irreversible behavior of the open system S and leads to decoherence (randomization of phases) and dissipation of energy into the surroundings.
The evolution of the open quantum system S can be described by using a master equation which defines the dynamics of the reduced density operator ρ S (t) = tr E ρ(t) obtained after tracing out the environment degrees of freedom (e.g. [72]). Quantum Markovian process represents the simplest case of the dynamics of an open systems. Within the Markov approximation, the memory effects of the system under influence of the environment are negligible. The correlation functions of the environment decay sufficiently fast over the time τ E which is small compared to the characteristic timescale of the relaxation processes τ S of the system S (e.g. [72,[79][80][81]).
The dynamics of open systems in the case of Markov processes can be described by a first-order linear differential equation for the reduced density matrix , which is known as a quantum Markovian master equation in Lindblad form [81][82][83]: where ρ S (t) is the reduced density matrix of the system S, and D[ρ S (t)] is the so-called dissipator: Summation over k extends over all processes of coupling with the environment. The first term on the right-hand side of the eq. (1) describes the unitary evolution of the system S under the action of a Hamiltonian H. The dissipator D[ρ S ] describes the environmental influence on the system state (e.g. [72]). The 'jump' operators L k (eq. (2)) describe a random evolution of the system which suddenly (at the time scale of the evolution) changes under the influence of the environment. Each L k ρ S L † k term induces one of the possible quantum jumps, while the remaining terms are needed to normalize properly the case in which no jump occurs. The simplest quantum system is a two-level system whose Hilbert space is spanned by two states, an excited state and a ground state. Such a two-level system is a very important basic model for an atom or a system of spins which is often used in quantum mechanics. The two-level system can be successfully used, provided that the transitions to other levels can be neglected.
In case of LSPs such picture can be related to the classical description based on Maxwell equations (see Section 2), where the problem is solved separately for each EM mode l of the SLEM field, and the final result is a sum over the solutions for consecutive modes with l = 1, 2, 3.... Basing on such analogy, we describe the plasmonic system S as a sum of S l of independent, open subsystems: Each subsystem S l is a two-level system: the excited states |ψ l and the ground non-oscillatory state |ψ 0 . The state of each subsystem S l can be described by the 2×2 matrix operator ρ S l (t). Each system S l is coupled to the environment independently and all assumptions about the coupling of the system S l to the environment remain fulfilled for subsystems S. The dynamics of the system S thus results from the independent dynamics of the systems S l : governed by the Lindblad eq. (1). The form of the Lindblad equation guarantees, that also dynamics of each matrix operator ρ S l (t) is governed by the equation in the Lindblad form. Therefore we can use the standard basis {|ψ l , |ψ 0 } represented by the two-dimensional column vectors: and represent the operator ρ S l (t) in matrix form: The diagonal elements of the matrix ρ S l represent the relative populations (more precisely, the population probability densities) of the relevant states. The off-diagonal elements represent quantum-mechanical coherences and are complex conjugates of each other, carrying the same information.

Dynamics of the radiative plasmon damping
An excited plasmon, similarly as an excited atom, decays to the state of lower energy spontaneously emitting a photon. In the theory of open quantum systems, such decay is assumed to be due to the coupling of the system to the EM vacuum (fluctuating) fields. The random vacuum fluctuations (which pervades all space even at zero temperature) cause random jumps influencing plasmon evolution. In the case when only radiative damping is present, the Hamiltonian in the Linblad eq. (1): H r l = H l (eq. (8)) and the eigenenergy E l =hω r l . A single jump operator describs random, sudden emission of a photon from the state |ψ l to the state |ψ 0 of the S l , with Γ r l being the (spontaneous) radiation rates. The Lindblad master equation (1) takes then the form: The evolution of the diagonal and off-diagonal elements of ρ S l d dt leads to the solutions (t 0 = 0): Therefore, the radiative damping rate Γ r pop l = Γ r l of populations ρ ll is twice as large as the radiative rate Γ r coh l = Γ r l /2 of coherences ρ l0 = ρ * 0l , so the relation between the corresponding radiative lifetimes of populations and coherences are 2T r pop l = T r coh l .

Dynamics of the total (radiative and collisional) plasmon damping
Electrons in a metal inevitably undergo collisions which lead to the dissipation of energy and release of heat. To account for nonradiative relaxation processes, the system S interacting with fluctuations of the vacuum is assumed to be immersed in a dissipative heat-bath in a thermal equilibrium state with an infinite number of degrees of freedom. The heat reservoir dynamics is assumed to be much faster than those of the open system S l , so the dynamics of S l (and those of S) is Markovian. The jump operators which describe the collisional transition from the state |ψ l to the state |ψ 0 are: where Γ nr l are the nonradiative rates describing the collisional processes leading to the heat release. Summing over radiative and nonradiative contributions in the dissipator D[ρ S ] (eq. (2)) we get the master equation: The random jumps in the evolution of ρ S l (t) with the rates Γ r l and Γ nr l are assumed to be uncorrelated. The relaxation of the diagonal and off-diagonal elements of ρ S l is governed by the total relaxation rate Γ l = Γ r l + Γ nr l . The Hamiltonian H l (eq. (8)) defines the eigenenergies E l =hω l of the excited states |ψ l in presence of the radiative and nonradiative damping processes.
Using the formalism recalled in the previous section and the equation 4, the evolution of the whole system S is described by the dynamics of the density matrix ρ S (t), (eqs. (4,6)): with the diagonal and off-diagonal temporal dependence given by: The populations of the oscillatory states are exponentially damped with the rates: while the oscillations of coherences are damped with two times lower rate: So the lifetime of coherences is twice as large as the lifetime of populations. The temporal evolution of coherences ρ S l proceeds according to exp(iω l − Γ l /2)t, similarly to the evolution of SLEM fields, which resulted from the dispersion relation (Section 2). So, the size dependence of ω l and Γ l in eq. (20) can be found by solving the dispersion relation for SLEM field for successive R + ∆R, as it was done in our previous papers [56][57][58][59].
The density matrix ρ S (t) (eq. (20)) fulfills all the properties which are of basic importance for quantum statistics: it is positive, self-adjoint and with the trace equal 1. The relative populations of the states |ψ n , n = 0, l: l N n (t)/N = ρ 00 (t) + l ρ ll (t) (23) are the same at any time t as the initial populations over the excited, oscillatory states |ψ l , l =1,2,3... at t = t 0 (see eq. (20)):

Pure dephasing
Pure dephasing takes place where the (external) reservoir is the source of fluctuations which do not change the average energy of the system. These fluctuations lead to a loss of coherence in the system resulting in the decay of the off-diagonal elements of the density matrix, without affecting the diagonal elements. So, if there are processes which lead to decoherence without affecting populations, the total damping rate is increased by the "pure" decoherence rate Γ * coh l : The population damping rates Γ pop l remain unchanged: In the case when the nanoparticle is illuminated by light, such 'pure' dephasing processes introduce additional broadening to the observed spectra in which LSP excitations manifest.

Discussion and conclusions
Localized Surface Plasmons in plasmonics are commonly described as the electromagnetic excitations coupled to coherent electron charge densities oscillations on a metal/dielectric interface. LSP's parameters such as damping rates of plasmon oscillations in function of MNP's size (related to the decoherence processes) are usually derived from the widths of scattering or absorption spectra for the dipole (l =1) plasmon mode only. In common practice, such damping rates (and the frequencies corresponding to the maxima in the farfield intensity signals) are suggested to describe LSP dynamics what may lead to narrow understanding of LSP damping, which in fact is the process embracing dephasing and depopulation of all the plasmon modes involved, including those not manifested as distinct maxima in the spectra. Moreover, it is known that in general LSPs resonances can manifest in different manner in divers spectra in the near-and far-field regions [84][85][86][87][88].
Solutions of the dispersion relation for SLEM fields (see the example in Figure 1) allow to find the intrinsic LSP dephasing rates and resonance frequencies of multipolar plasmons in absence of illumination and predict their size dependence in the large range of MNP radii [56][57][58][59].
In the applied quantum description an MNP is treated as "quasi-particle", similarly to an atom or molecule. Such picture enables studying the intrinsic decay dynamics of both: populations and coherences of the quasi-particles oscillatory states involving physical quantities and their relations which have not been discussed in the previous approaches. In particular, the damping rates of populations and of coherences of consecutive plasmon modes occur to be intrinsically different: Γ pop l = 2Γ coh l , regardless of the MNP size. Fig. 3 The correspondence of deriving the radius dependent quantities which describe the dynamics of the SLEM fields within classical EM modelling (left) and the corresponding scheme for the dynamics of coherences in the quantum description (right).
The present paper builds the bridge between the results of the classical electrodynamic description supplying size dependence of the dephasing rates and the quantum description of the intrinsic dynamics of plasmon decay processes in MNPs which are essentially not restricted in size (see Figure 3). Such dynamics is affected by several factors to consider such as electron dumping in bulk metal γ bulk , the electron-surface scattering γ R (the parameters of the dielectric function), radiation damping and prospectively the interface damping resulting in energetic charge carriers production. In case of solving the dispersion relation for SLEM fields using the Drude-like dielectric function (electron dumping in bulk metal with the rate γ bulk and the radius dependent electronsurface scattering rate γ R accounted), the damping rates Γ l (R) define the size dependence of the total population and and coherence damping rates resulting from the quantum modelling: Γ l (R) = Γ coh l (R) = Γ pop l (R)/2 (see Figure 3). Such link can be a useful tool in tailoring MNPs plasmonic performance in experiments and applications also in the size regions still practically unavailable. Suggested (Figure 3) decomposition of the contribution of the radiative rates from the total damping rates after appropriate modification of the input dielectric function of gold and silver (γ bulk =0) will be the subject of our next study.
The quantum picture offers the attractive prospects for designing the plasmonic devices that operate at the quantum level and exploit their lossy nature (e.g. [36,[89][90][91]), in addition to the broad range of applications based on the enhancement of EM fields at metal/dielectric interfaces. In such a practical context, accounting for various dissipative decay channels and quantification of their importance in the plasmonic systems built of diverse material, size, and shape, embedded in various matrices, is of basic importance.