Kohn-Sham decomposition in real-time time-dependent density-functional theory: An efficient tool for analyzing plasmonic excitations

The real-time-propagation formulation of time-dependent density-functional theory (RT-TDDFT) is an efficient method for modeling the optical response of molecules and nanoparticles. Compared to the widely adopted linear-response TDDFT approaches based on, e.g., the Casida equations, RT-TDDFT appears, however, lacking efficient analysis methods. This applies in particular to a decomposition of the response in the basis of the underlying single-electron states. In this work, we overcome this limitation by developing an analysis method for obtaining the Kohn-Sham electron-hole decomposition in RT-TDDFT. We demonstrate the equivalence between the developed method and the Casida approach by a benchmark on small benzene derivatives. Then, we use the method for analyzing the plasmonic response of icosahedral silver nanoparticles up to Ag$_{561}$. Based on the analysis, we conclude that in small nanoparticles individual single-electron transitions can split the plasmon into multiple resonances due to strong single-electron-plasmon coupling whereas in larger nanoparticles a distinct plasmon resonance is formed.

The structure of the article is as follows. In Sec. II we derive the linear response of the time-dependent density matrix in the KS electron-hole space. We review the formulation of the same quantity in the Casida ap-proach and describe the decomposition of the photoabsorption spectrum in KS electron-hole contributions. In Sec. III we benchmark the numerical accuracy of the implemented method by analyzing the KS decomposition of small benzene derivatives using both the real-timepropagation and the Casida method. This is followed by an analysis of the plasmonic response of large silver nanoparticles, which yields microscopic insight into the plasmon formation in nanoparticles. In Sec. IV we discuss the general features of the presented methodology. Our work is concluded in Sec. V.

II. METHODS
A. Linear response of the density matrix in the real-time propagation method The time-dependent Kohn-Sham equation is defined as where H KS (t) is the time-dependent KS Hamiltonian and ψ n (r, t) is a KS wave function. The density matrix operator is defined as where f n is an occupation factor of the nth KS state. In order to proceed with KS decomposition, we express the density matrix in the KS basis, spanned by the groundstate KS orbitals ψ (0) n (r), which fulfill the ground-state KS equation where H KS is the ground-state KS Hamiltonian and n the KS eigenvalue of nth state. The density matrix can be written in this KS basis as This equation establishes a link between a timedependent density matrix and the usual KS (electronhole) basis set used in linear-response calculations, see Sec. II B. Previously, similar or related quantities have been used within the real-time propagation method for analysing the response. [45][46][47][85][86][87][88] When the real-time propagation method is applied in the linear-response regime, the usual approach is to use a δ-pulse perturbation. 82,83 This corresponds to the Hamiltonian where the interaction with external electromagnetic radiation is taken within the dipole approximation. The electric field is assumed to be aligned along the z direction and the constant K z is proportional to the external electric field strength, which is assumed to be small enough to induce only negligible non-linear effects. After the perturbation by the δ-pulse at t = 0, Eq. (1) is propagated in time and the quantities of interest are recorded during the propagation. As a post-processing step time-domain quantities, e.g., ρ nn (t), can be Fourier transformed into the frequency domain.
It is important to note that the size of the density matrix ρ nn (t) can be significantly reduced since only its electron-hole part is required in linear-response theory. 75,77 It is thus sufficient to consider only ρ ia (t), where i and a represent occupied and unoccupied KS states, respectively. Then, we obtain the linear-response of the density matrix in electron-hole space as where ρ ia (0 − ) is the initial density matrix before the δpulse perturbation and the superscript z indicates the direction of the perturbation.
In common TDDFT implementations, there is no mechanism for energy dissipation and the lifetime of excitations is infinite. A customary way to restore a finite lifetime is to apply the substitution ω → ω + iη, where the parameter η is small. This leads to an exponentially decaying term in the integrand in Eq. (6), i.e., e iωt → e iωt e −ηt , and to the Lorentzian line shapes in the frequency domain. The decaying integrand also means that a finite propagation time is sufficient in practical calculations. The Gaussian line shapes can be obtained by replacing the Lorentzian decay e −ηt with the Gaussian decay function e −(σt) 2 /2 , where the parameter σ determines the spectral line width.

