Optical spectra of plasmonic Au clusters and nanoparticles obtained using the TDDFT+U method

Abstract

The description of localized surface-plasmon resonances (LSPRs) in gold clusters has been a long-standing problem for quantum calculations, complicated by the presence of the filled d shell close to the Fermi energy and the resulting strong influence of interband transitions in the visible spectral range. A full quantum-mechanical description that represents the atomistic structure and treats all relevant electrons at the same footing is needed. In this study, we demonstrate that a Hubbard U correction to the 5d electrons in the real-time TDDFT+U approach is able to model the optical response of plasmonic Au clusters and nanoparticles. The corrected LSPR energies show good agreement with “semi-quantal” calculations known to provide reliable plasmon energies, as well as with measurements when the matrix-induced/surface-induced red shifts are taken into account. Our findings establish the RT-TDDFT+U method as a valuable approach for modeling the optical properties of plasmonic Au clusters and nanoparticles, enabling calculations of clusters of up to 923 atoms, representing a 3 nm diameter nanoparticle.

---

1. Introduction

Noble-metal clusters and nanoparticles (NPs) exhibit outstanding optical properties driven by localized surface-plasmon resonances (LSPRs) that dominate their response from the visible to the ultraviolet region, enabling a wide range of applications. They are essential components of emerging plasmonic technology including surface-enhanced Raman spectroscopy (SERS),^1,2 biosensing,^3,4 photothermal cancer therapy,^5,6 plasmon-enhanced catalysis⁷ and many more. Thus, having theoretical means to predict their properties and size-evolution could accelerate rational design, guide experiments, and optimize device performance across these applications.

Historically, the quantum description of metal clusters and their optical properties was based on the jellium model, which neglects the atomistic structure of the clusters.^8–10 In the case of the noble metals gold, silver, and copper, the presence of the respective electronic d shells close to the Fermi surface influences the optical spectra strongly. The d-shell electrons are relatively strongly localized around the atoms as compared to the free s electrons, but their polarizability plays an essential role in the optical response of the materials, leading to the screening of the surface plasmons in noble-metal clusters.^11–14 Their localization leads to the model idea that introduces a shell of reduced screening at the cluster surfaces, which is essential for explaining the size dependence of surface-plasmon resonances.¹¹ A self-consistent “semi-quantal” model that uses this concept was elaborated by Lermé et al.,¹⁵ obtaining reliable plasmon energies for free and embedded spherical noble-metal clusters.^13,15

Today, the most widely used method for obtaining the spectra of noble-metal clusters is time-dependent density-functional theory (TDDFT) with some atomistic description of the material using mostly pseudopotentials.¹⁶ The simplest local or semi-local approximations for the XC functionals like local-density approximations (LDAs) and generalized-gradient approximations (GGAs) can fail dramatically to produce correct spectra notably due to the poor description of the d states.¹⁷ Hybrid functionals yield excellent spectra for small silver clusters but are limited in their application to sizes of about 150 to 200 atoms.^18,19 However, even they may fail to obtain a good description of small gold clusters.²⁰

Recently, we have shown that using real-time TDDFT with an empirical Hubbard-U correction (RT-TDDFT+U), one can obtain reliable and accurate spectra for both small, molecular and larger plasmonic Ag clusters.^17,21 In particular, we have demonstrated that a U of fixed value, corresponding to the one found adequate for bulk silver,²² could be applied for all the clusters studied in that work.¹⁷

Following the success of the RT-TDDFT+U method for Ag clusters across all size regimes—from molecular clusters to nanoparticles—we explore here the applicability of this approach to the closely related element Au. Despite the similarities between the valence electronic configurations of Ag and Au, their optical properties are very different. The key complication lies in the position of the d-band edge, which lies much closer to the Fermi level in Au than in Ag. This leads to strong dynamic screening of the LSPR by the d electrons and significant direct contributions to the optical response in the visible spectral range, the plasmon emerging on top of a rising background of interband transitions from d to sp states.^18,23,24

In the present work, we use the DFT+U method for the ground state and, subsequently, RT-TDDFT+U calculations to compute the optical spectra. We compare the energy and the size dependence of the LSPR with the above-mentioned “semi-quantal” jellium-based calculations, which, for spherical clusters, provide good plasmon energies.^13,15 The optical properties of larger gold clusters showing a surface plasmon in the spectra are well represented, whereas the description of the spectra of very small gold clusters with discrete electron–hole-type spectra remains poor. We demonstrate the practical applicability of the RT-TDDFT+U approach to clusters of up to 923 atoms, representing a nanoparticle of 3 nm diameter.

2. Methods

Atomistic calculations

We consider two types of Au clusters: (i) the usual series of icosahedral clusters (Au_n with n = 55–923) to investigate the size trend of the LSPR without any change of morphology and (ii) the lowest-energy quasi-molecular clusters taken from Anak et al.¹⁸

For the molecular structures, both the structural relaxation and the subsequent optical absorption calculations were carried out using the LDA (PZ81)²⁵ XC functional. For the investigation of the size dependence of the LSPR, we considered icosahedral Au clusters of five different sizes. These structures were relaxed in VASP^26,27 using the LDA (PZ81)²⁵ prior to the LDA+U optical calculations and using the GGA (PBE)²⁸ prior to the GGA+U optical calculations for comparison. Structural relaxations were performed with symmetry held fixed, using an SCF energy convergence of 10⁻⁶ eV and a force threshold of 0.02 eV Å⁻¹ per atom.

The reason for using LDA for most calculations is that it is known to yield better structural agreement for Au than GGA, as also demonstrated by López Lozano et al.²⁹—who demonstrated that for the experimentally determined structure of Au₁₄₆, LDA provides better agreement with measured Au–Au bond lengths. Moreover, the LDA has previously been used to compute the optical spectra of relatively large Au clusters,³⁰ largely due to its lower computational cost.

