Plasmon Hybridization Model for a Nanoparticle ... - ACS Publications

Aug 27, 2015 - Jean Lermé*. Institut Lumière Matière, UMR5306 Université Lyon 1-CNRS, Université de Lyon, Bât A. Kastler, 43 Bld du 11 Novembre 1918,...
2 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Plasmon Hybridization Model for a Nanoparticle above a Dielectric Interface: Dielectric Effects, Comparison with the Dimer System, Range of Applicability, and Limits Jean Lermé* Institut Lumière Matière, UMR5306 Université Lyon 1-CNRS, Université de Lyon, Bât A. Kastler, 43 Bld du 11 Novembre 1918, F-69622 Villeurbanne Cedex, France

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

S Supporting Information *

ABSTRACT: In the presence of a nearby highly reflecting dielectric substrate, the optical response of a metal nanoparticle exhibits a rich multipolar plasmonic structure. For quantifying the multimodal plasmonic response of such hybrid particle/interface systems, the hydrodynamic plasmon hybridization (PH) method has been applied in the simple case of a solid metal sphere (PI-PH model). In this work I extend the formalism described in my previous paper (J. Phys. Chem. C 2014, 118, 28118−) to the more general case where underlying dielectric media are present in the whole space (polarizable ionic background and matrix). In the quasistatic limit the method of images allows the PI-PH model to be closely related to the two-sphere dimer PH model. However, because the dynamics of the electron fluid and of the polarization charges in the mirror particle are enslaved to those in the real metal particle, strong differences between both PH models result from the reduction by a factor of two of the number of dynamical variables describing the system. As compared to jellium spheres in vacuum, the PH formalism is much more complex because the induced surface polarization charges in the underlying media contribute to the energetics of the system. In this work a very detailed analytical description of both PH models, allowing a direct numerical implementation, is provided. In the frame of the original and physically transparent method developed in my previous paper, particular emphasis is given to modeling the coupling of the nanosystem with an external applied field. In the presence of underlying dielectric media, I show that the taking into account of the field created by the additional polarization charge distributions induced by the bare applied field is essential for defining the effective external field that is coupled to the nanosystem. The predictions of both PH models are also compared, through absorption and extinction cross-section spectra, to exact results computed in the frame of the multipole expansion method, allowing the range of applicability and limits of both models to be determined. Finally, the PH formalism is shown to be generalizable to deal with underlying dielectric media characterized by complex (that is, absorbing) frequency-dependent permittivities.

1. INTRODUCTION In the exponentially growing nanoscience and nanotechnology fields, single nanosized metal particles or aggregates of metal particles sustaining collective oscillations of conduction electrons play a central role.1−8 When driven by light, these collective excitations, referred to as localized surface plasmon resonances (LSPR), give rise to large absorption/scattering cross sections and concomitant enhanced local field at specific wavelengths. These appealing properties, which can be spectrally tuned through size, shape, and local particle environment effects,9−11 as well as through the geometrical arrangements in the case of aggregates of particles,3,12−27 are at the origin of the innumerable optical applications offered by metal particle-based devices of subwavelength dimensions. Whatever the application planned, optimizing the performances of nanodevices comprising metal particles requires a preliminary detailed characterization of the intrinsic properties of the elemental building blocks. Various sensitive far-field optical techniques, suitably designed for probing single © XXXX American Chemical Society

nanoobjects, often combined with electron-microscopy imaging, have been developed for ensuring such a thorough characterization.28−36 However, gaining reliable quantitative information about the intrinsic optical properties of the optically addressed nanoobject requires, in the theoretical analysis, taking into account the specific experimental conditions prevailing in such experiments. In particular, the strong convergence of the excitation light beam, focused close to the diffraction limit,37−39 and the substrate on which is deposited the nanoobject,40 are ingredients to be imperatively included in the modeling. Especially, a thorough analysis of the substrate effects is required in the context of LSPR-based sensing applications involving supported nanosensors of high figure of merit, since very tiny LSPR shifts are induced by adsorbed molecules.4 This past decade a lot of works have been devoted to these effects, emphasizing that they are of main Received: June 26, 2015

A

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C

analytical cross-section formula allows analysis and quantification of how the oscillator strengths of the dipole excitations of the isolated particles are progressively transferred to the whole set of plasmon excitations of the global system when the interparticle distance is decreased (emergence of the multimodal patterns in the spectra in the strong coupling regime).44,45 As in paper II the range of applicability and limits of both models are also determined through comparison with exact multipole expansion method (MEM) calculations.46,47 For avoiding an excessive manuscript length, most of the illustrative examples and spectra are reported on in the SI. I will also discuss a crucial point that, to my knowledge, has not been addressed in the previously published PH-related works. This important issue concerns the coupling of any nanoobject to an external applied electromagnetic field, in the presence of underlying polarizable media. Actually, in addition to the surface polarization charge distributions on both sides of each dielectric interface involved in the system, that are directly48 induced by the free charge densities, additional surface polarization charge distributions are induced by the external applied field on both sides of each dielectric interface. So, in the presence of underlying polarizable media, the total field generated by all these additional polarization charge distributions has to be added to the bare external applied field to define the ef fective external field that is coupled to the nanoobject. As a matter of fact, this point, which is necessary for ensuring quantitative predictions, has already been addressed in previous works aiming at introducing quantum finite-size effects in the classical quasistatic Mie theory, in the case of multilayered metal/dielectric spheres.49,50 It should be pointed out that, in the case of multiparticle systems such as those investigated in this work, the situation is more complex, as compared to spherical systems, because each particle experiences the strongly inhomogeneous field generated by the bare fieldinduced polarization charge distributions located on the dielectric interfaces of the other nearby particles. Finally I will show that the PH formalism can be generalized to deal with underlying dielectric media characterized by arbitrary complex (that is, absorbing) frequency-dependent electric permittivities. In this respect, I will provide explicit modified PH formulas for dealing with absorbing dielectric media, in the case of several nanosystems of interest. This paper is organized as follows. The D-PH and PI-PH models and their corresponding illustrative examples are reported in sections 2 and 3, respectively. A short comparison with previous works is given in section 4, and the extension of the “standard” PH formalism for dealing with frequencydependent absorbing underlying media is discussed in section 5. Finally, this work is summarized in section 6. Atomic units are used in this work. For convenience and for comparison purposes with previous works, the parameters e (minus the electron charge) and me (free electron mass) will be kept in the equations.

importance in many cases. Lists (nonexhaustive) of relevant works on these issues can be found in our previous detailed works: refs 37 and 40, respectively. This paper is a direct continuation of two previous works40,41 (referred to as paper I and paper II, respectively), aiming at studying the effects of a nearby dielectric interface on the optical properties of a metal nanoparticle. In paper I, calculations involving spherical noble metal particles have been reported on for illustrating the specific effects resulting from each relevant parameter characterizing the particle/ dielectric interface system, namely, the particle size, the particle−interface distance, the dielectric index of the substrate, the polarization, and the angle of incidence of the excitation light beam. These investigations have shown that, in the presence of a highly reflecting dielectric substrate, the optical response differs dramatically from that of an isolated particle. In particular, the spectra exhibit a rich multimodal plasmonic structure, strongly dependent on the particle/substrate distance and the irradiation conditions (polarization and angle of incidence), very similar to the one observed in the case of metal sphere dimers in the strong coupling regime.13,16 In paper II, I have shown that the multimodal plasmonic pattern can be rationalized and accurately quantified in the quasistatic limit (retardation effects neglected), in the framework of the appealing hydrodynamic plasmon hybridization (PH) approach, initially developed by the Nordlander group in the case of concentric nanoshells42,43 and jellium sphere dimers.12 The tight relationship between the PH models for the metal sphere dimer and particle/dielectric interface systems (referred to as the D-PH and PI-PH models, respectively) has been brought into light in noting that, in the quasistastic limit, the method of images allows the electrostatic metal sphere−plane dielectric interface interaction to be expressed as the interaction of the metal sphere with a fictitious metal sphere of identical radius located below the interface. In paper II the D-PH and PI-PH models have been described in great detail for the simplest generic systems, namely: (1) two jellium spheres in vacuum and (2) a jellium sphere in vacuum located above a dielectric interface of refractive index Nsub. In the present work, I extend the PH modeling to the more realistic case where homogeneous polarizable media are present in the whole space (polarizable underlying backgrounds inside the jellium spheres and polarizable embedding matrix). For the axially symmetric dimer and particle/interface systems, the introduction of underlying polarizable media leads to a much more complex PH formalism, in stark contrast with the corresponding PH modeling appropriate to spherically symmetric layered particles (the so-called nanomatryushka structures), for which the PH formalism remains of reasonable complexity.43 Since the PI-PH model in the presence of underlying dielectric media is tightly related to the corresponding D-PH model, a detailed description of the D-PH model, including explicit formulas allowing a direct numerical implementation, is reported on in the Supporting Information (SI). Only the ingredients and main equations of the D-PH model, that are necessary for facilitating the reading, and useful for building the PI-PH model, are given in the main text. Second, the coupling of both systems to an external homogeneous electromagnetic field will be handled according to the original procedure described in paper II.41 This method allows the absorption cross section σ(ω) to be expressed in a simple way as the sum of the individual contributions arising from all the eigenmodes of the dimer. In particular the

2. PLASMON HYBRIDIZATION (PH) MODEL FOR A SPHERE DIMER In the presence of underlying homogeneous linear dielectric media, the formalism is indeed much more complex but the general strategy for elaborating the PH model is not significantly modified as compared to the case of jellium sphere dimers.41 Before reporting on the detailed equations and technical aspects of the PH formalism, I summarize and B

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C emphasize here the main points to be kept in mind when dielectric media are present. First, the polarization of the underlying dielectric media will act through screening effects, reducing the plasmon eigenfrequencies, as expected, and, second, will contribute to the field-induced dipole and dynamical polarizability when the coupling to an external field will be introduced. In the Maxwell equations, expressed in terms of the free charge distributions only, the dynamics of the bound charges (polarization in the media) is actually enslaved to that of the free charge distributions. This comment can be inferred in just noting that, for given free charge distributions, the solution of the Maxwell equations and the constitutive equations allow determining the polarization of the underlying media.51 So, no extra dynamical variables have to be introduced for characterizing the polarization charge distributions that could develop in the ionic backgrounds and the matrix: only the equations of motion of the dynamical variables associated with the electron fluid deformations (giving rise to the free-charge densities) are required in the PH modeling. Besides, in the case of homogeneous linear dielectric media, the electrostatic potential created by the polarization of the bound charges can be replaced by that created by surface charge densities located on both sides of each dielectric interface present in the nanosystem.51 As will be seen in the following subsection, a linear relation can be established between the free and polarization charge distributions. In the dimer system, in the absence of an external field, two independent polarization charge densities add up on both sides of the two spherical interfaces. The first distribution, which develops on the surfaces of both spheres, is directly48 induced by the electron fluid deformations in sphere no. 1, and the second distribution by those in sphere no. 2 (as a consequence of the linearity of the Maxwell equations). Actually, it is not necessary to determine each individual density contributing to the overall surface charge distributions for elaborating the PH modeling, because, to compute the electrostatic energy of the system (V), only the free-charge densities and the total electrostatic potential are required. Four surface charge distributions are thus involved in the modeling: the surface free-charge densities σ(i)(Ωi) (with i = 1 or 2) and the surface total-charge densities σt(i)(Ωi) (with i = 1 or 2), with each of the latter including the corresponding free-charge density σ(i)(Ωi) and the polarization charge densities on both sides of the interface that are induced by the two free-charge densities. Thanks to the linear relation linking the two sets {σ(i)(Ωi), i = 1,2} and {σt(i)(Ωi), i = 1,2}, the electrostatic potential V and, thus, the Lagrangian are finally expressed in terms of the dynamical variables describing the two free charge densities. The free problem, namely the determination of the eigenfrequencies and eigenvectors, is then solved as in the case of jellium sphere dimers. When the dimer is subject to an external electromagnetic field, the handling of the dimer/field interaction problem, and the successive steps leading to the cross-section formula, are quite similar to those prevailing in the jellium sphere dimer case.41 The novel crucial point in the PH formalism consists in determining the additional field generated by the surface polarization charge distributions that are directly48 induced by the bare external field. 2.1. Lagrangian, Equations of Motion, and Eigenmodes. Figure 1 displays the dimer geometry and the main ingredients involved in the modeling. The dimensionless

Figure 1. Geometry of the dimer, and notations and coordinate systems used. d is the interparticle distance. σ(i) and σt(i) are the surface free-charge and total-charge densities, respectively.

relative electric permittivities of the three underlying dielectric media, assumed homogeneous and isotropic, are ε1, ε2, and εm. The surface free-charge densities in both spheres are expanded upon the complete set of scalar spherical harmonics (standard definitions of these complex angular functions are assumed)52 S



σ (i)(Ωi , t ) =

∑∑

σS(, mi)(t )Y Sm(Ωi)

S = 1 m =−S

(1)

⎛ S ⎞1/2 = −nie⎜ 3 ⎟ SS(,im) (t ) ⎝ Ri ⎠

(2)

with σS(, mi)(t )