Implementation
We have implemented the density matrix formalism outlined above in the RT-TDDFT code 71 that is part of the open source gpaw package. [93][94][95] Our implementation uses the LCAO basis set 92 and the projector-augmented wave (PAW) 96 method. In the LCAO method the wave function ψ n (r, t) is expanded in localized basis functions φ µ (r) centered at atomic coordinates with expansion coefficients C µn (t). The density matrix is reads in the LCAO basis set as Then, Eq. (4) can be written in LCAO formalism as (using implied summation over repeated indices) where S µµ = φ * µ (r)φ µ (r)dr is the overlap integral of the basis functions. A detailed derivation of Eq. (9) is given in Appendix, in which it is shown that the PAW transformation affects only the evaluation of the overlap integral.
The emphasis in our implementation is to minimize the computational footprint of the analysis methods. Thus, instead of calculating Eq. (9) at every time step during the time propagation, we only store the alreadycalculated matrix C z µn (t) at every time step. Then, as a post-processing step, we calculate ρ z µν (t) with Eq. (8) and Fourier transform the result to obtain δρ z µν (ω). The latter quantity can be subsequently transformed to δρ z ia (ω) via Eq. (9) keeping only the electron-hole part. Thus, in practical implementation, the linearity of the equations allows exchanging the order of Fourier transformation and matrix multiplications.
Finally, we note that in our experience it is advantageous to store the whole time-dependent evolution of the system, i.e., C z µn (t), as done in the present implementation. While alternative on-the-fly Fourier or other transformations would reduce the amount of required storage space, they would restrict the analysis to the set of parameters specified at the outset of the calculation.
B. Linear response of the density matrix in the Casida method In Casida's linear-response formulation of TDDFT 75,77 the response is obtained by solving the matrix eigenvalue equation yielding excitation energies ω I and corresponding Casida eigenvectors F I . The matrix Ω is constructed in the KS electron-hole space. Using a double-index ia (jb) to denote a KS excitation from an occupied state i (j) to an unoccupied state a (b), the elements of the matrix can be written as where (3), and the matrix K ia,jb represents the coupling between the excitations i → a and j → b. 75 The linear response of the density matrix at frequency ω can be obtained as 75 where the summation runs over electron-hole pairs (eh) and involves the dipole matrix elements Using the spectral , allows us to write Eq. (12) as The term G I (ω) is divergent at excitation energies ω I in the common TDDFT implementations due to the infinite lifetime of the excitations. Analogously to the time domain, a finite lifetime for the excitations can be restored by the substitution ω → ω + iη, where the arbitrary parameter η determines the lifetime. This leads to a Lorentzian line shape and the imaginary part is given by where is the Lorentzian function. Alternatively, the Gaussian line shape can be obtained by using the Gaussian function g

C. Kohn-Sham decomposition
The linear response of the density matrix in the KS electron-hole space, δρ z ia (ω), can be calculated equivalently using both the real-time propagation [Eq. (6)] and the Casida approach [Eq. (13)]. While this quantity would already allow the analysis of the response at frequency ω in terms of its components in the KS electron-hole space, a more intuitive analysis can be obtained by connecting δρ z ia (ω) to an observable photoabsorption cross-section describing the resonances of the system. First, the dynamical polarizability is given by 75 and the photo-absorption is described the dipole strength function which is normalized to integrate to the number of electrons in the system N e , i.e., ∞ 0 S z (ω)dω = N e . This is similar to the sum rule 2 is the oscillator strength of the discrete excitation I. 75 By comparing Eqs. (15) and (16), we can now define the KS decomposition of the absorption spectrum as This quantity is used to analyze the response of silver nanoparticles in Sec. III B below. Previously, similar photo-absorption decompositions have been used in the electron-hole space 88 and based on, e.g., spatial location 74,89 or angular momentum. 74