For computing the spectra, the real-space code Octopus was used,^31,32 with the grid spacing and the radius of the simulation set to 0.20 Å and 5.0 Å, respectively. This radius corresponds to spheres centered on each atom, whose superposition determines the final simulation domain used in the calculation. The interactions between the electrons and the ions were described using norm-conserving Troullier–Martins pseudopotentials,³³ treating 11 valence electrons corresponding to the ten 5d electrons and a single 6s electron per Au atom explicitly. The absorption spectra were calculated using the RT-TDDFT+U method, which is based on the Yabana–Bertsch time-evolution formalism,³⁴ and involves real-time propagation of the wavefunctions after a δ-kick at t = 0. We have used the approximated enforced time-reversal symmetry (AETRS) propagator, with a time step of ≈0.0016 fs (0.0024 ħ eV⁻¹) and a total propagation time of ≈13 fs (20 ħ eV⁻¹) to obtain a broadened spectrum useful for identifying the plasmonic peak; alternatively, 26 fs (40 ħ eV⁻¹) was used to improve the resolution of the spectral peaks.

Our investigation involves the use of the rotationally invariant formulation³⁵ of the DFT+U method, where U_eff = U − J is used but generally referred to as just U. It acts locally on the 5d orbitals at all Au atomic sites. We used a Hubbard-U correction of 2.0 eV for all calculations. Subsequently, the extension RT-TDDFT+U is used for the spectra.³²

Jellium calculations

In addition to the full quantum-mechanical RT-TDDFT+U calculations reported in this work, the absorption spectra of free Au_n clusters have been computed in the framework of a mixed classical/quantal jellium-type model. In this simple approach only the 6s electrons are quantum-mechanically treated, whereas the d-shell electrons are implicitly taken into account through their absorption/screening properties in superimposing a homogeneous dielectric medium on the positively charged background consisting of the Au⁺ ions. In the jellium model¹⁰ the atomistic discrete ionic background, of total positive charge Q = nq (q is the elementary charge), is replaced by a step-walled uniformly positively-charged sphere of radius R_N = r_sn^1/3, where r_s = 3.01 a₀ is the Wigner–Seitz radius per conduction electron (6s) in bulk gold (note that, in the present work, the icosahedral structures involved are indeed quasi-spherical; see Fig. 1 in the main text). For handling noble metals, such as gold, as well as matrix embedded-particles, the pioneering jellium model¹⁰ has been successfully extended in order to include polarizable dielectric backgrounds self-consistently.^13,15,36 In the present calculations, the inner homogeneous polarizable/absorbing dielectric medium extends up to R₁ = R_N − d, where d is the skin thickness of ineffective ion polarizability¹¹ (for gold the appropriate value, determined with reference to experimental results on alumina-embedded gold NPs, is d = 3.5 a₀). The inner dielectric medium is characterized by a frequency-dependent complex dielectric function ε_ib(ω) (dimensionless input data of the model) corresponding to the interband transitions, assumed to be bulk-like. The dielectric function ε_ib(ω) has been extracted from the experimental gold dielectric function using Kramers–Kronig analysis (a brief description of the procedure is provided in the SI of ref. 36). The properties of the ground state (GS), computed thanks to the DFT, are obtained from the electron density associated with the n occupied DFT-Kohn–Sham (KS) wavefunctions minimizing the energy. The absorption spectra are then computed using the time-dependent local-density-approximation (TDLDA),³⁷ appropriately generalized in order to include underlying dielectric media self-consistently.^13,15,36 For both the GS-KS calculations and the TDLDA optical response, the local Gunnarsson–Lundqvist exchange–correlation energy functional,³⁸ of widespread use in the cluster physics literature, has been selected. Moreover, the electron self-interaction correction (SIC) of the Amaldi type has been included, as in ref. 39. In this approximate SIC formulation all the occupied DFT-KS wavefunctions are treated on an equal footing, in both the Hartree and exchange–correlation contributions to the total energy of the NPs. Let us specify that the SIC applies obviously only to the 6s electrons (quantum-mechanically treated). Nevertheless, it should be emphasized that, for neutral Au_n clusters under vacuum, the optical responses computed in the jellium model, either with or without the SIC, are almost identical, even though the effective GS-KS confining potentials differ greatly. The free electron mass m_e = 1 is used in DFT-KS calculations. This value is consistent with the effective optical electron mass values determined by various authors from a careful analysis of the measured bulk gold dielectric function. An additional parameter must be briefly discussed, namely the infinitesimal δ parameter entering the Green's functions G₀(r, r′, ω) in the TDLDA formalism. In practical calculations a finite value is used, and this parameter turns out to act as an effective smoothing parameter, which may be used to mimic phenomenologically line broadening effects (essentially dissipative) that are not included in the modelling. The role played by this parameter has been extensively discussed in the previous work (see, for instance, the SI of ref. 36). It is sufficient here to say that this parameter amounts to attributing an intrinsic width 2ħδ to each bound–bound one-electron excitation line (Lorentzian-shaped curve peak).

Fig. 1 Absorption spectra of Au icosahedral clusters in vacuum, calculated using (a) TDDFT+U with a U correction of 2 eV, and (b) pure TDDFT. The inset shows the absorption spectra of Au clusters embedded in an alumina matrix, which produces an estimated red shift of 0.3 to 0.4 eV, taken from Cottancin et al.¹³ A second red shift is induced in these experiments by the size distribution of clusters and the fact that the biggest clusters will contribute most to the measurement. A rough agreement of the pure-TDDFT spectra with these measurements shows that the calculated energies are substantially too low. This error is corrected in the TDDFT+U calculations; the shift is of the order of 0.4 to 0.5 eV. In addition, the total intensity of the plasmon peaks increases. The LDA (PZ81)²⁵ XC functional was used for the relaxation of the structures (without the U correction) as well as in the spectral calculations. The curves are obtained using an evolution time of 13 fs after the δ-kick.

3. Results

Optical spectra