In eq 2, −nie (e > 0) is the homogeneous charge density of the incompressible electron fluid in the sphere i. Specific normalization constants have been introduced for convenience, as in refs 12 and 43, and paper II. The deformation amplitudes SS(,im) (unit: length5/2) are the dynamical variables in the PH model. The electron kinetic energy [T(i)(t)] in each sphere can be expressed in a compact form in terms of these dynamical variables T (i)(t ) =

nime 2



S

∑∑

(i)

(i)

( −1)m SṠ , m(t )SṠ , −m(t )

S = 1 m =−S

(3)

where me is the electron mass. The above expression is quite general for solid spheres and holds in the absence or in the presence of underlying dielectric media. The total electrostatic energy of the dimer (adiabatic building-up of the free-charge densities) is computed thanks to the general equation of electrostatics, suitable on the condition that the underlying dielectric media are linear,51 namely V = 1/2∫ ρ(r)Φt(r)d3r, where ρ(r) = σ(1)(r) + σ(2)(r) is the overall free charge density and Φt(r) the total electrostatic potential (generated by both the free and enslaved polarization charges). Since Φt(r) can be expressed in terms of the surface total-charge densities σt(i)(Ωi) only, and depends linearly on both densities, V is the sum of four terms 1 1 1 V = V (σ (1) , σ t(1)) + V (σ (1) , σ t(2)) + V (σ (2) , σ t(1)) 2 2 2 1 + V (σ (2) , σ t(2)) (4) 2 where (4πε0 = 1) V (σ (i) , σ t(j)) = C



σ (i)(ri, t )σ t(j)(rj, t ) | M i M j|

d 2Sid 2Sj

(5)

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C In eq 5 Mi and Mj are the two points on the surfaces of the spheres i and j specified by the coordinate-vectors ri and rj in their respective Cartesian coordinate frames. Each term can be expressed in terms of the multipolar components of both involved surface densities, thanks to the formulas listed in section A of the SI (eqs 11 and 31), namely V (σ (i) , σ t(i)) =

∑ (−1)m S, m

σ ′(2) S2, m =

σ ′(Sii,)m

(6)

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

m

1 2

(7)

V (σ (2) , σ t(1)) ⎛ R 3R 3 ⎞1/2 = 2π ∑ ( −1) ⎜ 1 2 ⎟ Va(S1, |m| , S2)σS(2) σ t(1) 2, m S1, −m S S ⎠ ⎝ 1 2 S ,S ,m m

1 2

εm

εm

∂Φ (r1) ∂r1 ∂Φt (r2) ∂r2

R1+

∂Φ (r1) ∂r1

− ε2 R 2+

∂Φt (r2) ∂r2

introduce linear combinations of the dynamical variables SS(,im) and SS(,i−) m. The resulting new variables will allow (i) obtaining

= −4πσ (2)(Ω 2)

associated densities with specific and well-defined symmetries relative to the Cartesian axes (see eq 8 in paper II) and (ii) facilitating the discussion when the interaction with a linearly polarized plane wave will be addressed. The new dynamical variables are

(10)

Since Φt is linear with respect to the total charge densities σt(i), eqs 9 and 10 provide directly the expressions of the free-charge densities as a function of the total-charge ones. This essential step in the PH modeling consists, first, in expressing Φt in the coordinate frames attached to spheres 1 and 2 in order to handle the left-hand sides of eq 9 and eq 10, respectively, and second, in introducing the multipolar expansions of the surface densities (see eq 1) for deriving the linear relations linking the multipolar expansion coefficients of the surface free- and totalcharge densities, for each m value (due to the axial symmetry of the dimer, the problem is diagonal with respect to the m value). This important step is detailed in the SI (section C.1). One obtains σ ′(1) S1, m =

[(S1 + 1)εm + S1ε1] t(1) σ ′S1, m 2S1 + 1 1 + (ε1 − εm) ∑ Va(S1, |m| , S2)σ ′tS2(2) ,m 2 S2

(14)

four terms and of the Lagrangian L = T − V in terms of these dynamical variables are given in the SI (sections C.2 and C.3, respectively)]. As in paper II, it is valuable at this stage to

(9)

R 2−

(13)

of the dynamical variables SS(,im) [the explicit expressions of the

= −4πσ (1)(Ω1) R1−

⎛ R 3 ⎞1/2 = ⎜ i ⎟ σSti,(mi) ⎝ Si ⎠

Note that, for a given S max value, the size of the four square matrixes [χij], as well as that of the column vectors, depends on m [Si varies between Smin = max(1,|m|) and Smax ]. Using eqs 2, 6−8, and 13−14 the four terms contributing to the total electrostatic energy V, eq 4, can be expressed in terms

t

− ε1

σ ′tSi(, mi)

⎡ σ ′t(1) ⎤ ⎡ 11 12 ⎤⎡ σ ′(1) ⎤ ⎢ S1, m ⎥ = ⎢ χ χ ⎥⎢ S1, m ⎥ ⎢⎣ σ ′t(2) ⎥⎦ ⎣ χ 21 χ 22 ⎦⎢⎣ σ ′(2) ⎥⎦ S2, m S2, m

(8)

In eqs 7 and 8, Va is the dimensionless coupling integral (depends on R1, R2, and d) that was already introduced in paper II (eq 25 in the present SI). The central task is now to establish the linear relation, giving the total-charge densities in terms of the free-charge ones. As done in ref 43 for nanoshells, the relationship can be determined thanks to the boundary condition expressing the discontinuity of the normal component of the electric displacement across both dielectric interfaces due to the presence of the surface free-charge densities,51 namely t

⎛ R 3 ⎞1/2 = ⎜ i ⎟ σS(i,im) ; ⎝ Si ⎠

Assuming a maximum S value (S max) in the multipolar expansions of the densities, eqs 11 and 12 define a matrix equation for each m value, namely σ′ = [E]σ′t (the matrix [E] does not depend on the sign of m). Defining the |m|-dependent matrix [χ] as the inverse of [E], the required linear relationship is finally obtained (for lightening the notations, the |m|dependence is omitted in [χ])

V (σ (1) , σ t(2)) ⎛ R 3R 3 ⎞1/2 = 2π ∑ ( −1) ⎜ 1 2 ⎟ Va(S1, |m| , S2)σS(1) σ t(2) 1, m S2, −m S S ⎠ ⎝ 1 2 S ,S ,m

(12)

where the multipolar components σ′ are defined by

4π R i3σS(, mi)σSt,(−i)m 2S + 1

(i = 1 and 2)

[(S2 + 1)εm + S2ε2] t(2) σ ′S2, m 2S2 + 1 1 + (ε2 − εm) ∑ Va(S1, |m| , S2)σ ′tS1(1) ,m 2 S1

ZS(,0i)(t ) = SS(,0i)(t )

XS(,im) (t ) = −

(m = 0)

(15)

1 [SS(,im) (t ) + ( −1)m SS(,i−) m(t )] 2

(m > 0) (16)

Y S(,im) (t ) = −

1 i[SS(,im) (t ) − ( −1)m SS(,i−) m(t )] 2

(m > 0) (17)

In terms of these new variables, the Lagrangian L = T(1) + T(2) − V of the dimer writes as (with the following convention: for m = 0, X is replaced by Z, and Y = 0)

(11) D

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C L=

+





n1me 2 n2me 2

Smax

S

Smax

(1)2

∑ ∑ m = 0 S1= Smin Smax

[Ẋ S(1)2 + YṠ 1, m ] 1, m

Smax

(2)2

∑ ∑

[Ẋ S(2)2 + YṠ 2, m ] 2, m

m = 0 S2 = Smin

men1 (1)2 ωP 2 men2 (2)2 ωP 2

Smax

Ẍ ′(2) S2, m +

Smax

Smax

∑ ∑



m = 0 S1= Smin S1′= Smin Smax

Smax

∑ ∑



Smax

m = 0 S2 = Smin S2 ′= Smin

+

(X S(1) A S11 X (1) + Y S(1) Y (1) ) 1S1′ 1, m S1′ , m 1, m S1′ , m

Smax

Smax

∑ ∑

Smax



Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

m = 0 S1= Smin S2 = Smin

(X S(2) A S22 X (2) + Y S(2) Y (2) ) 2S2 ′ 2, m S2 ′ , m 2, m S2 ′, m

(A S12 )(X S(1) X (2) + Y S(1) Y (2) ) + A S21 1S2 2S1 1, m S2, m 1, m S2, m

(18)

+

ωP(1)ωP(2) 2

S1′= Smin

S2 = Smin

(23)

[Vk̈ , m] + [M ][C ][M ]−1 [Vk , m] = [Vk̈ , m] + [ωk2, m][Vk , m] = 0 (24)

[ω2k,m]

where is a diagonal square matrix (dimension 2Sd ), the diagonal coefficients of which (eigenvalues of the matrix problem) are the squares of the eigenmode frequencies of the dimer system. 2.2. Coupling to a Monochromatic Electromagnetic Field within the Quasistatic Approximation. The dimer is now subject to a spatially homogeneous monochromatic electric field E(t) = E0e−iωt (bare external applied field). In this simple case, the bare interaction can be derived from the scalar potential Φ E(r, t ) = −r·E(t ) = −xEx(t ) − yEy(t ) − zEz(t )

(25)

In the following the time dependence of the Cartesian components of the electric field will be omitted for lightening the equations. The electrostatic energy term corresponding to the dimer-external field interaction [building up of the surface free-charge density ρ(r) = σ(1)(r) + σ(2)(r) in the presence of the external field] expresses as

(1) (A S11 + A S11 ′1S1)X ′S1′ , m ] 1S1 ′

∫ ρ(r)[ΦE(r) + Φ pol ,E(r)]d3r bare = V ext + ∫ ρ(r)Φ pol , E(r)d3r

Vext =

Smax



(22)

The equations of motion of the variables Vk,m take the compact matrix form

Smax



(21)

⎡ X ′(1) ⎤ S1, m ⎥ [Vk , m] = [M ]⎢ (2) ⎢⎣ X ′ ⎥⎦ S2, m

S S χ ii = ωP(i)2 δS, S ′ 2S + 1 S, S ′ (S + 1)εm + Sεi