A. Benzene derivatives
To benchmark the presented methods and their computational implementation, we now analyze the optical response of the molecular systems benzene (C 6 H 6 ), naphthalene (C 10 H 8 ), and anthracene (C 14 H 10 ) using both the RT-TDDFT and Casida implementations in gpaw package. [93][94][95]97 These characteristic conjugated molecules are suited for the present benchmark as they have well-defined π → π * transitions that exhibit a systematic red-shift as the extent of the conjugated π-system increases. 98,99 As the real-time propagation uses the full timedependent Hamiltonian matrices, the end result includes contributions from all electron-hole pairs and the limit of the full KS space is automatically achieved by propagating only the occupied orbitals. This is in contrast to the gpaw implementation of the Casida approach, 97 which commonly requires setting an energy cut-off that determines the KS transitions included in the calculation of the Casida matrix. In order to ensure comparability of the results, we have included in the calculation of the Casida matrix all the transitions that are possible within the KS electron-hole space spanned by the orbitals.
Both the RT-TDDFT and Casida calculations were carried out using the default PAW data sets and the default double-ζ polarized (dzp) basis sets within the LCAO description. While these dzp basis sets might not be sufficient for yielding numerical values at the complete-basisset limit, 92,100 they are suitable for qualitative analyses and for the benchmarking study presented here. The Perdew-Burke-Ernzerhof (PBE) 101 exchange-correlation functional was employed in the adiabatic limit. A coarse grid spacing of 0.3 Å was chosen to represent densities and potentials and the molecules were surrounded by a vacuum region of at least 6 Å. The Hartree potential was evaluated with a multigrid Poisson solver using the monopole and dipole corrections for the potential.
For the RT-TDDFT calculations, we used a small time step of ∆t = 5 as in order to achieve high numerical accuracy. The total propagation time was T = 30 fs, which is sufficient for the used Gaussian broadening with σ = 0.07 eV corresponding to a full width at halfmaximum (FWHM) of 0.16 eV. The calculated photo-absorption spectra of the benzene derivatives are shown in Fig. 1. The Casida and RT-TDDFT methods yield virtually indistinguishable spectra. For conciseness, we only present an analysis for excitations along the long axis (x) of the molecules. Note, however, that the response in the other directions can be analyzed in similar fashion.

Casida approach
The response of each of the molecules is dominated by a single absorption peak (see Fig. 1), which results from discrete excitations. In Table I, we show the KS decomposition of these excitations as described by the components of the normalized Casida eigenvectors F I,ia . Due to the normalization, ia F 2 I,ia = 1 for each excitation I.
For benzene (C 6 H 6 , point group D 6h ) the excitation at 7.2 eV corresponds to the first E 1u transition from the doubly degenerate highest occupied molecular orbital (HOMO; E 1g ) to the doubly degenerate lowest unoccupied molecular orbital (LUMO; E 2u ). In the present calculations the symmetry of the molecule has not been enforced and the orbitals π −0/1 and π * +0/1 span the E 1g and E 2u symmetries, respectively. Implementationdependent numerical factors slightly lift their degeneracy and determine the exact unitary rotation between the states. Naphthalene (C 10 H 8 ) and anthracene (C 14 H 10 ) belong to the D 2h symmetry point group. In both molecules the most prominent excitation is the B 3u transition, which is mainly composed of transitions from HOMO to LUMO+1 and HOMO−1 to LUMO. While in naphthalene the other contributions amount to less than 1%, in anthracene, a minor contribution originates also from a transition from HOMO−4 to LUMO+2. I. Casida analysis of the most prominent excitations of benzene (C6H6), naphthalene (C10H8), and naphthalene (C14H10). Orbitals are enumerated with respect to HOMO (π−0) and LUMO (π * +0 ). The orbital characters are given in brackets based on the point groups D 6h (benzene) and D 2h (naphthalene, anthracene).

RT-TDDFT approach
The Casida eigenvector F I,ia considered in Table I is directly related to the linear response of the density matrix, see Eq. (13), and is employed here for benchmarking the RT-TDDFT methodology described Sec. II A. In order to proceed with comparison, consider a discrete excitation J that is energetically separated from other excitations. Since Im[G I (ω J )] in Eq. (14) is approximately zero when I = J, only the excitation J contributes in Eq. (13).
yields the components of the Casida eigenvector F J,ia . This connection allows us to calculate the Casida eigenvector also from the RT-TDDFT approach. This is demonstrated in Table II, in which we show the calculated KS decompositions at the peak energies of the photo-absorption spectrum (Fig. 1).
In the case of benzene (C 6 H 6 ), we inevitably obtain a superposition of the two underlying degenerate excitations (see Table I). We can, however, calculate the equivalent superimposed F x ia (ω) eigenvector also from the Casida approach (shown in the last column of Table II). For this quantity, we obtain an excellent match between the RT-TDDFT and Casida approaches.