In Fig. 1, we present the calculated absorption spectra of the series of icosahedral gold clusters Au_n (where n = 55–923), using TDDFT with the U correction (Fig. 1(a)) and without (Fig. 1(b)). A U value of 2 eV was selected based on the findings of Avakyan et al., who reported good agreement with the experimental dielectric function of bulk gold using this value.²²

Using the TDDFT+U method, we find that the spectra of Au₅₅ and Au₁₄₇ show no clear LSPR signature, in agreement with the finding that the LSPR emerges for sizes above 150 atoms.^40,41 In contrast, a well-defined LSPR peak appears for Au₃₀₉, emerging from the dense interband background around 2.9 eV (top panel of Fig. 1). One might argue that a rather flat damped bump at around 3.1 eV in Au₁₄₇ could represent an LSPR, but its weak oscillator strength likely renders it unobservable under experimental conditions. For the larger two clusters, Au₅₆₁ and Au₉₂₃, the LSPR becomes more pronounced and red-shifts to approximately 2.8 and 2.7 eV, respectively, following the expected trend with increasing size.

With “pure” RT-TDDFT without the Hubbard correction, here using the LDA functional, the overall spectra are noticeably red-shifted relative to TDDFT+U results. The distinct LSPR peaks appear at lower energies and with a significantly weaker intensity relative to the interband contributions, clearly visible for the larger clusters (Au₃₀₉, Au₅₆₁ and Au₉₂₃). Au₅₅ and Au₁₄₇ show no LSPR signature, while weak resonances begin to emerge at approximately 2.6, 2.4 and 2.3 eV for Au₃₀₉, Au₅₆₁ and Au₉₂₃, respectively (see Fig. 1(b)). The weakness of the resonances in the calculations without the Hubbard-U correction, as compared to the U-corrected results, clearly indicates the occurrence of stronger LSPR coupling (resulting in strong damping) with interband transitions, which can be directly attributed to the d-band edge lying closer to the Fermi level than it should.

Validation of TDDFT+U results

To the best of our knowledge, there are no experimental reports of plasmonic behavior in pure Au clusters studied via free-beam experiments under vacuum, unlike in the case of Ag clusters. The vast majority of optical absorption measurements have been performed either in solid matrices such as alumina or silica, in solution, or on soft-landed size-selected particles on thin carbon films—all of which can significantly influence the LSPR energy and intensity.^13,42–44

For Au particles embedded in an alumina matrix, Cottancin et al.¹³ observed the emergence of an LSPR peak between 2.3 and 2.4 eV for particles larger than 2 nm. These data are reproduced in the inset of the lower panel of Fig. 1. Interestingly, these energies roughly correspond to those of our pure DFT results without the Hubbard correction. However, this only shows that the energies of these calculations of clusters under vacuum are too low, similar to the result of Ag clusters¹⁷ since embedding matrices are known to red-shift the LSPR energy. An estimate of the matrix-induced red shift of approximately 0.3 to 0.4 eV is provided by the “semi-quantal”, jellium-based results of Lermé.¹⁵ Such a redshift brings our TDDFT+U calculations into good qualitative agreement with the experimental results. Moreover, the matrix-induced red shift reduces the overlap with interband transitions from the d band, thereby decreasing damping. This reduced damping should allow the LSPR to emerge more readily in smaller particles than would be possible under vacuum.

In the absence of experimental results under vacuum, we compare our TDDFT+U results with the semi-quantal model¹⁵ that has been successfully applied in previous studies^13,36 in the cases of free and embedded noble metal particles (see, in particular, the short description of the modeling in the SI of ref. 36). In these structureless jellium calculations, the screening of the LSPR by the d-shell electrons is ensured using a complex interband-related dielectric function that has been extracted from experimental data on bulk gold through Kramers–Kronig analysis. It should be emphasized that this semi-quantal modeling includes both the interband absorption and screening/polarization effects self-consistently. In addition, a shell of reduced d-electron-induced screening at the cluster surfaces is introduced,¹¹ the appropriate thickness of which has been determined with reference to experimental results. Thus, even while neglecting the atomistic structure of the clusters, these calculations obtain reliable plasmon energies for free and embedded spherical noble-metal clusters,^13,15 and they can likewise produce the size-dependence of plasmons in different environments well.³⁶ Owing to the success of the model calculations in reproducing plasmon energies and size effects in plasmonic clusters, we use these calculations to validate our atomistic TDDFT+U calculations.

In Fig. 2, we present a comparison between TDDFT+U calculations (the same as those presented in Fig. 1 upper panel) and the semi-quantal calculations. Both calculations produce similar energies and a clear red shift with an increase in cluster size. In general, as expected, the plasmon emerges out of the background of interband transitions, which is substantial in the spectra, and the assignment of the plasmon peak is tricky. In particular for the Au₅₆₁ and Au₉₂₃ clusters, the most prominent peak clearly consists of two main peaks, of similar energy between the two approaches. While this corroborates the values of the plasmon energies calculated using TDDFT+U, we need to point out the differences which remain, probably mostly due to the simplifying treatment of the d-shell electrons and of the positively-charged ionic background in the jellium approach. First, the overall intensity of the semi-quantal calculations is higher. Second, the energy locations of the bumps over the entire spectra, including those underlying the fragmented pattern of the LSPR resulting from the coupling of the collective excitation with electron–hole transitions (Landau damping mechanism), do not coincide in general. This is not surprising since these bumps are tightly related to the independent-electron energy spectrum of the GS-KS potential. Indeed, it is clear that atomistic and homogeneous jellium background descriptions will result in slightly different GS-KS energy spectra, and thus in their respective bump patterns in the optical responses.

Fig. 2 Comparison between (a) atomistic TDDFT+U (U = 2 eV) calculations using LDA (PZ81)²⁵ to approximate the XC functional, and (b) jellium calculations for Au clusters of five sizes (see legend). In (a), thick curves correspond to a propagation time of 13 fs, while thin, better-resolved curves correspond to 26 fs. In (b), broadened curves correspond to ħδ = 120 meV, whereas thin, better-resolved curves use ħδ = 60 meV.