The right-hand side of eq 19 is nothing else but the square S mode plasmon frequencies of the isolated solid spheres in the presence of inner and outer dielectric media, proving that, in the limit Va → 0, the Lagrangian eq 18 reduces to the sum of the Lagrangians associated with two isolated matrix-embedded spheres (see eqs 48 and 42 in the SI). The structure of the Lagrangian, eq 18, shows that the hybridization mechanism, i.e. the coupling of the plasmon modes of both spheres, is diagonal with respect to the m value (m ≥ 0) and, for m > 0, that the X and Y variables are uncoupled and will correspond to degenerate plasmon modes (same eigenvalues and eigenvectors in the resulting matrix problem described hereafter). These properties stem from the axial symmetry of the dimer and the equivalence of the x and y axes. Since the m = 0 (variables Z) and m > 0 (variables X or Y) problems are uncoupled and lead to identical formal equations, only the equations of motion corresponding to the m > 0 (variables X) problem are reported on. Introducing the convenient dynamical variables X′(i) = (ni)1/2X(i) (the resulting matrix [C] to be diagonalized will be symmetric), the Euler− Lagrange equations derived from the Lagrangian write as ωP(1)2 [ 2

S1= Smin

[(A S12 )X ′(1) + A S21 S1, m ] = 0 1S2 2S1

The square matrix [C] of dimension 2Sd is real and symmetric and can be diagonalized thanks to standard numerical algorithms. In particular, the 2Sd eigenvalues ω2k,m (k = 1, 2, ..., 2S d) are real, and the corresponding eigenvectors {Vk,m; k = 1, 2, ..., 2S d} (that is, the uncoupled dynamical variables) can be chosen real and orthonormalized. The uncoupled dynamical variables are therefore related to the initial dynamical variables (2) {X ′(1) S1, m } and {X ′S2, m } through a real orthogonal matrix [M] (mdependent square matrix of dimension 2Sd such as [M]−1 = [M]T) in which the k-th row consists of the expansion coefficients of the uncoupled dynamical variables Vk,m as a function of the 2Sd initial dynamical variables, namely [hereafter [Vk,m] is a column vector containing the 2Sd uncoupled dynamical variables]

(19)

Ẍ ′(1) S1, m +

Smax



⎡ Ẍ ′(1) ⎤ ⎡ X ′(1) ⎤ ⎢ S1, m ⎥ + [C ]⎢ S1, m ⎥ = 0 ⎢⎣ Ẍ ′(2) ⎥⎦ ⎢⎣ X ′(2) ⎥⎦ S2, m S2, m

2 1/2 where the ω(i) values (with i = 1 and 2) are the P = (4πnie /me) bulk plasma frequencies of both metals.7 The explicit expressions of the matrix elements of the four |m|-dependent square matrixes [Aij] are given in the SI (eqs 104−107). When no dielectric media is present (ε1 = ε2 = εm = 1; [χ] is the unit 21 matrix), one has ASii, S ′ = (S /(2S + 1))δS, S ′ and AS12 , S ′ = AS ′ , S = 1 /2Va(S ,|m|,S ′), and the jellium-sphere dimer PH formalism is recovered (paper II).41 It is easy to check also that the expected results for isolated spheres are ensured when the interparticle distance is very large relative to the sphere radii (that is, formally, in the limit Va → 0). In this case [E] and [χ] are diagonal, and A12 = A21 = 0. One obtains χSii, S ′ = δS, S ′(2S + 1)/[(S + 1)εm + S εi] and thus

ωP(i)2ASii, S ′ = ωP(i)2

ωP(1)ωP(2) 2

The above set of second-order differential eqs [2(S max − S min + 1) = 2S d equations and dynamical variables] is indeed an eigenvalue problem which can be solved thanks to matrix techniques. In matrix form, eqs 20−21 express as

m (n n )1/2 (1) (2) ωP ωP − e 1 2 2 ×

ωP(2)2 max (2) [ ∑ (A S22 + A S22 ′2S2 )X ′S2 ′ , m ] 2S′2 2 S2 ′= Smin

[(A S12 )X ′(2) + A S21 S2, m ] = 0 1S2 2S1 (20) E

(26)

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C where Φpol,E(r) is the electrostatic potential generated by the total surface polarization charge densities in the spheres 1 and 2 that are directly induced by the applied field E [denoted σE(i)(ri)]. These additional induced polarization charge densities are determined in solving the electrostatic problem (matrix-embedded metal spheres in the presence of the external applied field E) in the absence of the free-charge densities σ(i).48 They can be obtained straightforwardly thanks to the boundary condition expressing the continuity of the normal component of the electric displacement across both dielectric interfaces, through equations similar to eqs 9−10. In eqs 9 and 10, the right-hand sides are now equal to zero and the electrostatic potentials Φt(ri) are replaced by their counterparts ΦE(ri) + Φpol,E(ri) [expressions of the total electrostatic potential ΦE + Φpol,E in the coordinate frames attached to spheres 1 and 2]. Since the first derivative of ΦE eq 25 is continuous across both interfaces, the equations can be cast in the following form εm

∂Φ

pol , E

(r1)

∂r1

− ε1

∂Φ

R1+

pol , E

(r1)

∂r1

Smax

Vext = (a

S1= 1

∂Φ pol , E(r2) ∂r2

− ε2

+ (a(1)

R2

+

+ (a(1)

TS(1) = δ S1,1 − 1, m



TS(2) 2, m

= −4π Θ(2)(Ω 2)

(29)

multipolar components of the auxiliary components σ′SEi,(mi) and Θ′(Sii,)m (defined as in eq 13). All the ingredients required to compute explicitly Vext, eq 26, have been determined. The first term, Vbare ext , is identical to the one appropriate to jellium-sphere dimers (eq 20 in paper II) ⎛ 4π ⎞1/2 (i) (i) (i) ⎟ [EzZ1,0 ] + ExX1,1 + EyY1,1 3 ⎠

∑ nieR i3/2⎜⎝ 2

=

TS(1) X (1) + a(2) 1,1 S1,1

∑ S2 = 1

TS(2) X (2))Ex 2,1 S2,1

Smax

TS(1) Y (1) + a(2) 1,1 S1,1

∑ S2 = 1

TS(2) Y (2))Ey = − D·E 2,1 S2,1

1 2

Smax

∑ S2 = 1

⎛ ⎛ R 2 ⎞3/2 12 ⎞ S1 ⎜ 11 ε − ε χ + ε − ε ( ) ( ) ⎜ ⎟ χS ,1 ⎟⎟ 1 m 2 m S ,1 1 1 2S1 + 1 ⎜⎝ ⎝ R1 ⎠ ⎠

⎛ ⎛ R ⎞3/2 ⎞ Va(S1, m , S2)⎜⎜(ε1 − εm)χS21,1 + (ε2 − εm)⎜ 2 ⎟ χS22,1 ⎟⎟ 2 2 ⎝ R1 ⎠ ⎝ ⎠

⎛ ⎞ ⎛ R1 ⎞3/2 21 S2 ⎜ 22 ⎟ = δ S2,1 − ε − ε χ + ε − ε χ ( ) ( ) ⎜ ⎟ 1 m 2 m S2,1 ⎟ S2,1 2S2 + 1 ⎜⎝ ⎝ R2 ⎠ ⎠ 1 2

Smax

∑ S1= 1

⎞ ⎛ ⎛ R ⎞3/2 Va(S1, m , S2)⎜⎜(ε1 − εm)⎜ 1 ⎟ χS11,1 + (ε2 − εm)χS12,1 ⎟⎟ 1 1 ⎝ R2 ⎠ ⎠ ⎝

eq 31 defines the “dipole” D from which the absorption cross section σ(ω) is computed. One can remark that, due to the presence of the induced polarization charge densities σE(i), the effective dimer/field coupling term, Vext, involves, from the outset, dynamical variables associated with multipoles of high S orders. The relative contributions to absorption stemming from the various S order terms will depend strongly on the intersphere distance. The high S order terms will be promoted when the spheres approach each other, due to the increasing inhomogeneity of the effective external field experienced by each sphere, whereas only the dipolar terms will contribute noticeably for large enough distances. In the presence of the external applied field E, −Vext has to be added to the Lagrangian eq 18, modifying therefore the equations of motion of the m = 0 and m = 1 related dynamical variables. Because the coupling term is linear with respect to the dynamical variables involved, the problem remains diagonal with respect to the m value and the symmetry (X or Y character for m = 1). Moreover, the three independent sets consisting of the X, Y, and Z related modes will be specifically excited by electric fields linearly polarized along the three axes of the Cartesian reference frames, namely x, y, and z, respectively (see Figure 1). Only an additional time-dependent driving term has to be introduced in the equations of motion to be modified (eqs 141, 142 in the SI). For instance, for a polarization along the x axis (E0 = E0ex), the changes of the equations of motion, eqs 20−21, consist in just adding in the right-hand side of these equations the corresponding driving term: − (a(i)/(me√ni))TS(i,1i) Ex(t) (with i = 1 in eq 20 and i = 2 in eq 21). The absorption cross section σ(ω) can be easily derived by following the same procedure as that described in paper II.41 The method consists in expressing the dipole D (defined in eq 31) as a function of the time-dependent steady-state solutions of the equations of motion of the dynamical variables, Vk,m, which remain uncoupled in the presence of the external field E. The details

(θi,ϕi). One obtains Θ(S,im) = (εi − εm)(12π)−1/2US, mδS,1 with U1,0 = Ez, U1,1 = −(1/√2)(Ex − iEy), and U1,−1 = (1/√2)(Ex + iEy). Equations 11−12 can thus be directly applied to relate the

i=1

Smax

(33)



The multipolar components of these fictitious densities can be obtained in expressing the Cartesian coordinates xi, yi, and zi in Φ(E), eq 25, as a function of the scalar spherical harmonics Y Sm= 1

2

S2 = 1

TS(2) Z (2))Ez 2,0 S2,0

(32)

where the fictitious surface free-charge densities Θ are defined by

bare = V ext



with (m = 0 or 1, keeping in mind that the square matrixes χij depend on m)



ri = R i

+a

(31)

(28)

εm − εi ∂Φ(E)(ri) 4π ∂ri

∑ S1= 1

(i)

Θ(i)(Ωi) =

(2)

Smax

= −4π Θ(1)(Ω1)

R2

∑ S1= 1

R1−

∂Φ pol , E(r2) ∂r2

TS(1) Z (1) 1,0 S1,0

Smax

(27)

εm

Smax



(1)

(i) (i) (i) ] + ExX1,1 + EyY1,1 ∑ a(i)[EzZ1,0 i=1

(30)

The second term in eq 26 can be computed from eqs 4−8 (the factor 1/2 has to be suppressed in eq 4; see section C.5 in the SI for details), with the substitution σt(i) → σE(i) (i = 1 and 2). The final result is F

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

creased. Analysis of the eigenvectors (rows of the matrix [M] in eq 23) shows that the expansion coefficients of the uncoupled dynamical variables Vk,m as a function of the 2S d initial dynamical variables obey specific relationships, namely Mk, S1= S

of the derivation are reported in section C.6 of the SI. One obtains σ(ω) =

4πω Im[α(ω)] c εm

4πω c εm

=

2Smax



βk′

k=1

ωΓ − ω 2)2 + ω 2 Γ 2

(ωk2,1

= −( − 1)S + m Mk, S2 = S for the bonding branches and Mk, S1= S = ( − 1)S + m Mk, S2 = S for the antibonding branches. Since (34)

(2) TS(1) , m = TS, m for homodimers, it is concluded that the antibonding m = 0 and bonding m = 1 modes are not coupled to the electromagnetic field (Ck′ = 0 in eq 35; dark modes), whereas the bonding m = 0 and antibonding m = 1 modes are coupled to the electromagnetic field (Ck′ ≠ 0 in eq 35; bright modes), as for jellium-sphere homodimers. By way of illustration, Figures 2 and 3 display the evolutions of the

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

where α(ω) is the dynamical polarizability defined by Dx(t) = α(ω)Ex(t) and Γ is the damping parameter characterizing the friction force that has been added in the equations of motion for ensuring energy exchanges between the dimer and the applied electromagnetic field (see paper II and the SI). βk′ is the oscillator strength corresponding to the eigenmode k, namely βk′ =

⎡S 1 3 (1)2⎢ max (1) R1 ωP ∑ TS,1 Mk , S ⎢⎣ S = 1 3

ω (2) ⎛ R ⎞ + P(1) ⎜ 2 ⎟ ωP ⎝ R1 ⎠

3/2 Smax

∑ S=1

⎤2

⎥ TS(2) ,1 Mk , Smax + S ⎥⎦

=

1 3 (1)2 R1 ωP Ck′ 3 (35)

where the Mk,S values are the elements of the k-th row of the orthogonal matrix [M] (defined in eq 23). βk′ and the dimensionless coefficients Ck′ will be referred to as the “oscillator strength” and “normalized oscillator strength” of the mode k, respectively, throughout this work. Strictly identical results are obtained for the m = 1 (variables Y, polarization along the y axis) case, since the x and y axes are equivalent. The m = 0 (variables Z, polarization along the z axis) case leads to the same formal equations, with TS(,0i) instead of TS(,1i) and the involvement of different [M]-matrix coefficients. Since the matrix [M] is orthogonal, we have 2Smax

S=

∑ βk′ = k=1

S

1 (1)2 3 max (1) 2 ωP R1 ∑ (TS,1 ) 3 S=1

Figure 2. Evolution of the plasmon energies ℏωk±≡ S, m of a homodimer (Ri = 10 nm; εi = 4.64; ωP(i) = 9 eV; εm = 2.25) as a function of the intersphere distance d. For large d the plasmon energies tend toward the multipolar plasmon energies of the individual isolated spheres,

Smax

+

1 (2)2 3 2 ωP R 2 ∑ (TS(2) ,1 ) 3 S=1

(36)

2.3. Results, Range of Applicability, and Limits of the D-PH Model. Numerous calculations have been carried out for illustrating the main features of the D-PH formalism described in the previous subsections. In particular the range of applicability and limits of the modeling have been carefully scrutinized. One or several examples have been selected for each point investigated. Most of them have been reported on in the SI. Moreover, the reader will be referred in some places to paper II41 because the involvement of underlying dielectric media does not change fundamentally the main features of the PH modeling. (2) First we have checked that, for homodimers (ω(1) P = ωP , ε1 = ε2, and R1 = R2), the concepts of bonding and antibonding, and bright and dark modes hold true in the presence of underlying dielectric media. In this case, for a given m value, each S order multipolar plasmon frequency of the isolated spheres (the square frequencies are given on the right hand side of eq 19) gives rise to a lower steadily decreasing bonding frequency branch, ωk−≡ S, m(d), and an upper steadily increasing antibonding frequency-branch, ωk+≡ S, m(d), when the interparticle is de-

( ((

)

namely ℏωS = ℏωP(i) S / S + 1 εm + Sεi

1/2

))

. (a) Four lowest bright

(red curves) and dark (black curves) plasmon frequency branches ℏωk,m(d) for m = 0. (b) Same as (a) for m = 1.

eigenfrequencies ωk,m and associated normalized oscillator strengths Ck′ as a function of the interparticle distance d for the matrix-embedded homodimer characterized by the parameter set (Ri = 10 nm; εi = 4.64; ω(i) P = 9 eV; εm = 2.25). The d-evolutions of the plasmon frequencies and of the oscillator strengths, in particular the differences observed between the m = 0 and m = 1 d-evolutions, are similar to those noticed in the case of jellium-sphere dimers. So the physical analysis and comments reported on in the SI of paper II (pages 22−23) are directly applicable to the present case. Two main differences deserve to be pointed out. First, the spectral range of the plasmon excitations is lower [the ωP value, 9 eV, appropriate for silver and gold, is the same as that assumed in paper II for illustrating the jellium-sphere D-PH formalism)], as expected, due to the screening effects from the underlying polarizable backgrounds. Second, the sum of the G

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C

