Efficient modal analysis of plasmonic nanoparticles: from retardation to nonclassical regimes

2.1 Incorporation of retardation and nonclassical effects

Consider a plasmonic system (Figure 1(A) and (B)), which is driven by an incident electric field E _in. We incorporate the nonclassical effects of electrons by introducing a surface polarization P ^Q source on the background side of the particle boundary ∂Ω [16, 25, 31]:

Here, d ̄ ̄ ≡ d ⊥ n ̂ n ̂ + d ∥ ( I ̄ ̄ − n ̂ n ̂ ) . d _⊥,∥ are the so-called Feibelman d-parameters, which are the centroids of the induced charge and the normal derivative of the tangential polarization current, respectively. The use of the Feibelman d-parameters allows us to elegantly include three major nonclassical effects, that is, surface-enabled intraband electron–hole excitations through I m ( d ̄ ̄ ) , nonlocality of permittivity and electron spill-out through R e ( d ̄ ̄ ) . On the other hand, quantum tunneling [24] and size quantization [34] effects are outside the applicability of the d-parameters, which, nevertheless, only become important when feature sizes are below 1 nm.

The optical response of the system is represented by expanding the total electric field E _tot with the modal fields of the electrostatic SPRs [E _n ≡ − ∇ ϕ _n , see Eq. (1)]:

We emphasize that this expansion is approximate and only includes the curl-free component of E _tot, i.e., assuming that ∇ × E _tot = 0. Nevertheless, as shall be validated below and also has been concluded in Reference [33], this approximation gives rather accurate results for particle size up to hundreds of nanometers as long as the retardation corrections are properly added.

We derive the modal coefficients |a⟩, | a 〉 ≡ a 1 ; a 2 ; … , by plugging the modal expansion Eq. (3) into the retarded Maxwell’s wave equation that includes the nonclassical surface polarization P ^Q . This treatment automatically includes both the retardation and nonclassical effects. Next, we express the electric field E _tot with the Green’s function surface integral formalism, concentrating on the fields on the particle boundary ∂Ω (see Supplementary Section S1). Finally, the equation that determines |a⟩ is derived as (see Supplementary Section S2):

Here, |S⟩ ( | S 〉 ≡ S 1 ; S 2 ; … ) is the vector of the incident electric field projected into the SPR basis, with S n = ∮ ∂ Ω ϕ n ( r ) E in ( r ) ⋅ n ̂ d 2 r . H 0 = Λ n δ n m with SPR eigenvalue Λ _n defined in Eq. (1). H R , H Q and H R − Q account for the corrections from retardation (R), quantum (Q), and coupled retardation-quantum (R–Q) effects, respectively. To reveal the geometrical dependencies of H R , Q , R − Q , we introduce the following dimensionless parameters [12]

where L is a characteristic length of the particle and k bg ≡ ω ε bg / c is the wavenumber of light in the background medium. The normalized modal fields relate to L according to

where σ n ≡ lim δ → 0 + E n ( r − δ n ̂ ) ⋅ n ̂ Λ n represents the surface charge of the nth mode, and δ _nm = 1 for n = m and otherwise 0. The matrix components of H R , Q , R − Q are then given by

Here, f n m R , Q , R − Q are dimensionless geometrical factors, expressed by

and f ̄ ̄ n m ; k R − Q is given in Eq. (S.2.9) in the SI.

Equations (4)–(6) determine the modal excitation coefficients, with which the optical responses of the system can then be reconstructed. They are derived by solving the Green’s function surface integral formalism [see Eq. (S.1.15) in the Supplementary Material], which characterizes the optical response of the system, with the SPR basis. The derivations are given in Supplementary Section S2. Specifically, Eq. (6a) gives the retardation corrections. The associated terms, f n m ; k R (k = 2, 3, …), are expressed with double surface integrals. Numerically, we find that the retardation corrections can be accurately computed by truncating the series to k = 5. Note that, in Reference [33], the retardation corrections are alternatively given by double volume integrals. In this regard, our surface-integral formulations are more computationally efficient. Moreover, the expression of f n m ; k R only involves electric fields on the particle boundary, thus being compatible with most existing SPR solvers that are implemented based on the surface-integral/boundary-element methods and solve boundary fields.

