Theoretical and numerical investigation of the sizedependent optical effects in metal nanoparticles
Abstract
We further develop the theory of quantum finitesize 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 infinitelyhigh spherical potential well. This model (referred to as the HRFR model below) yields mathematical expressions for the linear and the thirdorder 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 thirdorder 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 (19282009) 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 surfaceenhanced 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 finitesize 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 thirdorder nonlinear polarizability involves a twelvefold 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 eightfold 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 fivefold 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 thirdorder nonlinear polarizability from eightfold to fivefold. 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 finitesize effects in conducting nanoparticles Rautian (1997). The physical system under consideration is a gas of noninteracting electrons placed inside a spherical, infinitely deep potential well of radius and subjected to harmonicallyoscillating, spatiallyuniform electric field
(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 electronfield interaction can be described in the dipole approximation by the timedependent operator
(2) 
where is the dipole moment operator. Under the additional assumption that is linearly polarized with a purely real amplitude , we can write
(3) 
where
(4) 
Rautian made use of the interaction representation in which the wave function is expanded in the basis of the unperturbed Hamiltonian eigenstates. The singleelectron unperturbed states are
(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
(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
(7) 
where
(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
(9) 
where
(10) 
are the Kronecker deltasymbols, and
(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
(12) 
The density matrix of the system, , can be written in the form , where are the transition frequencies and is the socalled slowvarying amplitude, which obeys the following master equation Rautian and Shalagin (1991):
(13) 
Here are the equilibrium state populations and are phenomenological relaxation constants. Following Rautian, we assume that
(14) 
Eq. (14) is the least complex assumption on , which still distinguishes the relaxation rates for the diagonal and offdiagonal 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
(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 wellknown analytical formula for the Fermi’s energy,
(16) 
holds with a good accuracy. Here is the characteristic atomic scale, defined by the relation
(17) 
where
(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
(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 nonrelativistic quantum mechanics.
The solution to (II) has the form of a Fourier series:
(20) 
The expansion coefficients obey the system of equations:
(21) 
where
(22) 
are Lorentzian spectral factors.
The optical response of the nanoparticle is determined by the quantummechanical expectation of its total dipole moment, which is given by
(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 fourwave mixing. Denoting the component of , which oscillates at the frequency , by , we can write
(24) 
where
(25) 
The coefficients and the amplitude in (24) can be expanded in powers of . Namely, we can write
(26) 
where
(27) 
Here we have introduced the characteristic atomic field
(28) 
and have used the assumption that is realvalued; 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:
(29a)  
(29b) 
where
(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 12fold 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 8fold 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 5fold 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:

Assume that there are two dominant contributions to the series (29). The offresonant (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
(31) 
where
(32) 
is the plasma frequency and , are dimensionless realvalued 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,
(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
(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 lowfrequency 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:
(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 sizedependent relaxation constant is given by
(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 thirdorder nonlinear susceptibility , we obtain, with the same accuracy as above,
(37) 
Here is the fine structure constant, is the wavelength at the plasma frequency ( in silver) and , is another set of dimensionless realvalued 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:
(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 ,
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 nonequilibrium 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 welldefined “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:

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:
(39a)  
(39b) 
to evaluate summations over all magnetic quantum numbers and over all orbital quantum numbers but one. This leaves us with a fivefold summation over four main quantum numbers and one orbital quantum number. After some rearrangements, we arrive at the following expression
(40) 
where
(41a)  
(41b) 
This expression is exact within the HRFR model. The twolevel 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, spatiallyuniform 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 densityfunctional 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 spatiallyuniform, 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 spatiallyuniform 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:
(42a)  
(42b) 
where
(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 realvalued simultaneously. In the theory presented above, we assume that is realvalued, 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
(44a)  
(44b) 
It is convenient to introduce the linear field enhancement factor according to
(45) 
where is the linear dielectric permittivity. Then the solutions to (44) have the form
(46) 
The factor , which relates the external and internal field amplitudes according to , is then found from
(47) 
Using (46), we find that, to first order in ,
(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 higherorder 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 lowfrequency limit for and does not affect the dependence of . Regarding the lowfrequency 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 sizecorrected 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 lowestenergy 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
(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 sizecorrected. 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 sizecorrected 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 finitesize 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, firstprinciple 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