resonances (ℏΓ). Figures 1−3(SI) in the SI illustrate this striking phenomenon in the case of a slightly asymmetric dimer. Indeed, the overall frequency branch pattern of heterodimers, as well as the transfer and sharing of the oscillator strengths induced by the hybridization mechanism, depend strongly on both metals involved. Figures 4(SI)−6(SI) in the SI illustrate the case of a heterodimer involving two noticeably different (2) metals (ω(1) P = 9 eV, ε1 = 4.64, R1 = 10 nm, ωP = 6 eV, ε2 = 3.5, R2 = 13 nm, εm = 1.5). In this case, for each m value, the overall frequency branch pattern consists of two well separated bunches of frequency branches, with each of them converging for large interparticle distances toward the multipolar excitation energies of a given isolated metal sphere. In this specific example the lower bunch corresponds to metal 2, and the upper bunch to metal 1. Moreover the respective d-evolutions of the frequency branches and normalized oscillator strengths corresponding to each bunch are strongly different. All the lower frequency branches exhibit a bonding behavior when the interparticle distance is decreased, and the oscillator strength transfers to high-order modes induced by the hybridization mechanism are large, sequential, and rather selective, for both the axial and transverse excitations (m = 0 and m = 1 modes, respectively). In stark contrast, the upper branches are flat and, concomitantly, the oscillator strength transfers are very weak, whatever the distance d is, and are less selective, and the increase of the oscillator strengths of the high-order hybridized modes with decreasing distance is slow and monotonous. These illustrative examples show that manyqualitatively different cases can be encountered in the case of heterodimers, making it difficult to draw general conclusions and predict general trends, contrary to the homodimer case. The range of applicability and limits of the electrostatic DPH model have been investigated in comparing, through absorption and extinction cross section spectra, the D-PH model predictions with “exact” MEM calculations (retardation effects fully included).46,47 This powerful method consists in expanding the involved electromagnetic fields in terms of vector spherical harmonics (VSH) and determining the unknown expansion coefficients in imposing the boundary conditions on the dielectric interfaces. Provided that the selected z-axis coincides with the dimer axis (this axis will specify the azimuthal “quantum number” m of the VSH’s), the MEM formalism is diagonal with respect to the m value. Consequently, the scattering, absorption, and extinction cross sections are the sums of the various m contributions for any incoming applied external electromagnetic field. Therefore, each m ≥0-eigenfrequency and associated oscillator strength set in the D-PH model can be directly compared to the corresponding |m| contribution to the exact MEM absorption spectrum (sum of the m and −m contributions for m ≠ 0). A home-written code developed in the course of a previous work 16 has been used, and the convergence of the computations has been carefully checked in varying the maximum S order value in the field-expansions. For consistency the MEM computations have been carried out in assuming full metal dielectric functions of the form ε(i)(ω) = εi − ω(i)2 P /[ω(ω + iΓ)] for describing the dielectric properties of both metal spheres. In this formula the second term is the electric susceptibility (Drude−Sommerfeld model), corresponding to the conduction electrons, namely the electron fluids of the PH model. A very small damping parameter ℏΓ (0.002 eV) has been selected in order to clearly separate the various hybridized plasmon excitation lines in the spectra.

Figure 3. Evolution of the normalized oscillator strengths Ck′ (defined in eq 35) of the bright plasmon modes as a function of the intersphere distance d (see the caption of Figure 2 for details). When the distance d is decreased, the oscillator strengths of the degenerate dipolar excitations (S = 1,m) of both isolated spheres are progressively transferred to high-order hybridized multipolar modes S = k > 1 of the dimer system (for each m value), through the electrostatic coupling of the charge distributions located on the surfaces of the two spheres. In each panel the blue curve is the sum of the normalized oscillator strengths Ck′.

oscillator strengths depends on the distance d (blue curves in Figure 3), more especially for the axial polarization (m = 0 related modes), contrary to the case of jellium-sphere dimers in vacuum (in this case TS(,1i) = δS,1 in eq 36). This feature is a direct consequence of the dependence on the distance of the surface polarization charge densities σE(i)(ri), in other terms of the effective external field to which the dimer is coupled. In Figure 2 the exact degeneracy of the plasmon modes at the frequency branch crossings is the direct signature of the symmetry of the nanosystem with respect to the transverse x−y plane cutting the z-axis at equal distance from the facing spherical surfaces. In heterodimers (different radii or/and metals involved) this symmetry is broken and the expansion coefficients of the eigenvectors do not obey the above-written relationships. As a consequence, all the modes can be coupled to the external electromagnetic field. In such asymmetric dimers, nearby plasmon frequency branches repel from each other and avoided crossings are exhibited in the overall energy branch pattern.12,41 This phenomenon is particularly evidenced when the symmetry is only very slightly broken. In such cases the physical characteristics of nearby frequency branches, in particular the bright and dark characters, are shared and exchanged at the avoided crossings, resulting in broad and asymmetric pseudoresonances (spectral overlap of two plasmon resonances) in absorption spectra (see section 2.4 in paper II). Obviously, the overall bandwidth and asymmetry of the pseudoresonances depend on the degree of asymmetry of the dimer and the intrinsic width of the individual plasmon H

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C This comparative analysis, carried out first for several nanosized matrix-embedded dimers, has proved that the electrostatic D-PH model described in subsections 2.1 and 2.2, although of different nature compared to more conventional electrodynamics formalisms, is an accurate quantitative approach within the quasistatic limit. Both the plasmon energies and oscillator strengths are found in perfect agreement with their MEM counterparts. Moreover, these computations show that involving the parameters (ωP(i),εi, Γ) in the “mechanical/electrostatic” D-PH model amounts to assuming full metal dielectric functions with identical parameters in more conventional electrodynamics methods. Figures 7−8(SI) in the SI illustrate this study in the case of a nanosized matrixembedded heterodimer. Besides, the perfect agreement found between the D-PH model and the MEM predictions proves that the taking into account of the polarization charge densities induced by the bare external field, σE(i), is essential for defining the correct dimer/external field coupling term eq 31. Second, calculations involving larger dimers have been carried out in order to estimate the dimer sizes beyond which the electrostatic D-PH model necessarily fails (the quasistatic approximation is assumed in the PH formalism). For addressing this objective, I have proceeded as in paper II, namely in testing the scaling properties of the D-PH model. A simple examination of the various formulas shows that, due to the scale invariance of the dimensionless coupling integrals Va(S 1,m,S 2) (eq 25 in the SI), the plasmon eigenfrequencies, the eigenvectors, and the quantities TS(i,im) , eqs 32−33, are unchanged if the geometrical parameters, namely R1, R2, and d, are multiplied by a common scaling factor α. On the other hand, the absorption cross section scales as α3 (see eq 35). So, keeping constant the ratios d/R1 and R2/R1, a series of MEM extinction cross section spectra have been computed for increasing radius R1. These investigations have shown that the predictions of the electrostatic D-PH model may lead to incorrect spectral analyses (quantitatively and also qualitatively) if applied to dimers of overall sizes exceeding a few tens of nanometers. All the discrepancies with the MEM results (redshift and broadening of the plasmon resonances, reordering of the eigenmodes on the frequency scale, emergence of the dark modes) stem from the size- and mode-dependencies of the retardation effects, which are completely disregarded in the PH approach. These drawbacks are magnified in the strong coupling regime (small d/R values) because, in this regime, many hybridized modes can contribute to the absorption/ extinction spectra in narrow spectral ranges, making therefore difficult the labeling of multimodal spectral structures. Figure 9(SI) in the SI illustrates this study in the case of a homodimer (radius R ranging from 1 nm up to 20 nm). For the two smallest R values, the extinction is exhausted by the absorption and the spectra scale as R3, as predicted by the D-PH model, but strong departures from the D-PH predictions are observed for the largest dimer sizes. We refer the reader to the subsection 2.3 in paper II for a detailed discussion of the effects observed (the comments are directly applicable to the present case).

Figure 4. Geometry of the metal sphere/dielectric interface system, and notations and coordinate systems used. The vertical red line (at equal distance between the real sphere and its image) symbolizes the plane interface. 2d is the distance between the real sphere (sphere 2 in the left) and its image (sphere 1 on the right). σ(i) and σt(i) are the surface free-charge and total-charge densities, respectively.

electrostatic metal sphere−dielectric interface interaction to be described as the interaction of the metal sphere with a fictitious mirror sphere of identical radius located below the interface, the enslaved surface total-charge density of which is proportional to that of the real metal sphere. In order to take advantage of the formulas and ingredients involved in the DPH model, the same notations are used: sphere 2 will play the role of the real metal sphere, and sphere 1, that of the fictitious mirror sphere. The dimensionless relative electric permittivities of the three underlying dielectric media, assumed homogeneous and isotropic, are ε2, εm, and εsub. The surface free- and total-charge densities in the real sphere, σ(2)(Ω2) and σt(2)(Ω2), are expanded according to eq 1. The surface total-charge density in the image sphere is related to that of the real sphere through the equations (see eqs 27 and 31 in paper II)41 σ t(1)(θ1 , ϕ1 , t ) = −ξσ t(2)(θ2 = π − θ1 , ϕ2 = ϕ1 , t )

(37)

S + m t (2) σSt,(1) σS, m (t ) m (t ) = − ξ( − 1)

(38)

with ξ = (εsub − εm)/(εsub + εm). The Lagrangian of the system writes as 1 1 L = T (2) − V (σ (2) , σ t(2)) − V (σ (2) , σ t(1)) (39) 2 2 where the two electrostatic terms are given by eq 6 (with i = 2) and eq 8, respectively. The linear relations linking the multipolar expansion coefficients of the surface free- and total-charge densities in the real metal sphere are determined thanks to the boundary condition expressing the discontinuity of the normal component of the electric displacement across the dielectric interface r2 = R2 due to the presence of the surface free charge density σ(2), eq 10. From eqs 12 and 38 one obtains (for each m and S2 values) σ ′(2) S2, m =

3. PLASMON HYBRIDIZATION MODEL FOR THE SPHERE/INTERFACE SYSTEM 3.1. Lagrangian, Equations of Motion, and Eigenmodes. Figure 4 displays the sphere/interface system geometry and the main ingredients entering the modeling. As emphasized in the Introduction, the method of images51 allows the

[(S2 + 1)εm + S2ε2] t(2) σ ′S2, m 2S2 + 1 1 − (ε2 − εm)ξ ∑ Va(S1, |m| , S2)( −1)S1+ m σ ′tS1(2) ,m 2 S1 (40)

where the multipolar components σ′ are defined as in eq 13. Assuming a maximum S value (Smax ) in the multipolar I

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C expansions of the densities, eq 40 defines a matrix equation for each |m| value (the matrix, of dimension Sd , does not depend on the sign of m), namely σ′(2) = [E]σ′t(2), where σ′(2) and σ′t(2) are column vectors. Let us note that the |m|-dependent symmetric (the symmetry is ensured by eq 26 in the SI) square matrix [E] defined above is different from the four submatrixes [Eij] involved in the D-PH model. As in section 2 we define the |m|dependent symmetric matrix [χ] as the inverse of [E]. Using eqs 2, 6, 8, and 38 and the definition of the matrix [χ], one obtains

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

1 V (σ (2) , σ t(2)) 2 mn S2 = e 2 ωP(2)2 ∑ ( −1)m χS , S′ S S(2) S (2) 2, m S2 ′ , −m 2 2 2 2 1 + S 2 m , S2, S2 ′

where the coefficients of the column vectors are ordered according to [X S(2) ] = [X S(2) , X S(2) , ..., X S(2) ]T 2, m min, m min + 1, m max, m