Equation (6b) gives the nonclassical correction in the absence of the retardation effects. Its diagonal terms involving f ̄ ̄ n n Q describes the leading-order nonclassical corrections, and has been comprehensively studied in Reference [32]. Interestingly, as will be discussed later, by recognizing the geometrical meaning of the Feibelman d-parameters, this nonclassical correction can be straightforwardly generalized to account for geometrical deformations (perturbations) of the particle, thereby rendering the present formalism an even broader applicability. Equation (6c) characterizes the coupling between the retardation and nonclassical effects. Despite its complicated expression, this term is numerically found to be tiny and is omitted through our numerical examples.

2.2 Independent-mode approximation

By retaining a few dominantly excited modes in Eq. (4), the excited modal coefficients can be computed by solving a small matrix problem. Nevertheless, a further simplification can be made by discarding modal couplings enabled by the retardation and nonclassical effects. More precisely, we neglect the off-diagonal terms in H R , Q , R − Q , i.e., assuming that the independence of different electrostatic SPRs in the presence of the retardation and nonclassical effects. The approximation is accordingly termed as independent-mode approximation. Under this approximation, the excitation coefficients are given by

with perturbed eigenvalue Λ n ′ expressed as

As shown in the results below, the independent-mode approximation gives rather accurate results for single nanoparticles (Figure 2). In relation with two closely related publications, References [32, 33], Eq. (7b) represents a straightforward generalization that its pole gives complex resonance frequencies including the leading-order corrections due to the retardation and nonclassical effects. For complicated plasmonic structures, such as dimer structures, this approximation loses its predictive accuracy, and the inclusion of the modal couplings due to the retardation and nonclassical effects is recommended (Figure 3).

Figure 2:

Modal analysis of extinction cross-section spectra for metal spheres (A) and rods (B), incorporating both retardation and nonclassical effects. The nanoparticles, made of Ag and Al and embedded in a dielectric medium with a refractive index of 1.33, are considered. The aspect ratio of the rods (diameter dived by height) in B is fixed to be 0.5. The results predicted by the modal method show quantitative agreements with the numerical solutions obtained with COMSOL Multiphysics. In the modal method, a single dipolar mode is retained with its excitation coefficient computed with Eq. (7a). The retardation series in Eq. (6a) is truncated to the power k = 5. The resonances show rich size-dependent shifts (dashed lines) and broadening due to the retardation and nonclassical effects. The insets sketch scattering of the nanoparticles by plane waves and profiles of dipolar modes.

Figure 3:

Modal analysis of extinction cross-section spectra for spherical dimers (A) and bow-tie structures (B), incorporating both retardation and nonclassical effects. The nanoparticles, made of Ag, are embedded in a dielectric medium with refractive index of 1.33. The gap varies from 10 nm (top) to 1 nm (bottom). The spheres in (A) have a diameter of 80 nm, while the triangles in (B) have a side length of 100 nm and a height of 25 nm. The results predicted by the modal method agree well with the numerical solutions obtained with COMSOL Multiphysics. In the modal method, four modes in both (A) and (B) are used, and the modal coefficients are obtained by solving Eq. (4). The retardation series in Eq. (6a) is truncated to the power k = 5. As the gap narrows progressively, plasmon hybridizations become stronger, which renders the resonances a redshift (visualized by dashed lines). To illustrate the effects of the nonclassical effects at small gap sizes, we include the extinction cross-section spectra under the local-response approximation (LRA) by setting d _⊥ = 0. The insets sketch scattering of the dimer structures by plane waves and profiles of the lowest-order bounding dipolar resonances.

2.3 Nonclassical effects and geometrical deformation

