A Functions g_{1}(\kappa) and g_{3}(\kappa)

# Theoretical and numerical investigation of the size-dependent optical effects in metal nanoparticles

## Abstract

We further develop the theory of quantum finite-size effects in metallic nanoparticles, which was originally formulated by Hache, Ricard and Flytzanis [J. Opt. Soc. Am. B 3, 1647 (1986)] and (in a somewhat corrected form) by Rautian [Sov. Phys. JETP 85, 451 (1997)]. These references consider a metal nanoparticle as a degenerate Fermi gas of conduction electrons in an infinitely-high spherical potential well. This model (referred to as the HRFR model below) yields mathematical expressions for the linear and the third-order nonlinear polarizabilities of a nanoparticle in terms of infinite nested series. These series have not been evaluated numerically so far and, in the case of nonlinear polarizability, they can not be evaluated with the use of conventional computers due to the high computational complexity involved. Rautian has derived a set of remarkable analytical approximations to the series but direct numerical verification of Rautian’s approximate formulas remained a formidable challenge. In this work, we derive an expression for the third-order nonlinear polarizability, which is exact within the HRFR model but amenable to numerical implementation. We then evaluate the expressions obtained by us numerically for both linear and nonlinear polarizabilities. We investigate the limits of applicability of Rautian’s approximations and find that they are surprizingly accurate in a wide range of physical parameters. We also discuss the limits of small frequencies (comparable or below the Drude relaxation constant) and of large particle sizes (the bulk limit) and show that these limits are problematic for the HRFR model, irrespectively of any additional approximations used. Finally, we compare the HRFR model to the purely classical theory of nonlinear polarization of metal nanoparticles developed by us earlier [Phys. Rev. Lett. 100, 47402 (2008)].

## I Introduction

This paper is dedicated to the memory of Sergey Glebovich Rautian (1928-2009) who was a teacher to some of us and inspiration to all.

Metal nanoparticles have received an extraordinary amount of attention recently because of their ability to greatly enhance local fields. The enhancement is attributed to the excitation of surface plasmons and it has a variety of applications in photovoltaics Huynh et al. (2002), sensing Sun and Xia (2002) and surface-enhanced Raman scattering Kelly et al. (2003); Nikoobakht et al. (2002); Hao and Schatz (2004). Currently, nanopartices of very small sizes, up to a few nanometers, are customarily used in experiments. The theoretical description of the optical properties of such nanoparticles is most frequently based on the macroscopic electrodynamics. At least, this is typical in the field of plasmonics. However, macroscopic electrodynamics can not capture certain effects of finite size. Hache, Ricard and Flytzanis Hache et al. (1986) and Rautian (in a somewhat modified form) Rautian (1997) have developed an elaborate theory of quantum finite-size effects in metal nanoparticles. In Refs. Hache et al., 1986; Rautian, 1997, a nanoparticle was modeled as degenerate Fermi gas confined in an infinite potential well of spherical shape (below, the HRFR model). Despite being fairly simple, the HRFR model results in very complicated formulas, which can not be evaluated numerically even with the aid of modern computers. For example, the expression for the third-order nonlinear polarizability involves a twelve-fold nested summation. Rautian has reduced the number of summations from twelve to eight by performing summation over the magnetic sublevels analytically; he then obtained a number of remarkable approximations to the resulting eight-fold summation Rautian (1997). However, these approximations have never been verified directly due to the overwhelming numerical complexity involved.

In this contribution, we develop the analytical theory of Rautian a step further by reducing the number of nested summations involved from eight to five without making any additional approximations. This turns out to be sufficient to render the formulas amenable to direct numerical implementation. We then compare the results of numerical evaluation of the five-fold series derived by us to the results, which follow from Rautian’s approximate formulas, and discuss various physical limits, including the limits of low frequency and large particle size.