The m-dependent real and symmetric square matrix [C] can be diagonalized thanks to standard numerical algorithms. The Sd eigenvalues ω2k,m (k = 1, 2, ..., S d) are real, and the corresponding eigenvectors {Vk,m; k = 1, 2, ..., S d} can be chosen real and orthonormalized. The uncoupled dynamical variables are } therefore related to the initial dynamical variables {X S(2) 2, m through a real orthogonal matrix [M] (m-dependent square matrix of dimension Sd such as [M]−1 = [M]T) in which the k-th row consists of the expansion coefficients of the uncoupled dynamical variables Vk,m as a function of the Sd initial dynamical variables, namely [hereafter [Vk,m] is a column vector containing the Sd uncoupled dynamical variables]

(41)

1 1 men2 (2)2 V (σ (2) , σ t(1)) = −ξ ωP 2 2 2



×

m , S2, S ′2 , S1

[Vk , m] = [M ][X S(2) ] 2, m

( −1)m ( −1)S1+ m Va(S1, |m| , S2)χS , S ′ S S(2) S (2) 2, m S2 ′ , −m 1 2

[Vk̈ , m] + [M ][C ][M ]−1 [Vk , m] = [Vk̈ , m] + [ωk2, m][Vk , m] = 0

In terms of the dynamical variables defined in eqs 15−17, the Lagrangian, eq 39, writes as (with the convention: for m = 0, X is replaced by Z, and Y = 0)

Smax

×

Smax

Smax

∑ ∑ m = 0 S2 = Smin

2 [Ẋ S22, m + YṠ 2, m] −

Smax

Smax

∑ ∑



m = 0 S2 = Smin S2 ′= Smin

(49) 2 where [ωk,m ] is a diagonal square matrix, the diagonal coefficients of which (eigenvalues of the matrix problem) are the square plasmon frequencies of the metal sphere/interface system. 3.2. Coupling to a Monochromatic Electromagnetic Field within the Quasistatic Approximation. The real metal sphere is now subject to an external spatially homogeneous monochromatic electric field Ei(t) = Ei0e−iωt (incoming plane wave). As discussed in paper II, the effects on the metal particle of the polarization charge distribution on the plane interface that is induced by the incoming applied electric field Ei(t) (polarization charge density in the absence of the real metal sphere) can be taken into account in a very simple way, namely in adding the reflected field ER(t) to the incoming field Ei(t) to define the external applied field irradiating the metal sphere. In the present context E(t) = Ei(t) + ER(t) = E0e−iωt is thus the bare external applied field. As in the D-PH model, the bare interaction can be derived from the scalar potential, eq 25. The electrostatic energy term corresponding to the particle/interface system−external field interaction (Vext) is given by eq 26, namely

n2me (2)2 ωP 2

A S2, S2 ′(X S(2) X (2) + Y S(2) Y (2) ) 2, m S2 ′ , m 2, m S2 ′ , m (43)

with A S2S2 ′ =

S2 χ 2S2 + 1 S2S2 ′ 1 + ( − ξ) 2

Smax

∑ S1= Smin

( −1)S1+ m Va(S1, |m| , S2)χS , S ′ 1 2

(44)

In the absence of underlying dielectric media, one has [χ] = [1] and the Lagrangian reduced to the one obtained in paper II for a jellium sphere in vacuum above a dielectric interface. Due to the axial symmetry of the system, the Lagrangian is diagonal with respect to the m value (m ≥ 0) and the symmetry (X or Y character for m > 0). In particular the equations of motion derived from the X and Y related Lagrangians will lead to strictly identical eigenvalue problems, as in the D-PH model. Since the m = 0 (variables Z) and m > 0 (variables X or Y) problems are uncoupled and lead to identical formal equations, only the equations of motion corresponding to the m > 0 (variables X) problem are reported on. The Euler−Lagrange equations derived from the Lagrangian, eq 43, expressed as (for each m and S2 values)

bare Vext = V ext +

∫ σ(2)(r2)Φ pol ,E(r2)d3r2

ωP(2)2 max [ ∑ (A S2, S2 ′ + A S2 ′ , S2 )X S(2) ]=0 2′,m 2 S ′= S 2

min

where Φ (r2) is the electrostatic potential generated by the surface polarization charge densities in the real and image53 spheres, that are directly induced by the bare external field E [denoted σE(i)(ri)]. We have therefore

∫ σ(2)(r2)Φ pol ,E(r2)d3r2 = V (σ(2), σ E(2)) + V (σ(2), σ E(1)) (51)

σSE, m(1)(t )

−ξ( −1)S + m σSE, m(2)(t ).

= with The polarization charge density σE(2)(r2) is determined in solving the electrostatic problem (matrix-embedded metal sphere above the plane dielectric interface in the presence of the external applied field E) in the absence of the free-charge density σ(2). As in the D-PH model, the density can be

(45)

In matrix form the above set of Sd second-order differential equations express as (for a given m) [Ẍ S(2) ] + [C ][X S(2) ]=0 2, m 2, m

(50)

pol,E

S

Ẍ S(2) + 2, m

(48)

The equations of motion of the variables Vk,m take the compact matrix form

(42)

nm L= 2 e 2

(47)

(46) J

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C motion for ensuring the energy particle/interface system and the field. In the βk expression, the Mk,l coefficients of Vk,1 as a function of

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

determined thanks to the boundary condition expressing the continuity of the normal component of the electric displacement across the dielectric interface r2 = R2, written according to the convenient form eq 28. Since the expression of Φpol,E(r2) as a function of the field-induced surface densities σE(i)(ri) is identical to that of Φt(r2) as a function of the surface totalcharge densities σt(i)(ri), we have (for each m value) the matrix equations Θ′(2) = [E]σ′E(2) and σ′E(2) = [χ]Θ′(2), where Θ′(2) and σ′E(2) are column vectors (the multipolar components are defined as in eq 13). The two square matrixes of dimension Sd have been defined previously in section 3.1, and the three nonvanishing dipolar components Θ(2) 1,m (m = 0, ±1) are given in section 2.2. All the ingredients required to compute the right-hand side of eq 51 have been therefore determined (Vbare ext is given by eq 30, with i = 2). The final result is (a(2) is defined in eq 30) Smax

variables XS(2) ,1 . After expressing the steady-state solution of eq 55, and using eq 48 and the orthogonality property of the matrix [M], the x-component of the dipole D defined in eq 52 can be expressed in the following generic form (holds for any nanosystem) Smax

Dx (t ) =

(56) S

max 1 1 βk ′ = R 23ωP(2)2[∑ TS,1Mk , S]2 = R 23ωP(2)2Ck′ 3 3 S=1

S2 = 1

Smax

+ a(2) ∑ TS2,1Y S(2) E = − D·E 2,1 y

S=

(52)

Smax

∑ (−)S + m Va(S1, m , S2)χS ,1 1

1

(53)

Smax

×

In the absence of dielectric media, TS2, m = δ S2,1, and the result eq 42 of paper II is recovered. In the presence of the external applied field E, −Vext has to be added to the Lagrangian eq 43, modifying therefore the equations of motion of the m = 0 and m = 1 related dynamical variables. As in the D-PH model, the problem remains diagonal with respect to the m value and the symmetry (X or Y character for m = 1), and the three independent sets consisting of the X, Y, and Z related modes will be specifically excited by electric fields linearly polarized along the three axes of the Cartesian reference frames. Assuming a polarization along the x axis (E0 = E0ex), the equations of motion corresponding to the m = 1 (variables X) problem become (for each S2 value)

=−

(54)

Using eqs 48−49, the equations of motion of the uncoupled variables are obtained straightforwardly Smax ⎡ (2) ⎤ a TS,1 2 ̈ ̇ Vk ,1 + ΓVk ,1 + ωk ,1Vk ,1 = −∑ ⎢ Mk , S ⎥Ex(t ) ⎢ ⎥⎦ S = 1 ⎣ me n 2

= βk Ex(t )

(58)

ωΓ (ωk2,1 − ω 2)2 + ω 2 Γ 2

(59)

Strictly identical results are obtained for the m = 1 (variables Y, polarization along the y-axis) case. The m = 0 (variables Z, polarization along the z-axis) case leads to the same formal equations, with TS,0 instead of TS,1 and the involvement of different [M]-matrix coefficients in βk′. 3.3. Results, Range of Applicability, and Limits of the PI-PH Model. By way of illustration, Figure 5 displays the m = 0 and m = 1 frequency branch pattern of a matrix-embedded metal sphere/interface system characterized by the parameter set (R2 = 10 nm; ε2 = 4.64; ω(2) P = 9 eV; εm = 2.25; εsub = 5.76). Figure 6 shows the 2d-evolutions of the normalized oscillator strengths C′k (defined in eq 57). The main features and trends are similar to those of the jellium-sphere/interface system. As emphasized in papers I and II, the differences with the D-PH model are basically rooted in the fact that the dynamics of the charge densities in the image sphere is enslaved to that of the densities in the real metal sphere, resulting in a reduction by a factor of 2 of the number of independent dynamical variables. As a consequence, for given m value and symmetry, each plasmon mode of the real sphere gives rise to a single bonding (that is, red-shifted) bright frequency branch which evolves steadily. Because of the reduction by a factor of 2 of the number of eigenmodes, and of the monotonous evolutions of the eigenfrequencies and oscillator strengths, the absorption spectra will evolve steadily as a function of the distance d, as in the case of nanosized homodimers. Actually this property will be observed too when the scaling property of the particle/interface system will be investigated (no drastic reordering of the modes on the frequency scale).

ωP(2)2 max [ ∑ (A S2, S2 ′ + A S2 ′ , S2 )X S(2) ] 2 ′ ,1 2 S2 ′= Smin a(2) TS ,1Ex(t ) men2 2

∑ βk′ k=1

S

Ẍ S(2) + 2,1

1 3 (2)2 max R 2 ωP [∑ (TS,1)2 ] 3 S=1

Finally, the absorption cross section is obtained in dividing the time-averaged absorbed power Pabs = 1/2ωIm[α(ω)]E20 by the electromagnetic field intensity in the matrix I0 = cE20(εm)1/2/ (8π) 4πω 4πω Im[α(ω)] = σ(ω) = c εm c εm

S2 1 (ε2 − εm)χS ,1 + ξ (ε2 − εm) 2 2S2 + 1 2

S1= 1

S

∑ βk′ = k=1

with (with m = 0 or 1, keeping in mind that the square matrix [χ] depends on m)

×

(57)

Since the matrix [M] is orthogonal, we have

Smax S2 = 1

1 Ex(t ) = α(ω)Ex(t ) (ωk2,1 − ω 2) − iω Γ

with

Vext = a(2) ∑ TS2,0Z S(2) E + a(2) ∑ TS2,1X S(2) E 2,0 z 2,1 x

TS2, m = δ S2,1 −

∑ βk′ k=1

Smax

S2 = 1

exchanges between the applied electromagnetic values are the expansion the Smax initial dynamical

(55) S

max with βk = −(e/me)(4π/3)1/2R3/2 2 ∑S = 1 TS,1Mk , S . As in the D-PH model, a “friction force”, characterized by the damping parameter Γ, has been added in the equation of

K

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C

Figure 6. Evolution of the normalized oscillator strengths Ck′ (defined in eq 57) of the plasmon modes of the metal sphere/interface system as a function of the real sphere-image sphere distance 2d (see the caption of Figure 5 for details). When the distance d is decreased, the oscillator strength of the dipolar excitation (S = 1,m) of the isolated metal sphere is progressively transferred to high-order hybridized multipolar modes S = k > 1 of the sphere/interface system (for each m value), through the electrostatic coupling of the charge distributions located on the surfaces of the real sphere and its image. In each panel the thick black curve is the sum of the normalized oscillator strengths Ck′.

Figure 5. Evolution of the plasmon energies ℏωk ≡ S, m of a matrixembedded metal sphere (R2 = 10 nm; ε2 = 4.64; ω(2) P = 9 eV; εm = 2.25) located above a dielectric interface (εsub = 5.76) as a function of the real sphere-image sphere distance 2d (d is the separation between the real sphere and the interface). For large distances, the plasmon energies tend toward the multipolar plasmon energies of the metal

( ((

)

sphere, namely ℏωS = ℏωP(2) S / S + 1 εm + Sε2

1/2

))

. (a) Five

lowest plasmon frequency branches ℏωk,m(d) for m = 0. (b) Same as (a) for m = 1.

electric field, E = Ei + ER [assumed homogeneous in the volume of the real sphere since the PI-PH model is grounded on the quasistatic approximation] divided by |Ei0|2, where Ei0 is the amplitude of the incoming electromagnetic field Ei. The weighting factors az and ax for the axial and transverse excitations are given in eqs 49−50 in paper II [the weighting factor cos2(θ) = sin2(θ) = 0.5 has been involved in the dimer case for comparing the MEM and D-PH model spectra displayed in Figures 7(SI) and 8(SI) in the SI]. Typical examples illustrating these investigations are provided in the SI [Figures 10−12(SI)]. First, the predictions of the PI-PH model described in this work are found to agree perfectly with the exact MEM results within the quasistatic approximation [Figures 10−11(SI)] with regard to both the plasmon frequencies and the peak intensities, proving again that the taking into account of the polarization charge densities induced by the bare external field E = Ei + ER, namely σE(2) and its enslaved image σE(1), is essential for defining the correct coupling term (eq 52). It is worthwhile emphasizing that, for many “simple” systems, such as sphere pairs,56−59 sphere/interface systems,60 and, obviously, single spheres and nanoshells,7,8 the electrostatic problems have been solved for a long time by several authors in the frame of classical electrodynamics (in solving the Laplace equation). Despite this fact, the hydrodynamic PH is of tremendous usefulness because, except for the simplest systems, the formalism is technically simpler61 and provides, in a simple way, new physical insights and concepts, such as, for instance,

We refer the reader to paper II for additional relevant comments regarding (i) the dependence of the magnitude of the interface effects on the electric permittivity of the substrate (the relevant parameter is the ratio εsub/εm) and (ii) the (rather slight in this example) differences observed between the devolutions of the eigenfrequencies and of the oscillator strengths of the m = 0 and m = 1 modes. Numerical tests have been carried out: first, in order to assess the accuracy of the PI-PH model within the quasistatic limit and, second, to determine its range of applicability for large spheres. As in section 2.3, absorption cross sections computed in the frame of the PI-PH model, eq 59, have been compared against results computed using a home-developed code (see the SI in paper I),40 based on an exact MEM formalism.54,55 A Drude−Sommerfeld dielectric function has been taken in the MEM calculations for describing the dielectric properties of the electron fluid in the real sphere, and a very small damping parameter Γ has been chosen to separate the plasmon excitation lines on the energy scale. In the MEM calculations, a p-polarized light excitation at the oblique incidence θ = 45° with respect to the z-axis has been assumed to ensure the simultaneous excitation of the m = 0 modes (axial modes; ez polarization) and m = 1 modes (transverse modes, e x polarization). As explained in paper II, for allowing a direct quantitative comparison with the m = 0 and |m| = 1 contributions to the MEM spectra, the PI-PH formula, eq 59 (evaluated for the ez and ex polarizations), has to be multiplied by the square z- or x-components of the total bare irradiating L

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C

It should be stressed that, if eq 61 is used, the external fieldinduced polarization charge density ρE is no longer involved in the PH modeling. Actually, the alternative intuitively attractive formula eq 61 leads to incorrect oscillator strengths, but surprisingly, although eqs 60 and 61 differ radically from each other, the incorrect values predicted by eq 61 have been found to be related to the correct ones through a mere modeindependent multiplicative factor. In order to address this issue, we have investigated analytically the relationship between both formulas for the three following nanosystems: (i) a single solid sphere, (ii) a solid sphere dimer, and (iii) the solid sphere/ interface system. 4.1. Single Solid Sphere. The PH formalism for this very simple system is described in detail in section B of the SI. Assuming eq 60, one obtains (eqs 60, 61, and 74 in the SI)

the bonding/antibonding and bright/dark characteristics of the modes. Moreover, the PH approach allows determining straightforwardly, for each collective excitation, (i) the electron fluid deformations of the global system in terms of those of the subunits (expansion coefficients of the uncoupled variables Vk,m) and (ii) the oscillator strengths of the modes (eqs 35 and 57; see the note ref 106 in Paper I). Actually, the PH picture is nowadays of almost systematic use by the “plasmonics community” for analyzing the optical properties of metal particle aggregates. The determination of the range of applicability of the PI-PH model has been addressed in testing its scaling property, in enlarging progressively the sphere radius R2 and the sphere− interface distance d, and in keeping constant the ratio d/R2 [Figure 12(SI); d/R2 = 0.1 with R2 ranging from 0.1 to 15 nm)]. Roughly, the same conclusions as those stated in the case of the dimer system apply. For very small radii, the absorption prevails and the extinction spectra scale as the real sphere volume. For larger sizes, a strong departure from this scaling law is observed, with regard to both the spectral location of the plasmon modes and the characteristics of the line-shapes (intensity and width). For this system the spectra are, however, less intricate, as compared to the dimer case, because the number of modes is two times smaller, and the evolution of the spectra for increasing system size is more regular, as noticed also in paper II.