We observe that the nonclassical surface polarization, Eq. (4), is formally equivalent to the effective surface polarization used to represent geometrical deformations (perturbations) of a nanoparticle. More specifically, consider that a nanoparticle is geometrically deformed with its boundary shifted by a vector h ( r ) n ̂ (see Figure 4(A)). To the leading order in the deformation, the optical responses of the deformed particle can be represented by placing a surface polarization P D ( ω ; r ) = lim δ → 0 + Δ ε ( ω ) h ( r ) E tot ( r − δ n ̂ ) on the particle boundary of the original particle [35]. Comparing P ^D with P ^Q [Eq. (4)], we recognize that the nonclassical polarization has the same mathematical form as the deformation polarization. In view of this, the nonclassical effects and the geometrical deformations can thus be treated under the same framework by defining a generalized d-parameter, d _eff, with

Apparently, this treatment can be applied for modeling statistical optical responses of particle ensembles that generally have tiny shape variations due to fabrication imperfections (Figure 4).

Figure 4:

Modal analysis of geometrically deformed nanoparticles.

(A) Optical response of a deformed nanoparticle can be effectively represented by introducing a generalized d-parameter, d ̄ ̄ eff . (B) Extinction cross-section spectra of Ag spheroids embedded in a dielectric medium with a refractive index of 1.33. In the modal method, a single dipolar mode (polarized in the direction of the incident electric field) of the sphere is used to predict the extinction spectra of spheroids with D ₁ = 100 nm and D ₂/D ₁ varied from 0.8 to 1.2. The retardation series in Eq. (6a) is truncated to the power k = 5.

The above observation offers a geometrical interpretation to the Feibelman d-parameters. Particularly, when d _⊥ = d _∥, the nonclassical effects are equivalent to an isotropic geometrical deformation with h = d _⊥. In general cases with d _⊥ ≠ d _∥, this strict equivalence breaks down. Nevertheless, the recognition of the unified treatment of the nonclassical effects and the geometrical deformation still offers an appealing perspective for better understanding the physical meaning of the Feibelman d-parameters.

3.1 Single nanoparticles

We first study scattering of a plane wave by metal spheres. Two different metals, Ag (silver) and Al (aluminum), are considered. The diameters of the spheres are varied from 80 nm (with noticeable retardation effects) to 2 nm (where the nonclassical effects become important). The synthesis of metal nanoparticles widely uses colloid chemistry in liquid environments with, e.g. water as solvent [33, 36]. We thus here choose water as the embedding medium with a refractive index of 1.33. For Ag, the nonclassical effects are included by setting d _⊥ = −0.5 + 0.3i nm. Specifically, the real part of d _⊥ is determined by fitting the experimental data of the size-dependent resonance shift [20] to the perturbation theory established in Reference [32], while the imaginary part is chosen to reproduce the well-known Kreibig damping [37] v _F/R (v _F, Fermi velocity; R, sphere radius). For Al, we use d _⊥ predicted with the time-dependent density-functional theory [38], which reasonably fits with the experimental observations [39]. The impact of d _∥ is known to be negligible for charge-neutral materials, and we thus set d _∥ = 0 [16]. Moreover, we note that d-parameters depend on specific dielectric environments surrounding metals. The used d _⊥’s here omit this dependence [40], and are thus approximate. Nevertheless, the main purpose of the present work focusing on clarifying the validity of the modal method is unaffected by this approximation.

Figure 2(A) plots the extinction cross-section spectra of the spheres, comparing the predictions from the modal method and the full numerical results obtained with COMSOL Multiphysics [25, 41]. Note that, in the COMSOL simulations, the effects of the d-parameters are incorporated by introducing an auxiliary potential, see Reference [25] for details. We focus on the spectral region wherein the lowest-order dipolar mode is dominant (see the inset in Figure 2(A) for the modal profile). Therefore, only a single dipolar mode is used in the modal method, and the independent-mode approximation [Eq. (7a)] is employed to compute the modal coefficient. The comparisons validate the accuracy of the modal method. In Supplementary Figure S1, we further validate the modal method by comparing its results wiwith the analytic solutions predicted by the mesoscopic Mie’s theory that incorporates d-parameters [42]. Moreover, as shown in Supplementary Figure S2(A), the modal method continually to provide reliable predictions for spheres with diameter as large as 200–300 nm.