Moreover, the LSPR emerges much earlier and more strongly for the smallest clusters than in the TDDFT+U results. For Au₅₅ and Au₁₄₇ clusters, the TDDFT+U calculations show no clear sign of the LSPR, but the jellium model calculations do. It is likely that some discrepancies could be reduced by relaxing the drastic bulk-like assumption concerning the interband dielectric function. Obviously, one could adjust the dielectric function to the smallest sizes, but developing an analytical, physically sound, size-dependent dielectric function converging toward bulk values at large sizes, while highly desirable, would be difficult. Ultimately, a challenge would be to introduce, within the jellium approach, an interband dielectric function involving some relevant, physically sound quantum-mechanical ingredients.

In Fig. 3, we compare our TDDFT, TDDFT+U, and jellium-model results with available measurements. Assigning a single LSPR energy to Au clusters is inherently challenging because of the strong coupling of the LSPR to the interband transitions that leads to fragmentation of the collective excitation. In the TDDFT+U (and TDDFT) calculations, we therefore identify the LSPR energy as the position of the first maximum in the absorption spectrum obtained with an evolution time of 13 fs. For jellium-model calculations, we define the LSPR energy as the centroid of the first two main peaks present in the curve (ħδ = 60 meV). The well-foundedness of these choices is corroborated in the following section, where we indicate the possibility of identifying the LSPR by modifying the parameters of the calculations.

Fig. 3 LSPR energies plotted as a function of inverse particle diameter. The TDDFT, TDDFT+U, and jellium-model results are shown as red, blue, and orange crosses, respectively. Experimental data are shown as squares: the black squares correspond to the measurements of Cottancin et al.¹³ and are shifted upward by 0.38 eV to correct for the alumina matrix, unshifted points are presented with thin dashed squares; the green squares correspond to the data of Lu et al.⁴⁴ The latter values are likewise expected to exhibit a red shift, but its magnitude cannot be reliably estimated due to the complex experimental environment, so the measured values are plotted as reported. The apparent agreement between the TDDFT results without the Hubbard correction (red crosses) and the measured points (green squares) is fortuitous, as the calculations assume perfect vacuum whereas the experimental measurements will include the influence of the substrate surface.

For comparison with experiment, we have used two types of experimental measurements: (a) size-selected Au clusters embedded in an alumina matrix,¹³ and (b) clusters soft-landed on a 3 nm-thick carbon film.⁴⁴ The measurements performed in the alumina matrix¹³ are rigidly shifted by 0.38 eV to account for the matrix-induced red shift. This value is based on the findings of Lermé et al. (EPJD, 2000),¹⁵ where the difference between the LSPR energies calculated using the jellium model for an Au₈₃₂ cluster with and without the matrix was found to be approximately 0.38 eV. Compared to these shifted values, the TDDFT+U results are in relatively good agreement, and certainly better than the TDDFT calculations without the Hubbard correction.

For the measurements on size-selected Au clusters deposited on a carbon film (green squares) by Lu et al.,⁴⁴ a red shift in the LSPR energy is also expected relative to clusters of the same size under ideal vacuum conditions, but the precise estimation of the red shift is complicated by the complex local environment. In addition, the systems may themselves evolve during the EELS measurement due to the effect of the electron beam, which may dynamically change the plasmon energies.³⁶ Therefore, the experimental points are taken directly as measured, and appear to be in good agreement with the TDDFT results without the Hubbard correction. This agreement with the TDDFT results is spurious, as the calculated and measured values correspond to two different systems: the TDDFT results refer to clusters/nanoparticles under perfect vacuum, whereas in the experiments the clusters are in contact with the carbon film, while a substantial portion remains exposed to vacuum-like conditions. The partial vacuum-like environment may influence the size-dependence of the LSPR in a similar way the porosity of the matrix does in matrix embedded experiments, as discussed by Campos et al.³⁶ Additionally, if the values measured by Lu et al. are extrapolated to the large-size limit, the resulting LSPR energy is unrealistically low.

In conclusion, the TDDFT+U spectra provide a corrected description compared to pure TDDFT with simple functionals, and the resulting LSPR energies are consistent with measurements on matrix-embedded clusters as well as those on clusters confined to surfaces, provided that a finite matrix/surface-induced red shift is taken into account. For measurements performed on surface-supported clusters, however, the complexity of the local environment prevents a reliable estimation of this red shift, making a direct comparison with the calculated results difficult. Finally, the validity of the computed spectra is further reinforced by the agreement with the ‘semi-quantal’ calculations, which are known to reproduce plasmon energies accurately, even though the detailed spectral shapes differ to some extent.

Analysis of the results: PDOS

For Au in both bulk and thin-film forms, the d-band edge is experimentally known to lie approximately 2 eV below the Fermi level.^45–47 This holds for clusters as well, as supported by the work of Taylor et al.,⁴⁸ who found an increase in the energy onset of both the 6s and 5d bands with a cluster size, while the difference between the two onsets, assuming a 6s onset as the Fermi level, remains approximately the same for medium to large clusters.

The blue-shift of the LSPR upon applying the U correction is due to increased localization of the d electrons (see Fig. 4(b)), which reduces the dynamic screening of the collective oscillations of the sp electrons responsible for the LSPR. The enhanced oscillator strength may result from the upward shift of interband transition energies, thereby decreasing coupling between LSPR and the transitions involving d states.

Fig. 4 (a) PDOS for the Au₅₅ cluster with and without an effective Hubbard-U correction of 2 eV. (b) Ground state total (5d + 6s) charge-density difference between DFT+U (2 eV) and DFT for Au₅₅, represented using a colour-coded slab passing through the central atom, confirming increased localization of the d electrons around atoms.

In Fig. 4(a), we present the projected density of states (PDOS) for the Au₅₅ cluster calculated using DFT (top panel), where the d-band edge appears approximately between 1.5 and 2.0 eV below the Fermi level. This result highlights the electron delocalization error inherent in local functionals, which stems from the insufficient cancellation of the spurious coulombic self-interaction. This issue, previously discussed in the context of Ag clusters,¹⁷ also affects Au.