Φtot , E(|r| = R , t ) =

3εm σ (R)(r)Φ E(|r| = R , t )d r 2εm + εi 3εm bare = V ext (t ) 2εm + εi

Vext(t ) =

D(t ) =

4. NANOSYSTEM/EXTERNAL FIELD COUPLING TERM: ADDITIONAL COMMENTS Throughout this work the electrostatic energy term Vext describing the nanosystem/external field interaction has been defined according to the general equation of electrostatics, namely Vext(t ) =

∫ ρ(r)Φtot ,E(r, t )d r

Vext(t ) =

(60)

∫ ρ (r)Φ (r, t )d r = − D ·E(t )

(63)

2 1 3 2⎛ 3εm ⎞ 1 R ωP ⎜ E(t ) ⎟ 2 2 3 2 ε + ε ( ) − iω Γ ω − ω ⎝ m i⎠ 1

(64)

where R is the sphere radius, εi the relative electric permittivity of the polarizable medium in the metal sphere, and ω1 = ωP/ (2εm + εi)1/2 the dipolar plasmon frequency. Since the (S = 1, m) components of the surface free- and total-charge densities are related by a simple multiplicative factor [σt(R) 1,m = (3/(2εm + εi))σ(R) 1,m; see eq 37 in the SI], the Vext and D(t) expressions obtained in using the alternative formula eq 61 underestimate the above ones by the factors εm and εm2, respectively. The absorption cross section is thus underestimated by the factor εm2. In this simple case, it is worthwhile to compare the induced dipole D(t) in eq 64 with the induced dipole DMie(t) that is obtained in the quasistatic Mie theory,7,8,63 assuming a full metal dielectric function of the form ε(ω) = εi − ω2P/[ω(ω + iΓ)]. DMie(t) expresses as (atomic units)

where ρ(r) is the free charge density and Φ (r,t) is the electrostatic potential generated by the total effective external field, that is, the sum of the bare applied external field E(t) = E0e−iωt and the field created by the surface polarization charge density ρE that is directly induced by the bare field (eqs 26 and 50 for the dimer and the sphere/interface systems, respectively). We have seen that the “dipole” D, more precisely the three linear combinations of the dynamical variables defined by the equation Vext(t) = −D·E(t) (eqs 31 and 52), leads, after solving the equations of motion of the uncoupled dynamical variables, to polarization-dependent dynamical polarizabilities αu(ω) [the induced dipole is D(t) = αu(ω)E(t); u = x, y, or z], the imaginary parts of which give the correct quasistatic absorption cross sections σabs,u(ω) (comparison with the exact MEM predictions). In ref 44 a different expression is assumed,62 namely eq 25 in ref 44 E

(62)



= α(ω)E(t )

tot,E

t

3εm Φ E(|r| = R , t ) 2εm + εi

DMie (t ) = εmR3

ε(ω) − εm E(t ) ε(ω) + 2εm

(65)

After some algebra DMie(t) can be cast in the following form DMie (t ) = D(t ) + εmR3

t

εi − εm E(t ) εi + 2εm

(66)

The second term on the right-hand side of eq 66, which does not contribute to the absorption, since the underlying dielectric media are assumed to be nonabsorbing, is patently the dipole corresponding to the polarization charge density σE(R) induced by the applied external field E(t). This analysis shows that the dipole D(t) that is computed in the PH formalism corresponds (for any nanosystem) to the contribution arising from the induced free charge densities and their corresponding enslaved induced polarization charge densities. 4.2. Solid Sphere Dimer. For this system eq 61 writes as

(61)

with Dt = ∫ rρt(r)dr (eq 61 defines three alternative linear combinations of the dynamical variables). In eq 61, ΦE(r,t) is the electrostatic potential created by the bare applied external field E(t) eq 25 and ρt(r) is the total charge density of the nanosystem in the absence of the external field [sum of the free charge density ρ(r) and the polarization charge density that is directly induced by ρ(r)]. At first sight, both formulas do not seem equivalent when polarizable dielectric media are present. M

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Vext(t ) =

∫ (σ t(1)(r) + σ t(2)(r))ΦE(r, t )d r

established are suitable for handling this case. In the presence of frequency-dependent permittivities, the price to be paid for computing, for instance, the absorption cross section of a dimer over a broad spectral range, is the necessity of computing all the required ingredients for each frequency ω, namely the matrixes [χ], [Aij], [C], and [M], the eigenvalues of the matrix [C] (the ω2k,m’s), and the coupling parameters TS(i,im) . The formula eq 34 is directly applicable, but the oscillator strengths βk′, the 2 ’s entering the eigenvalues of the matric [C] (the ωk,m denominators), and, eventually, εm depend now on the frequency ω. As a matter of fact, in the case of ω-dependent permittivities, the steps consisting in diagonalizing the large matrix [C] and solving the equations of motion of the uncoupled dynamical variables Vk,m are not the best strategies from a computational point of view, and even are not relevant from a physical point of view. Indeed, in the case of strongly ωdependent electric permittivities, the eigenvalues and associated eigenvectors resulting from the diagonalization of the matrix [C(ω)] (built for a given frequency ω) are, for most ω values, not related to the plasmon modes of the nanosystem in hand, but to those of a specific fictitious (moreover ω-dependent) system having constant electric permittivities [namely the values εi = εi(ω) and εm = εm(ω)]. In practice, the plasmon modes (frequencies and oscillator strengths) of the system can only be determined accurately through the examination of the entire α(ω) or σ(ω) curves [probably the dark modes of homodimers will be lacking and accurate MEM calculations would be necessary for determining their spectral locations (see, for instance, the Figure 3 in paper II)]. So, in such a case, it is simpler, and less computer time-consuming, to determine the dynamical polarizability from the steady-state solutions of the equations of motion of the coupled dynamical variables, as explained hereafter for the dimer and particle/interface systems. Second, is the PH formalism applicable to the case of complex ionic background electric permittivities? For metals of interest in plasmonics, such as silver and gold, the ionic background permittivity εi(ω) is complex above the interband transition threshold, and this issue is of particular relevance. Actually, it can be conjectured that the PH modeling is applicable to this case too, since no condition on the electric permittivities of the ionic backgrounds is prescribed in the whole formalism, as in standard electrodynamics approaches. Obviously, all the ingredients entering the PH formalism are now complex quantities. For the nanosystems studied in the course of this work, perfect agreement between the PH and MEM (or Mie) predictions has been obtained in the quasistatic limit when complex frequency-dependent permittivities are considered (including the nanoshell case). It should be emphasized that, if complex electric permittivities are involved, the dynamical polarizability αE,u(ω) (with u = x, y, or z, depending on the polarization of the applied external field) associated with the surface polarization charge density ρE that is directly induced by the applied external field E(t) = E0e−iωteu = E(t)eu has to be added to the PH dynamical polarizability αu(ω) associated with the induced surface total-charge density ρt (see eq 66 in the case of a single solid sphere). In this more general case the absorption cross sections of the nanosystem are thus given by

(67)

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

After some algebra, one obtains a formula similar to eq 31, with (quantities denoted T′) ⎛ R 2 ⎞3/2 21 11 TS(1) ′ = χ ( m ) + ⎜ ⎟ χ1, S1 (m) 1, S1 1, m ⎝ R1 ⎠

(68)

⎛ R1 ⎞3/2 12 22 TS(2) ′ = χ ( m ) + ⎜ ⎟ χ1, S2 (m) 1, S2 2, m ⎝ R2 ⎠

(69)

Computations involving several embedded dimers have shown that the above quantities are numerically related to the expressions eqs 32 and 33 by a mere multiplicative factor (T′ = T/εm). As before, the induced dipole D(t) and the absorption cross sections σabs,u(ω) computed in using eq 67 are underestimated by the factor εm2. Indeed, this simple (and rather unexpected) relationship has been analytically proved in using the equation [E][χ] = [1] and the structure of the matrix [E] defined in eqs 11 and 12 (see section D in the SI, where are outlined the main steps allowing the raw expressions eqs 32−33 to be greatly simplified; see the final results eqs D41 and D42). 4.3. Solid Sphere/Interface System. For this system I have conjectured that the counterpart of eq 61 writes as (only the real metal sphere is involved) Vext(t ) =

∫ σ t(2)(r)ΦE(r, t )d r

(70)

One obtains again TS2, m′ = χ1, S (m) = TS2, m/εm (see section D in 2 the SI; eqs D52 and D53). In the course of this work I have also found that a similar relationship holds for nanoshells. In view of these results, it can be conjectured that, for any nanosystem, eq 61 has to be multiplied by the electric permittivity of the matrix for ensuring correct quantitative predictions.

5. GENERALIZATION TO COMPLEX FREQUENCY-DEPENDENT ELECTRIC PERMITTIVITIES In the preceding sections, as well as in previous PH-related works (to my knowledge; see, for instance, refs 43-44), real, that is nonabsorbing, frequency-independent electric permittivities have been assumed for exemplifying the PH formalism. For plasmonic metals, such as silver or gold, this assumption is a reasonable approximation in the red and infrared spectral ranges, far below their respective interband transition threshold (about 4 and 2 eV, respectively). Actually the PH approach can be straightforwardly generalized to the more general case where frequency-dependent, and even complex, electric permittivities are involved. In this section, I discuss the important points to be kept in mind for dealing with complex frequency-dependent permittivities. I describe also the advisable modifications to be taken into account, and new ingredients to be included in the PH formalism, for ensuring perfect agreement with quasistatic MEM calculations. First, it is clear that, in the case of real frequency-dependent ionic background and matrix permittivities [εi(ω) and εm(ω), respectively], the PH models described in sections 2 and 3 for the dimer and particle/interface systems, as well as that reported on in sections B of the SI for a solid sphere, can be directly applied. All the analytical expressions that have been

σu(ω) = N

4πω c εm(ω)

Im[αu(ω) + αE , u(ω)] (71) DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C For the three nanosystems, αE,u(ω) can be computed from the general equation (inferred from the equation established in section 4, namely D = εmDt)

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

αE , u(ω) = εm(ω)

1 E (t )



αx(ω) =

S

S

1

2

max 1 (1)2 3 max (1) ̂ (1) ωP R1 [ ∑ TS1,1 X ′S1,1 + K ∑ TS(2) X̂ ′(2) S2,1] 2,1 3 S =1 S =1

(76)

uρ E (r, ω , t )d r

with

(72)

⎡ X̂ ′(1) ⎤ ⎡ T (1) ⎤ ⎢ S1,1 ⎥ = [D]−1 ⎢ S1,1 ⎥ ⎢⎣ X̂ ′(2) ⎥⎦ ⎢⎣ KT (2) ⎥⎦ S2,1 S2,1