In Sec. II we review Rautian’s theory. Here we use somewhat simplified notation and, in particular, avoid the use of irreducible spherical tensors and -symbols. In Sec. III, we develop the theory further by utilizing the orbital selection rules and reduce the nested summation involved in the definition of the third-order nonlinear polarizability from eight-fold to five-fold. In Sec. IV, we describe a simple method to relate the internal and applied fields, which is to the first order consistent with the approach proposed in Ref. Drachev et al., 2004a, but is more mathematically rigorous. In Sec. V, the results of numerical computations are reported. A summary of obtained results and a discussion are contained in Sec. VI.

## Ii Rautian’s theory

We start by reviewing Rautian’s theory of quantum finite-size effects in conducting nanoparticles Rautian (1997). The physical system under consideration is a gas of non-interacting electrons placed inside a spherical, infinitely deep potential well of radius and subjected to harmonically-oscillating, spatially-uniform electric field

 Ei(t)=Aiexp(−iωt)+c.c. (1)

Note that is the electric field inside the nanoparticle. It will be related to the external (applied) field in Sec. IV below.

Since the nanoparticle is assumed to be electrically small (that is, ), the electron-field interaction can be described in the dipole approximation by the time-dependent operator

 V(t)=−d⋅E(t) , (2)

where is the dipole moment operator. Under the additional assumption that is linearly polarized with a purely real amplitude , we can write

 V(t)=Gexp(−iωt)+c.c. , (3)

where

 G=−d⋅Ai . (4)

Rautian made use of the interaction representation in which the wave function is expanded in the basis of the unperturbed Hamiltonian eigenstates. The single-electron unperturbed states are

 ϕnlm(r)=1Znljl(ξnlr/a)Ylm(^r) , (5)

where are the spherical Bessel functions of the first kind and order ; is the -th positive root () of the equation , are the spherical functions (viewed here as functions of the polar and azimuthal angles of the unit vector ) and

 Znl=√a32jl+1(ξnl) (6)

are normalization factors. The energy eigenstates are labeled by the main quantum number , the orbital number and the magnetic number . The unperturbed energy levels are given by the formula

 Enl=E0ξ2nl , (7)

where

 E0=ℏ22mea2 (8)

and is the electron mass.

In what follows, we use the composite indices , , and to label the eigenstates. Each composite index corresponds to the triplet of quantum numbers . By convention, if , then . The matrix elements of the -projection of the dipole moment operator are given by

 ^z⋅dμμ′=eaΔμμ′ , (9)

where

 Δμμ′=δmm′Rn′l′nl(blmδl−1,l′+bl+1,mδl+1,l′) , (10)

are the Kronecker delta-symbols, and

 blm=√l2−m24l2−1 ,  Rn′l′nl=4ξnlξn′l′(ξ2nl−ξ2n′l′)2 . (11)

Note that the diagonal elements of are all equal to zero, as is the case for any system with a center of symmetry. Finally, the matrix elements of the operator are given by

 Gμμ′=−eaΔμμ′Ai . (12)

The density matrix of the system, , can be written in the form , where are the transition frequencies and is the so-called slow-varying amplitude, which obeys the following master equation Rautian and Shalagin (1991):

 (∂∂t+iωμν+Γμν)~ρμν=δμνΓμμNμ −iℏ∑η[Vμη(t)~ρην−~ρμηVην(t)] . (13)

Here are the equilibrium state populations and are phenomenological relaxation constants. Following Rautian, we assume that

 Γμν=Γ1δμν+Γ2(1−δμν) . (14)

Eq. (14) is the least complex assumption on , which still distinguishes the relaxation rates for the diagonal and off-diagonal elements of the density matrix.

It can be seen that, for the case of zero external field, . The Fermi statistics is introduced at this point by writing

 Nμ=2exp[(Eμ−EF)/(kBT)]+1 , (15)

where is the Fermi energy, is Boltzmann’s constant, is the temperature, and the factor of two in the numerator accounts for the electron spin. Conservation of particles reads . When , the well-known analytical formula for the Fermi’s energy,

 EF=E0(3π2)2/3(aℓ)2=(3π2)2/3ℏ22meℓ2 , (16)