As shown in the bottom panel of Fig. 4(a), the application of the U-correction with a value of 2 eV shifts the d band to a more reasonable position, approximately 2 eV below the Fermi level, bringing it closer to the experimental results.

Enhanced d-electron localization under the Hubbard-U correction is confirmed by the ground-state charge-density difference plot, Δρ(r) = ρ_DFT+U(r) − ρ_DFT(r), as shown in Fig. 4(b). The colour-coded slab visualization shows an increase in electron density around the Au atoms, exhibiting a distinct characteristic shape of the d orbitals, while in the intermediate region between the atoms, a decrease in the electron density is clearly visible. This increased localization is essential to correct the dynamic screening of the LSPR, an effect which is overestimated when using only DFT. Thus, the U correction brings the calculated ground state closer to the expected real ground state of the system.

Identifying the LSPR

Unlike in the case of silver, where identifying the LSPR peak in the calculated spectra is relatively straightforward due to the strong contrast between the plasmonic peak and other spectral features, this task is considerably more challenging for Au clusters. The electron–hole-like interband transitions give rise to numerous peaks, some of which may be mistakenly interpreted as an LSPR. Additionally, it should be stressed that because of the Landau damping mechanism, which results in fragmentation of the LSPR band, several bumps may indeed underlie the LSPR band when a large propagation time (TDDFT calculations) or equivalently a small ħδ parameter (jellium calculations) is used (see Fig. 2). A key advantage of the TDDFT+U and the jellium calculations is that they allow direct tuning of the d-electron influence on the LSPR. By employing a higher U value, the intensity of the LSPR peak can be artificially enhanced, enabling a controlled study through a series of calculations with progressively smaller U values to track the evolution of the plasmonic response. Similarly, in the jellium calculations the thickness d of the surface shell of ineffective interband screening can be used to enhance the plasmonic peak (the d parameter is known to phenomenologically mimic the degree of increased localization of the d-electron orbitals).^11,15,36

The two approaches are illustrated in Fig. 5(a) and (b), where we performed a series of TDDFT+U calculations with U values ranging from 0 eV to 3 eV, and jellium calculations with an ineffective-screening shell of different thicknesses. With an increase in U, the overall screening is reduced and the LSPR peak exhibits a clear and consistent trend—shifting to higher energies and becoming more pronounced. From this analysis, we conclude that the peak observed at around 2.9 eV in the TDDFT+U calculations in Fig. 1 for Ag₃₀₉ indeed corresponds to the LSPR. This also holds for the peak that appears in the pure TDDFT spectrum at a (too low) energy of 2.6 eV. In the jellium calculations, lower values of thickness d increase the overall screening of the LSPR and hence shift it to lower energies, away from the interband transitions.

Fig. 5 (a) Variation in the TDDFT+U absorption spectra of Au₃₀₉ as a function of the U parameter (0 to 3 eV), with the LDA (PZ81)²⁵ to approximate the XC functional (time evolution of 13 fs) and (b) variation in the absorption spectra calculated using the jellium model, with a reduced screening layer of different thicknesses (0 to 3.5 Bohr), as indicated in the legend (ħδ = 60 meV).

The procedures used (in TDDFT+U and jellium calculations) in Fig. 5 clearly show that, in all cases, the LSPR lies at the lower-energy edge of the absorption spectra, consistent with how it was assigned in the results discussed earlier in Fig. 1 and 2.

Influence of the XC functional

Until now, we have shown that TDDFT+U calculations using the LDA for the XC functional and a U value of 2 eV yield good agreement with the expected position of the LSPR in Au clusters for the relevant sizes.

The reason for this choice was that the LDA is known to yield better structural agreement with the experiment for Au than, for instance, the GGA. This well-known fact was corroborated for gold clusters by López-Lozano et al.,²⁹ who demonstrated that for the experimentally determined structure of Au₁₄₆, the LDA provides very good agreement with measured Au–Au bond lengths. Moreover, LDA has previously been used to compute the optical spectra of relatively large Au clusters,³⁰ largely due to its low computational cost. However, the question of how the choice of the XC functional influences the TDDFT+U calculations must be addressed. For that reason, we studied the performance of the GGA functional (PBE)²⁸ in the TDDFT+U calculations.

In this section, when we write ‘LDA+U’ (or, respectively, ‘GGA+U’), the LDA (or, respectively, GGA) XC functional is used in both steps: in the ground-state calculation and in the subsequent TDDFT+U calculation.

In the first two panels of Fig. 6, we show the calculated absorption spectra of the Au icosahedral cluster series with the LDA (the same as that in Fig. 1(a)) and the GGA, obtained using TDDFT+U with a fixed U correction of 2 eV. For the GGA+U calculation, the LSPR emerges again for the correct size and red-shifts with an increase in size but the energies are consistently higher than what we have obtained with the LDA shown in Fig. 6(a). Additionally the integrated intensity under the LSPR peak (area under the curve in the visible region) is also larger in the latter case.

Fig. 6 Influence of the XC functional used in the TDDFT+U calculations. The panels (a) and (b) show absorption spectra computed using the TDDFT+U method using LDA (PZ81)²⁵ and GGA (PBE),²⁸ with the same XC functional also used to relax the structures. All the calculations shown in the first two panels are consistent, i.e., the structures were relaxed using the same functional, the GGA or the LDA, that was used in the TDDFT+U calculations. Thick colored curves are obtained with an evolution time of 13 fs, while the thin black curves overlaid on them correspond to a longer evolution time of 26 fs, showing a higher resolution of the same spectra. To distinguish the influence of the structural differences in the calculations from the effect of using the different functionals in the TDDFT+U part, in panel (c) we present a comparison of the TDDFT+U spectra computed with the LDA and GGA for Au₃₀₉, using the structure relaxed with the LDA. An effective-U correction of 2 eV was used in all the calculations shown.