I summarize the basic equations allowing the absorption cross section to be computed, for the three nanosystems studied in this work. 5.1. Single Solid Sphere. All the equations reported on in section B of the SI are suitable in the case of arbitrary complex ω-dependent permittivity ε i(ω) and real ω-dependent permittivity εm(ω). The analytical expression of αu(ω) is given in eq 64 [note that the “dipolar frequency” ω1(ω) is also a complex ω-dependent quantity in this expression] and the expression of αE,u(ω) is given by the second term on the righthand side of eq 66. Both polarizabilities do not depend on u because of the spherical symmetry of the system. 5.2. Solid Sphere Dimer. In section 2 an identical damping parameter Γ had been implicitly taken for characterizing the absorption properties of the incompressible electron fluids in both spheres (see ref 64). This assumption was crucial for preserving the uncoupling of the dynamical variables Vk,m resulting from the diagonalization of the matrix [C], when Vext is added to the Lagrangian. Here I allow the Γ parameter to have different values in spheres 1 and 2 (Γ1 and Γ2, respectively). As previously, only the equations of motion corresponding to the m = 1 (variables X) problem are reported on. In matrix form the equations of motion of the X′(i) = (ni)1/2X(i) variables express as (see eq 22, and eqs 141−142 in the SI)

(77)