holds with a good accuracy. Here is the characteristic atomic scale, defined by the relation

 ℓ3=Ω/N , (17)

where

 Ω=4πa3/3 (18)

is the nanoparticle volume. Thus, is the specific volume per conduction electron. We note that is, generally, different from the lattice constant . Many metals of interest in plasmonics have an fcc lattice structure with four conduction electrons per unit cell. In this case . For example, in silver, , and . At room temperature (), , so that is a good approximation. In this case, if and otherwise. Most numerical results shown below have been obtained in this limit. However, to illustrate the effects of finite temperature, we have also performed some computations at . Finally, the Fermi velocity is given by the equation

 vF=√2EFme=(3π2)1/3ℏmeℓ . (19)

In silver, and, correspondingly, . Electron velocities in excited states are expected to be no larger than a few times Fermi velocity, still much smaller than . This justifies the use of non-relativistic quantum mechanics.

The solution to (II) has the form of a Fourier series:

 ~ρμν(t)=∞∑s=−∞~ρ(s)μνexp(−isωt) . (20)

The expansion coefficients obey the system of equations:

 ~ρ(s)μν=δμνδs0Nμ −Λ(s)μν(ω)ℏω∑η[Gμη(~ρ(s−1)ην+~ρ(s+1)ην) −Gην(~ρ(s−1)μη+~ρ(s+1)μη)] , (21)

where

 Λ(s)μν(ω)=ωωμν−sω−iΓμν (22)

are Lorentzian spectral factors.

The optical response of the nanoparticle is determined by the quantum-mechanical expectation of its total dipole moment, which is given by

 ⟨d(t)⟩=ea∑μνΔμν~ρμν(t) . (23)

Upon substitution of (20) into (23), we obtain the expansion of into temporal Fourier harmonics. We now consider the optical response at the fundamental frequency , which describes degenerate nonlinear phenomena such as the four-wave mixing. Denoting the component of , which oscillates at the frequency , by , we can write

 ⟨dω(t)⟩=Dexp(−iωt)+c.c. , (24)

where

 D=ea∑μνΔμν~ρ(1)μν . (25)

The coefficients and the amplitude in (24) can be expanded in powers of . Namely, we can write

 D=Ωχ(Ai)Ai , (26)

where

 χ(Ai)=χ1+χ3(Ai/Aat)2+χ5(Ai/Aat)4… (27)

Here we have introduced the characteristic atomic field

 Aat=e/ℓ2 (28)

and have used the assumption that is real-valued; in the more general case, the expansion contains the terms of the form , etc. Note that the definition of in (26),(27) is somewhat unconventional. The nonlinear susceptibility , as defined in standard expositions of the subject Boyd (1992), has the dimensionality of the inverse square of the electric field, that is, of in the Gaussian system of units. Here we find it more expedient to define dimensionless coefficients , , , etc., and expand the dipole moment amplitude in powers of the dimensionless variable .

The first two coefficients in the expansion (27), and , have been computed by Rautian explicitly and are given by the following series:

 χ1 =(ea)2(ℏω)Ω∑μνNμνΛ(1)μνΔμνΔνμ , (29a) χ3 =(ea)4A2at(ℏω)3Ω∑μνηζBζημνΔμζΔζηΔηνΔνμ , (29b)

where

 Bζημν =Λ(1)μνNζμ[Λ(0)μη(Λ(1)μζ+Λ(−1)μζ)+Λ(2)μηΛ(1)μζ] +Λ(1)μνNνη[Λ(0)ζν(Λ(1)ην+Λ(−1)ην)+Λ(2)ζνΛ(1)ην] −Λ(1)μνNηζ(Λ(1)ζη+Λ(−1)ζη)(Λ(0)ζν+Λ(0)μη) −Λ(1)μνNηζΛ(1)ζη(Λ(2)ζν+Λ(2)μη) . (30)