The difference between the two can be attributed to two factors influenced by the choice of the functional: first, structural differences, as the GGA is known to yield longer Au–Au bond lengths, which can affect the optical properties; and second, the functional's direct impact on electronic interactions, which also modifies the optical response. To isolate the influence of structural changes, we performed a TDDFT+U calculation using the GGA on an LDA-relaxed structure and found that the LSPR energy remained higher than in the LDA-based calculation for the Au₃₀₉ cluster (see Fig. 6(c)). This confirms that the trend of higher LSPR energies with the GGA is only weakly dependent on how the structure is relaxed.

It is important to note that the U correction is used empirically in our work and, in principle, can also depend on the choice of the functional. Therefore, using a smaller U value with the GGA may lead to results that agree with those obtained using LDA+U. It may be helpful in future investigations to determine U from first principles, as this could provide a more coherent combination with the chosen functional.

The second key difference, namely the larger area under the absorption curve in the visible range for the GGA+U calculation, can be explained by the density of states (DOS) plot. In Fig. 7, we show the DOS—specifically the superatomic projected DOS (SAC-PDOS)—for the Au₃₀₉ cluster calculated using the LDA, GGA, LDA+U, and GGA+U. The d band is readily identifiable as a dense region of states, starting below approximately −1.5 eV in the DFT calculations, and shifting further below −2.0 eV with the inclusion of the U correction, displaying no distinct superatomic character (mixture of many colours indicating no single angular momentum).

Fig. 7 Superatomic projected density of states (SAC-PDOS) of the Au₃₀₉ cluster, calculated using the LDA, GGA, LDA+U, and GGA+U, with a Hubbard-U correction of 2 eV. The LDA and LDA+U structures were relaxed using the LDA, while the GGA and GGA+U structures were relaxed using GGA. The projection onto spherical harmonics was carried out up to L = 10, ensuring that the superatomic states near the HOMO with higher angular momentum character are appropriately captured. The atomic 5d band lies below approximately −1.5 eV in the DFT calculations, and below −2.0 eV when the U correction is applied, showing no superatomic character, as expected. The lowest-lying delocalized states, originating from the atomic 6s electron, exhibit clear superatomic character, following a sequence of S, P, S, D, … (see ref. 49).

Upon closer inspection, it becomes evident that the d band is significantly narrower in the GGA and GGA+U calculations compared to the LDA and LDA+U cases. In the context of absorption spectra, this implies that the GGA-based calculations include a greater contribution from d states in the energy region (LSPR) of our interest in the absorption spectra, explaining the larger area under the absorption curves in the visible (enhanced magnitude and broadening of the LSPR).

Superatomic states

One characteristic property of metal clusters, often invoked to explain their structural stability, is the presence of superatomic states.^50,51 The superatomic projected density of states for the Au₃₀₉ cluster shown in Fig. 7 shows that the Hubbard correction primarily affects the d states, while the superatomic character near the HOMO remains essentially unchanged. A notable effect is that the U-induced downward shift of the d band onset exposes additional superatomic states, previously concealed and possibly hybridized within the d band, which may influence the LSPR response. Confirming this influence would require further analysis, such as transition maps using the Kohn–Sham decomposition method,⁵² which is beyond the scope of this work.

Below the d band, 6s-derived states are also present, forming deeper-lying SAC states, following the characteristic superatomic filling sequence 1s², 1p⁶, 2s², 1d¹⁰, and so on.⁴⁹ These SAC states appear at much lower energies relative to the HOMO when the LDA is used compared to the GGA. The application of the Hubbard-U correction results in partial overlap between the lower edge of the d band and some 6s-derived SAC states, though the lowest SAC states remain intact. As these states lie far below the HOMO, they are unlikely to contribute to the plasmonic response in the visible region.

Spectra of very small clusters

Up to this point, we have demonstrated that applying a DFT+U correction to the 5d electrons in Au clusters allows for tuning of the screening and the coupling between interband d transitions and the LSPR, thereby enhancing its visibility and improving the resonance position. We now turn our attention to very small clusters, where the optical response is dominated by discrete electron–hole transitions. Structural models for these calculations were taken from Anak et al.,¹⁸ and the geometries were relaxed using LDA before performing the optical-response calculations.

In Fig. 8, we present a comparison between the spectra calculated using the TDDFT+U method with a Hubbard-U correction of 2 eV, with the LDA for the XC functional, and the experimental measurements performed on size-selected small Au clusters embedded in a neon matrix at a temperature of 7 K, as reported by Lecoultre et al.⁵³ From all the plots, it is evident that the agreement between theory and experiment is poor, and a simple Hubbard U correction is insufficient to capture all the features observed in the measurements, unlike in the cases of larger Au clusters and Ag clusters of all sizes.¹⁷ This disagreement is not surprising since many high level calculations have not been able to reproduce the experimental measurements until now.^18,19,54

Fig. 8 Comparison between the experimentally measured spectra of size-selected clusters in a neon matrix⁵³ and the calculated response using the TDDFT+U method with a Hubbard-U correction of 2 eV, employing the LDA to approximate the XC functional.

The large disagreement between theory and experiment for very small clusters could be rooted, to a large extent, in the assumption of a fixed value of the U correction applied to all the atoms (in very small clusters, all atoms probe the NP surface). In principle, however, the appropriate correction required to recover the correct ground state may depend on both the cluster size and the specific atomic sites. In the case of Ag clusters, the relevant transitions in the energy region of interest (LSPR) primarily involve 5s-derived states, with the 4d electrons contributing mainly through screening. Thus, shifting the interband d transitions to higher energies led to good agreement with experiment, as the GGA already provides a reasonable description of the s electrons.¹⁷ In contrast, for Au clusters, the d band not only contributes to screening but also directly participates through d-interband transitions within the relevant energy range (visible range). As a result, a simple energy shift may not be sufficient to achieve agreement with experiment for very small clusters, which could be another reason for the limitations of this approach.