3/2 (1) and K = (ω(2) P /ωP )(R2/R1) . αE,x(ω) is computed from eq 72 E E(1) E(2) with ρ = σ + σ and u = x. The multipolar components of the surface densities σE(i)(r) are given by eqs 126−127 in the SI. One obtains (omitting the ω-dependence of the various quantities for lightening the formulas) ε 11 21 αE , x(ω) = m (ε1 − εm)[R13χ1,1 + (R1R 2)3/2 χ1,1 ] 3 ε 22 12 ] + m (ε2 − εm)[R 23χ1,1 + (R1R 2)3/2 χ1,1 (78) 3

keeping in mind that the elements of the matrix [χ] depend on m (m = 1 for the x- and y-component, m = 0 for the zcomponent). By way of illustration Figure 7 shows the results obtained in the case of a tiny matrix-embedded silver/gold heterodimer, using parameters appropriate for both metals. The interbandtransition electric permittivities [εi(ω) with i = 1 and 2] have

⎡ Ẍ ′(1) ⎤ ⎡ ⎡ X ′(1) ⎤ ⎤⎡ Ẋ ′(1) ⎤ ⎢ S1,1 ⎥ + ⎢ Γ1 0 ⎥⎢ S1,1 ⎥ + [C ]⎢ S1,1 ⎥ ⎢⎣ Ẍ ′(2) ⎥⎦ ⎣ 0 Γ2 ⎦⎢⎣ Ẋ ′(2) ⎥⎦ ⎢⎣ X ′(2) ⎥⎦ S2,1 S2,1 S2,1 ⎡ a(1)T (1) ⎤T (2) a(1)TS(1) a(2)TS(2) a(2)T1,1 1,1 max,1 max,1 ⎥ ⎢ = − , ..., − ,− , ..., − Ex(t ) ⎢⎣ me n1 me n1 me n2 me n2 ⎥⎦

(73)

In the first line of eq 73 the “damping-matrix” is diagonal with diagonal elements equal to either Γ1 (Smax first elements) or Γ2 (Smax last elements). The steady-state time evolution of the dynamical variables is of the form X ′(S i,1) (t) = X̃ ′(S i,1) e−iωt [E(t) =

E0e−iωtex = Ex(t)ex], and one obtains

i

i

⎡ X̃ ′(1) ⎤ S1,1 ⎥ [D]⎢ (2) ⎢⎣ X̃ ′ ⎥⎦ S2,1 ⎡ a(1)T (1) ⎤T (2) a(1)TS(1) a(2)TS(2) a(2)T1,1 1,1 max,1 max,1 ⎥ = ⎢− , ..., − ,− , ..., − E0 ⎢⎣ me n1 me n1 me n2 me n2 ⎥⎦

(74)

Figure 7. m = 1 and m = 0 contributions to the absorption cross section σabs of a matrix-embedded Ag/Au heterodimer characterized by the parameter set (silver sphere: R1 = 1 nm, ω(1) P = 9 eV, Γ1 = 0.1 eV, interband transition electric permittivity ε1(ω); gold sphere: R2 = 1 nm, ω(2) P = 9 eV, Γ2 = 0.07 eV, interband transition electric permittivity ε2(ω); intersphere distance d = 0.1 nm; εm = 1.44), for p-polarized light excitation at the oblique incidence θ = 45° with respect to the dimer axis, computed using the exact multipole expansion method [MEM; thin black curves] and within the D-PH model [thick gray dashed curves]. The black and gray curves are perfectly superimposed. The red and blue curves are the absorption in silver and gold, respectively.

The square complex matrix [D] (dimension 2Smax ) is given by ⎡ −iω Γ1 0 ⎤ [D] = [C ] − ω 2[1] + ⎢ ⎥ ⎣ 0 −iω Γ2 ⎦

(75)

where [1] is the unit matrix. Multiplying eq 74 by the inverse of [D] allows determining the amplitudes X̃ ′(Sii,1) and thus the induced dipole Dx(t) = αx(ω)Ex(t) (defined in eq 31). After some algebra the dynamical polarizability αx(ω) can be expressed in the following convenient form O

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

created by the additional polarization charge distributions, that are directly induced by the bare external applied field, has to be added to the bare external field for defining the effective external field experienced by the nanosystem and ensuring correct quantitative predictions. Both PH models have been described in great detail, and explicit analytical formulas are provided in the main text and the SI, allowing direct use and computer implementation. The predictions of the D-PH and PI-PH models have also been compared, through absorption and extinction cross section spectra, to results computed in the frame of the multipole expansion method (retardation effects fully included). This comparative analysis has proved that both PH models, although of different nature compared to more conventional electrodynamics formalisms, are perfect quantitative approaches within the quasistatic limit. However, for large overall dimer sizes, quantitative and qualitative discrepancies are noticed, due to retardation effects. These effects, disregarded in the PH approach, induce large modifications of the line-shapes (intensity and width) and red-shifts of the modes, with an eventual spectral mixing and reordering of the mode spectrum on the frequency scale (case of the dimer system). It should be emphasized that the limit of applicability of the PH approach depends strongly on the system in hand, in particular the strength of the interparticle electrostatic couplings. Roughly, beyond an overall size of the system exceeding 50−60 nm, caution has to be taken for analyzing and labeling spectral structures on the basis of the PH model predictions only, especially when modes of different symmetries are simultaneously excited (several modes may contribute to broad structures; see for instance Figure 9 in the SI). In a large part of this work, as well as in all previous PH-related works, real frequency-independent electric permittivities have been assumed for simplicity’s sake. I have shown that the PH approach can be indeed generalized to deal with underlying dielectric media characterized by arbitrary complex frequencydependent electric permittivities. This work and paper II show that the PH formalism gives, in the quasistatic limit, the exact absorption cross section (and thus eigenmodes) of any nanosystem including arbitrary underlying homogeneous linear isotropic dielectric media and metal subunits involving incompressible free electron fluids, the optical properties of which are appropriately described by a Drude−Sommerfeld model susceptibility.

been extracted from the experimental complex refractive indexes of bulk silver and gold through a procedure described in a previous work.65 For the axial excitation (lower panel) one can note a clear Fano profile in the gold absorption contribution due to the coupling between the narrow silver LSPR (the “discrete state” in the Fano theory) and the continuum of gold interband transitions, as predicted theoretically.66 The perfect superimposition of the MEM and D-PH model curves (black and gray curves) proves that the PH formalism is a perfect numerical tool in the quasistatic limit, for any involved underlying dielectric media. 5.3. Solid Sphere/Interface System. The PH modeling [m = 1 (X variables) problem] is elaborated in following the same procedure as previously. In matrix form (the square matrixes are of dimension Smax ) the amplitudes of the variables Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

X S(2) express as (see eqs 45, 46, and 54) 2,1 [X̃ S(2) ] 2,1

T ⎡ a(2)T a(2)TSmax,1 ⎤ 1,1 ⎥ = [D] − , ..., − E0 ⎢⎣ men2 men2 ⎥⎦ −1 ⎢

(79)

where the matrix [D] is defined as in eq 75, namely [D] = [C] − (ω2 + iωΓ2)[1]. One obtains from the expression of the component Dx defined in eq 52 αx(ω) =

S

S

2

2

max 1 (2)2 3 max ωP R 2 ∑ TS2,1[ ∑ DS−2,1S2 ′TS2 ′ ,1] 3 S =1 S ′= 1

(80)

αE,x(ω) is computed from eq 72 with ρ = σ and u = x. One obtains (omitting the ω-dependence of the various quantities) ε αE , x(ω) = m (ε2 − εm)R 23χ1,1 (m = 1) (81) 3 E

E(2)

6. CONCLUSION In this work the D-PH and PI-PH models reported on in paper II41 have been extended to the case where homogeneous underlying polarizable dielectric media are present in the full space. In the presence of underlying polarizable media, the PH formalism is indeed much more complex as compared to jellium spheres in vacuum because the surface polarization charge densities on both sides of each dielectric interface involved in the dimer and particle/interface systems, that are induced by the deformations of the incompressible electron fluids, contribute to the energetics of the systems, and thus to the Lagrangian and the resulting equations of motion. A central step in the formalism consists in relating the enslaved surface polarization charge densities to the surface free charge densities resulting from the deformations of the electron fluids. From the equations of motion deduced from the Lagrangian, expressed in terms of the dynamical variables selected for describing the electron fluid deformations, result an eigenvalue matrix problem (similar to that corresponding to a large set of onedimensional coupled oscillators), the diagonalization of which yields, in the case of frequency-independent electric permittivities, the eigen-frequencies and eigenmodes (uncoupled oscillators) of the dimer or particle/interface system. The coupling of both systems to an external applied electromagnetic field has been handled in using the method described in paper II, allowing the absorption cross section to be expressed in a compact analytical form, as the sum of the individual contributions from all the eigenmodes of the dimer or particle/interface system. I have also shown that the field



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b06129. Basic formulas; plasmon hybridization model for a single sphere; plasmon hybridization model for a sphere dimer; simplification of the raw formulas; illustrative examples (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Kelsall, R.; Hamley, I. W.; Geoghegan, M. Nanoscale Science and Technology; Wiley: New York, 2006. P

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C (2) Giannini, V.; Fernandez-Dominguez, A. I.; Heck, S. C.; Maier, S. A. Plasmonic Nanoantennas: Fundamentals and Their Use in Controlling the Radiative Properties of Nanoemitters. Chem. Rev. 2011, 111, 3888−3912. (3) Halas, N. J.; Lal, S.; Chang, W. S.; Link, S.; Nordlander, P. Plasmons in Strongly Coupled Metallic Nanostructures. Chem. Rev. 2011, 111, 3913−3961. (4) Mayer, K. M.; Hafner, J. H. Localized Surface Plasmon Resonance Sensors. Chem. Rev. 2011, 111, 3828−3857. (5) Stewart, M. E.; Anderton, C. R.; Thompson, L. B.; Maria, J.; Gray, S. K.; Rogers, J. A.; Nuzzo, R. G. Nanostructured Plasmonic Sensors. Chem. Rev. 2008, 108, 494−521. (6) Jain, P. K.; El-Sayed, I. H.; El-Sayed, M. A. Au nanoparticles target cancer. Nano Today 2007, 2, 18−29. (7) Kreibig, U.; Vollmer, M. Optical Properties of Metal Clusters; Springer: Berlin, 1995. (8) Bohren, C. F.; Huffman, D. P. Absorption and scattering of light by small particles; Wiley: New York, 1983. (9) Kelly, K. L.; Coronado, E.; Zhao, L. L.; Schatz, G. C. The Optical Properties of Metal Nanoparticles: The influence of Size, Shape, and Dielectric Environment. J. Phys. Chem. B 2003, 107, 668−677. (10) Wiley, B. J.; Im, S. H.; Li, Z. Y.; McLellan, J.; Siekkinen, A.; Xia, Y. Maneuvering the Surface Plasmon Resonance of Silver Nanostructures through Shape-Controlled Synthesis. J. Phys. Chem. B 2006, 110, 15666−15675. (11) Liz-Marzan, L. M. Nanometals: formation and color. Mater. Today 2004, 7, 26−31. (12) Nordlander, P.; Oubre, C.; Prodan, E.; Li, K.; Stockman, M. I. Plasmon Hybridization in Nanoparticle Dimers. Nano Lett. 2004, 4, 899−903. (13) Romero, I.; Aizpurua, J.; Bryant, G. W.; Garcia de Abajo, F. J. Plasmons in nearly touching metallic nanoparticles: singular response in the limit of touching dimers. Opt. Express 2006, 14, 9988−9999. (14) Jain, P. K.; Huang, W.; El-Sayed, M. A. On the Universal Scaling Behavior of the Distance Decay of Plasmon Coupling in Metal Nanoparticle Pairs: A Plasmon Ruler Equation. Nano Lett. 2007, 7, 2080−2088. (15) Grillet, N.; Manchon, D.; Bertorelle, F.; Bonnet, C.; Broyer, M.; Cottancin, E.; Lermé, J.; Hillenkamp, M.; Pellarin, M. Plasmon Coupling in Silver Nanocube Dimers: Resonance Splitting Induced by Edge Rounding. ACS Nano 2011, 5, 9450−9462. (16) Marhaba, S.; Bachelier, G.; Bonnet, C.; Broyer, M.; Cottancin, E.; Grillet, N.; Lermé, J.; Vialle, J. L.; Pellarin, M. Surface Plasmon Resonance of Single Gold Nanodimers near the Conductive Contact Limit. J. Phys. Chem. C 2009, 113, 4349−4356. (17) Brandl, D. W.; Mirin, N. A.; Nordlander, P. Plasmon Modes of Nanosphere Trimers and Quadrumers. J. Phys. Chem. B 2006, 110, 12302−12310. (18) Mirin, N. A.; Bao, K.; Nordlander, P. Fano Resonances in Plasmonic Nanoparticle Aggregates. J. Phys. Chem. A 2009, 113, 4028−4034. (19) Wang, H.; Brandl, D. W.; Nordlander, P.; Halas, N. J. Plasmonic Nanostructures: Artificial Molecules. Acc. Chem. Res. 2007, 40, 53−62. (20) Encina, E. R.; Coronado, E. A. Plasmon Coupling in Silver Nanosphere Pairs. J. Phys. Chem. C 2010, 114, 3918−3923. (21) Lassiter, J. B.; Sobhani, H.; Fan, J. A.; Kundu, J.; Capasso, F.; Nordlander, P.; Halas, N. J. Fano Resonances in Plasmonic Nanoclusters: Geometrical and Chemical Tunability. Nano Lett. 2010, 10, 3184−3189. (22) Funston, A. M.; Novo, C.; Davis, T. J.; Mulvaney, P. Plasmon Coupling of Gold Nanorods at Short Distances and in Different Geometries. Nano Lett. 2009, 9, 1651−1658. (23) Slaughter, L. S.; Wu, Y.; Willingham, B. A.; Nordlander, P.; Link, S. Effects of Symmetry Breaking and Conductive Contact on the Plasmon Coupling in Gold Nanorod Dimers. ACS Nano 2010, 4, 4657−4666. (24) Atay, T.; Song, J. H.; Nurmikko, A. V. Strongly Interacting Plasmon Nanoparticle Pairs: From Dipole-Dipole Interaction to Conductively Coupled Regime. Nano Lett. 2004, 4, 1627−1631.

(25) Jain, P. K.; El-Sayed, M. A. Universal Scaling of Plasmon Coupling in Metal Nanostructures: Extension from Particle Pairs to Nanoshells. Nano Lett. 2007, 7, 2854−2858. (26) Barrow, S. J.; Funston, A. M.; Gomez, D. E.; Davis, T. J.; Mulvaney, P. Surface Plasmon Resonances in Strongly Coupled Gold Nanosphere Chains from Monomer to Hexamer. Nano Lett. 2011, 11, 4180−4187. (27) Gunnarsson, L.; Rindzevicius, T.; Prikulis, J.; Kasemmo, B.; Käll, M.; Zou, S.; Schatz, G. C. Confined Plasmons in Nanofabricated Single Silver Particle Pairs: Experimental Observations of Strong Interparticle Interactions. J. Phys. Chem. B 2005, 109, 1079−1087. (28) Munechika, K.; Smith, J. M.; Chen, Y.; Ginger, D. S. Plasmon Line Widths of Single Silver Nanoprisms as a Function of Particle Size and Plasmon Peak Position. J. Phys. Chem. C 2007, 111, 18906−18911. (29) Sönnichsen, C.; Franzl, T.; Wilk, T.; von Plessen, G.; Feldmann, J. Plasmon resonances in large noble-metal clusters. New J. Phys. 2002, 4, 93.1−93.8. (30) Mock, J. J.; Barbic, M.; Smith, D. R.; Schultz, D. A.; Schultz, S. Shape effects in plasmon resonance of individual silver nanoparticles. J. Chem. Phys. 2002, 116, 6755−6759. (31) Lindfors, K.; Kalkbrenner, T.; Stoller, P.; Sandoghdar, V. Detection and Spectroscopy of Gold Nanoparticles Using Supercontinuum White Light Confocal Microscopy. Phys. Rev. Lett. 2004, 93, 037401. (32) Billaud, P.; Marhaba, S.; Cottancin, E.; Arnaud, L.; Bachelier, G.; Bonnet, C.; Del Fatti, N.; Lermé, J.; Vallée, F.; Vialle, J. L.; Broyer, M.; Pellarin, M. Correlation between the Extinction Spectrum of a Single Metal Nanoparticle and its Electron Microscopy Image. J. Phys. Chem. C 2008, 112, 978−982. (33) Sönnichsen, C.; Geier, S.; Hecker, N. E.; von Plessen, G.; Feldmann, J.; Ditlbacher, H.; Lamprecht, B.; Krenn, J. R.; Aussenegg, F. R.; Chan, V. Z. H.; Spatz, J. P.; Möller, M. Spectroscopy of single metallic nanoparticles using total internal reflection microscopy. Appl. Phys. Lett. 2000, 77, 2949−2951. (34) Kalkbrenner, T.; Hakanson, U.; Sandoghdar, V. Tomographic Plasmon Spectroscopy of a Single Gold Nanoparticle. Nano Lett. 2004, 4, 2309−2314. (35) Billaud, P.; Marhaba, S.; Grillet, N.; Cottancin, E.; Bonnet, C.; Lermé, J.; Vialle, J. L.; Broyer, M.; Pellarin, M. Absolute optical extinction measurements of single nano-objects by spatial modulation spectroscopy using a white lamp. Rev. Sci. Instrum. 2010, 81, 043101. (36) Crut, A.; Maioli, P.; Del Fatti, N.; Vallée, F. Optical absorption and scattering spectroscopies of single nano-objects. Chem. Soc. Rev. 2014, 43, 3921−3956. (37) Lermé, J.; Bachelier, G.; Billaud, P.; Bonnet, C.; Broyer, M.; Cottancin, E.; Marhaba, S.; Pellarin, M. Optical response of a single spherical particle in a tightly focused light beam: application to the spatial modulation spectroscopy technique. J. Opt. Soc. Am. A 2008, 25, 493−514. (38) Lermé, J.; Bonnet, C.; Broyer, M.; Cottancin, E.; Marhaba, S.; Pellarin, M. Optical response of metal or dielectric nano-objects in strongly convergent light beams. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 245406. (39) Mojarad, N. M.; Sandoghdar, V.; Agio, M. Plasmon spectra of nanospheres under a tightly focused beam. J. Opt. Soc. Am. B 2008, 25, 651−658. (40) Lermé, J.; Bonnet, C.; Broyer, M.; Cottancin, E.; Manchon, D.; Pellarin, M. Optical Properties of a Particle above a Dielectric Interface: Cross Sections, Benchmark Calculations, and Analysis of the Intrinsic Substrate Effects. J. Phys. Chem. C 2013, 117, 6383−6398. (41) Lermé, J. Nanoparticle above a Dielectric Interface: Plasmon Hybridization Model, Comparison with the Dimer System, and against Exact Electrodynamics Calculations. J. Phys. Chem. C 2014, 118, 28118−28133. (42) Prodan, E.; Radloff, C.; Halas, N. J.; Nordlander, P. A Hybridization Model for the Plasmon Response of Complex Nanostructures. Science 2003, 302, 419−422. (43) Prodan, E.; Nordlander, P. Plasmon hybridization in spherical nanoparticles. J. Chem. Phys. 2004, 120, 5444−5454. Q

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

Downloaded by UNIV OF MICHIGAN ANN ARBOR on September 4, 2015 | http://pubs.acs.org Publication Date (Web): August 27, 2015 | doi: 10.1021/acs.jpcc.5b06129

The Journal of Physical Chemistry C

(63) In the quasistatic Mie theory the field outside the metal sphere is the superposition of the applied external field and the field created by an ideal point-like dipole located at the origin [the corresponding electrostatic potential is DMie·r/(εmr3) (atomic units)]. (64) This restriction is implicitly assumed in previous works (refs 43−44) where the absorption is ensured in adding an imaginary part to the real frequency in the expression of the dynamical polarizability [αu(ω) → αu(ω + iδ)]. For small Γ values this procedure leads to absorption cross sections identical to those computed in the frame of the PH formalism described in this work, when 2δ = Γ. However, noticeable differences are observed when large damping parameters are involved because of the presence of a spurious additional term (δ2) in the energy denominators of eqs 34 and 59. Since the present approach leads to absorption spectra in perfect agreement with the (exact) MEM predictions, whatever the Γ value is, this approach, consisting in adding a frictional term in the equations of motion of the dynamical variables, is numerically more efficient. Moreover, the method allows attributing different parameters Γi in each metallic subunit. (65) Lermé, J.; Palpant, B.; Prével, B.; Cottancin, E.; Pellarin, M.; Treillleux, M.; Vialle, J. L.; Perez, A.; Broyer, M. Optical properties of gold metal clusters: A time-dependent local-density-approximation investigation. Eur. Phys. J. D 1998, 4, 95−108. (66) Bachelier, G.; Russier-Antoine, I.; Bénichou, E.; Jonin, C.; Del Fatti, N.; Vallée, F.; Brevet, P. F. Fano Profiles Induced by Near-Field Coupling in Heterogeneous Dimers of Gold and Silver Nanoparticles. Phys. Rev. Lett. 2008, 101, 197401.

(44) Brandl, D. W.; Oubre, C.; Nordlander, P. Plasmon hybridization in nanoshell dimers. J. Chem. Phys. 2005, 123, 024701. (45) In previous PH-related works, the absorption cross section is computed formally through matrix product formulas where all the physical aspects are buried. See, for instance, eq 72 in ref 43 and eq 30 in ref 44. Even in the simple case of a single jellium nanoshell, where the hybridization mechanism involves only two plasmon modes for each multipolarity (eq 72 in ref 43), the explicit expression is not physically transparent. (46) Mishchenko, M. I.; Hovenier, J. W.; Travis, L. D. Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications; Academic Press: San Diego, 2000. (47) Mishchenko, M. I.; Travis, L. D.; Lacis, A. A. Scattering, Absorption, and Emission of Light by Small Particles; Cambridge University Press: Cambridge, 2002. (48) The physical meaning of the terminology “the polarization charge density σpol (or σpol,E) is directly induced by the free charge density σ (resp. by the bare external field E)” is the following: σpol (respectively σpol,E) is the polarization charge density that is obtained when the Maxwell equations are self-consistently solved disregarding all the other free charge densities σ′ (respectively all the free charge densities). (49) Lermé, J.; Palpant, B.; Cottancin, E.; Pellarin, M.; Prével, B.; Vialle, J. L.; Broyer, M. Quantum extension of Mie’s theory in the dipolar approximation. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 60, 16151−16156. (50) Lermé, J. Introduction of quantum finite-size effects in the Mie’s theory for a multilayered metal sphere in the dipolar approximation: Application to free and matrix-embedded noble metal clusters. Eur. Phys. J. D 2000, 10, 265−277. (51) Jackson, J. D. Classical Electrodynamics, 3rd ed.; John Wiley & Sons, Inc.: New York, 1998. (52) Arfken, G. B.; Weber, H. J. Mathematical Methods for Physicists, 4th ed.; Academic Press: San Diego, 1995. (53) The induced polarization charge density σE(2), which develops on the surface of the real metal sphere polarizes the matrix and the substrate, giving rise to a complex surface polarization charge distribution on the plane matrix/substrate interface. The method of images [illustrated in the case of point charges in textbooks] states that the electrostatic potential in the matrix, which is created by the polarization charge distribution on the plane interface, is identical to that created by a charge density located on the surface of a fictitious mirror sphere located below the interface. (54) Fucile, E.; Denti, P.; Borghese, F.; Saija, R.; Sindoni, O. I. Optical Properties of a Sphere in the Vicinity of a Plane Interface. J. Opt. Soc. Am. A 1997, 14, 1505−1514. (55) Wriedt, T.; Doicu, A. Light Scattering from a Particle on or near a Surface. Opt. Commun. 1998, 152, 376−384. (56) Goyette, A.; Navon, A. Two dielectric Spheres in an Electric Field. Phys. Rev. B 1976, 13, 4320−4327. (57) Poladian, L. General theory of electrical images in sphere pairs. Q. J. Mech. Appl. Math. 1988, 41, 395−417. (58) Ruppin, R. Surface Modes of Two Spheres. Phys. Rev. B: Condens. Matter Mater. Phys. 1982, 26, 3440−3444. (59) Chu, P.; Mills, D. L. Electromagnetic Response of Nanosphere Pairs: Collective Plasmon Resonances, Enhanced Fields, and LaserInduced Forces. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 045416. (60) Ruppin, R. Surface Modes and Optical Absorption of a Small Sphere above a Substrate. Surf. Sci. 1983, 127, 108−118. (61) For instance, for the two-sphere and sphere/interface systems, the introduction of bispherical coordinates is required to derive the recursion relations allowing the expansion coefficients of the analytical solution of the Laplace’s equation to be determined. (62) In ref 43, devoted to single spherical layered particles (with or without underlying polarizable media), the Vext formula that is given ̀ E (r,t)dr, is suitable only to systems (eq 65), namely Vext(t) = ∫ ρ(r)Φ not including polarizable backgrounds. R

DOI: 10.1021/acs.jpcc.5b06129 J. Phys. Chem. C XXXX, XXX, XXX−XXX