Here . Note that all quantities inside the summation symbols are dimensionless and so are the factors in front of the summation symbols.

The expression (29b) involves a staggering 12-fold summation (recall that each composite index , , and consists of three integer indices). Rautian used the mathematical formalism of irreducible spherical tensors and -symbols to perform summation over magnetic sublevels analytically and to reduce the expression to an 8-fold summation. However, this approach does not make use of the orbital selection rules, which are explicit in (10). In Sec. III, we will use the orbital selection rules to analytically reduce (29b) to a 5-fold summation. The resultant formula is amenable to direct numerical implementation, as will be illustrated in Sec. V.

Having performed the summation over the magnetic sublevels, Rautian has evaluated the resulting series by exploiting the following two approximations:

1. Adopt the two-level approximation Hache et al. (1986). In this approximation, only the terms with and are retained in the right-hand side of (29b).

2. Assume that there are two dominant contributions to the series (29). The off-resonant (Drudean) contribution is obtained by keeping only the terms with in the Lorentzian factors . The resonant contribution is obtained by keeping only the terms with . Each contribution is then evaluated separately by replacing summation with integration.

Using the same approximations, we have reproduced Rautian’s analytical results. For , we obtained

 χ1=−14πω2pω+iΓ2[F1ω+iΓ2−ig1vF/aω2] , (31)

where

 ωp=√4πe2meℓ3 (32)

is the plasma frequency and , are dimensionless real-valued functions, which weakly depend on the parameters of the problem and are of the order of unity. More specifically, is very close to unity for all reasonable particle sizes and approaches unity asymptotically when (we have verified this numerically). In what follows, we assume that . The function depends most profoundly on the ratio . We can write, approximately,

 g1≈1κ∫11−κx3/2(x+κ)1/2dx ,  κ=ℏωEF . (33)

An analytical expression for this integral and a plot are given in the Appendix.

Equation (31) is equivalent to combining equations 3.16 and 3.23 of Ref. Rautian, 1997. On physical grounds, one can argue that these expressions are applicable only if . Indeed, in the classical Drude model, we have

 χDrude1=−14πω2pω(ω+iγ) , (34)

where is a relaxation constant. We expect the Drude model to be accurate in the limit , when the second term in the square brackets in (31) vanishes. Thus, (31) has an incorrect low-frequency asymptote. We argue that the asymptote is incorrect because the HRFR model disregards the Hartree interaction potential. This will be discussed in more detail in Sec. IV below. At this point, we assume that and expand (31) in , neglect the correction to the real part of the resulting expression, and obtain:

 χ1≈−14π(ωpω)2[1−i2Γ2+g1vF/aω] . (35)

This expression corresponds to equation 3.28 of Ref. Rautian, 1997. Note that neglecting the correction to the real part but retaining the correction to the imaginary part in the above equation is mathematically justified because the terms and can be of the same order of magnitude, as we will see below.

Comparing (35) to a similar expansion of (34), we conclude that the size-dependent relaxation constant is given by

 γ≈γ∞+g1vFa ,  γ∞=2Γ2 , (36)

where is the relaxation constant in bulk. It can be seen that the ratio plays the role of the collision frequency. The analytical result (36) is widely known and frequently used; it will be confirmed by direct numerical evaluation of (29a) below.

For the third-order nonlinear susceptibility , we obtain, with the same accuracy as above,

 χ3= Γ2Γ1α210π3(aℓ)2(λpℓ)2(ωpω)4 × [F3−i(F3γ∞ω+g3(vF/a)5ω3γ2∞)] . (37)

Here is the fine structure constant, is the wavelength at the plasma frequency ( in silver) and , is another set of dimensionless real-valued functions of the order of unity. For realistic parameters, the function varies only slightly Rautian (1997); Drachev et al. (2004a) between and ; we have taken in the numerical computations of Sec. V. The function can be approximated by the following integral:

 g3≈1κ∫11−κx5/2(x+κ)3/2dx ,  κ=ℏωEF . (38)