Notably, as the sphere diameter increases in the range (roughly) above 10 nm, the dipolar resonances show red shifts along with line broadening. This phenomenon is due to the retardation effects and is precisely reproduced by the modal method. On the other hand, as the sphere diameter decreases in the deep-nanometer scale below 10 nm, the resonance shows 1/D (D, diameter) size-dependent shifts and broadening. This is due to the nonclassical effects. Specifically, for Ag with Re(d _⊥) < 0, the resonant blueshifts are observed [19, 20], while for Al with Re(d _⊥) > 0, the opposite redshifts instead occur [21]. Moreover, the 1/D-dependent broadening is known as the Kreibig damping [37], which in our formalism is characterized by the imaginary part of d _⊥.

To further validate the predictive accuracy of the modal method, we additionally test nanoparticles of other shapes. Figure 2(B) and Supplementary Figure S2(B) illustrate the extinction cross-section spectra for rod structures. The comparisons between the modal method and the full numerical solutions present quantitative agreements.

3.2 Multi-scale plasmonic architectures

Plasmonic dimers, film-coupled plasmonic nanoresonators and other similar structures, with particle dimensions D comparable to light wavelength and gap sizes g close to electronic length scales (Figure 1(B)), are well known for supporting extreme light–matter interactions inside tiny-gap regions [17, 23, 35]. In these multi-scale structures, the retardation and nonclassical effects coexist, interplay with plasmon oscillations, thereby significantly affecting optical responses. In this regard, their simulations require electrodynamic solvers that incorporate both the retardation and nonclassical effects, which just fits the application scope of the present modal method.

We first consider a spherical dimer. The dimer hosts strong hybridizations of plasmon resonances of individual spheres, leading to rich resonance features. To include multiple resonance features, we here retain the four lowest-order electrostatic SPRs in the modal method (see Supplementary Figure S3 for modal profiles). The dimer is made of Ag with spherical diameter 80 nm. We vary the gap distance between two spheres from 10 nm to 1 nm. Figure 3(A) compares the extinction cross-section spectra computed with the modal method with the numerical solutions obtained with COMSOL Multiphysics. Referring to the modal method, we note that the modal couplings due to the retardation and nonclassical effects cannot be neglected in this case (see Supplementary Figure S5). Hence, we compute the modal coefficients by solving the full matrix of Eq. (4), i.e., retaining the off-diagonal terms in the Hamiltonian, instead of using the simplified independent-mode approximation. The excellent agreements between the modal results and the fully numerical solutions prove the soundness of the modal method for describing optical responses of complex plasmonic structures. Besides, to highlight the significant impact of the nonclassical effects at small gap sizes, we also plot the spectra under the LRA by setting d _⊥ = 0 for 1 nm and 2 nm gaps. The comparisons between the nonclassical and LRA results show that the nonclassical effects lead to considerable resonance shifts and damping broadening, which have been observed in recent experiments [23, 25]. Next, we consider a bow-tie structure. Figure 3(B) plots the extinction cross-section spectra. In the modal method, we use the four lowest-order static modes in the modal method (see Supplementary Figure S4 for modal profiles). The comparisons between the modal method and the full numerical solutions again show quantitative agreements.

3.3 Geometrically deformed nanoparticles

As discussed above, optical responses of a set of nanoparticles, deformed from an original nanoparticle, can be conveniently treated (to the leading order of the deformation) by defining a generalized d-parameter, d ̄ ̄ eff , by simply adding the deformation contribution on top of the nonclassical d-parameter [see Eq. (8) and Figure 4(A)]. To validate this idea, we consider a sphere and take its electrostatic plasmon modes to predict the extinction cross-section spectra of spheroids. A single dipolar mode (polarized in the direction of the incident electric field) is retained in the modal formalism. In Figure 4(B), we compare the spectra predicted with the modal method to the full numerical solutions obtained with COMSOL Multiphysics. The qualitative agreements between two approaches validate the effectiveness of the modal method for predicting optical responses of geometrically deformed particles. Note that, the numerical results show two spectral peaks, while the modal method only predicts the long-wavelength one due to the use of a single dipolar mode. Moreover, since d ̄ ̄ eff includes both perpendicular and parallel components (with respect to the particle boundary ∂Ω), the above results also validate that the present method can handle general cases with both nonzero d _⊥ and d _∥ equally well.