B. Silver nanoparticles
TDDFT calculations of noble metal nanoparticles up to diameters of several nanometers are computationally demanding, but the have become feasible with recent developments. [70][71][72][73][74] Here, we focus on silver nanoparticles as prototypical nanoplasmonic systems with a strong plasmonic response in the visible-ultraviolet light regime. 58,59 Using the methodology described above in conjunction with a recent RT-TDDFT implementation, 71 we can analyze the response of silver nanoparticles with reasonable computational resources. For illustration, a full real-time propagation of 3000 time steps for Ag 561 can be realized in 110 hours using 144 cores on an Intel Haswell based architecture. 102 Kuisma et al. have previously studied icosahedral silver nanoparticles composed of 55, 147, 309, and 561 atoms corresponding to diameters ranging from 1.1 nm to 2.7 nm. 71 Here, we consider the same nanoparticle series and use the same geometries and computational parameters as in Ref. 71. We employ optimized LCAO basis sets 71  The calculated photo-absorption spectra of the nanoparticles are shown in Fig. 2. The non-interactingelectron spectra calculated from the KS eigenvalue differences ω ia and transition dipole matrix elements µ x ia

FIG. 2.
Photo-absorption spectra of icosahedral silver nanoparticles. The non-interacting-electron spectra shown for comparison are vertically shifted and scaled by a factor of 0.2.
are also shown to facilitate the discussion below. In Ref. 71 it was observed that the plasmon resonance is well-formed in Ag 147 and in larger nanoparticles, whereas the response of Ag 55 consists of multiple peaks, the origin of which cannot be readily resolved. In the following, we analyze the response of nanoparticles in terms of the KS decomposition, which enables us to shed light on the response of the Ag 55 nanoparticle.

Transition contribution maps
In order to analyze the response in terms of the Kohn-Sham decomposition, we present the decomposition as a transition contribution map (TCM; see Fig. 3 below), 40,107 which is an especially useful representation for plasmonic systems in which resonances are typically superpositions of many electron-hole excitations. The TCM represents the KS decomposition weight w ia (ω) at a fixed ω in the two-dimensional (2D) plane spanned by the energy axes for occupied and unoccupied states. More specifically, the 2D plot is defined by where g ia is a 2D broadening function of the discrete KS states. Here, we employ the Gaussian function with σ = 0.07 eV. The same σ parameter is also used in the spectral broadening. For the weight w ia (ω), we use the absorption decomposition of Eq. (17) normalized by the total absorption, i.e., Due to the icosahedral symmetry of the nanoparticles their response is isotropic, S x (ω) = S y (ω) = S z (ω), and the decomposition is degenerate (compare the case of benzene in Sec. III A). Alternatively, instead of Eq. (20) one could use, e.g., the normalized transition density matrix (w ia (ω) = |ρ x ia (ω)| 2 ) as the weight. Equation (20), however, has the advantage that it retains the information about the sign of the response in the KS decomposition and has a physically sound interpretation as the photo-absorption decomposition.
TCMs of the nanoparticles at different resonance energies are shown in Fig. 3 along with the density of states (DOS), which has been colored to indicate the sp and d character of the states. The latter decomposition is based on the angular momentum quantum number l µ of the LCAO basis functions indexed by µ. For example, the d character of the nth state is estimated as µ:lµ=2 |C (0) µn | 2 , where the coefficients are normalized such that µ |C (0) µn | 2 = 1.