The approximate formula (38) applies only for . However, we are interested in the spectral region . In silver, and . This leaves us with the spectral range in which (38) is not applicable. The integral (38) can be evaluated analytically; the resulting expression and plot are given in the Appendix.

Expression (II) contains several dimensionless parameters. For a silver nanoparticle of the radius ,

 α210π3(aℓ)2(λpℓ)2≈71.6 .

The ratio is more puzzling. While can be related to the Drude relaxation constant through (36), does not enter the analytical approximations (31), (35) or the exact expression (29a). Therefore, can not be directly related to any measurement of the linear optical response. It was previously suggested Drachev et al. (2004a) that, based on the available experimental studies of non-equilibrium electron kinetics in silver Groeneveld et al. (1995); Del Fatti et al. (1998); Lehmann et al. (2000), . This ratio will be employed below.

Another interesting question is the dependence of the results on the particle radius, . It follows from the analytical approximation (35) that approaches a well-defined “bulk” limit when . The characteristic length scale is ( in silver). Of course, direct numerical evaluation of according to (29a) is expected to reveal some dependence of on , which is not contained in the analytical approximation (35), and this fact will be demonstrated below in Sec. V.2. However, it will also be demonstrated that (35) becomes very accurate in the spectral range of interest when . Thus, the HRFR model yields a result for , which is consistent with the macroscopic limit.

The situation is dramatically different in the case of the nonlinear susceptibility . It follows from (II) that . Therefore, there is no “bulk” limit for . This is an unexpected result. While some studies suggest that a positive correlation between and in a limited range of is consistent with experimental measurements Drachev et al. (2004a), we can not expect this correlation to hold for arbitrarily large values of , as this would, essentially, entail an infinite value of in bulk. Such a prediction appears to be unphysical. Of course, Rautian’s theory is not expected to apply to arbitrarily large values of because the interaction potential (2) is written in the dipole approximation and, moreover, it assumes that the electric field inside the nanoparticle is potential, that is, is a good approximation. Still, the absence of a “bulk” limit for is troublesome. We, therefore, wish to understand whether the quadratic dependence of on is a property of the HRFR model itself or an artifact of the additional approximations made in deriving the analytical expression (II). More specifically, we can state the following two hypotheses:

1. The quadratic dependence of on is an artifact of the approximations made in deriving the analytical expression (II) from (29b) [these approximations are listed explicitly between Eqs. (30) and (31)]. In this case, we can expect that direct evaluation of (29b) will not exhibit the quadratic growth.

2. The quadratic dependence of on is a property of the HRFR model itself. In particular, the absence of a “bulk” limit for can be caused by the following reasons: (i) The HRFR model neglects the retardation effects in large particles. (ii) The HRFR model does not account for the Hartree interaction potential. (In reality, interaction of the conduction electrons with the induced charge density may be important, especially, for computing nonlinear corrections.) (iii) The HRFR model makes use of a phenomenological boundary condition at the nanoparticle surface.

Verification of these hypotheses was previously hindered by the computational complexity of the problem. In what follows, we render Rautian’s theory amenable to direct numerical validation. Then we show that the analytical approximation (II) is surprisingly good. Therefore, the second hypothesis must be correct.

## Iii Rautian’s theory further developed

It is possible to simplify (29b) without adopting any approximations. To this end, we deviate from Rautian’s approach of using irreducible spherical tensors and -symbols. Instead, we directly substitute the expressions (9),(10),(11) into (29b). We use the selection rules in (10) and the following results:

 Zl≡l∑m=−lb4lm=l(4l2+1)15(4l2−1) , (39a) Sl≡l∑m=−lb2lmb2l+1,m=2l(l+1)15(2l+1) (39b)

to evaluate summations over all magnetic quantum numbers and over all orbital quantum numbers but one. This leaves us with a five-fold summation over four main quantum numbers and one orbital quantum number. After some rearrangements, we arrive at the following expression

 χ3=(ea)4A2at(ℏω)3Ω∞∑l=1(ZlPl+SlQl) , (40)