Furthermore, relativistic effects are known to be significant in small Au clusters.⁵⁵ In our calculations, we employed a norm-conserving Troullier–Martins pseudopotential,³³ which includes relativistic corrections within the scalar-relativistic approximation. However, spin–orbit coupling, which is known to significantly affect the spectra, was not included.⁵⁵ At the same time, Anak et al.¹⁸ showed that spin–orbit coupling becomes less significant for Au clusters containing 20 atoms or more. Therefore, neglecting it is justified to a certain degree for the larger clusters, particularly those large enough to support localized surface-plasmon resonances.

4. Conclusion

In this study, we investigated the ability of the TDDFT+U approach within the real-time formalism to model the optical response of Au clusters and nanoparticles, building upon its proven success in the case of Ag.

Our results show that the inclusion of a Hubbard U correction applied to the 5d electrons, combined with the LDA for the XC functional, significantly improves the position of the LSPR in Au clusters. The corrected LSPR energies show good agreement with “semi-quantal” calculations known to provide reliable plasmon energies, as well as with measurements of Au clusters embedded in alumina, assuming a plausible matrix-induced red shift of 0.3 to 0.4 eV. Additionally, the oscillator strength of the LSPR in the large quasi-spherical clusters increases compared to a pure TDDFT calculation result. The U correction effectively shifts the d band to lower energies. This shift reduces its coupling with the LSPR and, consequently, improves both the energetic position and the intensity of the LSPR.

However, our investigation also highlights the limitations of this approach when applied to very small clusters, where the optical response is dominated by discrete electron–hole transitions of a finite system. For these spectra, a simple U correction is not sufficient to produce good agreement with experimental results. The possible reasons include the increased importance of relativistic effects such as spin–orbit coupling that are not captured by the scalar-relativistic pseudopotentials used, and the possible need for size- and site-dependent U parameters to describe individual d states correctly, not just the overall screening.

In addition, our comparative study between LDA and GGA functionals revealed that the GGA consistently provided higher LSPR energies than the LDA for an equal effective Hubbard U correction of 2 eV. The DOS further suggested that the GGA results in a narrower d band, resulting in more d interband transitions contributing to the absorption spectra in the visible region.

Finally, the analysis of superatomic states confirms that the DFT+U correction preserves the overall structure of the delocalized states near the Fermi level while exposing previously obscured states within the d band, potentially influencing the plasmonic response.

In conclusion, our findings establish the TDDFT+U method and, in particular, its applications within the real-time formalism as a valuable tool for modeling the optical properties of plasmonic Au clusters and nanoparticles. We have, here, an atomistic quantum description including the d electrons appropriately, at a numerical cost that is only slightly larger than that pure DFT calculations. Nonetheless, our findings underscore the need for more sophisticated treatments—such as inclusion of spin–orbit effects and refined U parametrization—for accurately describing the spectra of small, non-plasmonic clusters. Future work will aim to integrate these features to improve agreement with experimental observations.

Author contributions

HCW and MC devised the study. MC carried out most of the (TD)DFT+U calculations. JL carried out the jellium-based calculations. MC and HCW prepared the manuscript, which was discussed and finalized by all authors. HCW supervised the project.

Conflicts of interest

The authors declare no competing interests.

Data availability

No relevant data beyond the ones presented in the figures of the manuscript are associated with this work. All data will be made available from the corresponding authors upon reasonable request.

Acknowledgements

The authors thank Matthias Hillenkamp, Gérard, Franck Rabilloud, and Lucia Reining for enlightening discussions. We further thank Franck Rabilloud for providing the geometries of the small clusters. We acknowledge support from the French National Research Agency (Agence Nationale de Recherche, ANR) in the frame of the projects “SchNAPSS” (ANR-21-CE09-0021) and “GNOME” (ANR-24-CE09-2403). The work has used HPC resources from GENCI-IDRIS (Grant 2023-0906829). Moreover, the authors acknowledge the contribution of the International Research Network IRN Nanoalloys (CNRS). Mohit Chaudhary thanks ED352 of Aix-Marseille University for the PhD scholarship.

References