Analysis of Ag147, Ag309, and Ag561
First, we consider the largest nanoparticles Ag 147 , Ag 309 , and Ag 561 , the TCMs of which are shown in Figs. 3(d-f). The TCMs highlight two major features in their response. First, there is a strong positive constructive contribution 41 (red features in Fig. 3) from the KS transitions whose eigenvalue differences are significantly lower than the plasmon resonance energy ω. The same low-energy sp transitions are responsible for the strong peaks in the non-interacting-electron spectra (see Fig. 2), which are indicated in Fig. 3 by dashed lines. Thus, TCM shows how the resonance energy is blue-shifted as the interaction is turned on from the non-interacting case (λ = 0) to the fully interacting one (λ = 1). This demonstrates the plasmonic nature of the excitation in the socalled λ-scaling approach for plasmon identification, 39,108 and illustrates the importance of low-energy transitions for plasmon formation. 47 Another prominent feature in the response is the damping due to d electrons, which is seen in the TCMs as large negative contributions from occupied d states into unoccupied states (blue features at ε o ≈ −4 eV in Fig. 3). Interestingly, the plasmon peak appears close to the onset of d electron transitions, corresponding to the intersection of the line ε u − ε o = ω and the horizontal Fermi level line. Generally, with increasing nanoparticle size the DOS becomes increasingly continuous, which is also visible in the increasing uniformity of the TCMs.
In Ref. 73, TCMs for charged silver nanoparticles up to Ag 309 have been studied. The two main features in Fig. 3, the low-energy sp transitions and the d electron damping, are in agreement with these TCMs reported earlier. In contrast to Fig. 3, the TCMs in Ref. 73 show, however, also a significant contribution from sp transitions close to the ε u − ε o = ω line. We consider this to be due to the different choice of the TCM weight w ia (ω) in Ref. 73. In the absorption decomposition we used in Fig. 3 [Eqs. (17) and (20)] the KS components are essentially weighted with the dipole matrix element µ x ia , which affects the relative magnitudes observed in TCM.

Analysis of Ag55
Next, we consider the Ag 55 nanoparticle that exhibits multiple strong peaks in the absorption spectrum, resulting in difficulties in identifying the plasmon resonance. The TCM analyses for the three prominent peak energies are shown in Figs. 3(a-c). Due to its small size, Ag 55 has well defined, discrete KS states as visible in DOS. The overall features in TCMs are similar to those of the larger nanoparticles, i.e., the low-energy sp transitions and the d electron transitions yield positive and negative contributions, respectively, though the low-energy transitions that form the plasmon are energetically clearly separated.
In contrast to the larger nanoparticles, in the Ag 55 nanoparticle some of the strongly contributing sp transitions are located close to the peak frequencies, i.e., close to the ε u − ε o = ω lines in the TCMs. These excitations are marked in Figs. 3(a-c) by green circles numbered as 1 and 2. By examining these KS transitions as a function of frequency ω (TCMs with the 0.01 eV resolution are provided as Supplemental Material 109 ), we note that the first transition changes its sign at ω = 3.85 eV, close to the minimum between the peak maxima at 3.71 eV [ Fig. 3(a)] and 4.00 eV (b). Similarly, the second transition changes its sign at ω = 4.06 eV between the maxima at 4.00 eV (b) and 4.20 eV (c). At the same time, the lowenergy transitions forming the plasmon remain mainly unchanged over this frequency window. Thus, the presence of multiple peaks in the Ag 55 spectrum seems to correspond to a strong coupling between the marked KS transitions and the plasmon. This is seen as the splitting of the plasmon into multiple resonances with antisym-metric and symmetric combinations of the KS transition and the plasmonic transitions. In the larger nanoparticles, the interaction between the plasmon and the nearby KS transitions is weak and the coupling is merely seen as a broadening of the plasmon peak.
A detailed inspection reveals that some d electron transitions also change their sign in the frequency range where the peak splitting occurs. The changes in their sign, however, do not match the maxima and minima of the absorption spectrum like in the case of the marked KS transitions. Thus, we expect the marked sp transitions to be the major cause for the plasmon splitting.
In the literature, Ag 55 has been reported to have slightly varying spectra depending, e.g., on the exact geometry, the exchange-correlation functional, and the numerical parameters used. 14,47,73,74,100,110,111 Correspondingly, the electronic structures are different and the Ag 55 spectra have single or multiple peaks. We expect, however, that the splitting behavior observed here can be a useful general concept for understanding the response of small plasmonic nanoparticles.