where

 Pl=∑n1,n2,n3,n4[Bn1,l,n2,l−1n3,l−1,n4,l+Bn2,l−1,n4,ln1,l,n3,l−1] ×Rn1,ln3,l−1Rn2,l−1n1,lRn4,ln2,l−1Rn3,l−1n4,l , (41a) Ql=∑n1,n2,n3,n4[Bn1,l,n2,l+1n3,l−1,n4,l+Bn2,l+1,n4,ln1,l,n3,l−1 +Bn4,l,n3,l−1n2,l+1,n1,l+Bn3,l−1,n1,ln4,l,n2,l+1] ×Rn1,ln3,l−1Rn2,l+1n1,lRn4,ln2,l+1Rn3,l−1n4,l . (41b)

This expression is exact within the HRFR model. The two-level approximation corresponds to keeping only the first term in the brackets in (40) and, further, keeping only the terms with and in (41a).

## Iv Relating the internal and the applied fields

In the HRFR model, electrons move in a given, spatially-uniform internal field (1). In practice, one is interested in the optical response of the nanoparticle to the external (applied) field. We denote the amplitude of the external field by . The two fields differ because of a charge density induced in the nanoparticle. The interaction of the conduction electrons with the induced charge density is described by the Hartree potential. However, rigorous introduction of the Hartree interaction into the HRFR model is problematic. Doing so would require the mathematical apparatus of density-functional theory. We can, however, apply here the classical concept of the depolarizing field, although this approach is less fundamental.

In the macroscopic theory, a sphere (either dielectric or conducting), when placed in a spatially-uniform, quasistatic electric field of frequency and amplitude , is polarized and acquires a dipole moment of an amplitude . The electric field inside the sphere is also spatially-uniform and has the amplitude . The induced charge accumulates at the sphere surface in a layer whose width is neglected. Under these conditions, . Note that a linear dependence between and is not assumed here. The form of the depolarizing field, , follows only from the assumption of spatial uniformity of the internal field and from the usual boundary conditions at the sphere surface. Then the Hartree interaction can be taken into account as follows.

Let us introduce the dimensionless variables and . Then we can expand in both variables:

 D=ΩAat(χ1x+χ3x|x|2+χ5x|x|4+…) , (42a) D=ΩAat(α1y+α3y|y|2+α5y|y|4+…) , (42b)

where

 x=y−DAata3=y−4π3(α1y+α3y|y|2+…) . (43)

Here we have accounted for the fact that there can be a phase shift between the internal and the external fields; therefore, and can not be real-valued simultaneously. In the theory presented above, we assume that is real-valued, and this can always be guaranteed by appropriately choosing the time origin. In this case, is expected to be complex.

The coefficients in (42a) can be found from Rautian’s theory; our task is to find the coefficients in (42b) given the constraint (43). To this end, we substitute (43) into (42a) and obtain a series in the variable . We then require that the coefficients in this series and in (42b) coincide. This yields an infinite set of equations for , the first two of which read

 χ1(1−4π3α1)=α1 , (44a) χ3(1−4π3α1)∣∣∣1−4π3α1∣∣∣2−χ14π3α3=α3 . (44b)

It is convenient to introduce the linear field enhancement factor according to

 f1=11+(4π/3)χ1=3ϵ1+2 , (45)

where is the linear dielectric permittivity. Then the solutions to (44) have the form

 α1=f1χ1 ,  α3=f21|f1|2χ3 . (46)

The factor , which relates the external and internal field amplitudes according to , is then found from

 f=xy=1−4π3(α1+α3|y|2+…) . (47)

Using (46), we find that, to first order in ,

 f=f1−4π3f21|f1|2χ3IIat . (48)

Here we have introduced the intensity of the incident beam, , and the “atomic” intensity .