1. K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan, R. R. Dasari and M. S. Feld, Phys. Rev. Lett., 1997, 78, 1667–1670.
2. S. Nie and S. R. Emory, Science, 1997, 275, 1102–1106.
3. R. Elghanian, J. J. Storhoff, R. C. Mucic, R. L. Letsinger and C. A. Mirkin, Science, 1997, 277, 1078–1081.
4. J. R. Mejía-Salazar and O. N. J. Oliveira, Chem. Rev., 2018, 118, 10617–10625.
5. X. Huang, I. H. El-Sayed, W. Qian and M. A. El-Sayed, J. Am. Chem. Soc., 2006, 128, 2115–2120.
6. C. Li, Y. Zhang, Z. Li, E. Mei, J. Lin, F. Li, C. Chen, X. Qing, L. Hou, L. Xiong, H. Hao, Y. Yang and P. Huang, Adv. Mater., 2018, 30, 1706150.
7. S. Linic, P. Christopher and D. B. Ingram, Nat. Mater., 2011, 10, 911–921.
8. W. Ekardt, Phys. Rev. B:Condens. Matter Mater. Phys., 1984, 29, 1558–1564.
9. W. Ekardt, Phys. Rev. Lett., 1984, 52, 1925–1928.
10. M. Brack, Rev. Mod. Phys., 1993, 65, 677–732.
11. A. Liebsch, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 11317–11328.
12. L. Serra and A. Rubio, Phys. Rev. Lett., 1997, 78, 1428–1431.
13. E. Cottancin, G. Celep, J. Lermé, M. Pellarin, J. R. Huntzinger, J. L. Vialle and M. Broyer, Theor. Chem. Acc., 2006, 116, 514–523.
14. H.-C. Weissker and X. López-Lozano, Phys. Chem. Chem. Phys., 2015, 17, 28379–28386.
15. J. Lermé, Eur. Phys. J. D, 2000, 10, 265–277.
16. E. Luppi, H.-C. Weissker, S. Bottaro, F. Sottile, V. Veniard, L. Reining and G. Onida, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 78, 245124.
17. M. Chaudhary and H.-C. Weissker, Nat. Commun., 2024, 15, 9225.
18. B. Anak, M. Bencharif and F. Rabilloud, RSC Adv., 2014, 4, 13001–13011.
19. J. V. Koppen, M. Hapka, M. M. Szczȩśniak and G. Chałasiñski, J. Chem. Phys., 2012, 137, 114302.
20. C. Yu, W. Harbich, L. Sementa, L. Ghiringhelli, E. Aprá, M. Stener, A. Fortunelli and H. Brune, J. Chem. Phys., 2017, 147, 074301.
21. I. Smirnov, Z. Kaszkur, M. Chaudhary, H.-Ch. Weissker, R. Ferrando and P.J. Kulesza, Vacancy-driven twinning in metal nanoparticles: from bulk morphology transformation to optical and electrochemical effects, Nanoscale, 2026 10.1039/D5NR04865F.
22. L. Avakyan, V. Durimanov, D. Nemesh, V. Srabionyan, J. Ihlemann and L. Bugaev, Opt. Mater., 2020, 109, 110264.
23. G.-T. Bae and C. M. Aikens, J. Phys. Chem. C, 2012, 116, 10356–10367.
24. N. Durante, A. Fortunelli, M. Broyer and M. Stener, J. Phys. Chem. C, 2011, 115, 6277–6282.
25. J. P. Perdew and A. Zunger, Phys. Rev. B:Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
26. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
27. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
28. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
29. X. López-Lozano, G. Plascencia-Villa, G. Calero, R. L. Whetten and H.-C. Weissker, Nanoscale, 2017, 9, 18629–18634.
30. K. Iida, M. Noda, K. Ishimura and K. Nobusada, J. Phys. Chem. A, 2014, 118, 11317–11322.
31. N. Tancogne-Dejean, M. J. T. Oliveira and A. Rubio, Phys. Rev. B, 2017, 96, 245133.
32. N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. Le Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, U. De Giovannini, A. Delgado, F. G. Eich, J. Flick, G. Gil, A. Gomez, N. Helbig, H. Hübener, R. Jestädt, J. Jornet-Somoza, A. H. Larsen, I. V. Lebedeva, M. Lüders, M. A. L. Marques, S. T. Ohlmann, S. Pipolo, M. Rampp, C. A. Rozzi, D. A. Strubbe, S. A. Sato, C. Schäfer, I. Theophilou, A. Welden and A. Rubio, J. Chem. Phys., 2020, 152, 124119.
33. N. Troullier and J. L. Martins, Phys. Rev. B:Condens. Matter Mater. Phys., 1991, 43, 1993.
34. K. Yabana and G. F. Bertsch, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 4484–4487.
35. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
36. A. Campos, N. Troc, E. Cottancin, M. Pellarin, H.-C. Weissker, J. Lermé, M. Kociak and M. Hillenkamp, Nat. Phys., 2019, 15, 275–280.
37. A. Zangwill and P. Soven, Phys. Rev. A, 1980, 21, 1561–1572.
38. O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B, 1976, 13, 4274–4298.
39. C. Yannouleas and U. Landman, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 8376–8387.
40. H. Qian, Y. Zhu and R. Jin, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 696–700.
41. S. Malola, L. Lehtovaara, J. Enkovaara and H. Hakkinen, ACS Nano, 2013, 7, 10263–10270.
42. Y. Negishi, T. Nakazaki, S. Malola, S. Takano, Y. Niihori, W. Kurashige, S. Yamazoe, T. Tsukuda and H. Häkkinen, J. Am. Chem. Soc., 2015, 137, 1206–1212.
43. M. M. Alvarez, J. T. Khoury, T. G. Schaaff, M. N. Shafigullin, I. Vezmar and R. L. Whetten, J. Phys. Chem. B, 1997, 101, 3706–3712.
44. S. Lu, L. Xie, K. Lai, R. Chen, L. Cao, K. Hu, X. Wang, J. Han, X. Wan, J. Wan, Q. Dai, F. Song, J. He, J. Dai, J. Chen, Z. Wang and G. Wang, Natl. Sci. Rev., 2020, 8, nwaa282.
45. A. Sekiyama, J. Yamaguchi, A. Higashiya, M. Obara, H. Sugiyama, M. Kimura, S. Suga, S. Imada, I. Nekrasov and M. Yabashi, et al., New J. Phys., 2010, 12, 043045.
46. P. B. Johnson and R.-W. Christy, Phys. Rev. B, 1972, 6, 4370.
47. P. Sheverdyaeva, R. Requist, P. Moras, S. Mahatha, M. Papagno, L. Ferrari, E. Tosatti and C. Carbone, Phys. Rev. B, 2016, 93, 035113.
48. K. Taylor, C. Pettiette-Hall, O. Cheshnovsky and R. Smalley, J. Chem. Phys., 1992, 96, 3319–3329.
49. H. Häkkinen, Chem. Soc. Rev., 2008, 37, 1847–1859.
50. M. Walter, J. Akola, O. Lopez-Acevedo, P. D. Jadzinsky, G. Calero, C. J. Ackerson, R. L. Whetten, H. Grönbeck and H. Häkkinen, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 9157–9162.
51. P. Jena and Q. Sun, Chem. Rev., 2018, 118, 5755–5870.
52. T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen and P. Erhart, J. Chem. Theory Comput., 2017, 13, 4779–4790.
53. S. Lecoultre, A. Rydlo, C. Félix, J. Buttet, S. Gilb and W. Harbich, J. Chem. Phys., 2011, 134, 074302.
54. C. Yu, W. Harbich, L. Sementa, L. Ghiringhelli, E. Aprá, M. Stener, A. Fortunelli and H. Brune, J. Chem. Phys., 2017, 147, 074301.
55. A. Castro, M. A. L. Marques, A. H. Romero, M. J. T. Oliveira and A. Rubio, J. Chem. Phys., 2008, 129, 144110.

---