IV. DISCUSSION
The RT-TDDFT approach provides a more favorable scaling with the system size than the Casida approach. The latter, however, achieves a smaller pre-factor, especially when using non-local (e.g., hybrid exchangecorrelation functionals), 84 which renders it computationally more efficient for small and moderately-sized systems. In contrast, the RT-TDDFT approach becomes very attractive for systems comprising thousands of electrons (and typically hundreds of atoms) such as the silver nanoparticles considered in the present work. Previously, the lack of a decomposition scheme on par with the Casida method has been identified as a drawback of the RT-TDDFT approach. 84 Here, we have introduced and demonstrated the performance of a method that overcomes this limitation and represents an efficient tool for analyzing electronic excitations within RT-TDDFT in general, and plasmonic response in particular.
It should be noted that in the RT-TDDFT approach the observable response is sensitive to the external perturbation used to initialize the time propagation. If the perturbation is chosen to be, say, a dipole perturbation along the x direction, only excitations with a dipole component parallel to x are observable in the response. By combining at most three separate time-propagation calculations (possibly even less in the cases of higher symmetry) with dipole perturbations along the x, y, and z axes, one can recover the full dynamical polarizability tensor. However, for obtaining optically dark (dipole-forbidden) excitations from RT-TDDFT calculations, one would need to run the time propagation with different initial perturbations. This is in contrast to the Casida approach, where also dipole-forbidden excitations are obtained by diagonalizing the Ω matrix.
It was illustrated in Sec. III A that the RT-TDDFT method does not yield direct access to the discrete spectrum, but rather allows an analysis at chosen frequencies yielding the combined response coming from all the contributing discrete excitations. Usually, this is not a significant restriction as in experimental measurements the energy resolution is limited by instrumental broadening and the excitation lifetimes. Computationally, the energy resolution is determined by the broadening parameter, which can be always reduced by increasing the propagation time. Furthermore, for larger systems that are the primary application area for RT-TDDFT, the electronic spectrum becomes increasingly dense and the distinction of individual excitations is less relevant.

V. CONCLUSIONS
In this work we have presented that the linear response of the density matrix in the Kohn-Sham electron-hole basis can be obtained from real-time propagation TDDFT via a basis transformation. The methodology has been implemented in a recent RT-TDDFT code 71 and is to be made publicly available as a part of the open source electronic structure code gpaw. [93][94][95] The present approach provides access to the same information via RT-TDDFT that is usually available only with the Casida approach. This was specifically demonstrated by a careful comparison of the results for benzene derivatives, which were shown to be numerically almost identical for the Casida and RT-TDDFT calculations.
Using the presented methodology, we analyzed the plasmonic response of icosahedral silver nanoparticles in the Kohn-Sham electron-hole space. The Ag 55 nanoparticle was considered in detail and the multiple resonances in its response were shown to reflect the splitting of the plasmon due to the strong coupling between the plasmon and individual single-electron transitions. In the larger Ag 147 , Ag 309 , and Ag 561 nanoparticles, the interaction between plasmon and individual single-electron transitions close to the resonance is weaker and a distinct plasmon resonance emerges from the constructive superposition of the low-energy Kohn-Sham transitions 39,47,108 accompanied by the damping due to d electron transitions.
In summary, the present work raises the analysis capabilities of the RT-TDDFT to the same level as with the Casida approach, without compromising the computational benefits of RT-TDDFT. and Kalle Väisälä Foundation of the Finnish Academy of Science and Letters, and Finnish Cultural Foundation for support. We also thank the Swedish Research Council, the Knut and Alice Wallenberg Foundation, and the Swedish Foundation for Strategic Research for support. We acknowledge computational resources provided by CSC -IT Center for Science (Finland), the Aalto Science-IT project (Aalto University School of Science), the Swedish National Infrastructure for Computing at NSC (Linköping) and at PDC (Stockholm). where ψ n (r, t) is a pseudo wave function and T denotes the PAW transformation operator. 96 In the LCAO method, the pseudo wave function ψ n (r, t) is expanded in localized basis functions φ µ (r) centered at atomic coordinates where the all-electron basis functions have been defined as φ µ = T φ µ .
The transformation of the real-space density matrix to the basis defined by the ground-state KS orbitals ψ Here, we have isolated the overlap integrals S µµ used regularly in LCAO calculations, i.e., S µµ = drφ * µ (r)φ µ (r) = dr φ * µ (r)T † T φ µ (r).
(A.7) After simplifying the overlap integrals in Eq. (A.6), we obtain Eq. (9). We note that the PAW transformation affects only the evaluation of the overlap integrals S µµ , see Eq. (A.7).