Note that our approach to finding the field enhancement factor is somewhat different from that adopted in Ref. Drachev et al., 2004a, where the expansions (42a) has been truncated at the third order and the truncated expression was assumed to be exact. The results obtained in the two approaches coincide to first order in . In Ref. Drachev et al. (2004a), higher order corrections to this result have also been obtained. In our approach, these corrections depend on the higher-order coefficients , , etc., which have not been computed by Rautian.

We finally note that the phenomenological accounting for the Hartree interaction described in this section, while is necessary for comparison with the experiment, does not remove the two main difficulties of the HRFR model. Specifically, it does not fix the low-frequency limit for and does not affect the dependence of . Regarding the low-frequency limit, we note that and the internal field in the nanoparticle tends to zero in this limit. The induced macroscopic charge is localized at the sphere surface where the electric field jumps abruptly. In a more accurate microscopic picture, the width of this surface layer is finite and the electric field changes smoothly over the width of this layer. Unfortunately, the classical concept of depolarizing field can not capture surface phenomena of this kind.

## V Numerical results

### v.1 Convergence

We have computed the Bessel function zeros using the method of bisection and achieved a numerical discrepancy of the equation of less than for all values of indices. Since the function is approximately linear near its roots, we believe that we have computed with sufficiently high precision.

The summation over in (40) was truncated so that and the quadruple summation in (41) was truncated so that . A typical set of energy levels used in the summation is shown in Fig. 1 for the case , . Here the energy levels (normalized to ) are shown by dots and the horizontal axis corresponds to the orbital number, . Referring to Fig. 1, we note that has been chosen so that all states with are above the Fermi surface. Since the electron transitions occur between two states with and such that , the factors for any transition involving the states with are zero (or exponentially small at finite temperatures). It can be seen that convergence with is very fast – contribution of the terms in (40) with is either zero (at ) or exponentially small.

The choice of is a more subtle matter. Since there are no selection rules on , transitions can occur between two states (one below and one above the Fermi surface) with very different values of and, correspondingly, very different energies. However, transitions with energy gaps much larger than are suppressed by the Lorentzian factors (22). In most numerical examples, we have chosen so as to account for, at least, all transitions with the energy gaps of . Many (but not all) transitions with larger energy gaps were also accounted for. This approach yields a result with seven significant figures. However, it results in too many terms in the summation when and . For these values of parameters, we have used a smaller so as to account for, at least, all transitions with . We estimate that the relative error incurred by this truncation is .

### v.2 Linear response

We begin by considering the linear susceptibility . In computations, we use the commonly accepted parameters for silver, () and . Here the relaxation constant , which enters (29a), is determined from (see (36)). The frequencies used satisfy the condition . More specifically, the ratio varies in the range . We do not consider the frequencies above because silver exhibits strong interband absorption in that spectral range. Except when noted otherwise, all computations have been carried out at .

Fig. 2 displays the quantity computed numerically by direct evaluation of (29a) and by the Drude formula (34) with the size-corrected relaxation constant (36). The factor in (36) has been computed using the analytical formula (A). At sufficiently high frequencies, the Drude model predicts that and this behavior is reproduced for all radiuses considered with good precision. However, at smaller frequencies, there are differences between the analytical approximation and the numerical results. These differences are especially apparent for . In this case, the optical response of the sphere is, effectively, dielectric rather than metallic for . A similar behavior has been observed at (data not shown). The emergence of a dielectric response in metal nanoparticles of sufficiently small size at sufficiently low frequencies has been overlooked in the past. It occurs due to discreteness of the energy states. Consider a particle with at zero temperature. In this case, the lowest-energy electronic transition, which is allowed by Fermi statistics (that is, a transition with ), occurs between the states and . The corresponding transition frequency is . It can be seen from Fig. 2(a) that the particle becomes dielectric for .

We now turn to consideration of the relaxation phenomena. To this end, we plot in Fig. 3 the quantity

 Z=−14πωpωIm1χ1 (49)

as a function of frequency. We note that , as defined in (49), is positive for all passive materials and, in the Drude model, ; here is size-corrected. It can be seen that the analytical formula (36) captures the relaxation phenomena in the nanoparticle surprisingly well. However, as in Fig. 2(a), the analytical approximation breaks down when and . A similar breakdown was observed for (data not shown). For all other values of parameters, the numerically computed is reasonably close to the size-corrected value of and exhibits the same overall behavior. The small systematic error at higher frequencies is, most likely, caused by the approximation (33) for . It was, in fact, mentioned by Rautian that (33) is hardly accurate when .

The fine structure visible in Fig. 3(a,b) is due to discreteness of electron states. The allowed transition frequencies can be “grouped”, which results in the appearance of somewhat broader peaks, clearly seen in Fig. 3(b,c). While spectral signatures of discrete states in metal nanoparticles have been observed experimentally (including the effect of “grouping”) Drachev et al. (2004b), the positions of individual spectral peaks should not be invested with too much significance. In any realistic system, these peaks will be smoothed out by particle polydispersity, variations in shape, and by nonradiative relaxation and energy transfer to the surrounding medium.

The finite-size correction (36) to the Drude relaxation constant is widely known and used. However, the derivations of (36) have been, so far, either heuristic or relied on poorly controlled approximations. In Fig. 3, we have provided, to the best of our knowledge, the first direct, first-principle numerical verification of (36) and of its limits of applicability.

### v.3 Nonlinear response

We next turn to the nonlinear susceptibility . The same parameters for silver as before will be used. In addition, the calculations require the relaxation constant . As was mentioned above, the experimental value of can not be inferred by observing the linear optical response. It was previously suggested Drachev et al. (2004a) that . This value will be used below.

In Fig. 4, we plot the absolute value of as a function of the particle radius for and for the Frohlich frequency , and compare the results of direct numerical evaluation of (40) to the analytical approximation (II). In the case , the analytical approximation is very accurate for and gives the correct overall trend for . A systematic discrepancy of unknown origin between the approximate and the numerical results is observed for . In the case , the analytical approximation gives the correct trend in the whole range of considered. Note that, in the case , the absolute value of is dominated by and for large and small values of , respectively. When , the real part of is dominating for all values of used in the figure.

Consider first particles with . As expected, the discreteness of energy levels plays an important role in this case and results in a series of sharp maxima and minima of . As is shown in the inset of Fig. 4(a), the function is discontinuous. These discontinuities are artifacts of the zero temperature approximation. The introduction of a finite temperature () removes the discontinuities (see the inset) but does not eliminate the fine structure of the curve. Note, however, that the computations have been carried out with a very fine step in , which is, arguably, unphysical: the parameter in a real nanosphere can change only in quantized steps of the order of the lattice constant ( in silver). Moreover, the fine structure of is unlikely to be observable experimentally due to the unavoidable effects of particle polydispersity. Therefore, the general trend given by the analytical approximation (II) can be a more realistic estimate of for .

Next consider the large- behavior. For , the analytical and the “exact” formulas yield results, which are scarcely distinguishable. In particular, the quadratic growth of with has been confirmed up to in the case – the largest radius for which numerical evaluation of (40) is still feasible. This confirms Hypothesis 2 stated above, namely, that the quadratic growth of is a property of the HRFR model itself rather than of the additional approximations, which were made to derive the analytical results.

In Fig. 5, we study the dependence of on the frequency for fixed values of . It can be seen that the accuracy of Rautian’s approximation improves for larger particles and higher frequencies. At , the approximation is nearly perfect in the full spectral range considered.

### v.4 Magnitude of the nonlinear effect and comparison with the classical theory of electron confinement

In the previous subsection, we have plotted the coefficient , which appears in the expansion (27). The dimensionless parameter of this expansion, , contains the amplitude of the internal electric field, . However, it is the amplitude of the external (applied) field, , which can be directly controlled in an experiment. The incident beam intensity is given by . We can use the results of Sec. IV to write